Three-dimensional Lattice Model for the Water/Ice System BY G. M. BELL* Mathematics Department, Chelsea College (University of London), London SW3 AND D. W. SALT Mathematics Department, Portsmouth Polytechnic, Hampshire Terrace, Portsmouth Received 12th May, 1975 The sites of a body-centred cubic lattice are either vacant or occupied by molecules with four hydrogen-bonded “ arms ”, two positive and two negative, directed tetrahedrally. Appropriately oriented nearest-neighbour molecules can form bonds by linking bonding “arms ” of opposite polarity. There are also interaction energies for unbonded nearest neighbours and second-nearest neighbours. As well as disordered states, two long range ordered structures are possible, one an open bonded network on a diamond sub-lattice, like ice I(c), and the other a close-packed structure with intertwined bonded networks, like ice VII.A free energy, depending on three order variables, is derived using a zeroth-order statistical approximation modified to allow for asymmetrical bonding. For suitable energy parameter values, there exist four phases, vapour, liquid, open ice and close- packed ice. There is a critical point, an open ice/liquid/vapour triple point and a close-packed ice/open ice/vapour triple point at a lower temperature. Between the triple points increase of pres-sure causes the open ice to transform into a liquid of higher density than the ice while below the lower critical point it causes it to transform to close-packed ice. 1. INTRODUCTION The structure of both fluid water and ice at low pressures depends on the water molecule’s capacity for forming tetrahedrally directed hydrogen-bonds [for more recent discussions see ref.(2) and (3)]. The resulting open four-coordinated network becomes less stable than closer-packed forms as the pressure increases. The present paper follows a number of studies of lattice fluids which show water-like properties resulting from the competition between regions of open structure and regions of high density and energy. These have been ~ne-dimensional,~* two-dimensional 6*’* and three-dimensional,8* the three-dimensional work being most relevant here. Several more recent lattice fluid models of water have, like that of been based on the body-centred cubic (b.c.c.) lattice.Lavis l3 has introduced long-range order into the two-dimensional model of Bell and Lavk7 The aim of the present paper is to introduce long-range order into a three-dimensional model similar to Bell’s and hence consider ice phases as well as fluid phases. The molecule is regarded as having four “ bonding arms ”, two of positive and two of negative polarity, and a bond is formed between nearest neighbours when a positive arm links with a negative arm. The lattice is of the b.c.c. type on which two forms of ice structure are possible. These are ice I(c), an open structure with half the b.c.c. lattice sites vacant, and ice VII, a close-packed structure consisting of two intertwining ice I(c) networks. The use of the b.c.c.lattice thus implies the approximations of replacing hexagonal ice I(h) with cubic ice I(c), which is close to it in local structure and free energyY3 and of replacing all the various close-packed 76 G. M. BELL AND D. W. SALT forms of ice by ice VII. A limitation is thus, in effect, imposed on the configuration space over which statistical averages are taken. However, the results of both the previous work on the fluid phases of water and that presented here indicate that enough configuration space is retained to represent some of the most important characteristic phenomena displayed by water. Since the possibility of long range ordering makes the analysis much more com- plicated than for the fluid phases alone,* the statistical treatment is simplified by using the zeroth-order (molecular field) approximation instead of the first-order approxima- tion.The zeroth-order method is adjusted to allow for the formation of bonds between arms of opposite polarity only and it is found that the entropy of a completely bonded network is then equal to that given by Pauling’s approximation. The zeroth- order approximation smoothes out the anomalous features of the liquid phase but does enable us to consider long-range ordered (ice) phases. Another modification here of the fluid phase theory concerns the “ close-packing ” energy parameter which is necessary to ensure that open ice is more stable at low pressures than the intertwined form. Bell imposed this close-packing “ penalty ” by assigning a posi-tive energy to the smallest fully occupied triads of sites on the b.c.c.lattice. Here better results were obtained by assigning a positive (i.e., repulsive) energy to second- neighbour pairs on the b.c.c. lattice. There may be some connection between this and the change of statistical approximation. After minimising the free energy with respect to order parameters we plot the equilibrium free energy against number density for various temperatures. For suit- able sets of parameter values we then derive phase diagrams involving the four phases, vapour, liquid, open ice and close-packed ice. The well-known experimental property of the low pressure form of ice melting into a denser liquid is reproduced theoretically. 2. THE MODEL The b.c.c.lattice can be divided into four facecentred cubic (f.c.c.) sub-lattices labelled 1, l’, 2 and 2’ respectively (see fig. 1). A pair of diamond sub-lattices can be constructed by grouping 1with 1’ and 2 with 2’ or alternatively by grouping 1 with 2’ and 2 with 1’. When preferential distribution on diamond sub-lattices is considered we shall choose, for definiteness, the first of these groupings. The tetrahedrally directed bonding arms of a molecule on any site are regarded as pointing towards four of the eight nearest neighbour sites on the b.c.c. lattice. Since two of its arms are positive and two negative a molecule on, say, a site of sub-lattice 1 has twelve distinct orientations in six of which its arms point towards 1’-sites (as shown in fig.1) and in the other six to 2’4tes. We denote the fractions of occupied sites in the four f.c.c. sub-lattices by pl, PI., p2 and p2. respectively and the overall number density p is then given by P = $(Pl+Pl,+P,+P,f). We denote the proportion of 1-sites occupied by molecules oriented towards 1’-sites by p11. and the proportion occupied by molecules oriented towards 2’-sites by p12.. Similarly defining p21p, p22J, plt1,plg2,p2t1and p2f2we have P1 = Pll’+P12‘1 Pz = P21’+P22’, P1* = P19+P1‘2’ P2‘ = P2’1+P2‘2. The basic parameters will now be introduced, starting with the volume uo per b.c.c. lattice site. This is regarded as a constant determined by the distance of closest approach of two molecular centres, which in our model is the distance between LATTICE MODEL FOR THE WATER/ICE SYSTEM nearest-neighbour lattice sites.An unbonded nearest-neighbour pair of molecules has energy -Eand a bonded pair -(E+w) (E > 0,w > 0), w being the bonding energy. Experimentally, forms of ice, like ice VII, with intertwined bonded networks only occur at high pressures and hence must involve high energies. In our model this requirement is met by assigning a positive energy u2 to each second-neighbour pair of molecules on the b.c.c. lattice. We now consider possible configurations on the lattice at zero absolute temperature (T= 0). .* FIG.1.-Sites on the body-centred cubic lattice :a molecule is placed on site 1 with bonding arms directed towards 1’-sites. OPEN ICE In perfect ice I(c) all sites of one diamond sub-lattice (taken as 11’) are occupied by a completely bonded network while all sites of the remaining diamond sub-lattice (22’) are empty.Hence, p1 = p11. = p1. = p1.1 = 1, p2 = p2’ = 0. (2.3) Eqn (2.3) implies that all molecules on I-sites are oriented towards 1’-sites and vice-versa. However, for complete bonding it is necessary, in addition, that each positive arm be directed towards a negative arm. The latter condition still allows a large number of orientational states, giving rise to the well-known zero-point en- tropy.14*l5 The value of the latter given by the approximate statistical theory of the present paper is discussed in section 3. There are no occupied second-neighbour pairs of sites (e.g., 12 or 1’2‘in fig.l), each molecule is bonded to four others and the volume per molecule is 2v0. Hence the configurational enthalpy of an assembly of A4 mole-cules is given by H,= ~M(~v~-E-w) (2.4) The perfect crystal to which eqn (2.4) applies exists only at T= 0 and in section 4 we introduce ordering parameters which enable us to consider ices (long-range ordered states) at T> 0. CLOSE-PACKED ICE In perfect ice VII there are two diamond sub-lattices occupied by completely bonded networks. Hence, p1 = PI11 = p1t = p101 = 1, p2 = p22. = p2. = p2.2 = 1. (2.5) The zero-point entropy per molecule is the same as in perfect open ice. As in the latter, each molecule is bonded to four others but it now has four unbonded G.M. BELL AND D. W. SALT 79 nearest-neighbour and six second-nearest neighbour molecules on the other diamond sub-lattice (e.g., in fig. 1, the molecule on site 1 is bonded to nearest neighbours on 1’-sites, has unbonded nearest neighbours on 2’-sites and second nearest neighbours on sites such as 2). Since the volume per molecule is now u0, the enthalpy of a close- packed assembly of M molecules is given by H, = M(Pv~-2~ -4~+3U2). (2.6) Like eqn (2.4), eqn (2.6) applies only at T = 0. SEPARATION PRESSURE Since enthalpy and Gibbs free energy are equal at T = 0, open ice is more stable than the close-packed form if the expression (2.4) is less than the expression (2.6). Now, they are equal when p = po,po being given by pouo = 3U,-2&. (2.7) Hence, if 3u2 > 28 there is a range of pressures 0 < p <po where open ice is more stable than close-packed ice at T = 0.The close-packed form is more stable at T = 0 for p > po. In our model it would also be possible, with an extra energy parameter, to consider ice VIII. The latter has a structure like VII but with the zero-point entropy eliminated by long-range ordering of the proton positions.2* However, we do not introduce this type of ordering, both for simplicity and because it is not likely to be important in the region of the triple points involving ice I. In the present model there is a possible ordered interstitial state, with, say, diamond sub-lattice 11’ occupied by a bonded network and f.c.c. sub-lattices 2 and 2’ respectively empty and occupied by unbonded molecules.It can be shown that there is a pressure range at T = 0 in which this form is stable if --E > w. This is undesirable since all forms of ice known experimentally are fully bonded, apart from imperfections. In fact, we take E =-0 in all calculations so that an ordered phase of this type is unlikely to occur. 3. FREE ENERGY Since a molecule on, say, f.c.c. sub-lattice 1 has twelve distinct orientational states, six with its bonding arms directed towards 1’-sites and six with them towards 2’-sites, the zeroth-order entropy of distribution on the lattice is -+p12.In P12’-+(l-pl) In (l-p,)+ 6 N being the total number of sites on the b.c.c. lattice. It is now convenient to introduce the term “ oriented pair ” which denotes a pair of nearest-neighbour mole- cules each with a bonding arm directed towards the other.Now AS is not the entire entropy since there is also a contribution due to the fact that the arms which an “oriented pair ” direct at each other are more likely than not to be of opposite polarity. Each 1-site, say, has four nearest-neighbour 1‘-sites and in the zeroth order approximation there are aN(4plplf)= Nplplf nearest-neighbour 1-1’ pairs of mole- LATTICE MODEL FOR THE WATER/ICE SYSTEM cules, of which Npl1~p,~,are oriented pairs. Denoting the total number of nearest- neighbour pairs by NAg and the total number of oriented pairs by No,, Nk2 = N(P1P 1 + P 1P2 + P2P1 + P2 P2J (3.2)and Nor = N(p11'PI'1 + ~12'~2'1+ + ~22'~2'2)* (3.3)~21'~1'2 Since each site has six second-neighbour sites, the number NA: of second-neighbour pairs of molecules is NIn",!= W(PlP2 + PlJP24. (3.4) The configurational energy contains terms --N~~E+N~~U~,due to the first and second neighbour energies, as well as a bonding energy contribution. To evaluate the latter, note that each oriented pair has 36 orientational states in half of which the bond- ing arms have opposite polarity.Since there is a bond energy -w we assign a free energy 18+ 18 exp(w/kT) 1+ exp(w/kT)-+(T) = -Win 36 = -kT In 2 (3 5) to each oriented pair. It is not difficult to show that, when kT/w < 1, $(T) = w -kT In 2 + 0{kTexp(-w/kT)) (3.6) while when kT/w % 1, 4(T)= +w+O(w2/kT).(3.7) We may now write, for the configurational Helmholtz free energy F,, Fc = -N(I)&+p)Umm mm z-Nor4(T)-TAS (3.8) where the terms on the right-hand side are given by equations (3.1-5). It is useful to define a free energyf, per site by fc(p, T)= FCm (3.9) Then the pressure p and the configurational part p, of the chemical potential are given by luc = (?fm)T, PUO = PPc-fc. (3.10) In the fluid state there is no long-range order and hence p11. = p12' = p21. = p22. = p1.1 = p1.2 = p2'1 = p292 = +p. (3.11) Hence, substituting in (3.1-4)and (3.8), fc(P, T) = -{d)(T)+4~-3~~)p~+kT Y2p In -+(l-p) ln(1-p)} (3.12){and then, from (3.10) pvo = -(~(T)+4~-33~~)p~-kTIn (3.13)(1-p). The critical-point conditions are apiap = 0, a2PpP = 0, (3.14) which, from (3.13) are satisfied by the critical values pc,pc and T,, given by pc = +, pcc0 = kT,(ln 2 -+), kT, = +(4(T,)+ 4~-3u2).(3.15) The first two relations of eqn (3.15) are the same as in the standard zeroth-order lattice fluid model but the last relation, which is an implicit equation for T,, is char- acteristic for the present model. From (3.7) the entropy contribution from +(T) G. M. BELL AND D. W. SALT disappears when kT/w $ 1 and hence, for high T, the configurational entropy S, is given by S, = AS = kN p In -+(l-p) In (1-p) .{pz I This is correct since it corresponds to a random distribution of orientations and justifies the factor 2 in the last expression for 4(T)of eqn (3.5). In the ice structure at very low temperatures it can be seen from eqn (3.6) that there is an entropy contribution of -k In 2 from &T) for each of the 3M oriented pairs.Hence, from eqn (2.3) and (3.1) the configurational entropy is S, = AS-2Mk In 2 = Mk(1n 6-2 In 2) = Mk In (3/2). This is just the Pauling approximation l4 which is, in fact, a reinarkably good estimate. 4. ORDER VARIABLES AND EQUILIBRIUM RELATIONS Choosing 11’ and 22’ as diamond sub-lattices of the b.c.c. we define order para- meters a, m, and m2 by the following relations, which satisfy eqn (2.2) : P1 = P(’+O) = P1’7 p2 = P(l-= p2.7 (4.1) pllt = +p(l+o)(l+m,) = plP17 plZf = +p(l+a)(l-nzl) = plI2, (4.3) p221= &1(1--0)(l+n2,) = p2p2, p2,’ = 3p(1-a)(l-nz2) = p2r1. (4.3) Both the open and close-packed forms of long-range order can be treated with the aid of these parameters.Non-zero values of a imply preferential concentration of molecules on one sub-lattice and positive values of m, and m2 that oriented pairs are more likely when both molecules are on the same sub-lattice. Substituting into eqn (3.1-3.4), (3.8) and (3.9), the free energy per site becomes fc = -p2[4& -h2(I -a2)+$$(T)(4+ (in, +rn2)2+20(m1-nz2)(nz1+inl +2) + ~~(m;-2mlm2 +4m1 +4m2)>]++in; P3kT[2p In -+p((l +a) In (1 +a)+(1-a) In (1-0))+ l2 (1 -p-pa) In (I -p-po)+(l -p+pa) In (1 -p+po)+ +p(l +a){(l +ml)In (1 +m,)+(l -ml)In (1 -inl)> + +p(l -){(I +m,) In (1 +m2)+(l-mJ In (1 -m,)> . (4.4)1 The free energy has a stationary value with respect to the order variables if the equilibrium relations are satisfied. After some manipulation these relations can be expressed in the form 6u2p0 = 0, (4.6) kT Ifin, -2 In 1-m, --p{!??1(1 +a)+m,(l -o)+2o)$(T) = 0, (4.7) -1n--kT l+m, p(1111(1 +0)+1722( 1-0) -2O)#( T) = 0.2 3 -11112 a2fIaa2 ayC/aaam, a2y~aQam, A = a2f,/aaam, a2f/am; a2f/dmidm2 > 0. (4.14) a2f,/aadm2 a2f,lam,am, a2ffclam; A = AF = (kT-2p4(T)){ -$u,kTp+-(W2 -P~~~(T)}.(4.15)4(1 -P) The boundary of stability between states F and C is given by equating the first factor in AF to zero, which yields eqn (4.12). The boundary of stability between states F and 0 is given by equating the second factor in AF to zero. At high values of p the first factor reaches zero before the second as the temperature is decreased, but for lower values of p the situation is reversed.It should be understood that even though, at a given point (p,T)a certain type of solution gives the lowest free energy it does not follow that the corresponding state will exist there since it may well be possible to further reduce the free energy by phase separation. 5. RESULTS AND DISCUSSION To perform numerical work values of the energy parameter, ratios E/W and u2/w must first be selected. The ranges of suitable parameter values are fully discussed by G. M. BELL AND D. W. SALT Salt? The main criteria are that the separation pressurep,, given by eqn (2.7),shall be positive and that the liquid/vapour critical temperature T,, given by eqn (3.15), shall be greater than the highest temperature at which equilibrium relation solutions of type 0 can exist.0:I 0,2 0:3 0.4 0:5 0,6 0,7 0,.8 0.9 1.00-I I I I I 13 I I I !z-0.1 -P FIG.2.-Curves of reduced Helmholtz free energy per site,f,/w, against number density p for u2/w = 1.0, E/W = 1.45. pm denotes the boundary between close-packed ice and the liquid phase. The broken lines are double tangents. (a)kT/w = 1.1, (b) kT/w = 1.15, (c) kT/w = 11275,-(d)kT/w = 1.45. To investigate the system at a given value of T,it is first necessary to find the ranges of p in which solutions of types F, C and 0 respectively give the lowest free energy LATTICE MODEL FOR THE WATERIICE SYSTEM per site fc. The next step is to plot this value of fc against p.Finally the range or ranges, if any, of p in which phase separation occurs are found by drawing double tangents to the (&I) curves. The full phase diagram results from repeating this procedure at different temperatures. Fig. 2 shows the type of (f,,p) curve occurring at various temperatures if a suitable choice of parameters is made. In fig. 2(a) the left-hand trough of the (f,,p) curve occurs in a low density range where only a disordered state exists and thus corresponds to the vapour phase. The middle trough occurs in a range round p = 4where solu- tions of type 0 give the lowest free energy and thus corresponds to open ice. Hence the left-hand double tangent represents a vapour/ice I phase equilibrium. On the right the free energyf, in the range (pm,l) is for a solution of type C, but for p < pm it is for a disordered state F.Hence the right-hand double tangent represents an 1.8-(0) . --.-.' * '-VAPO U R 1.4-i 3. /// ---___5 1.2-ii li 1.0-// CLOSE -PACKEDOPEN ICE ICE 0.8' I I l > I 5, 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.9 0.9 1.3 '-I1.0I I \ OPEN ICEL11 ,I0.80 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.e 0.9 1.0 e FIG.3.-Density against temperature phase diagram for (a)uz/w= 1.0, E/W = 1.45 ; (b) u2/w= 0.633, E/W = 0.9. The horizontal lines connect coexistent phases. open ice-close packed ice (ice VII) equilibrium. Fig. 2(b) corresponds to a higher temperature and the important change from 2(a)is that the right-hand double tangent now touches the high-density trough where p < pm.Hence there is equilibrium between open ice and a denser liquid phase. The ice I/ice VII/liquid triple point has thus occurred in the temperature range between fig. 2(a)and 2(b). Fig. 2(c) is drawn for a higher temperature at which the two double tangents coincide ;this is the ice I/ liquid/vapour triple point. In a temperature range above the triple point solutions of type 0 still give lower values offc than F or C in a range around p = 3. However fig. 2(d)shows that the corresponding trough in the (f,,p) curve lies above the double- tangent drawn between the liquid and vapour troughs. Hence phase-separation into liquid and vapour gives a lower free energy than formation of ice I.Fig. 3(a) is a density against temperature phase diagram for the parameter values used in the free energy curves of fig. 2 while fig. 4(a) is the corresponding pressure against temperature phase diagram. Fig. 3(b)and 4(b)are respectively density against G. M. BELL AND D. W. SALT 85 temperature and pressure against temperature phase diagrams for a different set of parameter values. The phase diagrams show a liquid/vapour critical point, a liquid/ vapour/open ice triple point and a liquid/open ice/close-packed ice triple point. Between the two triple points there is, at each temperature, an open ice phase in equilibrium with a denser liquid phase. The model thus reproduces the phenomenon of “ice floating on water ”. At temperatures lower than the lower triple-point tem- perature, increasing the pressure sends the open ice phase over into a close-packed ice phase.It is of interest to note that, in contrast to the “ second-neighbour exclu- sion ” quadratic lattice melting model discussed by Runnels l7 the existence of a zero-point entropy does not prevent liquid/solid transitions. The difference is pre-sumably that in the present model the structure of the solid “ices” is as much determined by the attractive directioiial bonding forces as by the “ hard-core ” repulsion. I 0.4 0.4 II CLOSE -’I 0.35 0.35 PACKED’ CLOSE, PACKED 0.3 ICE 0.3 I il 0.25 0.25 LIQUID 3 I -1 0 f$ 0.2 @. 2 ”./c.-j/ UR ”OPEN, VAPOUR 0.15 0.15 ICE/ i/ 0. OPEN-0./ /I I / 0.05 /” ,3.05 / I , .1---0.8 1.0 1.2 1.4 1.6 1.8 0.9 1.0 1.1 1.2 1.3 1.4 (a> kTlw (b) FIG.4.-Pressure against temperature phase diagram for (a) u2/w = 1.0, EIW= 1.45; (b) UZ/W = 0.633, E/W = 0.9. As well as the satisfactory features just mentioned there are some unsatisfactory features in the results. The transition between close-packed ice and liquid is second- order, which results from these two states not giving separate ‘‘troughs ” on the (f,,p) curve as seen in fig. 2 above. Second-order transitions are quite common in lattice melting models when the ‘‘hard core ” of the molecule occupies only a small number of sites,l7-I9 the small number being unity in the present case. Although the open icelliquid coexistence pressure decreases with temperature near the lower triple point the curve turns upwards near the higher triple point. Where this coexistence LATTICE MODEL FOR THE WATER~ICE SYSTEM curve has a positive slope the dense liquid must have a lower molar entropy than the open ice.The positive slope is very small in the case shown in fig. 4(a), but we have not been able to find parameter values where it completely disappears. It is hard to tell whether this results from deficiencies in the lattice model or the zeroth-order statistical approximation and we are currently attempting to apply a first-order approximation to the problem. However, in spite of these difficulties this simple model agrees well enough with experiment to confirm that the water phase diagram depends on the possibility, due to hydrogen bonding, of an open four-coordinated structure occurring.The non-bonding energy parameters E and 2i2 which have been used rep- resent, in a schematic way, a wide range of real interactions due to electrostatic, induction and dispersion forces as well as steric effects for unbonded molecules.2 The fact that the ranges of these parameters must be carefully chosen to give reasonable results indicates that many characteristic phenomena in water depend on a rather delicate balance between different types of interaction energy as well as on the essential structural factors. D. W. S. thanks Chelsea College for a research studentship. J. D. Bernal and R. H. Fowler, J. Chem. Phys., 1933, 1, 515.D. Eisenberg and W. Kauzmann, The Structures and Properties of Water (Clarendon, Oxford. 1969). N. H. Fletcher, The Chemical Physics of Ice (Cambridge University Press, London, 1970). G. M, Bell, J. Math. Phys., 1969, 10, 1753. G. M. Bell and D. W. Salt, MoZ. Phys., 1973, 26, 387. G. M. Bell and D. A. Lavis, J. Phys. A, 1970,3,427. G. M. Bell and D. A. Lavis, J. Phys. A, 1970,3,568. G. M. Bell, J. Phys. C, 1972, 5, 389. G. M. Bell and H. Sallouta, MoZ. Phys., 1975, 29, 1621. lo 0. Weres and S. A. Rice, J. Amer. Chem. SOC., 1972, 94, 8983. l1 D. E. O’Reilly, Phys. Rev. A, 1973, 7,1659. l2 P. D, Fleming and J. N. Gibbs, J. Stat. Phys., 1974, 10, 157, 351. l3 D. A. Lavis, J. Phys. C, 1973, 6, 1530. l4 L. Pauling, J. Amer. Chem. SOC.,1935, 57, 2680. l5 E. H. Lieb and F. Y.Wu, in Phase Transitions and Critical Phenomena, ed. C.Domb and M. S. Green (Academic Press, London, 1972), vol. 1. I6 D. W. Salt, Ph.D. Thesis (London, 1974). l7 L. K. Runnels, in Phase Transitions and Critical Phenomena, ed. C. Domb and M. S. Green (Academic Press, London, 1972), vol. 2. D. M. Burley, in Phase Transitions and Critical Phenomena, ed. C. Domb and M. S. Green (Academic Press, London, 1972), vol. 2. l9 R. D. Kaye and D. M. Burley, J. Phys. A, 1974, 7,843. (PAPER 5/889)