Modeling the morphogenesis of brine channels in sea ice
Brine channels are formed in sea ice under certain constraints and represent a habitat of different microorganisms. The complex system depends on a number of various quantities as salinity, density, pH-value or temperature. Each quantity governs the process of brine channel formation. There exists a strong link between bulk salinity and the presence of brine drainage channels in growing ice with respect to both the horizontal and vertical planes. We develop a suitable phenomenological model for the formation of brine channels both referring to the Ginzburg-Landau-theory of phase transitions as well as to the chemical basis of morphogenesis according to Turing. It is possible to conclude from the critical wavenumber on the size of the structure and the critical parameters. The theoretically deduced transition rates have the same magnitude as the experimental values. The model creates channels of similar size as observed experimentally. An extension of the model towards channels with different sizes is possible. The microstructure of ice determines the albedo feedback and plays therefore an important role for large-scale global circulation models (GCMs).
pacs:92.05.Hj, 82.40.Ck, 89.75.Kd, 47.54.-r
Formation and decay of complex structures depend on changes in entropy. In the long run structures tend to decay since the entropy of universe leads to a maximum and evolves into a ’dead’ steady state Clausius (1867). On the other hand not only living cells avoid the global thermodynamic equilibrium. A. M. TuringTuring (1952) showed in his paper about the chemical basis of morphogenesis which additional conditions are necessary to develop a pattern or structure. For instance, cells can be formed due to an instability of the homogeneous equilibrium which is triggered by random disturbances. In this sense it should be possible, that the habitat of microorganisms in polar areas, the brine channels in sea ice, can be described through a Turing structure.
The internal surface structure of ice changes dramatically when the ice cools below -23 or warms above -5 and has a crucial influence on the species composition and distribution within sea ice Light et al. (2003); Krembs et al. (2000). This observation correlates with the change of the coverage of organisms in brine channels between -2 and -6 Krembs et al. (2000). Golden et al Golden et al. (1998) found a critical brine volume fraction of 5 percent, or a temperature of -5 for salinity of 5 parts per thousand where the ice distinguishes between permeable and impermeable behavior concerning energy and nutrient transport. According to Perovich et al Perovich and Gow (1991) the brine volume increases from 2 to 37 and the correlation length increases from 0.14 to 0.22 mm if the temperature rises from -20 to -1. The permeability varies over more than six orders of magnitude Eicken (2003). Whereas Golden et al Golden et al. (1998) used a percolation model we will demonstrate how the brine channel distribution can be modeled by a reaction-diffusion equation similar to the Ginzburg-Landau treatment of phase transitions. A molecular dynamics simulation shows the change between the hexagonal arranged ice structure and the more disordered liquid water structure Nada et al. (2004).
After a short introduction into the key issue of the structure formation we describe the brine channel structure in sea ice and propose a phenomenological description. For the interpretation of the order parameter we discuss some microscopic properties of water using molecular dynamics simulation in the next chapter II. In chapter III we consider the phase transition and the conditions which allow a structure formation. We verify the model on the basis of measured values in chapter IV and give finally an outlook on further investigations in chapter V.
Ii Microscopic properties of water
ii.1 Formation of brine channels
Various publications report on the life condition for different groups of organisms in the polar areas in brine-filled holes, which arise under certain boundary conditions in sea ice as base- or brine channels (lacuna) Bartsch (1989); Trinks et al. (2005); Weissenberger (1992). They are characterized by the simultaneous existence of different phases, water and ice in a saline environment. Because already marginal temperature variations can disturb this sensitive system, direct measurements of the salinity, temperature, pH-value or ice crystal are morphologically difficult Weissenberger (1992). Weissenberger et al Weissenberger et al. (1992) developed a cast technique in order to examine the channel structure. Freeze-drying eliminates the ice by sublimation, and the hardened casts illustrate the channels as negative pattern. Figure 1 shows a typical granular texture without prevalent orientation.
Sometimes, both columnar and mixed textures occur. Using an imaging system Light et al Light et al. (2003) found brine tubes, brine pockets, bubbles, drained inclusions, transparent areas, and poorly defined inclusions. Air bubbles are much larger than brine pockets. Bubbles possess a mean major axis length of some millimeters and brine pockets are hundred times smaller Perovich and Gow (1996). Cox et al Cox (1983); Cox and Weeks (1983); Cox (1988) presented a quantitative model approach investigating the brine channel volume, salinity profile or heat expansion but without pattern formation. They also described the texture and genetic classification of the sea ice structure experimentally. A crucial factor for the brine channel structure formation is the spatial variability of salinity Cottier et al. (1999).
Different mechanisms are employing the mobility of brine channels which can be used to measure the salinity profile Cottier et al. (1999). Advanced micro-scale photography has been developed to observe in situ the distribution of bottom ice algae Mundy et al. (2007) which allows to determine the variability of the brine channel diameter from bottom to top of the ice. By mesocosm studies the hypothesis was established that the vertical brine stability is the crucial factor for ice algae growth Krembs et al. (2001). Therefore the channel formation during solidification and its dependence on the salinity is of great interest both experimentally and theoretically Worster and Wettlaufer (1997). Experimentally Cottier et al Cottier et al. (1999) presented images, which show the linkages between salinity and brine channel distribution in an ice sample.
To describe different phases in sea ice dependent on temperature and salinity one possible approach is based on the reaction diffusion system
where is the vector of reactants, the three dimensional space vector, the nonlinear reaction kinetics and the matrix of diffusitivities, where is the diffusion coefficient of water and is the diffusion coefficient of salt. For the one- and two-dimensional case we set respectively . The reaction kinetics described by can include the theory of phase transitions by Ginzburg-Landau for the order parameters. Referring to this M. Fabrizio Fabrizio (2008) presented an ice-water phase transition. Under certain conditions, spatial patterns evolve in the so-called Turing space. These patterns can reproduce the distribution ranging from sea water with high salinity to sea ice with low salinity. The brine channel system exists below a critical temperature in a thermodynamical non-equilibrium. It is driven via the desalination of ice during the freezing process that leads to a salinity increase in the brine channels. The higher salt concentration in the remaining liquid phase leads to a freezing point depression and triggers the ocean currents.
ii.2 Different states of water
Already Wilhelm Conrad Röntgen had described the anomalous properties of water with molecules of the first kind, which he called ice molecules and molecules of the second kind Röntgen (1892) and which represent the liquid aggregate state. Dennison Dennison (1921) determined the ordinary hexagonal ice-I modification from X-ray pattern methodically verified by Bragg Bragg (1922). This so called -ice is formed by four oxygen atoms which build a tetrahedron as illustrated in figure 2. In -ice each oxygen atom is tetrahedrically coordinated by four neighbouring oxygen atoms, each accompagnied by a hydrogen bridge. The arrangement is isomorphous to the wurtzite form of zinc sulphide or to the silicon atoms in the tridymite form of silicon dioxide. Bjerrum Bjerrum (1951) and Eisenberg Eisenberg and Kauzmann (2005) have provided a survey about the structure differences between the different polymorphic forms of ice and liquid water.
Molecular dynamics simulations with the TIP3P-model of water using the NAMD-software by the Theoretical and Computational Biophysics Group of the university of Illinois show the change from a regular hexagonal lattice structure to irregular bonds after the melting (figure 2). Nada et al Nada et al. (2004) developed a better six-site potential model of for a crystal growth of ice from water using molecular dynamics and Monte Carlo methods. They computed both the free energy and an order parameter for the description of the water structure. Also Medvedev et al Medvedev and Naberukhin (1987) introduced a ”tetrahedricity measure” for the ordering degree of water. It is possible to discriminate between ice- and water molecules via a two-state function ( if and if ). This tetrahedricity is computed using the sum
where are the lengths of the six edges of the tetrahedron formed by the four nearest neighbors of the considered water molecule. For an ideal tetrahedron one has and the random structure yields . The tetrahedricity can be used in order to define an order parameter according to the Landau - de Gennes model for liquid crystals, which refers to the Clausius-Mosotti-relation. Other simulations such as the percolation models of Stanley et al Stanley (1979); Stanley and Teixeira (1980) use a two-state model, in which a critical correlation length determines the phase transition. A mesoscopic model for the sea ice crystal growth is developed by Kawano and Ohashi Kawano and Ohashi (2009) who used a Voronoi dynamics.
Iii Reaction - diffusion model
iii.1 1+1-dimensional model equations
We consider the reaction diffusion system
in one space dimension. The order parameter according to the Ginzburg-Landau-theory is with and proportional to the tetrahedricity . If the variable is smaller than (), the phase changes from water to ice and vice versa. Thus, changes in reflect temperature variations. The variable is a measure of the salinity. The coefficient depends on the temperature as with the critical temperature . The salt exchange between ice and water is realized by the gain term and the loss term . The positive terms and are the temperature- and salt-concentration-dependent ”driving forces” of the system.
In order to realize the -dependent phase transition, one can expand the order parameter in a power series corresponding to Ginzburg and Landau in equation 3. In order to describe properly a temperature-induced phase transition of second order an expression is necessary. The first-order phase transition is dependent on . Supercooled or superheated phases can coexist, i.e. an hysteresis behavior is possible. Without the term we can also realize a brine channel formation. But this second-order phase transition does not allow us to consider the specific heat as a jump in the order parameter (figures 3 and 4). We write the equations (3) and (4) in dimensionless form (see appendix) by setting , , , and get
with , , and as well as the reaction kinetics
where is the dimensionless order parameter of the water/ice system and the dimensionless salinity. Thus the dynamics only depends on four parameters , , and . Without the salinity in equation (3) respectively (5) the above equation system is reduced to a Ginzburg-Landau equation for the first order phase transition.
iii.1.1 First Order Phase Transitions
When we neglect the salinity , the integration of the kinetic function (7) yields the Landau function for the order parameter of water/ice
as plotted in figure 3. It possesses three minima
When several different minima of equal depth exist, then there is a discontinuity in due to Maxwell construction and one has a first-order phase transition Haug (1997). This is the case if
in consequence of
Thus the critical parameter is determined by the temperature-dependent critical value . The jump in figures 3 and 4 is a measure for the latent heat of the phase transition from water to ice. Feistel and Hagen Feistel and Hagen (1998) have deduced theoretically the latent heat of sea ice for various salinities.
iii.2 Linear Stability Analysis
as outlined in appendix.
There are five fixed points for the kinetics (7) and (8) which satisfy and . In order to get a stable non-oscillating pattern we need a stable spiral point as fixed point. Moreover the associated eigenvalues have to possess a positive real part for a positive wavenumber, i.e. they have to allow to create unstable modes. Not each fixed point satisfies both conditions. Therefore for the following discussion we choose the steady state, and , which corresponds to the observable brine channel structures measured by a casting experiment Mundy et al. (2007); Cottier et al. (1999).
Short-time experiments may also record structures that are formed under non-steady conditions. Since those structures are beyond the scope of the present paper, we proceed with the steady state which leads from (15) to
and which is readily solved
iii.3 Turing space
Let us discuss the equations (16) and (17) concerning the conditions for the occurrance of Turing structures in detail. First we concentrate on the situation where a homogeneous phase is formed. Then the fixed points and are stable according to the eigenvalues (18) if
are negative. Otherwise we would have a globally unstable situation which we rule out. Also homogeneously oscillations don’t describe a brine channel formation. It is easy to see that the solution (LABEL:hom) gives only two negative values if
The trajectories of salinity and the order parameter converge to the steady state value zero by damped oscillations. Therefore the structure formation does not follow from the initial condition in time but from the range of the interaction in space.
Next we discuss the spatial inhomogeneous case, , where some spatial fluctuations may be amplified and form macroscopic structures, i.e. the Turing structure. Therefore we search for such modes of (18) which grow in time, i.e. . Time oscillating structures appear if which can be seen from the solution of (18) to be the case if . In this region we have
Demanding to be positive means . Due to (20) this cannot be fulfilled since the diffusion constant is positive and real. Therefore for a time-growing mode we do not have an imaginary part of in our model. In other words we do not have oscillating and time-growing structures. The restriction for the only allowed region is
In this region we search now for the condition . The term before the square in (18) is negative as can be seen from (20). Therefore we can only have positive if the square of this term is less than the content of the root. This leads to
which restricts the region to the interval
Moreover, the square root must be real, i. e.
which leads to
This has to be in agreement with (25) and discussing the different cases results finally into
The Turing space as phase diagram is determined by condition I and II and is plotted in figure 6.
One can see that the Turing space starts at the minimal (tricritical point)
which means that we have only a Turing space for sufficient large diffusivity . For the Turing space we obtain the possible wave numbers according to condition III as plotted in figure 7.
iii.4 Critical modes
The critical wavenumber can be found from the largest modes. These are given by the minimum of equation (17) from which we find the wavenumbers
and the minima
For and we find again the corresponding inequalities (25) and (26). The formation of a spatial Turing structure, a non-oscillating pattern, requires a negative for . In this case there is a range of wavenumbers which are linearly unstable as seen in figure 5. In figure 8 we illustrate the behaviour of for different diffusion constants. Only those which lead to negative are forming the Turing structure as discussed in the previous chapter. This critical range can be obtained, if the diffusion coefficient is greater than the critical diffusion coefficient of condition II, (LABEL:cond)
which we get from with the critical wavenumber
The size of the structure can be estimated from . The pattern size depends on the both parameters and . The parameters determine the brine channel size and vice versa. With the parameters chosen in fig. 5 we obtain a pattern size of . In the next chapter we compare this value with experimental quantities. With a small initial random perturbation we plot snapshots of the time evolution of the order parameter and the salinity in figure 9. The quantities and are opposite to each other; domains with low salinity correspond to domains with ice and domains with high salinity to water domains. We see the formation of a mean mode given by the wave length .
A positive , respectivly a negative , for guarantees that and converge to the stable fixed point . Therefore the structure formation does not follow from the initial oscillations in time (fig. 10). In order to obtain a new spatial structure there must exist at least a negative for , respectively a positive eigenvalue , for as was discussed in the last chapter.
iii.5 2+1-dimensional model
From the characteristic equation in the spatially two-dimensional case
we find the corresponding dispersion relation
which is illustrated in figure 11.
The Turing space is bounded by the sectional plane . The evolution of the order parameter and the salinity is illustrated in figures 12 and 13. Their behavior is inversely proportional and corresponds to the fact, that a high salinity occurs in the water phase and a low salinity in the ice phase. Similiar as in the one-dimensional case we see the dominant formaton of one wavelength. The model kinetics generate brine channels of similar size. In order to obtain a hierarchical net of brine channels of different size, the kinetics in the basic equations can be altered accordingly Murray (1993).
iii.6 Note on second and thirth order kinetics
If we replace the kinetics (7) by
we can carry out the same linear stability analysis for the fixed point and . Then we obtain the same characteristic equation as (16) and (17). Therfore we get the same Turing space for the structure formation. Consequently, both kinetics allow us to realize a brine channel formation but the thirth order kinetics describes a second order phase transition only. In this connection it is possible to discuss second order phase transitions with spin models, too.
Iv Connection to experimental data
The critical domain size is determined by the equation (34). Due to this relation we can infer other parameters in the model equations (3) and (4). From the relation between the dimensionless wavenumber and the dimensional wavenumber
The observed diameters of the brine channels range from to scale Weissenberger (1992). For a size of and a diffusion coefficient for -molecules we obtain the product and a transition rate . The rate is proportional to reorientations of the molecules per second, (Eisenberg) Eisenberg and Kauzmann (2005) and to the scaled temperature
where is the melting point depending on the salinity. The mean salinity in sea ice of corresponds to 1 -molecule per 100 -molecules, i.e. 1 -ion and 1 -ion per 100 -molecules in a diluted solution after the dissociation or a ratio of . From this facts we obtain according to Clausius-Clapeyron
a freezing point depression from to , where is the latent heat of the phase transition from water to ice, is the universal gas constant and . Thus we obtain correctly . For an environmental temperature of according to (41), a transition rate of follows which nearly corresponds to , which we estimated from the domain size (39).
Furthermore we find the transition rate and the diffusion coefficient . Due to the transformation of equation (3) into the dimensionless form (5) there exists a fixed relation (11) because of . We obtain for the equation (3) the rate with , of
which is proportional to . A transition rate yields a critical rate . From the knowledge of the diffusion coefficient and the size of brine channels we can deduce the two rates and . Both rates possess the same order of magnitude and are inside the Touring space of structure formation. If the experiments would lead to other parameters and , especially to rather different values, a brine channel could not arise because of the limitation of the Turing space in figure 6 and 7. In other words the model here seems to describe the experimental finding of brine channel formation.
Due to the small difference between the time constants and we obtain a dynamic interference between the reorientation of the water molecules and the desalinization. Both are evolving on nearly the same time scale. In particular, we cannot simplify the kinetics by separating time scales using the Tichonov theorem Tichonov (1952) in order to reduce our reaction diffusion system (3),(4) but have to consider both dynamics as demonstrated.
In this paper it has been shown that a reaction diffusion system which connects the basic ideas both of Ginzburg and of Turing can describe the formation of brine channels with realistic parameters. For the chosen parameters patterns of similar size emerged. Eicken Eicken (1991) and Weissenberger Weissenberger (1992) distinguished between six various texture classes of sea ice dependent on the crystal morphology, brine inclusions and the genesis. The different ice crystal growth depends on snow depostion, flooding, turbulent mixing, quiescent growth rate or supercooling. Each condition determines the character of the kinetics. Non-linear heat and salt dissipation for example leads to dendritic growth (snowflakes) whereas one observes in sea ice mostly lamellar or cellular structures rather than complete dendrits Eicken (2003). Hence, the morphology of sea ice is one criterion for the choice of an appropriate kinetics for the genesis of sea ice. Therefore, in order to simulate different structure sizes and textures, we can modify the dispersion relation by varying the parameters , and or by a modified kinetics Turk (1991). The crucial point is the shape of the dispersion function. If there are multiple different positive unstable regions for the wavenumbers with positive real part of eigenvalues we could expect that differently large channels evolve. For instance Worster et al Worster and Wettlaufer (1997) has presented a general theory for convection with mushy layers. The two different minima of the neutral curve, determined by the linear stability analysis, correspond to two different modes of convection, which affect the kinetics and determines the size distribution of the brine channels. We note that the initial conditions are decisive for the appearance of specific pattern Murray (1993). Hence one should investigate how dislocations or antifreeze proteins influence the formation of the brine channel distribution.
Appendix A Dimensionless Quantities
Appendix B Linear Stability Analysis
If we consider only linear terms
where is the Jacobian
Using the Fourier ansatz and in (15) we find
With and the eigenvalue equation
follows with the characteristic equation
Acknowledgements.The authors are deeply indebted Jürgen Weissenberger for the many exposures of casts brine channels. Dr. L. Dünkel and Dr. P. Augustin are thanked for their system theoretic discussion and Dr. U. Menzel at the Rudbeck laboratory, Uppsala, Sweden for his support. This work was supported by DFG Priority Program 1157 via GE1202/06, 444BRA-113/57/0-1, the BMBF, the DAAD and by European ESF program NES. The financial support by the Brazilian Ministry of Science and Technology is acknowledged.
- Clausius (1867) R. Clausius, Über den zweiten Hauptsatz der mechanischen Wärmetheorie (Friedrich Vieweg und Sohn, Braunschweig, 1867).
- Turing (1952) A. M. Turing, Phil. Trans. of the R. Soc. B: Biological Sciences 237, 37 (1952).
- Light et al. (2003) B. Light, G. A. Maykut, and T. C. Grenfell, J. Geophys. Res. 108, 3051 (2003).
- Krembs et al. (2000) C. Krembs, G. R, and M. Spindler, J. Exp. Mar. Biol. Ecol. 243, 55 (2000).
- Golden et al. (1998) K. M. Golden, S. F. Ackley, and V.I.Lytle, Science 282, 2238 (1998).
- Perovich and Gow (1991) D. K. Perovich and A. J. Gow, J. Geophys. Res. 96, 16943 (1991).
- Eicken (2003) H. Eicken, in Sea Ice: an introduction to its physics, chemistry, biology and geology, edited by D. N. Thomas and G. S. Dieckmann (Blackwell Science Ltd, Oxford, 2003), chap. 2.
- Nada et al. (2004) H. Nada, J. P. van der Eerden, and Y. Furukawa, J. Cryst. Growth 266, 297 (2004).
- Weissenberger (1992) J. Weissenberger, Ber. Polarforschung 111, 1 (1992).
- Bartsch (1989) A. Bartsch, Ber. Polarforschung 63, 1 (1989).
- Trinks et al. (2005) H. Trinks, W. Schröder, and C. K. Biebricher, Origins Life Evol. Biosphere 35, 429 (2005).
- Weissenberger et al. (1992) J. Weissenberger, G. G. R. Dieckmann, and M. Spindler, Limnol. Oceanogr. 37, 179 (1992).
- Perovich and Gow (1996) D. K. Perovich and A. J. Gow, J. Geophys. Res. 101, 16943 (1996).
- Cox (1983) G. F. N. Cox, J. Glaciol. 29, 425 (1983).
- Cox and Weeks (1983) G. F. N. Cox and W. F. Weeks, J. Glaciol. 29, 306 (1983).
- Cox (1988) G. F. N. Cox, J. Geophys. Res. 93, 449 (1988).
- Cottier et al. (1999) F. Cottier, H. Eicken, and P. Wadhams, J. Geophys. Res. 104, 15859 (1999).
- Mundy et al. (2007) C. J. Mundy, D. G. Barber, C. Michel, and R. F. Marsden, Polar. Biol. 30, 1099 (2007).
- Krembs et al. (2001) C. Krembs, T. Mock, and R. Gradinger, Polar. Biol. 24, 356 (2001).
- Worster and Wettlaufer (1997) M. G. Worster and J. S. Wettlaufer, J. Phys. Chem. B 101, 6132 (1997).
- Fabrizio (2008) M. Fabrizio, J. Math. Phys. 49, 102902 (2008).
- Röntgen (1892) W. C. Röntgen, Ann. Phys. (Leipzig) 45, 91 (1892).
- Dennison (1921) D. M. Dennison, Phys. Rev. 17, 20 (1921).
- Bragg (1922) W. H. Bragg, Proc. Roy. Soc. London 34, 98 (1922).
- Bjerrum (1951) N. Bjerrum, Det Kongelige Danske Videnskabernes Selskab, Matematisk-fysiske Meddelelser 27, 1 (1951).
- Eisenberg and Kauzmann (2005) D. Eisenberg and W. Kauzmann, The structure and properties of water (Clarendon Press, Oxford, 2005).
- Medvedev and Naberukhin (1987) N. N. Medvedev and Y. I. Naberukhin, J. Non-Cryst. Solids 94, 402 (1987).
- Stanley (1979) H. E. Stanley, L. Phys. A: Math. Gen. 12, L329 (1979).
- Stanley and Teixeira (1980) H. E. Stanley and J. Teixeira, J. Chem. Phys. 73, 3404 (1980).
- Kawano and Ohashi (2009) Y. Kawano and T. Ohashi, Cold Reg. Sci. Technol. (Accepted date: 2 February 2009) (2009).
- Haug (1997) H. Haug, Statistische Physik (Friedr. Vieweg & Sohn Verlagsgesellschaft mbH, Braunschweig/Wiesbaden, 1997).
- Feistel and Hagen (1998) R. Feistel and E. Hagen, Cold Reg. Sci. Technol. 28, 83 (1998).
- Murray (1993) J. D. Murray, Mathematical Biology (Springer Verlag, Berlin, 1993).
- Tichonov (1952) A. Tichonov, Mat. Sbornik 31, 575 (1952).
- Eicken (1991) H. Eicken, Ber. Polarforschung 82, 1 (1991).
- Turk (1991) G. Turk, Computer Graphics 25, 289 (1991).