1. Introduction

When water freezes in confined spaces, it can generate large stresses, often resulting in material damage. This is important across fields ranging from glaciology to geomorphology, food science, civil engineering and cryopreservation (Walder and Hallet, Reference Walder and Hallet1985; Karlsson and Toner, Reference Karlsson and Toner1996; Dash and others, Reference Dash, Rempel and Wettlaufer2006; Petzold and Aguilera, Reference Petzold and Aguilera2009; Rempel, Reference Rempel2010; Vlahou and Worster, Reference Vlahou and Worster2015; Jha and others, Reference Jha, Xanthakis, Chevallier, Jury and Le-Bail2019). Broadly speaking, ice can generate stresses via two different mechanisms (Wettlaufer and Worster, Reference Wettlaufer and Worster2006; Peppin and Style, Reference Peppin and Style2013). The first is due to the expansion of water as it freezes: in a closed cavity, freezing will generate pressure (Fig. 1a). The second is unrelated to the expansion of water and often dominates in porous materials – for example, in the process of frost heave (Peppin and Style, Reference Peppin and Style2013). Here, ice can form in open pores of a wet material (Fig. 1b), but no pressure builds up during the initial ice-formation process (any pressure is relieved by water flow away from the growing ice). However, after the initial ice formation, unfrozen water is sucked back toward the ice crystals. When this water freezes onto the existing ice, it can cause the ice to expand its confining pore. This cryosuction process is aided by the presence of thin, mobile layers of water at the surface of ice (known as premelted films) (Slater and Michaelides, Reference Slater and Michaelides2019). These allow growth of the ice, not just at the pore throat, but also along the pore/ice interface. In both cases, ice will continue to grow, building up pressure, until the pressure reaches a maximum value given by a temperature-dependent stall pressure, P _st (Peppin and Style, Reference Peppin and Style2013; Gerber and others, Reference Gerber, Wilen, Poydenot, Dufresne and Style2022). P _st is very similar to the concept of crystallization pressure, found when confined crystals grow from supersaturated solutions (Flatt, Reference Flatt2002; Steiger, Reference Steiger2005; Desarnaud and others, Reference Desarnaud, Bonn and Shahidzadeh2016) and to the concept of condensation pressure, when phase separation occurs in confinement (Style and others, Reference Style2018; Fernández-Rico and others, Reference Fernández-Rico, Sai, Sicher, Style and Dufresne2021).

Figure 1. (a) Ice generates pressure as it grows in a closed cavity, due to the expansion of water upon freezing. (b) Ice growing in an open pore is fed by nearby water, and this growth wedges open the cavity, generating stresses.

Theoretical descriptions of these stress-generation mechanisms are rooted in the (generalized) Clapeyron equation, a fundamental equation that describes static equilibrium between a solid (ice) at pressure P _s, and a reservoir of liquid (water) at a different pressure P _l, but at the same temperature, T (Black, Reference Black1995; Henry, Reference Henry2000; Wettlaufer and Worster, Reference Wettlaufer and Worster2006):

(1)$$( P_{\rm s}-P_{\rm l}) + ( P_{\rm l}-P_0) \left(1-{\rho_{\rm s}\over \rho_{\rm l}}\right) = \rho_{\rm s} q_{\rm m}{( T_{\rm m}-T) \over T_{\rm m}}.$$

Here, ρ _l and ρ _s are the densities of water and ice respectively, q _m is the specific latent heat of freezing of ice and T _m is the melting temperature at a reference pressure, P ₀ (often taken as atmospheric pressure). For the freezing mechanisms described above, this equation can be used to predict P _st as a function of the temperature, T. For case (i) with ice growing in a closed cavity, the ice and water are both at the same pressure (P _s = P _l = P _st) (ignoring capillary effects), so

(2)$$( P_{{\rm st}}-P_0) \left({1\over \rho_{\rm s}}-{1\over \rho_{\rm l}}\right) = {q_{\rm m}( T_{\rm m}-T) \over T_{\rm m}}.$$

Using values from Table 1, we find that ice can exert pressures of ~11 MPa per degree of undercooling (T _m − T).

Table 1. Ice/water parameter values at atmospheric pressure and 273.15 K (Hobbs, Reference Hobbs2010)

For case (ii), the ice and water need no longer have the same pressure. If the water reservoir is held at the reference pressure P _l = P ₀, then P _st = P _s and

(3)$${( P_{{\rm st}}-P_0) \over \rho_{\rm s}} = {q_{\rm m}( T_{\rm m}-T) \over T_{\rm m}}.$$

In this case, ice can exert pressures of ~1 MPa per degree of undercooling.

Even when ice is not in equilibrium (e.g. it is growing), the Clapeyron equation gives us useful information. During growth, there is no macroscopic equilibrium, but water immediately adjacent to an ice surface can often be considered to be in equilibrium with the ice (Wettlaufer and Worster, Reference Wettlaufer and Worster2006). Then, the Clapeyron equation relates the local hydrodynamic pressure in the water, P _l, to the local pressure that has been built up in the ice (P _l is the pressure that would exist in a bulk reservoir of water that was connected to, and in thermodynamic equilibrium with water at the ice interface – note that this definition works even for water in premelted films). Water flows along non-hydrostatic gradients in P _l, so the Clapeyron equation allows us to predict how water is transported toward (or away from) ice, and thus gives ice growth/melting rates (Derjaguin and Churaev, Reference Derjaguin and Churaev1986; Wettlaufer and Worster, Reference Wettlaufer and Worster1995; Rempel and others, Reference Rempel, Wettlaufer and Worster2004; Style and Worster, Reference Style and Worster2005; Wettlaufer and Worster, Reference Wettlaufer and Worster2006).

The various applications of the Clapeyron equation make it a key tool for understanding freezing processes (e.g. Dash and others, Reference Dash, Rempel and Wettlaufer2006; Wettlaufer and Worster, Reference Wettlaufer and Worster2006; Vlahou and Worster, Reference Vlahou and Worster2015; Gerber and others, Reference Gerber, Wilen, Poydenot, Dufresne and Style2022). However, it makes a number of assumptions. For example, it assumes that ice can be described by an isotropic pressure, whereas ice is often characterized by an anisotropic stress state, σ _ij (Budd and Jacka, Reference Budd and Jacka1989) – for example, in glacier flow (Glen, Reference Glen1955), or because anisotropic stresses arise spontaneously in a temperature gradient (Gerber and others, Reference Gerber, Wilen, Poydenot, Dufresne and Style2022). It also uses linear approximations that are valid only near the bulk melting point of ice (see later). Thus, several key questions arise. In particular: What is the appropriate extension of the Clapeyron equation for anisotropically stressed ice? Over what range of conditions should the Clapeyron equation be applicable?

Surprisingly, we are not aware of a systematic derivation of the Clapeyron equation that would allow us to address these questions. However, there are several related works. For example, several authors have established the thermodynamic relations that govern the dissolution of anisotropically stressed solids into adjacent fluids (Gibbs, Reference Gibbs1879; Kamb, Reference Kamb1959, Reference Kamb1961), with notable applications to recrystallization and pressure solution processes (e.g. Paterson, Reference Paterson1973). Although melting was not a focus of these works, some of the consequences for ice melting were recognized by Nye (Reference Nye1967). He argued that the phenomenon of wire regelation requires a generalization of Eqn (2). For this case, P _s should be replaced by the normal stress −σ _nn, and not by the mean of the principal stresses −Tr(σ)/3, as had been argued by others. Finally, Sekerka and Cahn (Reference Sekerka and Cahn2004) examined the special case of a solid with σ _nn = −P _l, to show that anisotropically stressed solids in equilibrium with their melt will recrystallize to form an isotropically stressed state.

Here, we provide a first-principles derivation of the generalized Clapeyron equation, along similar lines to Paterson (Reference Paterson1973). We clearly lay out all the underlying assumptions, and present the appropriate extension for the melting behavior of anisotropically stressed ice.

2. Deriving the generalized Clapeyron equation

We consider thermodynamic equilibrium for the two scenarios shown in Figure 1, in both of which the temperature is held fixed at T < T _m. In case (i), water freezes in a closed cavity, so that the ice and and water both have the same pressure, P _s = P _l. In case (ii), ice has frozen in an open cavity, and is in equilibrium with neighboring bulk water, which has pressure P _l. At the same time, the ice exerts a normal stress, −σ _nn, on the walls of the cavity, but is assumed to exert negligible shear forces, due to the presence of premelted films which lubricate the ice/cavity interface (Gerber and others, Reference Gerber, Wilen, Poydenot, Dufresne and Style2022). The ice cannot grow through the small, connecting pore throat into the neighboring water due to capillarity (i.e. the Gibbs–Thomson effect (Hardy, Reference Hardy1977; Schollick and others, Reference Schollick2016)).

For each scenario, we establish equilibrium behavior by minimizing the relevant free energy of the ice/water system. The relevant free energy, G _sys satisfies ΔG _sys = ΔU _sys − TΔS _sys + W, where U _sys is the internal energy of the ice/water system, S _sys is its entropy and W is the work done by the system on its surroundings. For case (i), ΔG _sys = ΔU _sys − TΔS _sys + P _l(ΔV _s + ΔV _l), while for case (ii), ΔG _sys = ΔU _sys − TΔS _sys − σ _nnΔV _s + P _lΔV _l where V _s and V _l are the volumes of ice and water, respectively. The first case is just a specialized version of the second, where −σ _nn = P _l. Thus, without loss of generality, we can proceed with the case (ii) expression, and the result will describe both cases.

We consider a perturbation to the system in Figure 1b, where a small mass of ice, Δm, melts and flows into the reservoir. Thus, the volumes of ice and water change as ΔV _s = −v _sΔm, and ΔV _l = v _lΔm, where v _s(σ_ij,  T) and v _l(P _l,  T) are the specific volumes of the ice and water, respectively. At equilibrium, this perturbation must not change the free energy, so ΔG _sys = 0, which becomes

(4)$$u_{\rm l}\Delta m - u_{\rm s}\Delta m -T( s_{\rm l}-s_{\rm s}) \Delta m + \sigma_{nn} v_{\rm s}\Delta m + P_{\rm l} v_{\rm l}\Delta m = 0.$$

Here, u _s(σ_ij,  T) and u _l(P _l,  T) are the specific internal energies of the ice and water, respectively, and s _s(σ_ij,  T) and s _l(P _l,  T) are the respective specific entropies. Dividing through by Δm, we obtain

(5)$$-( \sigma_{nn} v_{\rm s} + P_{\rm l} v_{\rm l}) = ( u_{\rm l} - u_{\rm s}) -T( s_{\rm l} - s_{\rm s}) .$$

In principle, Eqn (5) completely describes equilibrium between ice and water – i.e. one could use tabulated values of u,  v, and s to find −σ _nn(P _l,  T). However, a more convenient form is found by expressing the equation relative to the pressure and temperature under bulk melting reference conditions, (P ₀,  T _m). With −σ _nn = P _l = P ₀ Eqn (5) becomes

(6)$$P_0( v_{\rm s}^{\rm o} - v_{\rm l}^{\rm o}) = ( u_{\rm l}^{\rm o} - u_{\rm s}^{\rm o}) -T_{\rm m}( s_{\rm l}^{\rm o} - s_{\rm s}^{\rm o}) ,\; $$

where the superscript ° indicates reference conditions. Subtracting Eqns (5) and (6), we find

(7)$$g_{\rm l}-g_{\rm l}^{\rm o} = g_{\rm s}-g_{\rm s}^{\rm o},\; $$

where the specific free energies g _l(T,  P _l) = u _l − Ts _l + P _lv _l, and g _s(T,  σ _ij) = u _s − Ts _s − σ _nn v _s. These can be Taylor-expanded to obtain the Clapeyron equation (e.g. Dash and others, Reference Dash, Rempel and Wettlaufer2006; Hütter and Tervoort, Reference Hütter and Tervoort2008):

(8)$$\eqalign{g_{\rm l}( T,\; v) = g_{\rm l}^{\rm o}( T_{\rm m},\; P_0) + \left({\partial g_{\rm l}\over \partial T}\right)_{P_{\rm l}}( T-T_{\rm m}) + \left({\partial g_{\rm l}\over \partial P_{\rm l}}\right)_T( P_{\rm l}-P_0) }$$

and

(9)

where δ _ij is the identity matrix. To evaluate the derivatives, we note that Δg _l = −s _lΔT + v _lΔP _l. Thus, at reference conditions,

(10)$$\left({\partial g_{\rm l}\over \partial T}\right)_{P_{\rm l}} = -s_{\rm l}^{\rm o},\; \quad \left({\partial g_{\rm l}\over \partial P_{\rm l}}\right)_T = v_{\rm l}^{\rm o},\; $$

and similarly in the solid at reference conditions (σ _ij = −P ₀δ _ij)

(11)$$\left({\partial g_{\rm s}\over \partial T}\right)_{\sigma_{ij}} = -s_{\rm s}^o.$$

To calculate the final derivative, we notice that $g_{\rm s} = ( f_{\rm s} + P_0v_{\rm s}) -\bar {\sigma }_{nn}v_{\rm s}$, where f _s is the specific Helmholtz free energy of the solid, and $\bar {\sigma }_{ij} = \sigma _{ij} + P_0\delta _{ij}$. Here, $( f_{\rm s} + P_0v_{\rm s}) /v_{\rm s}^{\rm o}$ is the free-energy per unit volume for deformations in an atmosphere at constant pressure, P ₀, and thus is the elastic energy per volume of ice in the reference state. Assuming that ice has linear-elastic behavior, we can write

(12)$$g_{\rm s} = {1\over 2}\bar{\sigma}_{ij}\epsilon_{ij}v_{\rm s}^{\rm o} - \bar{\sigma}_{nn} v_{\rm s}^{\rm o} \left(1 + {{\rm Tr}( \bar{\sigma}) \over 3K_{\rm s}}\right),\; $$

where K _s is now the bulk modulus of the solid. $\epsilon _{ij}$ is the strain of the ice relative to its shape in the reference state (T _m,  P ₀), and satisfies the linear-elastic constitutive relationship:

(13)$$\epsilon_{ij} = {1\over E_{\rm s}}\left[( 1 + \nu_{\rm s}) \bar{\sigma}_{ij}-\nu_{\rm s} \delta_{ij}{\rm Tr}( \bar{\sigma}) \right],\; $$

where E _s = 3K _s(1 − 2ν _s) is the Young's modulus of the ice and ν _s is its Poisson ratio. For small strains, $v_{\rm s} = v_{\rm s}^{\rm o}( 1 + {\rm Tr}( \epsilon ) )$, and we use this in the second term of Eqn (12).

With the two equations above, we can evaluate the remaining derivative at (P ₀,  T _m):

(14)$$\left({\partial g_{\rm s}\over \partial \sigma_{ij}}\right)_T( \bar{\sigma}_{ij} = 0) = -n_i n_j v_{\rm s}^{\rm o}.$$

Here, n _i is the normal vector to the surface of the ice, so that $\bar {\sigma }_{nn} = n_i\bar {\sigma }_{ij}n_j$.

Finally, we can insert these first-derivative expressions into Eqns (7)–(9) to obtain the Clapeyron equation for anisotropically stressed solids:

(15)$$-{( \sigma_{nn} + P_0) \over \rho_{\rm s}^{\rm o}} - {( P_{\rm l}-P_0) \over \rho_{\rm l}^{\rm o}} = {q_{\rm m}( T_{\rm m}-T) \over T_{\rm m}}.$$

Here, $\rho _{\rm l}^{\rm o} = 1/v_{\rm l}^{\rm o}$, and $\rho _{\rm s}^{\rm o} = 1/v_{\rm s}^{\rm o}$ are the densities of water and ice respectively at the bulk melting point, and $q_{\rm m}\equiv ( s_{\rm l}^{\rm o}-s_{\rm s}^{\rm o}) T_{\rm m}$. Consistent with the regelation analysis of Nye (Reference Nye1967), this version of the Clapeyron equation is identical to Eqn (1), but with P _s replaced by −σ _nn, and not −Tr(σ)/3, as some might assume (Verhoogen, Reference Verhoogen1951).

3. Field data supporting the anisotropic Clapeyron equation

While Nye (Reference Nye1967) has presented arguments supporting the form of Eqn (15) in the context of regelation, further evidence comes from simultaneous measurements of temperatures and liquid pressures in glacier boreholes. These measurements show that temperatures increase when changes in the hydrologic system cause borehole pressures, P _l, to decrease (e.g. Andrews and others, Reference Andrews2014).

The anisotropic Clapeyron equation indeed recovers this correlation. Along borehole walls, σ _nn = −P _l. Inserting this into Eqn (15), we find that changes in temperature are correlated with changes in borehole pressure by:

(16)$$\Delta T = -{T_{\rm m}\over q_{\rm m}}\left({1\over \rho_{\rm s}^0}-{1\over \rho_{\rm l}^0}\right)\Delta P_{\rm l}\approx \left(-7.16\times 10^{-8}\, {\rm K\, Pa^{-1}}\right)\Delta P_{\rm l}\; ,\; $$

in agreement with the field data.

By contrast, extending the isotropic Clapeyron Eqn (1), by replacing −P _s = Tr(σ)/3, does not match the experimental data. The classic analysis of Nye (Reference Nye1953) gives the complete stress tensor at the surface of an idealized cylindrical borehole containing liquid at pressure P _l. Far from the borehole, the ice has a far-field isotropic ice pressure P _∞, and creeps according to Glen's flow law with exponent n = 3 (Glen, Reference Glen1955; Hewitt and Creyts, Reference Hewitt and Creyts2019). In this case, −Tr(σ)/3 = P _l + (P _∞ − P _l)/n. Substituting P _s = −Tr(σ)/3 into the isotropic Clapeyron Eqn (1) and treating the far-field ice pressure as constant leads to

(17)$$\Delta T = -{T_{\rm m}\over q_{\rm m}}\left({1\over \rho_{\rm s}^0}-{1\over \rho_{\rm l}^0}-{1\over n\rho_{\rm s}^0}\right)\Delta P_{\rm l}\approx\left(2.26\times 10^{-7}\, {\rm K\, Pa^{-1}}\right)\Delta P_{\rm l}\; .$$

This predicts the opposite of the correlation seen in the field data.

4. Errors in the Clapeyron equation

In deriving this version of the Clapeyron equation, we have had to make two main assumptions. Firstly, strains in the ice are small, so we can use linear elasticity (Sekerka and Cahn, Reference Sekerka and Cahn2004). This is reasonable as the stresses in the ice (which are O(MPa) – see introduction) are much less than the ice's elastic moduli E _s,  K _s = O(GPa), so strains will be small.

Secondly, we assume that higher-order terms in the expansions of g _l and g _s are negligible. We can test this by reverting to the case of isotropically stressed ice (σ _ij = −P _sδ _ij). Then, we Taylor-expand Eqn (7) in T,  P _l and P _s to obtain the second-order version of the Clapeyron equation:

(18)$$\eqalign{& {( P_{\rm s}-P_0) \over \rho_{\rm s}^{\rm o}} - {( P_{\rm l}-P_0) \over \rho_{\rm l}^{\rm o}} = {q_{\rm m}( T_{\rm m}-T) \over T_{\rm m}}\cr & \quad-{c_{\rm l}^{\rm p}-c_{\rm s}^{\rm p}\over 2T_{\rm m}}( T_{\rm m}-T) ^2-{1\over 2\rho_{\rm l}^{\rm o} K_{\rm l}}( P_{\rm l}-P_0) ^2 + {1\over 2\rho_{\rm s}^{\rm o} K_{\rm s}}( P_{\rm s}-P_0) ^2 \cr & \quad-{\alpha_{\rm l}\over \rho_{\rm l}^{\rm o}}( T_{\rm m}-T) ( P_{\rm l}-P_0) + {\alpha_{\rm s}\over \rho_{\rm s}^{\rm o}}( T_{\rm m}-T) ( P_{\rm s}-P_0) . }$$

Here, we use the following identities (e.g. Venerus and Öttinger, Reference Venerus and Öttinger2018): ∂²g/∂T ² = c ^p/T _m, ∂²g/∂P ² = 1/(Kρ ^o) and ∂²g/∂T∂P = α/ρ ^o, where c ^p is the heat capacity at constant pressure, K is again the isothermal bulk modulus and α is the coefficient of thermal expansion.

We can now predict the pressure-melting curve for different freezing scenarios. For bulk ice/water equilibrium (Fig. 1a), P _s = P _l, and we take atmospheric pressure, P _a, as the reference pressure, and T _m = 273.15 K. Figure 2a compares the isotropic Clapeyron Eqn (1) (red, dashed) with experimental data (black, dotted) (Dunaeva and others, Reference Dunaeva, Antsyshkin and Kuskov2010). There is a significant error between the two results for an undercooling of more than $\sim\!\! 3\, ^\circ$C. However, when we use the full, second-order Clapeyron Eqn (18) (blue), we find good agreement down to an undercooling of at least 15$\, ^\circ$C. In this situation, the terms that are quadratic in pressure dominate the error, and to excellent approximation (Fig. 1a, orange dash-dotted):

(19)$$\left({1\over \rho_{\rm s}^{\rm o}}-{1\over \rho_{\rm l}^{\rm o}}\right)( P_{\rm s}-P_0) + \left({1\over 2\rho_{\rm l}^{\rm o} K_{\rm l}}-{1\over 2\rho_{\rm s}^{\rm o} K_{\rm s}}\right)( P_{\rm s}-P_0) ^2 = {q_{\rm m}( T_{\rm m}-T) \over T_{\rm m}}.$$

Note this equation offers a way to extract information about material properties from a pressure/temperature phase diagram, as the curvature of the liquidus is controlled by the quadratic term's prefactor. Comparing the first two terms in the equation, we see that the linear Eqn (15) is only appropriate when:

(20)$$\vert P_{\rm s}-P_{\rm a}\vert \ll \Delta P^\ast = \left\vert \left({1\over \rho_{\rm s}^{\rm o}}-{1\over \rho_l^{\rm o}}\right)\left({1\over 2\rho_{\rm s}^{\rm o}K_{\rm s}}-{1\over 2\rho_{\rm l}^{\rm o}K_{\rm l}}\right)^{-1}\right\vert \approx 420\, {\rm MPa}.$$

Typically, a factor of 10 suffices for such inequalities to hold. Thus, we expect the linear theory to hold when |P _s − P _a| < ΔP*/10 = 42 MPa: in good agreement with the data.

Figure 2. Evaluating the accuracy of the Clapeyron Equation. (a) The pressure of ice in bulk ice/water equilibrium in a closed cavity (P _l = P _s = −σ _nn), as a function of undercooling. The black, dotted curve shows experimental data (Dunaeva and others, Reference Dunaeva, Antsyshkin and Kuskov2010). (b) The stress exerted by ice in an open pore, as a function of undercooling. Both figures show the linear Clapeyron equation (dashed red), full second-order theory (Eqn (18), blue) and simplified second-order theory (Eqns (19) and (21), orange dash-dotted).

We can perform a similar analysis for freezing in an open system (Fig. 1b). We let P _l = P ₀ = P _a, and assume that the ice is in an isotropic state of stress, with pressure, P _s. Figure 2b compares the prediction of the Clapeyron Eqn (1) (red, dashed) with that obtained when we keep the extra quadratic terms (18) (blue). We are not aware of any experimental data precise enough to validate the theory (Gerber and others, Reference Gerber, Wilen, Poydenot, Dufresne and Style2022). However, here, the higher-order theory agrees well with the linear Clapeyron equation down to large undercoolings. The difference is dominated by the term in Eqn (18) that is quadratic in undercooling. Thus, to excellent approximation (Fig. 1a, orange dash-dotted):

(21)$$P_{\rm s}-P_{\rm a} = {\rho_{\rm s}^{\rm o} q_{\rm m}( T_{\rm m}-T) \over T_{\rm m}} - {\rho_{\rm s}^{\rm o}( c_{\rm l}^{\rm p}-c_{\rm s}^{\rm p}) \over 2T_{\rm m}}( T_{\rm m}-T) ^2.$$

Comparing terms on the right-hand side shows that we only recover the linear Clapeyron Eqn (1) if

(22)$$\vert T_{\rm m}-T\vert \ll\Delta T^\ast = \left\vert {2q_{\rm m}\over c_{\rm l}^{\rm p}-c_{\rm s}^{\rm p}} \right\vert \approx 320 \, {\rm K}.$$

Again, assuming a factor of 10 for the inequality to hold, we find that linear theory should work when |T _m − T| ≪ ΔT*/10 = 32 K. This requirement is certainly reasonable for many terrestrial temperatures. Thus, there is some justification for use of the linear Clapeyron equation down to relatively large undercoolings to model this type of freezing scenario.

To summarize, our results suggest that the linearized Clapeyron equation will be valid, provided that |P _l − P ₀| and |P _s − P ₀| are both less than ΔP*/10, while |T _m − T| < ΔT*/10. At larger pressures/undercoolings, the quadratic terms in Eqn (18) should be included.