Phase-Field Modeling of Freeze Concentration of Protein Solutions

Bulk solutions of therapeutic proteins are often frozen for long-term storage. During the freezing process, proteins in liquid solution redistribute and segregate in the interstitial space between ice crystals. This is due to solute exclusion from ice crystals, higher viscosity of the concentrated solution, and space confinement between crystals. Such segregation may have a negative impact on the native conformation of protein molecules. To better understand the mechanisms, we developed a phase-field model to describe the growth of ice crystals and the dynamics of freeze concentration at the mesoscale based on mean field approximation of solute concentration and the underlying heat, mass and momentum transport phenomena. The model focuses on evolution of the interfaces between liquid solution and ice crystals, and the degree of solute concentration due to partition, diffusive, and convective effects. The growth of crystals is driven by cooling of the bulk solution, but suppressed by a higher solute concentration due to increase of solution viscosity, decrease of freezing point, and the release of latent heat. The results demonstrate the interplay of solute exclusion, space confinement, heat transfer, coalescence of crystals, and the dynamic formation of narrow gaps between crystals and Plateau border areas along with correlations of thermophysical properties in the supercooled regime.


Introduction
Understanding the mechanisms and process impacts of freezing on protein solutions is critical in the pharmaceutical industry for the manufacturing of high quality biologics. During the freezing process, the heat, mass, and momentum transport influence the size, shape, and growth rate of ice crystals as well as the distribution of protein molecules in between. When a frozen product is freeze-dried, the ice crystal size and distribution resulting from freezing determine the sublimation kinetics during primary drying, and eventually the stability of the biologics [1][2][3][4]. Freeze concentration occurs due to the exclusion of protein and other solutes from ice crystals and space confinement between ice/freeze concentrate interfaces. Overall, the seeding or nucleation pattern and the freezing rate control the growth of ice crystals. Current understanding is that at fast cooling rates, dendritic crystals branch out quickly, trapping proteins among side branches of crystals to yield a more uniform protein distribution. In contrast, at slow cooling rates, protein molecules redistribute while cellular (instead of dendritic) crystal structure forms, resulting in a more significant concentration polarization as crystals continue to grow and coalesce. Freeze concentration and its impact on protein stability is more prominent during large scale freezing. Although rapid freezing is likely advantageous because of more uniform protein distribution in the frozen product, it can be difficult to achieve at the manufacturing scale. Experimentally it has been observed that upon freezing the local thermal and convective effects induced by phase transition, cooling, gravity, control of ice nucleation for freeze-drying, and protein residence time in the freeze concentrates all have impacts on the process performance and the distribution and stability of proteins [5][6][7][8].
Qualitative and quantitative analyses of the relevant dynamics in freeze concentration are complicated by the coupling of thermal and multi-component mass transport, phase transition physics, protein interactions with co-solutes, morphological change and the space confinement effect, thermophysical properties at the supercooled regime, and molecular crowding and protein-protein interactions. Although casting in conventional thermal manufacturing process is a well-established technology, where the relevant physics such as dendritic growth, patern formation, and microsegregation of pure and alloy systems has been studied extensively [9][10][11][12][13][14][15], the existing theoretical models and experimental results are not readily applicable to describe freezing of protein biologics. Even at the laboratory scale, only a few investigations on the phase transition of protein solutions have focused on fundamental heat transfer, phase change, and computational fluid dynamics simulations. Nakagawa et al. [16] proposed a thermal conduction model using the classical enthalpy method to trace phase change and temperature distribution during freezing of mannitol and BSA based solutions in a vial. Without further description of the growth dynamics of ice crystals, the estimated sizes of crystals and the freeze-dried layer permeability correlate well with the freezing front rate and temperature gradient in the frozen zone. Radmanovic et al. [17] modeled the heat transfer involved in freezing of monoclonal antibody solution and evaluated the product quality based on the size and polydispersity of the protein aggregates. Roessl et al. applied volume-of-fluid computational fluid dynamics to simulate the freezing process and compared the results with measurements of temperature and cryoconcentration fields [18]. On a smaller scale, Butler observed the solute profile at the ice/freeze concentrate interface due to freeze concentration using optical interferometry [19]. Kaempfer and Plapp developed phase-field modeling of sublimation dynamics and relevant vapor and thermal transport on the microstructure of dry snow and compared with microtomography images [20]. van der Sman developed a phase-field model to simulate growth and suppression of ice crystals in sugar solutions [21], and successfully incorporated polymer mean field theory into the phase transition dynamics with simplified thermal effect. Huang et al. [22] revealed directional growth of anisotropic ice crystals in ceramic colloidal suspensions, with the resulting porous material having potential for biomedical applications such as functional materials for implants. These relevant studies have only focused on either thermal or composition with phase transition in the freezing process, however, at the mesoscale in particular, multiple physics are often convoluted in the phase transition dynamics, requiring the inclusion of the effects of thermal, fluid flow, thermal mechanical response, composition and changes of thermophysical properties into the analysis [23,24]. Here we propose a new model to integrate these effects and provide mathematical descriptions and simulation of the freeze concentration effect.
The phase-field method is an Eulerian approach that naturally resolves fusion, merging, splitting, topological and morphological changes of dynamic interfaces and thermodynamic states of the materials involved. The underlying concept of a smooth transition between phases and the contribution of non-local or gradient energy was initially introduced by van der Waals in 1893 [25] to describe liquid-vapor phase transition over a diffuse interface. A general phase transition process is driven by the increase of entropy or reduction of free energy. The gradient effect due to spatial variation of the phase-field function or order parameter eventually leads to a continuous phase-field equation that describes the local state of the material, and can clearly distinguish the location of the transition interface between phases. This approach was further developed in many research areas in physics and material sciences, known as Ginzburg-Landau free energy theory to describe superconductivity, Allen-Cahn model to descibe transition dynamics of non-conserved phase field or order parameter, Cahn-Hilliard model for conserved parameters such as density or concentration, and Model A/B for isothermal and H for nonisothermal fluids near a critical state [26,27]. The phase-field function is uniform within a homogeneous phase, but has a narrow and smooth transition profile across the interface between phases. Rigorous derivation of phase-field equations for non-isothermal systems is based on the principles of irreversible thermodynamics [28,29], in which the entropy functional of the system, entropy transport equation, and the 2nd law of thermodynamics are the starting point of derivation. The phase-field method has been widely applied to model dendritic pattern formation in material sciences since the works of Kobayashi [9], Warren and Boettinger [10], Murray, Wheeler, and Glicksman [11], and Karma and Rappel [12,13] in the 1990s. The transient evolution of the solidification interfaces is driven by thermal gradient, interfacial instability, Gibbs-Thompson kinetics, and anisotropic interfacial energy along with hydrodynamics effect, which all play an important role in pattern formation, solute distribution, and coarsening and remelting of dendritic microstructures [14,15,30,31]. The method has been further extended to address multi-component problems with and without convective effect, and for a variety of applications in thermal fluid sciences and applied mechanics including geohydraulics, fracture mechanics, multiphase flows, and reaction engineering etc. [32][33][34][35]. In this paper we apply phase-field method to model freeze concentration of protein solutions and observe concentration polarization without considering anisotropic dendritic formation. Instead of anisotropic interfacial energy, the focus here is to couple the evolving dynamics of crystals with fluid flow, heat transfer, protein diffusion, and the convective transport induced by gravity effect and density variation across the moving ice/freeze-concentrate interface. In the mesoscale domain, the cooling effect is simplified by using a representative bulk cooling rate, which reduces the computational cost and potentially can be extended and integrated with multiscale analysis including container configuration and experimental settings at a larger scale.

Theoretical Analysis
A challenge in the theoretical analysis is to compile relevant thermophysical properties of protein solutions as in crystal or supercooled liquid state, including density, specific heat, thermal conductivity, mass diffusivity, and dynamic viscosity. These physicochemical properties in general depend on temperature and solute concentration, and the available measurements are very limited. Here we consider sucrose solution to demonstrate the modeling results based on available data. Although sucrose has very different properties than proteins, it is a very common excipient used as a cryoprotectant in protein formulations, freeze concentration of sucrose may also impact protein stability, and thus the investigation is very relevant to the biopharmaceutical industry. The temperature-composition phase diagram of a binary model protein solution (e.g., sucrose and water) at ambient pressure is shown in Figure 1a [6,36]. The freezing curve (liquidus) shows freezing point depression, and the extension of the freezing curve intercepts the glass transition curve at the T g point, where the solute concentration reaches the maximum value. Freezing of protein solutions often starts from room temperature to a supercooled liquid state ( Figure 1b) with a few to tens of degree Kelvin below the equilibrium freezing temperature at a given composition. The temperature raises due to heat release from freezing or solidification, possibly near the equilibrium freezing temperature and eventually drops due to the cooling rate applied to the material volume. In this study, assuming that the heterogeneous nucleation can be controlled by random seeding, ultrasound, electrical, or other freezing methods [3], the supercooled temperature and the size and locations of seeding are predetermined as the initial conditions, and then the simulation proceeds as a non-equilibrium freezing process. The continuous increase in viscosity and change of local freezing point from freeze concentration are considered. The following assumptions are made to facilitate and simplify the theoretical analysis and computation: (i) ice/freeze concentrate interfacial energy is assumed half of the free surface energy of water, (ii) the interfacial energy is assumed isotropic or no preferential directions for the formation of dendritic pattern, which can be considered in further investigations, (iii) thermal expansion, thermal stress, or elasticity in the crystal phase are neglected, and in fact ice crystals are treated as a highly viscous fluid with dynamic viscosity at least four orders of magnitude larger than fluid, (iv) entanglement or crystallization of the solute molecules at high concentrations is not considered as the corresponding mean field theory is much more complicated, (v) adsorption and molecular interactions of proteins with the ice/freeze concentrate interfaces are not considered, (vi) thermal radiation and van der Waals force between nearby interfaces may have effects on heat transfer and interfacial dynamics at the mesoscale, but are assumed negligible in this work.

Thermodynamics Approach
Irreversible thermodynamics provides a rigorous route to derive the phase-field governing equations for a non-isothermal system involved in thermal manufactuing processes, in which the coupling of phase-field functions with other transport equations is important. Here we consider phase field φ s (r, t) as a function of location r and time t to distinguish the solid (ice crystal, φ s = 1) and liquid (liquid protein solution, φ s = −1) phases, and the ice/freeze concentrate interface has a narrrow but smooth transition of φ s between −1 to 1. The 2nd phase-field function is the protein volume fraction φ c (r, t) as a conserved parameter at the range of 0-1. To derive the thermodynamically consistent formulation we follow the approach proposed by Penrose and Fife [28] and Wang et al. [29]. Starting from the entropy functional of the material volume that takes the contributions of the entropy density in the bulk phase and the non-local or gradient effects due to spatial variation of the phase fields φ s and φ c into account, we have where ρ is density, s is specific entropy as a function of specific internal energy e and the two phase-field functions, and ξ s and ξ c are the assumed constant coefficients corresponding to the gradient effects [15,29] . The gradient coefficient ξ s is connected with reference density ρ 0 , solid-liquid energy barrier coefficient h s , and the characteristic interfacial thickness W s as The solid-liquid interfacial energy γ s or the equivalent surface excess energy across the smooth interface is further associated with the above parameters through a one-dimensional approximation, expressed as where the factor 2 √ 2/3 comes from the hyperbolic tangent function to approximate φ s in equilibrium, which has a defined value from −1 to 1, and the reference temperature T 0 here represents the equilibrium freezing temperature of pure water. We assume γ s (1/2)γ g at the ice/freeze concentrate interface, where γ g is the surface tension of the free surface of water at T 0 , and all intensive thermodynamic properties are defined as per unit mass quantity.
By applying the Reynolds transport theorem to the entropy of the system, the general differential entropy transport equation can be written as where t is time, D/Dt ≡ ∂/∂t + v · ∇ indicates material derivative, J s represents entropy flux,Γ is the rate of entropy production, andΩ is a heat sink to mimic the cooling from the bulk fluid flow around the mesoscale volume. Considering internal energy e = e(s, φ s , φ c ), de = Tds + (∂e/∂φ s )dφ s + (∂e/∂φ c )dφ c , and thus ds = de/T − (∂e/∂φ s )dφ s /T − (∂e/∂φ c )dφ c /T, the time rate of change of entropy Ds/Dt in the transport equation (Equation (4)) can be replaced by the material derivatives of internal energy and the two phase fields, expressed as where the three material derivatives on the right-hand side eventually lead to the thermal energy equation and two phase-field equations. Substituting Equation (5) into (4) and by separating the contributions of entropy flux and entropy production terms due to heat conduction and the evolution of phase fields, that is, and whereq = −k T ∇T is the conduction heat flux and k T is the thermal conductivity. The resulting positive entropy production rate can be obtained and expressed as in which the viscous dissipation and capillary work have been neglected due to relatively low dissipated energy compared to heat conduction (small Brinkman number). To accommodate the positive entropy production, one can conclude the non-conserved phase-field equation for the growth of ice crystals as and the conserved phase-field equation for tracking the single-component protein volume fraction: where J φ c is the species flux, δS /δφ c is the first variation of the entropy functional, and M s and M c are assumed positive interfacial mobility coefficients. Equations (10) and (11) are essentially the Allen-Cahn (2nd-order) and Cahn-Hilliard (4th-order) equations for a non-isothermal system, respectively. The mobility M s may be acquired empirically or from scaling and asymptotic estimation of the the evolution kinetics of the diffuse interface [12], whereas M c is associated with the classical Fickian diffusion coefficient.

Internal Energy and Free Energy
The specific internal energy including solid and liquid phases and mixing energy of proteins in solutions can be approximated by where e s represents the internal energy of solid to liquid phases, R is gas constant, T is temperature, and χ is the Flory's interaction parameter of regular solutions, and along with the G function the mixing term indicates the increase of internal energy due to mixing of water and proteins in either liquid or crystal phases. Here we assume χ > 0 (net repulsion between interacting species). The typical G function for the mixing effect can be replaced by a double-well potential for keeping φ c within the range of 0 to 1 in the phase-field computation, written as Furthermore, the internal energy e s can be formulated as where the subscript s and indicate homogeneous solid and liquid phases, respectively, and L a is the assumed constant latent heat based on the reference freezing temperature T 0 of pure water at equilibrium. Developed by Wang et al. [29], the 5th-order interpolation polynomial P(φ s ) describes a smooth transition of the internal energy between the liquid (P = 1, φ s = −1) and solid (P = 0, φ s = 1) phases. Here a P function that satisfies P = P = 0 at φ s = ±1 is which gives P(1) = 0 and P(−1) = 1 so that e s (T, 1) = e s (T) and e s (T, −1) = e (T). The P function can be widely used for other thermodynamic properties for a smooth transition between solid and liquid phases. For example, we define Flory's parameter as with the interaction parameters χ s > χ , which adjusts the partition or exclusion effect for the proteins at the freezing front. Proteins are mostly excluded from ice crystals due to relatively high mixing energy, and are soluble (within the solubility limit) in the liquid phase (Figure 1a). Similar to the model proposed for binary alloy systems [15,32], in a regular protein solution the free energy density can be obtained by superposing the contributions of pure solid and liquid water, proteins, and the mixing entropy and enthalpy effects based on Flory-Huggins mean field theory for polymer solutions [21,37], applicable for polymer solutions with temperature above T g . Here we incorporate the G function into to the free energy as where N accommodates the size effect as the protein-to-water ratio of partial molar volume, and the Flory's χ parameter, assumed independent of solute volume fraction, controls the energy barrier in the mixing enthalpy term to ensure exclusion of proteins from the ice crystals. The free energy of ice and liquid water f s drives the phase transition based on the difference of local temperature to the equlibrium freezing point T eq , which follows the Gibbs-Thompson effect. Here we assume the interfacial curvature effect is negligible in the mesoscale model, and the equilibrium freezing temperature, T eq = T eq (φ c ), follows the liquidus line. From the Gibbs-Helmholtz relation, the free energy that describes the latent heat effect for the solid-liquid phase transition dynamics [29] can be expressed as where the first term is the driving force for freezing as temperature is different from the reference point T eq , whereas the 2nd term is often described by a double-well potential with free energy minima for the the ice or liquid phases at T = T eq , written as where h s is the energy barrier coefficient for the phase transition. The free energy has two minima at φ s = −1 (liquid) and φ s = 1 (solid) at the equilibrium freezing temperature T eq . From the internal energy and free energy above, the derivative of internal energy in Equation (10) thus can be expanded and described by the following expression: Therefore, Equation (10) becomes where the transient evolution of the phase field φ s is now determined by the balance of diffusion, latent heat, double-well potential, and the effect of freezing point depression at higher protein volume fractions. The diffusion effect tends to smear out the interface, the double-well potential keeps the bulk phases separated, whereas the temperature difference along with the latent heat effect provides a driving force for freezing. The increase of protein concentration at the interface suppresses (slows down) the growth of ice crystals due to freezing point depression and the increase of local viscosity at higher solute concentration. Similarly, the derivative of internal energy in Equation (11) can be expanded as where the free energy difference between the ground states f c and f s has been neglected. Therefore the phase-field Equation (11) that governs the transient distribution of protein volume fraction reduces to The above Cahn-Hilliard type model recovers to a conserved, Fickian diffusion equation by defining the mobility as The diffusivity D in liquid and crystal phases can be scaled (represented by tilde over a variable) and interpolated as where D is a concentration-and temperature-dependent protein diffusivity in liquid solution and D s is a constant diffusivity in the crystal phase, and both are scaled by reference value D 0 . The subscripts s and hereafter denote the solid and liquid phases, respectively. We assume D s D and the reference diffusivity D 0 = D (φ c → 0, T = T 0 ).

Thermal Energy and Momentum Equations
From Equations (12)- (14), the rate of change of specific internal energy can be partitioned using the P function and expressed as and with the Fourier law of conduction the thermal energy equation (Equation (6)) can be developed and written as where c p is specific heat, k T is thermal conductivity, and both are defined as phase φ s and temperature dependent properties: and Therefore the interpolation P function defined previously is applied here for the transition of properties between solid and liquid phases. The corresponding reference values are defined as c p 0 = c p (T → T 0 ) and k T 0 = k T (T → T 0 ) based on pure water. The volumetric heat sinkΩ drives the overall freezing process for the simplified mesoscale model without considering actual heat transfer surfaces or wall conditions. Here we apply a linear cooling model: with a constant cooling rate β f to represent a uniform influence on the mesoscale domain based on cooling condition for the bulk environment. The parameter β f may be modified for multiscale simulation that couples the representative local model with a full scale CFD simulation at the vial or bottle level. In principle, directional cooling can also be applied to inspect the anisotropic growth of ice crystals. The fluid flow within the liquid domain or in the confined interstitial space between ice crystals is primarily driven by the blowing mass flux due to density change upon freezing (or suction upon melting). The continuity equation based on local mass conservation is coupled with the phase field φ s and can be expressed as where v is velocity, ρ is density, and the material derivative D/Dt takes convective effect on the evolution of the phase field into account. The local density is defined as where the scaled density of ice crystals ρ s and supercooled water ρ are temperature dependent, and ρ 0 = ρ (φ c → 0, T → T 0 ). The Navier-Stokes momentum equation with Boussinesq approximation for the buoyancy effect is given by where g is gravity acceleration, and σ v is the viscous stress for the assumed Newtonian fluid, written as where the dynamic viscosity η is defined as with an assumed constant viscosity for the crystal phase η s η and a reference value η 0 = η (φ c → 0, T → T 0 ). At the mesoscale, the buoyancy effect introduces a small updraft and shows tendency of collective motion of ice crystals. Furthermore, in the proposed phase-field model, instead of elastic solid material the ice crystals are treated as a highly viscous fluid with viscosity at least four orders of magnitude higher than liquid viscosity. To further simplify the computation, the fluid flow is assumed quasi-incompressible, so that the density variation is decoupled from pressure field except for the phase change interface and Boussinesq approximation. The modified pressure that incorporates the interfacial dynamics can be approximated by the pressure Poisson equation: where the last term is considered an additional contribution to pressure due the non-solenoidal velocity field.

Scaling and Computation
The system has two length scales involved, the physical domain size 2πL and the apparent interfacial thickness W s for computing phase transition. This yields six characteristic time scales based on solid-liquid phase transition, protein diffusion, thermal diffusion, viscous diffusion, convection, and the assumed freezing rate. These time scales are estimated by the following expressions: respectively, where U is the characteristic velocity and T represents the characteristic subcooled temperature. All governing equations are scaled by length L and phase transition time τ s . The phase field φ s , protein volume fraction φ c , χ parameter, G and P functions are already normalized. Temperature is scaled by the subcooled temperature as and pressure and stress are scaled by the viscous effect based on the reference viscosity. From length scale L and time scale τ s , the resulting scaled phase-field equations reduce to and where the Peclet number P e , phase-change number Λ s , and the Cahn-Hilliard number C h are defined as Specifically the Peclet number measures the convective to diffusive effects on the redistribution of proteins, the phase-change number is the ratio of latent heat and interfacial energy, and the Cahn-Hilliard number controls the characteristic interfacial thickness to the length scale.
Furthermore, the scaled thermal energy equation can be expressed as where Peclet number P e has been defined above, interfacial Lewis number L e measures the ratio of interfacial evolution to thermal diffusion time scales, two Stefan numbers S te1 and S te2 compare the sensible heat and partition effect to the latent heat, respectively, and Ω is the scaled cooling rate, defined as respectively. The scaled continuity equation is The scaled Navier-Stokes momentum equation can be expressed as whereê g indicates the downward direction of gravity acceleration. The Schmidt number S c compares the phase transition to viscous time scales, Reynolds number R e characterizes the inertial to viscous effects and the local Grashof number G r measures the buoyancy to viscous forces adjusted by the local density, defined as respectively. The pressure can be obtained from the Poission equation (Equation (36)) and because of the small Reynolds number, the nonlinear inertial effect can be neglected. In summary, the phase-field equations describe the transient evolution of the freezing front and redistribution of solutes in the liquid solution, the thermal energy equation determines the transient temperature distribution, and the continuity, momentum, and auxiliary pressure Poisson equations govern the interplay of fluid flow with phase fields and thermal energy. These equations are discretized and computed by our in-house Matlab codes with algorithms developed based on Fourier spectral method [38]. The two-dimensional computational domain includes 1600 × 1600 uniform collocation points and periodic boundary conditions. The built-in Matlab functions of fast Fourier and inverse Fourier transforms are applied in the computation. The temporal discretization is based on forward Euler scheme with a uniform and scaled time step 2 × 10 −4 throughout the transient process. A pseudo spectral scheme is applied to the nonlinear terms of the governing equations. As the phase-field method naturally resolves deformation and coalescence of the crystals, there is no inter/extrapolation or any smooth or adaptive schemes applied to the interface or spatial derivatives in the computation. The initial conditions are based on predetermined nucleation sites and seeding size, and a uniform supercooled temperature for the onset of crystaliztion. The model is not limited to this simplified initial condition and can be further developed to incorporate various nucleation models. A semi-implicit scheme is applied to spatial or spectral domain for variable transport coefficients primarily associated with the 2nd-order derivatives [39]. The two-dimensional simulation results for the cases presented here are manageable by a conventional desktop computer with 16 gb ram, and for each case the transient simulation can be completed overnight.

Material Properties
First, water expands upon freezing and due to a smaller mass density of ice crystals than liquid solution, mass flux is generated at the freezing front toward the liquid side. We describe this effect by incorporating the temperature-dependent densities for the supercooled water and ice crystals based on the compiled dataset from the CRC Handbook [40], and formulated in terms of dimensional values in MKS (meter, kilogram, second) units as The sucrose density ρ c is assumed a constant 1587 kg/m 3 . Second, the specific heat of ice crystals is assumed independent of solute volume fraction due to the strong exclusion effect, and decreases as temperature decreases, with data [40] correlated as The specific heat of supercooled sucrose solution [8] is given by Third, the thermal conductivity of ice crystals increases as temperature decreases [40], approximated by The thermal conductivity of supercooled sucrose solution [8] is approximated by Several reference properties at temperature T 0 are denoted by the subscript 0 with test values listed in Table 1.
Fourth, the dynamic viscosity of water increases when temperature decreases. The temperaturedependent viscosity for the supercooled water can be correlated by the Vogel-Fulcher-Tamman (VTF) model [41]: in which the prefactor 4.442 × 10 −5 is the viscosity at temperature 168.9 K (>T g 136 K). Note that the scaling factor for viscosity η 0 (∼0.0018 Pa·s) is defined at the freezing temperature T 0 . At higher solute concentration, the mobility or self diffusivity of solutes can be significantly reduced due to molecular crowding effect and protein-protein interactions. Empirical model for the concentration-dependent diffusivity in the supercooled regime is required for a specific protein of interest. However, due to the lack of such data in the supercooled regime, here we compose the VTF model and the concentration-dependent viscosity of sucrose solutions measured at T 0 = 273.15 K [42]. With a conversion from weight percentage to volume fraction, the viscosity is correlated as where the exponential term follows the Mooney's viscosity model. In the crystal domain we assume η s 10 4 η 0 .
Finally, the reference solute diffusivity D 0 in sucrose solution at T 0 is about 2.1 × 10 −10 [19], from which the temperature and concentration dependency can be estimated by the proportionality to temperature and viscosity based on the Stokes-Einstein relation and described as In the crystal domain we assume D s 10 −4 D 0 . The liquidus temperature that accommodates the freezing-point depression effect is correlated from the sucrose-water phase diagram [43] as A summary of the parameters and dimensionless groups used for the test cases are listed in Tables 1 and 2. Table 1. Parameters, characteristic time scales, and reference properties for ice crystal and supercooled water used in the test cases. Parameter N is approximated by sucrose in water solution.   Figure 2 demonstrates the computational results at a scaled time instant t = 1.5, showing the fully coupled thermal fluid and phase transition dynamics driven by a supercooled temperature T = −1.0 initially, and then by a uniform volumetric cooling rate with β f = 100 K/s. The initial volume fraction is set to φ c = 0.05 in the liquid mixture and φ c = 0 for the seeding crystals. The phase field φ s (Figure 2a) has a sharp interface that clearly distinguishes the crystal domains from the liquid mixture. The solute volume fraction φ c (Figure 2b) shows that solutes are pushed away by the interfaces as they grow into the liquid mixture. This is due to solute exclusion from the crystals, and the solute concentration is further enhanced between close-approached crystals as the space is confined. The small dark-blue spot near the center of each crystal implies the seeding site where φ c = 0 initially. The narrow gap between ice/freeze concentrate interfaces further hinders diffusive mass transport with lower mass diffusivity due to the increase of viscosity at higher solute concentration. This is an additional influence on mass transfer due to temperature adjustment. As crystals approaches each other, the moving interfaces slow down and flatten, which creates a narrow gap in between. The slow down is due to higher local temperature from the exothermic effect upon freezing and the freezing point depression as local concentration increases. The solutes are mostly trapped in the narrow channels and the Plateau border areas between the interfaces. The temperature variation primarily due to heat conduction is shown in Figure 2c and the distribution is relatively diffuse compared with the phase and concentration fields. Near the freezing front the latent heat is released and this exothermic effect increases the local temperature significantly. Although thermal diffusion is much faster than mass transfer and phase transition, temperature variation over the computational domain is significant and thus an isothermal assumption often applied in solidification process simulation is not suitable here. The thermal convection effect is less important in the mesoscale regime due to relatively low Reynolds number and small Peclet and Grashof numbers. Because of smaller mass density of crystal phase than liquid solution, the volumetric flow rate (Figure 2d) creates a blowing mass flux from the ice/freeze concentrate interface towards the liquid solution. The local flux slightly enhances the mass transfer and the removal of latent heat.

Results and Discussion
Furthermore, the scaled density ( Figure 3a) and liquid viscosity (Figure 3b) corresponding to Figure 2 have taken local temperature and solute volume fraction into account. This is important for describing reduced mobility of solutes at higher concentration and lower temperature. Again the interfaces can be slowed down significantly due to higher solution viscosity, freezing point depression, and higher local temperature due to the release of latent heat. Although the increase of temperature slightly reduces the viscosity and enhances the mass diffusivity, the higher temperature lowers the driving force for the phase transition significantly. From the velocity field ( Figure 3c) one can observe a tendency of collective upward motion of crystals. This is due to buoyant effect in addition to the local mass flux. Finally, the vorticity field (Figure 3d) shows stronger circulation around the interstitial liquid space, which correlates well with the phase transition and higher temperature region. Figure 4 shows transient evolution of crystal-solution interfaces and freeze concentration between these crystals. The total number and locations of nuclei and the onset temperature ( T = −1.0 or dimensional temperature T = T 0 − 10 K) of crystallization are predefined, same as in Figures 2 and 3. The dynamic space confinement and freeze concentration are clearly resolved along with the fluid flow surrounding crystals, and for the movement of the crystals. In principle, as cooling continues, the non-equilibrium process may lead to fusion or coalescence of crystals. Solute exclusion from the crystal phase is controlled by the Flory's interaction parameter χ, which can be an empirical property from experimental analysis. The topological structure of the growth of crystals is quite similar to those found in bubble or foam materials. The drainage of the fluid within the thin film and the Plateau border area is however not as significant as in foam materials, in which both the continuous and dispersed phases have fluid-like behaviors. At t = 0.5, the scaled temperature spans a wide range from −1.8 to −0.2, this is because of the rapid cooling that competes with the release of latent heat. The raised temperature between two nearby crystals is clearly observed, with decelerated phase transition. At t = 2.5, highest local concentration increases about seven-fold within the narrow gap between crystals, with maximum φ c about 0.35. In principle, the concentration could further increase as cooling continues till T g , but the relevant physics involving strong protein-protein interactions are beyond the scope of this model. In Figure 5 we demonstrate a case with five times more nucleation sites (25 seeds in the computational domain) and lower initial protein volume fraction (φ c = 0.02), but under the same cooling condition and onset temperature of crystallization. By arranging higher seeding density presumable on the left-hand side of the domain and with the same seeding size, the growth of crystals is relatively uneven as expected. At the early stage t = 0.5 we observe a greatly increased temperature because of the larger amount of phase transition that happens simultaneously so that the released latent heat can not be removed efficiently. On the 3rd row, as many crystals grow spontaneously the collective updraft appears due to the buoyant effect, which also induces a locally sinking flow on the right-hand side of the domain that has more liquid solution and less crystals. At the later stage t = 2.8, although freeze concentration in the interstitial space would suppress freezing, many crystals still coalesce with surrounding crystals. This creates several spots that have almost sixfold increase of solute volume fraction. Coalescence creates spotty areas with solutes trapped inside. Overall, the crystal phase grows faster around high-density nucleation sites, but such growth is limited by the large amount of latent heat released, whereas the loosely seeding site has more room for the crystals to grow in the later stage, showing an interesting growth pattern, solute composition, temperature map, and the surrounding fluid flows. All these effects, along with the initial conditions determine the growth of ice crystals and the degree of freeze concentration. We observe that the onset temperature of crystallization is relatively less important because even the initial temperature is near the homogeneous nucleation temperature about −40 • C, the amount of sensible heat c p ∆T is about an order of magnitude smaller than the latent heat to be released. Therefore, without rapid cooling the sensible heat itself is insufficient to grow the crystals. The onset temperature may have influence on the sporadic precipitation or growth of dendritic microstructure, which will be investigated in future works.
In summary, freezing is a complicated exothermic process and the design and control of the process has a large parameter space. This preliminary and qualitative investigation has demonstrated a few important mechanisms involved in the freezing process. Dynamic freeze concentration due to volume exclusion and space confinement effect is demonstrated for a mesoscale system with domain size about tens of microns. Within the interstitial space between ice crystals, the increase of solute concentration suppresses the freezing point, reduces the solute mobility, and hinders the coalescence or further growth of ice crystals. Parameter and sensitivity tests with experimental validation are important for future investigation so that a lab-scale model can be further developed. The theoretical framework for a binary solution can be further extended to multiple components including co-solutes, and for a system under directional cooling or cooling from multiple heat transfer surfaces. Multiscale modeling, precision measurement of the transport properties especially in the supercooled regime, and the quantitative analysis of process dynamics for a larger scale experiment with configuration such as a vial, bag, or bottle that has volume from a few milliliters to tens of liters are all important for industrial applications. Finally the dimensionless groups list in Table 1 can be a reference for the design-of-experiment (DoE) approach for developing operation parameters or new manufacturing strategies to mitigate freeze concentration and enhance the quality of biologic products.

Conclusions
Freeze concentration is an important aspect of the freezing process in manufacturing of biopharmaceutical products. Freeze concentration alters the solution composition in which pharmaceutical proteins have been stabilized, often destabilizing the proteins. The freeze concentration affects ice crystal growth, which can also influence protein stability. We initiated a theoretical analysis on the dynamic evolution of ice crystals, crystal-crystal interactions, and the interplay of thermal, mass, and momentum transport phenomena with the dynamics of freeze concentration. The model is based on an asymmetric binary solution with thermophysical properties in the supercooled regime. A priori specified nucleation sites and a global and rapid cooling are applied to the computational domain. We found that the phase transition and freeze-concentration effect are largely controlled by the pattern of initial seeding and thermal transport. Other competing effects including mass diffusion, density variation, interstitial fluid flows, and the significantly increased local viscosity also play roles at different stages of the freezing process. The phase-field approach provides a great opportunity to model and investigate these process details.