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.