A non-equilibrium model for soil heating and moisture transport during extreme surface heating: The soil (heat-moisture-vapor) HMV-Model Version

Increased use of prescribed fire by land managers and the increasing likelihood of wildfires due to climate change require an improved modeling capability of extreme heating of soils during fires. This issue is addressed here by developing and testing the soil (heat–moisture–vapor) HMVmodel, a 1-D (one-dimensional) non-equilibrium (liquid– vapor phase change) model of soil evaporation that simulates the coupled simultaneous transport of heat, soil moisture, and water vapor. This model is intended for use with surface forcing ranging from daily solar cycles to extreme conditions encountered during fires. It employs a linearized Crank–Nicolson scheme for the conservation equations of energy and mass and its performance is evaluated against dynamic soil temperature and moisture observations, which were obtained during laboratory experiments on soil samples exposed to surface heat fluxes ranging between 10 000 and 50 000 W m. The Hertz–Knudsen equation is the basis for constructing the model’s non-equilibrium evaporative source term. Some unusual aspects of the model that were found to be extremely important to the model’s performance include (1) a dynamic (temperature and moisture potential dependent) condensation coefficient associated with the evaporative source term, (2) an infrared radiation component to the soil’s thermal conductivity, and (3) a dynamic residual soil moisture. This last term, which is parameterized as a function of temperature and soil water potential, is incorporated into the water retention curve and hydraulic conductivity functions in order to improve the model’s ability to capture the evaporative dynamics of the strongly bound soil moisture, which requires temperatures well beyond 150 C to fully evaporate. The model also includes film flow, although this phenomenon did not contribute much to the model’s overall performance. In general, the model simulates the laboratoryobserved temperature dynamics quite well, but is less precise (but still good) at capturing the moisture dynamics. The model emulates the observed increase in soil moisture ahead of the drying front and the hiatus in the soil temperature rise during the strongly evaporative stage of drying. It also captures the observed rapid evaporation of soil moisture that occurs at relatively low temperatures (50–90 C), and can provide quite accurate predictions of the total amount of soil moisture evaporated during the laboratory experiments. The model’s solution for water vapor density (and vapor pressure), which can exceed 1 standard atmosphere, cannot be experimentally verified, but they are supported by results from (earlier and very different) models developed for somewhat different purposes and for different porous media. Overall, this non-equilibrium model provides a much more physically realistic simulation over a previous equilibrium model developed for the same purpose. Current model performance strongly suggests that it is now ready for testing under field conditions.

Abstract.Increased use of prescribed fire by land managers and the increasing likelihood of wildfires due to climate change require an improved modeling capability of extreme heating of soils during fires.This issue is addressed here by developing and testing the soil (heat-moisture-vapor) HMVmodel, a 1-D (one-dimensional) non-equilibrium (liquidvapor phase change) model of soil evaporation that simulates the coupled simultaneous transport of heat, soil moisture, and water vapor.This model is intended for use with surface forcing ranging from daily solar cycles to extreme conditions encountered during fires.It employs a linearized Crank-Nicolson scheme for the conservation equations of energy and mass and its performance is evaluated against dynamic soil temperature and moisture observations, which were obtained during laboratory experiments on soil samples exposed to surface heat fluxes ranging between 10 000 and 50 000 W m −2 .The Hertz-Knudsen equation is the basis for constructing the model's non-equilibrium evaporative source term.Some unusual aspects of the model that were found to be extremely important to the model's performance include (1) a dynamic (temperature and moisture potential dependent) condensation coefficient associated with the evaporative source term, (2) an infrared radiation component to the soil's thermal conductivity, and (3) a dynamic residual soil moisture.This last term, which is parameterized as a function of temperature and soil water potential, is incorporated into the water retention curve and hydraulic conductivity functions in order to improve the model's ability to capture the evaporative dynamics of the strongly bound soil moisture, which requires temperatures well beyond 150 • C to fully evaporate.The model also includes film flow, although this phenomenon did not contribute much to the model's overall performance.In general, the model simulates the laboratoryobserved temperature dynamics quite well, but is less precise (but still good) at capturing the moisture dynamics.The model emulates the observed increase in soil moisture ahead of the drying front and the hiatus in the soil temperature rise during the strongly evaporative stage of drying.It also captures the observed rapid evaporation of soil moisture that occurs at relatively low temperatures (50-90 • C), and can provide quite accurate predictions of the total amount of soil moisture evaporated during the laboratory experiments.The model's solution for water vapor density (and vapor pressure), which can exceed 1 standard atmosphere, cannot be experimentally verified, but they are supported by results from (earlier and very different) models developed for somewhat different purposes and for different porous media.Overall, this non-equilibrium model provides a much more physically realistic simulation over a previous equilibrium model developed for the same purpose.Current model performance strongly suggests that it is now ready for testing under field conditions.

Introduction
Since the development of the theory of Philip and de Vries (PdV model) almost 60 years ago (Philip and de Vries, 1957;de Vries, 1958), virtually all models of evaporation and condensation in unsaturated soils have assumed that soil water vapor at any particular depth into the soil is in equilibrium with the liquid soil water (or soil moisture) at the same depth.(Note that such soil evaporation models also assume thermal equilibrium, so that at any given depth the mineral soil, Published by Copernicus Publications on behalf of the European Geosciences Union. the soil moisture, and the soil air and water vapor within the pore space are also at the same temperature.)In essence, this local equilibrium assumption means that whenever the soil moisture changes phase it does so instantaneously.This assumption is quite apropos for its original application, which was to describe the coupled heat and moisture transport in soils (and soil evaporation in particular) under environmental forcings associated with the daily and seasonal variations in radiation, temperature, precipitation, etc. (e.g., Milly, 1982;Novak, 2010;Smits et al., 2011).Under these conditions, assuming local equilibrium is reasonable because the time required to achieve equilibrium after a change of phase is "instantaneous" (short) relative to the timescale associated with normal environmental forcing.The great benefit to the equilibrium assumption is that for modeling purposes it is a significant simplification to the equations that describe heat and moisture flow in soils because it eliminates the need to include soil water vapor density, ρ v , as an independent model variable.More formally, under equilibrium ρ v is directly equated to the equilibrium vapor density, a function only of local soil temperature and soil water content (or more specifically the soil water potential).
Subsequent to the development of the original PdV model the equilibrium assumption has also been incorporated into models of heat and moisture transport (evaporation and condensation) in soils and other porous media under more extreme forcings associated with high temperatures and heat fluxes.For example, it has been applied to (i) soils during wildfires and prescribed burns (Aston and Gill, 1976;Campbell et al., 1995;Durany et al., 2010;Massman, 2012), (ii) drying of wood (Whitaker, 1977;di Blasi, 1997), (iii) drying and fracturing of concrete under high temperatures (Dayan, 1982;Dal Pont et al., 2011), (iv) high temperature sand-water-steam systems (e.g., Udell, 1983;Bridge et al., 2003), and (v) evaporation of wet porous thermal barriers under high heat fluxes (Costa et al., 2008).
Although the PdV model and the equilibrium assumption have certainly led to many insights into moisture and vapor transport and evaporation in porous media, they have, nonetheless, yielded somewhat disappointing simulations of the coupled soil moisture dynamics during fires (see Massman, 2012, for further details and general modeling review).Possibly the most interesting of these modeling "disappointments" is the soil/fire-heating model of Massman (2012), who found that as the soil moisture evaporated it just recondensed and accumulated ahead of the dry zone; consequently, no water actually escaped the soil at all, which, to say the least, seems physically implausible.He further traced the cause of this anomalous behavior to the inapplicability of the equilibrium evaporation assumption, which allowed the soil vapor gradient behind the drying front to become so small that the soil vapor could not escape (diffuse) out of the soil.Moreover, more fundamentally, the calculated vapor and its attendant gradient became largely meaningless because it is impossible for water vapor to be in equilibrium with liquid water if there is no liquid water.Of course, such extremely dry conditions are just about guaranteed during soil heating events like fires.Novak (2012) also recognized the inapplicability of the equilibrium assumption for very dry soils, but under more normal environmental forcing.On the other hand, even under normal (and much less extreme) soil moisture conditions both Smits et al. (2011) and Ouedraogo et al. (2013) suggest that non-equilibrium formulations of soil evaporation may actually improve model performance, which implies that the non-equilibrium assumption may really be a more appropriate description for soil evaporation and condensation than the equilibrium assumption.The present study is intended to provide the first test of the non-equilibrium hypothesis during extreme conditions.
Specifically, the present study develops and evaluates a non-equilibrium (liquid-vapor phase change) model for simulating coupled heat, moisture, and water vapor transport during extreme heating events.It also assumes thermal equilibrium between the soil solids, liquid, and vapor.It uses a systems-theoretic approach (e.g., Gupta and Nearing, 2014) focused more on physical processes than simply tuning model parameters; i.e., that whatever model or parameter "tuning" does occur, it is intended to keep the model numerically stable and as physically realistic as possible.
In addition, the present study (model) is a companion to Massman (2012).It uses much of the same notation as the earlier study.But, unlike its predecessor, this study allows for the possibility of liquid water movement (i.e., it includes a hydraulic conductivity function for capillary and film flow).It also improves on and corrects (where possible and as noted in the text) the mathematical expressions used in the previous paper to parameterize the high-temperature dependency of latent heat of vaporization, saturation vapor density, diffusivity of water vapor, soil thermal conductivity, water retention curve, etc.This is done in order to achieve the best representation of the physical properties of water (liquid and vapor) under high temperatures and pressures (see, e.g., Harvey and Friend, 2004).Lastly, in order to facilitate comparing the present model with the earlier companion model the present study displays all graphical results in a manner very similar to those of Massman (2012).

Model development
The present model is one-dimensional (1-D) (in the vertical) and is developed from three coupled partial differential equations.It allows for the possibility that the soil liquid and vapor concentrations are not necessarily in local equilibrium during evaporation/condensation, but it does assume local thermal equilibrium during any phase change.The present model has three simulation variables: soil temperature (≡ T [C] or T K [K]), soil water potential (≡ ψ [J kg −1 ] or ψ n ≡ normalized soil water potential [dimensionless]), and vapor density (≡ ρ v [kg m −3 ]).Here ψ n = ψ/ψ * and ψ * = −10 6 [J kg −1 ], which Campbell et al. (1995) identify as the water potential for oven-dry soil.This current model employs a linearized Crank-Nicolson (C-N) finite difference scheme, whereas the preceding (companion) model (Massman, 2012) used the Newton-Raphson method for solving the fully implicit finite difference equations.The present model further improves on its companion by including the possibility of soil water movement (hydraulic conductivity function driven by a gradient in soil moisture potential) and better parameterizations of thermophysical properties of water and water vapor.These latter parameterizations allow for the possibility of large variations in the amount of soil water vapor, which Massman (2012) suggests might approach or exceed 1 standard atmosphere and therefore could become the major component of the soil atmosphere during a heating event.This is quite unlike any other model of soil heat and moisture flow, which universally assume that dry air is the dominant component of the soil atmosphere and that water vapor is a relatively minor component.Finally, and also atypical of most other soil models, the model's water retention curve and hydraulic function include a dynamic residual soil moisture content as a function of soil temperature and soil water potential.

Conservation equations
The conservation of thermal energy is expressed as where is the volumetric heat capacity of soil, a function of both soil temperature and soil volumetric water content θ [m 3 m −3 ]; is the thermal conductivity of the soil, a function of soil temperature, soil moisture, and soil vapor density; η [m 3 m −3 ] is the total soil porosity from which it follows that (η − θ ) is the soil's air-filled porosity; is the mass density of the soil air, a function of temperature and soil vapor density; ] is specific heat capacity of ambient air, also a function of temperature and vapor density; u vl [m s −1 ] is the advective velocity induced by the change in volume associated with the rapid volatilization of soil moisture (detailed below); is the source term for water vapor.The conservation of mass for liquid water is where is the velocity of liquid water associated with surface diffusion of water, which may be significant at high temperatures (e.g., Kapoor et al., 1989;Medve ď and Černý, 2011).Note that switching variables from ψ < 0 to ψ n produces ψ n > 0 and K n < 0. This last equation can be simplified to dT varies by only 4 % between about 10 to 100 dT ∂T ∂z can be ignored.But the model does retain the temperature dependency ρ w = ρ w (T K ), except as noted in the section below on volumetric-specific heat capacity of soil, and it also specifically includes dρ w /dT for other components of the model.
The conservation of mass for water vapor is where molecular diffusivity associated with the diffusive transport of water vapor in the soil's air-filled pore space.As with Massman (2012), the present model also expresses Fick's first law of diffusion in terms of mass; i.e., the diffusive flux But there are other forms that could have been used.For example, Campbell et al. (1995) use a form discussed by Cowan (1977) and Jones (2014); i.e., , where e v [Pa] is the vapor pressure, M w = 0.01802 kg mol −1 is the molar mass of water vapor, and R = 8.314 Jmol −1 K −1 is the universal gas constant.Yet again, Jaynes and Rogowski (1983) suggest that J diff = J [Fraction]   diff = −D ve (ρ v + ρ d )∂ω v /∂z may be the more appropriate expression for J diff , where ρ d [kg m −3 ] is the dry air density (defined and discussed later) and is the mass fraction of water vapor within the soil pore space.The distinctions between these different formulations of Fick's first law is important to the present work because different forms of J diff can yield different numerical values for the fluxes (e.g., Solsvik and Jakobsen, 2012) that can diverge significantly with large temperature gradients.This issue is examined in more detail in a later section devoted to the model's sensitivity to different modeling assumptions and performance relative to different data sets and input parameters.
The final model equations are expressed in terms of the model variables (T , ψ n , ρ v ) and result from (a) expanding the spatial derivative ∂K H ∂z in terms of the spatial derivatives fying Eq. ( 4).These equations are which is the conservation of energy; which is the conservation of soil moisture; and which is the conservation of mass for water vapor.Apropos to these last three equations (i) D θψ = ∂θ/∂ψ n and D θT = ∂θ/∂T K are obtained from the water retention curve (WRC); (ii n [m 2 s −1 ] (which subsumes K n ) are related to V θ,surf and are defined in a later section; (iv) because ρ v ρ w the term (ρ w − ρ v ) ∂θ ∂t originally in Eq. ( 6) has been approximated by ρ w ∂θ ∂t ≡ ρ w D θψ ∂ψ n ∂t + ρ w D θT ∂T ∂t ; and (v) the total porosity η is assumed to be spatially uniform and temporally invariant.

Thermophysical properties of water, vapor, and moist Air
The algorithm for calculating water density, ρ w (T K ), is Eq.(2.6) of Wagner and Pruess (2002) and employed only within the temperature range 273.15K ≤ T K ≤ 383.15K (≡ T K,max ).At temperatures greater than T K,max , then ρ w (T K ) = ρ w (T K,max ).This approach yields a range for ρ w (T K ) of 950 kg m −3 < ρ w (T K ) < 1000 kg m −3 , which represents a compromise between the fact that the density of (free saturated liquid) water continues to decrease with increasing temperatures (Yaws, 1995) and the hypothetical possibility that in a bound state a mono-layer of liquid water ρ w (T K ) may reach values as high as ≈ 5000 kg m −3 (Danielewicz-Ferchmin and Mickiewicz, 1996).dρ w /dT is computed from the analytical expression derived by differentiating the expression for ρ w (T K ) and dρ w /dT = 0 for T K > T K,max .
The enthalpy of vaporization of water, Somayajulu (1988) augmented by the soil moisture potential, ψ, (Massman, 2012;Campbell et al., 1995) and is expressed as follows: where H 1 = 13.405538kJ mol −1 , H 2 = 54.188028kJ mol −1 , H 3 = −58.822461kJ mol −1 , and T crit = 647.096K is the critical temperature for water.Note that h v = h v (T K ) [J mol −1 ] will denote the enthalpy of vaporization without the additional −M w ψ term, i.e., . The present formulation differs from Massman (2012) because here h v (T K ≥ T crit ) = 0; whereas the equivalent from Massman (2012) (and Campbell et al., 1995) h v was a linear approximation of the present h v , which yielded h v (T K ≥ T crit ) 0. This distinction will become important when discussing the water vapor source term, S v .Note that because L v = H v /M w it also employs Eq. ( 8).
The formulations for thermal conductivity of water vapor, and liquid water, Huber et al. (2012).For water vapor their Eq.( 4) is used and for liquid water the product of their Eqs.(4) and ( 5) is used.The formulations for viscosity of water vapor and liquid water are taken from Huber et al. (2009) and are similar algorithmically to thermal conductivity.For water vapor, ≡ Pa s], their Eq.( 11) is used and for liquid water, µ w = µ w (T K , ρ w ) [Pa s], their Eq.( 36) is used.For liquid water both these formulations include a dependence on the density of water.Consequently, once soil temperature exceeds T K,max both λ w and µ w are assigned a fixed value determined at T K,max .On the other hand, λ v and µ v increase continually with increasing temperatures.
The formulation for the thermal conductivity of dry air, Kadoya et al. (1985) and for the viscosity of dry air, µ d = µ d (T K ) [Pa s], Eq. (3a) of Kadoya et al. (1985) is used.The model of the thermal conductivity of soil atmosphere, , is a nonlinear expression given by Eq. (28) of Tsilingiris (2008).The relative weights used in this formulation are determined using the mixing ratios for water vapor (χ v [dimensionless]) and dry air (χ d = 1 − χ v ): where χ v = e v /(P d + e v ), and P d [Pa] is the dry air pressure.Here P d will be held constant and equal to the external ambient atmospheric pressure, P atmos (i.e., 92 kPa), during the laboratory experiments (see Massman, 2012;Campbell et al., 1995).The vapor pressure, e v , is obtained from ρ v and T K using the ideal gas law.
The volumetric-specific heat for soil air, ρ a c pa [J m −3 K −1 ], is estimated for the soil atmosphere from ρ a c pa = c pv ρ v + c pd ρ d , where ρ d = M d P atmos /(RT K ) [kg m −3 ] is the dry air density; M d = 0.02896 kg mol −1 is the molar mass of dry air; and the isobaric-specific heats for water vapor, c pv [J kg −1 K −1 ], and dry air, c pd [J kg −1 K −1 ], use Eq. ( 6) of Bücker et al. (2003).
The saturation vapor pressure, e v,sat = e v,sat (T K ) [Pa], and its derivative, de v,sat /dT [Pa K −1 ], are modeled using Eqs.(2.5) and (2.5a) of Wagner and Pruess (2002).The saturation vapor density, ρ v,sat [kg m −3 ], is modeled using Eq.(2.7) of Wagner and Pruess (2002).Following Massman (2012), these saturation curves are restricted to temperatures below that temperature, T K,sat , at which e v,sat = P atmos .For the present case, T K,sat = 370.44K was determined from the saturation temperature equation or "the Backward Equation", Eq. (31) of IAPWS (2007).For T K ≥ T K,sat , the saturation quantities e v,sat and de v,sat /dT remain fixed at their values T K,sat , but ρ v,sat is allowed to decrease with increasing temperatures, i.e., ρ v,sat = ρ v,sat (T K,sat )[T K,sat /T K ][P atmos /P ST ], in accordance with Table 13.2 (p.497) of Wagner and Pruess (2002), where P ST = 101325 Pa is the standard pressure.The present treatment of ρ v,sat is different from Massman (2012), who assumed that ρ v,sat (T K ≥ T K,sat ) = ρ v,sat (T K,sat ).
Embedded in the hydraulic conductivities (K H and K n of Eq. 2) are the surface tension of water, σ w (T K ) [N m −1 ], and the static dielectric constant (or the relative permittivity) of water, w (T K ) [dimensionless].These physical properties of water are integral to the hydraulic conductivity model from Zhang (2011) associated with film flow, which is incorporated into the present model and detailed in a later section.The surface tension of water is modeled following Eq.(1) of Vargaftik et al. (1983) and w (T K ) is taken from Eq. (36) of Fernández et al. (1997).

Functions related to water vapor
where τ [m m −1 ] is the tortuosity of soil with τ = 0.66[(η − θ)/η] 3 after Moldrup et al. (1997); E [dimensionless] is the vapor flow enhancement factor and is discussed in Massman (2012); D v /(1 + e v /P atmos ) [m 2 s −1 ] is the molecular diffusivity of water vapor into the soil atmosphere, which will be taken as a mixture of both dry air and (potentially large amounts of) water vapor; and S F [dimensionless] is Stefan correction or mass flow factor.Externaliz-ing the 1/(1 + e v /P atmos ) term of the vapor diffusivity and combining it with S F allows for the following approximation for S F /(1 + e v /P atmos ) = 1/(1 − e 2 v /P 2 atmos ) ≈ 1 + e v /P atmos ; where the correct form for S F is 1/(1 − e v /P atmos ).The reason for this approximation is to avoid dividing by 0 as e v → P atmos .Massman (2012) and Campbell et al. (1995) avoided this issue by limiting S F to a maximum value of 10/3.This newer approximation for S F /(1 + e v /P atmos ) is an improvement over their approach.It is relatively slower to enhance the vapor transport by diffusion than the original approach, but this is only of any real significance when S F ≥ 1.On the other hand because the linear form is not limited to any preset maximum value, it compensates for these underestimations when S F > 10/3.
Mindful of the externalization of 1/(1 + e v /P atmos ), D v is estimated from the diffusivity of water vapor in dry air, D vd = D vd (T K ) [m 2 s −1 ] and the self-diffusivity of water vapor in water vapor, where and and α vv = 9/4.The parameters D vvST and α vv relate to the self-diffusion of water vapor and their numerical values were determined from a synthesis of results from Hellmann et al. (2009) and Yoshida et al. (2006Yoshida et al. ( , 2007)).The uncertainty associated with this value for D vvST is at least ±15 % and possibly more, e.g., Miles et al. (2012).Blanc's law (Marrero and Mason, 1972) combines D vd and D vv to yield the following expressions for The present model for the advective velocity associated with the volatilization of water, u vl , is taken from Ki et al. (2005) and is non-equilibrium equivalent to that used by Massman (2012) in his equilibrium model.Here where the basic assumptions are that both liquid water and vapor are Newtonian fluids and that only incompressible effects are being modeled.In essence Eq. ( 9) assumes that the vaporization of soil moisture acts as a steadystate (and rapidly expanding or "exploding") volume source term, which yields a 1-D advective velocity associated with volatilization of liquid water.For an equilibrium model of soil moisture evaporation that does not include water movement (i.e., K H ≡ 0, K n ≡ 0, and V θ,surf ≡ 0), then S v ≡ −ρ w ∂θ/∂t (from Eq. 2 above), which demonstrates the connection between present model of u vl with that used by Massman (2012).But unlike Massman (2012), the present model does not require any numerical adjustments to Eq. ( 9) in order to maintain numerical stability.
The functional parameterization of S v follows from the non-equilibrium assumption, i.e., S v ∝ (ρ ve − ρ v ), where ρ ve = a w ρ v,sat (T K ) [kg m −3 ] is the equilibrium vapor density and a w = e Mwψ * RT K ψ n [dimensionless] is the water activity, modeled here with the Kelvin equation.The difficult part is how to construct the proportionality coefficient.Nevertheless, there are at least a two ways to go about this: (a) largely empirically (e.g., Smits et al., 2011, and related approaches referenced therein); or (b) assume that S v = A wa J v (e.g., Skopp, 1985, or Novak, 2012), where A wa [m −1 ] is the volume-normalized soil water-air interfacial surface area and J v [kg m −2 s −1 ] is the flux to/from that interfacial surface.This second approach allows for a more physically based parameterization of the flux, viz., Novak (2012) proposed that the flux be driven be diffusion, so that R v = D v /r ep , where r ep [m] is the equivalent pore radius and D v is the diffusivity of water vapor in soil air.After a bit of algebra and some simple geometrically based assumptions concerning the relationships between r ep , a spherical pore volume, and A wa , one arrives at the Novak model of the source term: where S [N] * [dimensionless] is an adjustable model parameter.
But there is another way to model the vapor flux, J v , which is also used in the present study.This second approach is based on the Hertz-Knudsen equation, which has its origins in the kinetic theory of gases and describes the net flux of a gas that is simultaneously condensing on and evaporating from a surface.A general expression for the Hertz-Knudsen flux is where K e [dimensionless] is the mass accommodation (or evaporation) coefficient and is the thermal accommodation (or condensation) coefficient.For the present purposes K e ≡ 1 can be assumed.This model of J v yields the following model for S v : where [dimensionless] is an adjustable model parameter, to be determined by "tuning" it as necessary to ensure model stability.This model for S [M]   v is now more or less complete, but the model for S [N ]   v is neither quite complete nor precisely comparable to S [M]  v .This is now remedied by introducing K c into S [N]   v and subsuming a factor of A wa into S [N]   * , yielding: where S [N]    * [m −1 ] now has physical dimensions, but otherwise remains an adjustable parameter that will be scaled such that S [N]   v ≈ O(S [M]  v ).In this form these two models for S v can be used to test the sensitivity of the model's solution to different temperature forcing, because where T K,in [K] is the initial temperature of the laboratory experiments and E av −M w ψ [J mol −1 ] is an empirically determined surface condensation/evaporation activation energy.(Note that (a) T K ≥ T K,in , valid most of the time for any simulation, guarantees K c ≤ 1.(b) The enthalpy of vaporization, h v (T K ), is a logical choice for E av ; but, model performance was significantly enhanced by simply assigning a constant value for E av ≈ 30-40 kJ mol −1 rather than assigning E av ≡ h v .)The present formulation ensures that ∂K c /∂T K < 0, in accordance with experimental and theoretical studies (Tsuruta and Nagayama, 2004;Kon et al., 2014).Mathematically this present formulation of K c largely eliminates model instabilities by suppressing condensation relative to evaporation throughout the experiment and will be discussed in greater detail in a later section.
A wa is parameterized as a parabolic function to simulate the conceptual model of A wa proposed by Constanza-Robinson and Brusseau (2002; see their Fig.1b): where S w = θ/η is the soil water saturation and a 1 = 40, a 2 = 0.003, and a 3 = 1/8.This particular functional form ensures that A wa = 0 when the soil is completely dry, θ = 0, and when fully saturated, S w = 1.This particular parameter value for a 1 was chosen so that the maximum value of A wa occurs at S w = 0.025 (i.e., 1/a 1 ) in accordance with the model of Constanza-Robinson and Brusseau (2002).

Thermal transport properties: C s , λ s
The model for C s (T , θ ) is taken from Massman (2012): is the volumetric heat capacity of water; and the parameterization constants c s0 , c s1 , C w0 , C w1 , C w2 are given by Massman (2012).Note that the present C s (θ, T ) is the isobaricspecific heat capacity of water, and ρ w (T ST ) = 1000 kg m −3 .This substitution for ρ w (T K ) is only made for C w (T ).
The present formulation for isobaric heat capacity of water, c pw (T ), was developed from Yaws (1995) and confirmed by comparing to Wagner and Pruess (2002).In general c pw (T ) is also a function of pressure (e.g., Wagner and Pruess, 2002), but this dependency can be ignored for the present purposes.Other parameterization of c pw (T ) (i.e., Sato, 1990;Jovanović et al., 2009;Kozlowski, 2012) were also examined, but proved unsatisfactory.Finally, Kozlowski ( 2012) reported numerical values for the dry soil parameters c s0 and c s1 that are similar to those discussed in Massman (2012) and used with the present model.
The model of soil thermal conductivity, λ s , is the sum of two terms.The first, λ [1]  s (θ, T K , ρ v ), is taken principally from Campbell et al. (1994) and the second, λ [2]  s (θ, T K ), is taken from Bauer (1993).This second term incorporates the effects of high-temperature thermal (infrared) radiant energy transfer within the soil pore space, which may be significant for certain soils and high enough temperatures (e.g., Durany et al., 2010).λ [1]  s is summarized first and repeated here to emphasize the difference between the present model's functional parameterizations and those used in Massman (2012).
λ [1]  s is modeled as where is the apparent thermal conductivity of the soil air and is the sum of the thermal conductivity of moist soil air, λ a , and λ * v , which incorporates the effects of latent heat transfer; λ m is the thermal conductivity of the mineral component of the soil, which is assumed to be independent of temperature and soil moisture; and k w , k a , and k m [dimensionless] are generalized formulations of the de Vries weighting factors (de Vries, 1963).Campbell et al. (1994) formulated λ * v as proportional to the product of the enthalpy of vaporization (h v ), the vapor diffusivity (D v ), the Stefan factor (S F ), the slope of the saturation vapor pressure (de v,sat /dT ), and the parameter f w (θ, T ).For the present model where ρST = 44.65 mol m −3 is the molar density of the standard atmosphere.Equations ( 12) and ( 13) are the same as those used in Massman (2012), but numerically they yield quite different results due to the different formulations for h v , S F , D v , and e v,sat = e v,sat (T K ).Otherwise the de Vries (1963) shape factors, the parameter f w , and all related parameters are the same as in Massman (2012).
λ [2]  s is modeled as where σ = 5.670 is the soil's pore space volumetric radius; and the factor of 3.8 subsumes a numerical factor of 4, a [dimensionless] pore shape factor (i.e., 1 for spherical particles), and the [dimensionless] emissivity of the medium ≈ 0.95 (by assumption).In general R p is considered to be a model parameter and it will be varied to evaluate the model's sensitivity and performance relative to λ [2]  s and unless otherwise indicated, R p = 10 −3 m is the nominal or default value.In the present context variations in R p are purely model driven, but in reality such variations are likely to be most strongly associated with (or proportional to) changes in the soil particle dimensions and secondarily with other soil characteristics (e.g., Leij et al., 2002).

Water retention curve
In general a WRC is a functional relationship between soil moisture and soil moisture potential and temperature, i.e., θ = θ (ψ, T ), although the temperature dependency is often ignored and was of little consequence to the model of Massman (2012).The three WRCs tested in the present study have been adapted to include a residual soil moisture, θ r = θ r (ψ, T ) [m 3 m −3 ], which is an atypical parameterization for both θ r and the soil moisture's temperature dependency.Under more normal soil environmental conditions θ r is assumed to be bound so securely to soil mineral surfaces that it is normally taken as a fixed constant.For the present purposes, the principal WRC is adapted from Massman (2012) and is where and θ l [m 3 m −3 ] is the extrapolated value of the water content when ψ = ψ l = −1 J kg −1 ; α l = ln (ψ * /ψ l ) = 13.8155106; and θ h [m 3 m −3 ], α h [dimensionless] and p [dimensionless] are parameters obtained from Campbell and Shiozawa (1992); θ r * [m 3 m −3 ] is a constant soil-specific parameter, such that θ r * ≤ 0.03 is to be expected; and b 1 [dimensionless] and b 2 [dimensionless] are adjustable parameters, which are expected to satisfy b 1 > 0 and 0 ≤ b 2 < 1.Note that a further discussion concerning the original version of Eq. ( 15) can be found in Massman (2012).
There is a simple and physically intuitive argument for the parameterization of θ r (ψ n , T K ) in Eq. ( 16).First, under more www.geosci-model-dev.net/8/3659/2015/Geosci.Model Dev., 8, 3659-3680, 2015 normal soil environmental conditions, i.e., T K ≈ T K,in and at least ψ n < 1 (if not ψ n 1), then it is reasonable to expect that θ r ≈ θ r * and nearly constant throughout (what might be expected to be relatively small variations in) those conditions.But as the temperature increases, it is also reasonable to assume that the increasing amounts of thermal energy will begin to overcome the forces holding the bound water to the soil mineral surfaces and that θ r will decrease.Mathematically, one might therefore expect that ∂θ r /∂T K < 0. Massman (2012; see his discussion of ψ T ) made similar arguments when he included the temperature dependency in his version of the same basic WRC.Consequently, the Eqs.( 15) and ( 16) above offer a different approach to including temperature effects on the WRC that maintains the temperature dependent properties of WRCs outlined by Massman (2012).
Second, as the temperature increases and the soil moisture (including θ r ) begins to decrease, the soil moisture potential ψ n will begin to increase (or ψ < 0 will decrease in absolute terms while increasing in magnitude), which in turn (it is hypothesized) will tend to strengthen the forces holding the bound water.Therefore, one might expect that as the soil dries out ∂θ r /∂ψ n > 0, which will oppose, but not dominate the temperature effects (because b 2 < 1).Equation ( 16) is designed to capture these two opposing influences, assuming of course that T K > T K,in .But, it is itself not intended to be a fully physically based dynamical theory of the residual soil moisture.Such a theory is beyond the intent of the present study.The sole intent here is to test and evaluate whether a dynamical θ r can improve the model's performance.In so far as it may succeed at doing so, it will also indicate the value and need for a more detailed physically based dynamical model of θ r .The present study also includes similar adaptations to two other WRCs so as to test the model's sensitivity to different WRCs.These WRCs, which will not be shown here, are taken from Groenevelt and Grant (2004) and Fredlund and Xing (1994).

Functions related to liquid water transport
The hydraulic conductivity functions, K n (T K , ψ n , θ ) and K H (T K , ψ n , θ ), are given as follows: where g = 9.81 m s −2 is the acceleration due to gravity; K I [m 2 ] is the intrinsic permeability of the soil (a constant for any given soil); and is the hydraulic conductivity function (HCF), which for the present study is taken as the sum of K cap R (associated with capillary flow) and K film R (associated with film flow).The model for intrinsic permeability, which is taken from Bear (1972), is K I = (6.17× 10 −4 )d 2 g ; where d g [m] is the mean or "effective" soil particle diameter.For the soils used in the present work (Campbell et al., 1995;Massman, 2012) d g was estimated from Shiozawa and Campbell (1991) and Campbell and Shiozawa (1992) or simply assigned a reasonable value if no other information was available.
For the present study, five difference parameterizations for K cap R were tested: two were from Grant et al. (2010), i.e., their Eq.( 18) (Burdine) and Eq. ( 19) (Mualem); the Van Genuchten and Nielson (1985) model, their Eq.( 22) with the mathematical constraints imposed as suggested by Assouline and Or (2013); the Brooks and Corey (1964) model; and Eq. ( 18) of Assouline (2001).The reason for testing several models of the HCF is to determine how different formulations for the HCF might impact the model's performance when comparing to the laboratory observations.The following HCF is the model of Assouline (2001), which is a relatively simple formulation for the HCF and serves as the reference HCF for the model simulations.
where for the present application 0 < m < 1, and n > 1. Zhang (2011) and is expressed, in the present notation, as where (Sect.2.2.1) σ w (T K ) is the surface tension of water and w (T K ) is the static dielectric constant or relative permittivity of water; 0 = 8.85 × 10 −12 C 2 J −1 m −1 is the permittivity of free space; k B = 1.308568 × 10 −23 J K −1 is the Boltzmann constant; a = 1.6021773×10 −19 C is the electron charge; z [dimensionless] is the ion charge, for which z = 1 can be assumed; and F w [dimensionless] is a soil-specific parameter, for which Zhang (2011) found that (roughly) 10 < F w < 10 4 .The term ρ w V θ,surf in Eq. ( 2) represents the soil moisture movement caused by water molecules "hopping" or "skipping" along the surface of the water films due to a temperature gradient (e.g., Medve ď and Černý, 2011).The present model for V θ,surf is adapted from the model of Gawin et al. (1999) and is given as where is the surface diffusivity and is parameterized as where D θs0 ≈ 10 −10 m 2 s −1 ; θ b ≈ 0.02; and β = 1/4 when θ ≥ θ b or otherwise β = 1 when θ < θ b .By expressing V θ,surf in terms of the gradient of the "normalized" soil moisture potential, ψ n , in Eq. ( 20), K * n and K m , used in Eqs. ( 5) and ( 6), can be identified as K * n = K n + D θ s D θψ and K m = D θ s D θT , where D θs D θψ will be defined as K surf n , where K surf n ≤ 0 ∀ T K and θ .

Numerical implementation
The numerical model as outlined above and detailed in this section is coded as MATLAB (The MathWorks Inc., Natick, MA, Version R2013b) script files.

Crank-Nicolson method
The linearized Crank-Nicolson method is used to solve Eqs. ( 5), ( 6), and (7).For Eq. ( 5) this yields the following (canonical) linear finite difference equation T ψi are the linearized C-N coefficients, which will not be explicitly listed here, but they do largely follow conventions and notation similar to Massman (2012).Although containing more terms than Eq. ( 22), the finite difference equation corresponding to Eq. ( 6) is very similar.But to linearize Eq. ( 7), the Crank-Nicolson scheme requires linearizing the source term, S v (T , θ, ψ, ρ v ).This is done with a first-order Taylor series expansion of the C-N term S j +1 v as follows: ∂ψ n , which in turn yields the following linearized finite difference equation for Eq.(7):  7) and t [s] is the time step.Here t = 1.2 s and was chosen after testing the model at t = 0.3 s and t = 0.6 s to ensure no degradation in model performance or solution stability at the larger time step.

Upper-boundary conditions
The upper-boundary condition on heat and vapor transfer is formulated in terms of the surface energy balance and, except for the latent heat flux, is identical to the upper-boundary condition in Massman (2012).
where the "0" subscript refers to the surface and the terms from left to right are the incoming or down welling radiant energy, , absorbed by the surface, which is partitioned into the four terms (fluxes) on the right side of the equation, the infrared radiation lost by the surface, the surface sensible or convective heat, the surface latent heat, and the surface soil heat flux.Q ⇓ R (t) and T amb (t) are functions of time and are prescribed externally as discussed in Massman (2012).The soil surface emissivity, (θ 0 ) [dimensionless], is a function of soil moisture and is taken from Massman (2012), as is surface heat transfer coefficient C H [m s −1 ], and σ is the Stefan-Boltzmann constant.
The surface evaporation rate, E 0 [kg m −2 s −1 ], is parameterized as www.geosci-model-dev.net/8/3659/2015/Geosci.Model Dev., 8, 3659-3680, 2015 where h s0 ≡ a w0 = exp([M w ψ * ψ n0 ]/[RT K0 ]) [dimensionless] is the "surface humidity", here modeled as the water activity at the surface using the Kelvin equation; C E [m s −1 ] is the surface the transfer coefficient, an adjustable model parameter but one that can reasonably be assumed to be between about 10 −4 m s −1 (Jacobs and Verhoef, 1997) and 10 −3 m s −1 (Massman, 2012).Finally, in the case of the laboratory experiments of Campbell et al. (1995), ρ v amb (t), like T amb (t) and Q ⇓ R (t), is an external forcing function at the soil surface.The present formulation of E 0 results from combining and adapting the expressions for the potential evaporation rate for soils developed by Jacobs and Verhoef (1997) and Eq.(9.14) of Campbell (1985).For this formulation the surface relative humidity, h s0 , is the surface property that constrains or reduces the surface evaporation E 0 to less than the potential rate.
The upper-boundary condition on soil water is (∂θ/∂z) 0 = 0, which when employed with the WRC, Eq. ( 15), yields the following upper-boundary condition on the conservation of soil moisture, Eq. ( 6): The boundary forcing functions e v amb (t) [Pa] (the ambient vapor pressure), T amb (t), and Q ⇓ R (t) are taken from Massman (2012), which in turn were adapted to the laboratory data of Campbell et al. (1995).They take the following generic form: where V i is the value of the function at the beginning of the soil heating experiment, V f is the value of the function at the end of the experiment, and τ [s] is a time constant of the heating source, which varies with each individual soil heating experiment.ρ v amb (t) is obtained from e v amb (t) and T amb (t) using the ideal gas law.

Lower-boundary conditions and initial conditions
As with the companion model (Massman, 2012), a numerical (or extrapolative or "pass-through") lower-boundary condition (Thomas, 1995) is also used for the present model.Analytically this is equivalent to assuming that the second spatial derivative (∂ 2 /∂z 2 ) of all model variables is zero at the lower boundary.It is used here for the same reason as with the previous model: principally for convenience because it is likely to be nearly impossible to specify any other the lower-boundary condition during a real fire.The boundary condition on the advective velocity is u vl = 0 at the bottom boundary, which is also the same as with Massman (2012) and Campbell et al. (1995).Further discussion on the model's lower-boundary conditions can be found in Massman (2012).Except for the initial value of ψ in (or ψ n, in ), all initial conditions (soil temperature and moisture content), which are assumed to be uniform throughout the soil column for each soil type and heating experiment, are taken directly from Campbell et al. (1995).The initial value for ψ is obtained by inverting (solving for it using) the WRC after inputing the initial values for soil temperature and moisture content.Consequently, ψ n, in can vary with the specific WRC.

Recalibration of observed volumetric soil moisture
In the original soil heating experiments of Campbell et al. (1995) soil temperatures were measured with copperconstantan thermocouples at the sample surface and at 5, 15, 25, 35, 65, and 95 mm depth and changes in soil moisture were obtained by gamma ray attenuation at the same depths (except the surface).The moisture detecting system was linearly calibrated for each experimental run between (a) the initial soil moisture amounts, which were determined gravimetrically beforehand, and (b) the point at which the sample was oven-dried (also determined before the heating experiment) where θ = 0 is assumed.But oven-drying a soil will not necessarily remove all the liquid water from a soil, i.e., a soil can display residual water content, θ r , after ovendrying.Consequently, the soil moisture data obtained and reported by Campbell et al. (1995) show negative soil moisture at the time the soil dryness passes outside the oven-dry range.Massman (2012) commented on this issue.With the present study, all volumetric soil moisture data were first adjusted (using a linear transformation) to rescale the observed soil moisture, θ observed , so that the values of θ observed < 0 became θ observed ≈ 0. This rescaling had very little impact on any values of θ observed except those asymptotic data where θ observed < 0. Furthermore, this recalibration is reasonable so long as the original calibration was linear and based on a Beer's Law type extinction coefficient (which would be linearly related to the logarithm of the attenuation of gamma ray intensity).

Model performance
The present model is evaluated against five of the soil heating experiments from Campbell et al. (1995): (1) Quincy Sand, which has an initial soil moisture content = θ in = 0.14 m 3 m −3 ; (2) Dry Quincy Sand with θ in = 0.03 m 3 m −3 , (3) Dry Palousse B with θ in = 0.07 m 3 m −3 , (4) Wet Palousse B with θ in = 0.17 m 3 m −3 , and (5) Wet Bouldercreek with θ in = 0.22 m 3 m −3 .But here the focus is principally on Quincy Sand and Wet Palousse B. Most of the major conclusions regarding the present model can be drawn from these two experiments and including Quincy Sand here also benefits comparisons with Massman (2012), who also tested his equilibrium model against the Quincy Sand data.But in addition to a general assessment of model performance, these two experiments also serve as vehicles for a sensitivity analysis of the key physical processes and parameterizations dis- The solid lines are for a model simulation with R p = 1 mm; the dotted lines corresponds to a simulation for R p = 4 mm.Note as R p increases the infrared portion of the soil thermal conductivity, λ [2]  s , also increases, in accordance with Eq. ( 14).To compare with the equilibrium model, see Fig. 2 of Massman (2012).cussed in the previous two sections.Thus, the Quincy Sand experiment is also used to explore the importance of the infrared component of the soil thermal conductivity (λ [2]  s : Eq. 14) and the Wet Palousse B experiment is also used to evaluate the potential contribution of the dynamic residual soil moisture (θ r (ψ n , T K ): Eqs. 15, 16 and 17) to the model's performance.Finally, because the HCF (Eq.17), is central to the overall performance of the model it is discussed in more detail a separate section following the Quincy Sand results.

Quincy Sand
Figure 1 compares the measured (symbols) and modeled (lines) of soil temperature during the Quincy Sand heating experiment for two model runs with different values of R p .Increasing R p to 4 mm (dashed lines) over the default value of 1 mm (solid lines) increases the infrared component of the soil thermal conductivity, λ [2]  s .Because the model's performance was not enhanced by the inclusion of either the residual soil moisture, θ r , or film flow, K film R , neither are included as part of the Quincy Sand simulations.The colors indicate the depths (mm) of the experimental and model data.(Note that the same color is also used to denote to same depth for both the model and observed data in Figs. 2, 3, and 4 below.)The corresponding measured and modeled soil moisture is shown in Fig. 2.These two figures indicate that the present model produces results that are similar to both the original Campbell et al. (1995) model and the observations.But they also indicate that the R p = 4 mm simulation is supe- rior to R p = 1 mm.Furthermore, comparing these two figures with their counterparts in Massman (2012) clearly indicates that the non-equilibrium model is a substantial improvement over the (older) equilibrium model, a conclusion that is easily confirmed by comparing Figs. 4 through 7 below with their equivalents in Massman (2012).Figure 3 is a plot of the data trajectory (observed soil temperatures vs observed soil moisture for all the monitored depths).The model's solution trajectories (for the same depths and the two R p values) are shown in Fig. 4. Comparing these two figures suggest that the model does a reasonable job of capturing the rapid vaporization of soil water at temperatures between 70 and 90 • C (at least at the depths below about 10 mm); furthermore, the present configuration of the model does not do as well at capturing the amplitude (amount) of the recondensing moisture ahead of the rapid drying nor the duration (or width of the amplitude) of the recondensation.Furthermore, when compared with observations (Fig. 3) the model does not fully capture the amount of unevaporated soil moisture that remains at temperatures ≥ 150 C. In this regard neither version of the model is significantly different from the other, but some of this "lack of precision" with the soil moisture simulation results in part by how the HCF was calibrated and will be examined in more detail in the next section.
Figure 5 compares the vertical profiles of the soil temperatures at the end of the laboratory experiment with those at the end of the two numerical simulations and Fig. 6 makes a similar comparison for the volumetric soil moisture content.These figures also include the modeling results synchronized in time and space with the observations, which are in- cluded to make the model output more directly comparable to the observations.(Note that the final vertical profiles obtained from the laboratory experiment are not coincident in time with the measurements made at any other depth.This is a consequence of the experimental design, which required several minutes to complete one vertical scan for soil moisture.)Figure 5 clearly indicates that the R p = 4 mm simulation does a better job of capturing the final temperature profile than does the R p = 1 mm simulation.The difference between the these two model simulations is less obvious with Fig. 6.The final soil moisture profiles shown in this figure can be used to estimate the percentage of soil moisture evaporated and lost from the soil column at the end of the 90 min experiment, E loss .The laboratory observations suggests 31 % was lost.The (R p = 1 mm; red) model simulation indicates a 31.4 % loss and the corresponding (red) synchronized model yielded a 33.8 % loss.The (R p = 4 mm; blue) model simulation indicates a 34.6 % loss and the corresponding (blue) synchronized model yielded a 34.2 % loss.Because the fully sampled and synchronized model results give somewhat different percentage loss it is possible to conclude that the laboratory estimate of evaporative loss is likely imprecise because it is poorly resolved in time and space.So exact agreement between model and observations is, in general, unlikely.But it is possible to use the model itself to estimate the uncertainty in the observed moisture loss that is caused by this undersampling.This is done by comparing and synthesizing all fully sampled model estimates of E loss with all the corresponding synchronized model estimates.For the five experiments studied here this uncertainty, expressed as a fraction, is about ±0.03.
For the present Quincy Sand experiment this yields a fractional E loss = 0.31 ± 0.03, which in turn suggests that both the present simulations (R p = 1 and 4 mm) provide quite good estimates of the total water loss as well as the final profiles of soil moisture and may in essence really be indistinguishable.But combining these soil moisture results with the final temperature profiles suggest that the R p = 4 mm simulation is the preferred.In addition, the present model results are significantly better than the equilibrium model, which found that no water was lost during the experiment, a clearly implausible result.(Rather than actually transporting the evaporated water out of the soil column, the equilibrium model "pushed" the moisture deeper into the soil ahead of the evaporative front as discussed in the Introduction.)On the other hand, despite the fact that the present estimates of evaporative loss are clearly a major improvement over the equilibrium results, both the equilibrium and nonequilibrium model solutions produce a sharply delineated advancing drying front, which is reminiscent of a Stefan-like or moving-boundary condition problem (e.g., see Whitaker and Chou, 1983, 1984, or Liu et al., 2005).So neither simulation actually captures the final observed moisture profile behind the drying front.
Figure 7 shows the final (90 min) modeled profiles of (a) the (R p = 1 mm) soil vapor density (ρ v (z)), equilibrium vapor density (ρ ve (z)), and the condensation term (K c (z)ρ v (z), used with the non-equilibrium model source term, S v : Eqs. 10 and 11) and (b) the (R p = 4 mm) soil va- The maximum vapor density for these two simulations is between about 1.3 and 1.5 times the density of the standard atmosphere (i.e., 1.292 kg m −3 ) and is located near the position of the maximum in the vapor source term, S v .This figure can be compared with the equilibrium model result: Fig. 8 of Massman (2012).
por density [ρ v (z)].The solid lines are model simulations with R p = 1 mm; the dashed red line corresponds to the R p = 4 mm simulation.(For the sake of simplicity only one curve is shown for R p = 4 mm simulation.)The maximum soil vapor density occurs at about 40 mm where the evaporative source term is greatest, i.e., where ρ ve (z)−K c (z)ρ v (z) is maximal, and where the moisture gradient is steepest, which is just ahead of the drying front (Fig. 6).Furthermore, the ρ v profile suggests that there are both upward and downward diffusional fluxes of vapor away from the maximal evaporative source.The upward-directed flux escapes through the soil surface and into the ambient environment of the laboratory (the surface evaporative flux) and the downwarddirected flux eventually recondenses below of the dry front.
The equilibrium model, on the other hand, produced virtually no vapor gradient within the dry zone thereby contributing to the model's inability to allow any moisture to escape (evaporate) from the modeling domain.Unfortunately, there are no observations with which to check either models' predictions of vapor density.But, although both models' results tend to agree that within the dry zone where temperatures exceed the critical temperature for water (i.e., 373.95 • C) there should be a single phase fluid that is significantly denser than water vapor near standard temperature and pressure (e.g., Pakala and Plumb, 2012), the non-equilibrium model does predict a more realistic vapor gradient than the non-equilibrium model.The solid line is the model simulation with R p = 1 mm and the dotted line corresponds to the simulation with R p = 4 mm.In both cases the maximum vapor pressure occurs at the soil surface and/or near the level of the maximum S v .For these two scenarios the maximum vapor pressure is about 3.2 times the pressure of 1 standard atmosphere (i.e., P ST = 101.325kPa).
If there is an implausibility with the present model it might be the soil vapor pressure, e v , as shown in Fig. 8.With either model simulation e v at the top of the soil column is between about 3 standard atmospheres (≈ 300 kPa).This is a bit unexpected because pressure at the open end of the column might be expected (at least by this author) to be close to equilibrium with the ambient pressure (≈ 92 kPa).Although there are no data against which to check this result, there are other modeling results that lend some support to the present predictions for e v .First, the steady state model of a sand-water-steam system from (Fig. 5 of) Udell (1983) heated from above indicates that the environment within the modeling domain is likely to be supersaturated and that at a minimum e v is greater than P atmos by ≈ 5 %, but (depending on the algorithmic treatment of the saturation vapor pressure and the exact value of P atmos he used for his simulations) it is also plausible to expect that e v ≈ (2-5)P ST .(Note that for the simulations in Udell (1983) the maximum model temperature was about 180 • C and that he also modeled advective velocity using Darcy's law.)Second, two different models of heated cement (Dayan, 1982, andDal Pont et al., 2011) indicate that near the top surface of the model domain e v can display values of ≈ (2-15)P ST .The overall similarities between these three earlier models and the present non-equilibrium model make it impossible to completely invalidate the present model's predictions for e v .Furthermore, the non-equilibrium model imposes no particular constraint on e v -it is calculated using the ideal gas law and the profiles of vapor density and temperature, both of which appear plausible.Consequently, the somewhat surprising result shown in Fig. 8 appears to be a natural consequence of the physics underlying the basic model equations: the conservation of mass and thermal energy.

HCF -Quincy Sand
Figure 9 shows the hydraulic functions and |K surf n | as functions of soil moisture for the model simulation with R p = 1 mm.(Recall that the components of the hydraulic diffusivity are all negative, so this figure reflects only their absolute value, not their sign.)K film H (assuming F w = 10 3 , Eq. 19) was calculated after the model run using the model's solution for T K , θ , and ψ.Because K film H did not actually contribute anything to any of the model runs for any of the five soil heating experiments (even for F w = 10 4 ), the present approach for evaluating it is sufficient.Consequently although |K n | shown here is more properly termed |K cap n |, this distinction is rendered moot for the present study.The values for m and n of K cap R , Eq. ( 17), were m = 0.26 and n = 1.80 and were obtained by subjectively optimizing the model (with R p = 1 mm) to fit the data shown in Figs. 1, 2, 5, and 6.In other words, it is possible to improve the model's performance for either the temperature data or the soil moisture data individually, but not simultaneously.Any improvement in one set of observations (T K or θ ) comes at the expense of the model's performance for the other variable.As a consequence of this optimization the model's ability to fully capture the amplitude (amount) of soil moisture that evaporates and recondenses ahead of the drying front (Figs. 3 and  4) has been sacrificed to improve the simulation of the soil temperatures.So the present numerical solution is a compromise between trying to fit two sets of data with a single "best" parameterization of K cap R from Assouline (2001).Similar compromises were required for other heating experiments as well.These particular optimal values for m and n were also found valid for the Bouldercreek soil experiments; but for the Palousse B soil the optimal values were m = 0.29 and n = 1.82.
Unfortunately, there are no independent confirmations for any values of m and n because no soil hydraulic conductivity data were (or ever have been) obtained for any of the soil samples used by Campbell et al. (1995).Nor are any data likely at this point in time because the original samples were destroyed years ago (G. S. Campbell, personal communication, 2015).At the time of these laboratory experiments the basic assumption that the original researchers made was that the heating and evaporation rates would be so fast as to preclude any (presumably much slower) liquid water transport associated with gradients in soil water potential.The present results appear to invalidate that assumption.
Although it is undeniably true that the present model is an improvement over the equilibrium model, the inclusion of the HCF within this non-equilibrium model (and its lack of inclusion in the equilibrium model) makes it difficult to conclude unambiguously that the improvement over equilib-  17), and the model for film flow from Zhang (2011), Eq. ( 18).Numerically rium model is the sole consequence of the non-equilibrium assumption.But the non-equilibrium model was tested in a mode that basically "turned off" K R and reduced K surf n by several orders of magnitude and still it yielded a solution (not shown) that simulated the soil temperatures and soil moisture observations better than did the equilibrium model.Furthermore, the non-equilibrium assumption will always remain superior to the equilibrium assumption for dry or extremely dry soils.Nonetheless, it remains unknown how much, if any, improvement in the equilibrium model's performance is possible with the inclusion of a HCF.
Finally it should also be pointed out that, unlike K film R , K surf n does contribute positively and significantly to the present model's performance because without it the model can become unstable for very dry soils.Therefore, K surf n is a significant factor in the HCF's overall contribution to the performance of the non-equilibrium model.

Wet Palousse B
Figure 10 compares the measured (symbols) and modeled (lines) of soil temperatures during the Wet Palousse B heating experiment with the model run that includes the residual soil moisture, θ r (dashed), and that which does not (solid).The solid lines are for a model simulation that does not include the dynamic residual soil moisture, θ r ; the dotted lines correspond to a simulation that includes θ r .For this experiment the initial soil moisture, θ in , is 0.17 m 3 m −3 .much, but that the soil moisture dynamics are much better portrayed by the model with θ r than without it.In fact, the main difference between these two simulations is that the soil moisture that evaporates and recondenses ahead of the drying front is much less for the model that includes θ r than that which does not (Fig. 12).The significance of this aspect of the model's performance is demonstrated in Fig. 13,  The solid lines are for a model simulation that does not include the dynamic residual soil moisture, θ r ; the dotted lines correspond to a simulation that includes θ r .This is the solution space representation of the model's solutions.most exactly in agreement with the observed value of 28.8 % and the modeled moisture profile below the drying front is also in almost perfect agreement with the observed data.The model without θ r allows for much more evaporated moisture to diffuse downward and to recondense ahead of the drying front than does the θ r model; thereby underestimating the water loss by about half and significantly overestimating the amount of soil water in lower portion of the profile.

Further sensitivity analyses
The Quincy Sand and Palousse B results in general confirm that the non-equilibrium model's performance is enhanced (and quite significantly) with the incorporation of liquid water transport (HCF) and that its performance is sensitive to (and can be improved by) either or both (a) the infrared thermal conductivity, λ [2]  s , through the volumetric pore radius, R p , and (b) the dynamical residual soil moisture, θ r (ψ, T K ).Although these last two aspects of the present model may not apply equally to any given soil, there seems little doubt that they should be considered as potentially quite important for modeling soil moisture and temperature dynamics during heating of soil during fires.
The remainder of this section is a summary of various (secondary, but important) model sensitivity analyses performed with all soil heating experiments.The ultimate intent here is to shed light on which physical process are relatively more important and to provide some guidance for further research.

The source term, thermal conductivity, and surface evaporation
Central to the success of the present model, relative to the performance of the equilibrium model of Massman (2012), is the functional parameterization of the source term, S v , and the related condensation coefficient, K c (T K , ψ n ).Basically K c was required to maintain model stability especially at high temperatures; without it the model was unstable and the dynamic between moisture and vapor was non-physical.
Regarding K c , the model is weakly sensitive to the choice of the surface evaporation/condensation activation energy, E av , providing it does not vary much outside the range of 30 kJ mol −1 ≤ E av ≤ 40 kJ mol −1 .On the other hand, from a systems perspective it is very difficult to infer much about the details of the (high-temperature) physical processes associated with K c or of the generality/universality of E av , other than their apparent existence and utility to the present model.The most effective value for the scaling parameter, S [M]   * , was within the range of about 0.5 to 1.The Novak model of the source term, S [N ]  v , also required the same K c , but the additional temperature dependency of S [N ]   v over S [M]   v forced the soil moisture to evaporate at slightly lower temperatures (therefore sooner) than shown in Fig. 4 for the Quincy Sand S [M]  v -configured model.S [N ]   v also eliminated an initial transient/instability (not shown) that occurred with the Quincy Sand S [M]   v solution.Otherwise, the differences between S [N ]   v and S [M]   v were not significant.It should not be surprising that the model is sensitive to soil thermal conductivity, λ s ; but it was somewhat surprising that the model is as sensitive to λ [2]  s [R p ] as it is.Both Durany et al. (2010) and Massman (2102) found that λ [2]  s only contributed for relatively porous soils, i.e., R p > 10 −3 m.In the present study the model temperatures could often be "fine tuned" (improved) even for relatively small values of R p .Therefore, in general, it seems unwise to ignore λ [2]  s , at least when modeling soil heating during fires.
The most important parameter controlling surface evaporation rate is the surface transfer coefficient C E , to which the model is reasonably sensitive.In particular (and similar to the results from Massman, 2012), the best values of C E were universally about 10 −3 m s −1 and values much above this caused the model to become unstable.Values well below these values (and closer to the theoretical value of 10 −4 m s −1 ) did not produce results much different than those resulting from C E = 10 −3 m s −1 .Nevertheless, C E does play a weak role in determining the soil surface temperature, T K0 , and therefore can influence the magnitude of the surface convective heat flux -Eq.( 24) -and λ [2]  s0 (T K0 ), which can influence the soil temperatures to some depth below the soil surface.

Water retention curves and hydraulic conductivity functions
The two other WRCs tested for model performance were Groenevelt and Grant (2004) (GG04) and Fredlund and Xing (1994) (FX94).But prior to implementing them in the model they were both calibrated to be numerically similar near the dry end (θ ≤ 0.03) of Eq. ( 15).Their performance was initially tested with the Assouline (2001) HCF, Eq. ( 17).
In addition, and as listed in Sect.2.2.5, besides Eq. ( 17) four other HCFs were also tested.Although not all pairings of WRCs and HCFs were tested against every heating experiment, the following conclusions seem relatively robust: (1) once calibrated to the present data set, the Brooks and Corey (1964) HCF performed at least as well as Eq. ( 17), but with only one parameter rather than two, while the other three HCFs gave somewhat less satisfying results and sometimes would even produce an instability; (2) the FX94 WRC often produced an instability, but its performance was also somewhat dependent upon which HCF was used with it; and (3) the GG04 WRC and associated HCFs did perform better the FX94, but overall, it did not perform as well as the present model configuration with Eq. ( 15) for the WRC and Eq. ( 17) for the HCF.In general the only guidance offered here is that some care must be given to choice of WRCs and

W. J. Massman: The HMV model version 1
HCFs because the modeling results can be quite dependent upon the choices made.Universal to all of HCFs tested here is K surf n (i.e., V θ,surf and its scaling parameter D θ s0 ; see Eqs. 20 and 21), which as explained above was incorporated into the hydraulic diffusivity, K * n (see Eq. 6 and the related discussion).The model's performance was relatively insensitive to the exact value of D θ s0 and D θ s0 ≈ 10 −10 m 2 s −1 is as good a default value as any other value.Furthermore, the inclusion of V θ,surf (soil moisture movement associated with water molecules "hopping or skipping" along the soil and water surfaces) did provide model stability when the soil moisture reached extremely low values.Including film flow, K film r -Eqs.( 18) and ( 19), brought no discernible benefit to the model's performance.For less extreme conditions and for relatively coarsegrained sand Smits et al. (2012) reached a similar conclusion.Nonetheless, K film r should not be discounted as unimportant, so further investigation into it's utility for modeling still seems warranted.

Different soils with different initial conditions
Of the remaining three heating experiments only Wet Bouldercreek, which had an initial soil saturation level of about 50 %, showed anything unexpected.In general, the model was able to capture the observed soil temperatures and temperature dynamics extremely well, even better than shown in either Fig. 5 (Quincy Sand) or Fig. 14 (Palousse B).But, regardless of any adjustments to any of the model parameters the model consistently underestimated the total amount of water evaporated (by about half), thereby also overestimating the amount of recondensing water ahead of the drying front (see the Wet Palousse B model run without θ r , Fig. 15, as an example).On the other hand, the model was able to capture the complete drying (θ ≡ 0) of the top 30 mm of the soil column during the Wet Bouldercreek experiment, whereas the Quincy Sand (Fig. 6) and Palousse B (Fig. 15) experiments indicated some residual soil moisture (0 < θ < 0.025) within the model's predicted dry zone.The model's performance was noticeably degraded when θ r was included and it was also quite sensitive to the choice of R p , so much so that best model fit to temperatures required that R p = 0, i.e., that λ [2]  s be excluded from the model.Nonetheless, the cause for this unexpected divergence between the modeling and observed soil moisture during this heating experiment is not understood.For the present it can only be surmised that the model's description of evaporative dynamics of soil moisture, and possibly the transport of water (both liquid and vapor), is still incomplete.

The advective velocity, u vl
Unlike with the companion model (Massman, 2012), the present model did not require reducing the magnitude of u vl in order to maintain model stability, which again rein-forces the impression that the present non-equilibrium model is an improvement over equilibrium model.Nonetheless, the present model can produce extraordinary gradients of vapor density and vapor pressure, which begs the question of whether such gradients could induce other types of mass transport than that captured by the present formulation of u vl , Eq. ( 9).This was tested by using a model for Darcy's Law type formulation based on the assumption that the advective velocity is proportional to the vapor pressure gradient (u vl ∝ −∂e v /∂z).This formulation was tested by incorporating it into the model (excluding θ r ).But the model became unstable because mathematically it was strongly hyperbolic, rather than predominantly parabolic.Further modeling development and parameterization of u vl and vapor transport in general are well beyond the intent of the present study.But it is still possible to conclude that such exploration is warranted and could help improve model performance.Finally, it should be noted that the model's performance was degraded when u vl was excluded from the model.

Different forms of Fick's first law for the diffusive flux
The point was made earlier that the present model was developed assuming the mass form for the diffusive flux, i.e., J [Mass]   diff = −D ve ∂ρ v /∂z, but that there are other forms that could have been used.Most notable among them is likely to be the mass fraction form for Fick's first Law: J [Fraction]   diff = −D ve (ρ v + ρ d )∂ω v /∂z.But J [Pressure]   diff = −D ve [M w /(RT K )]∂e v /∂z was also previously mentioned and was used by Campbell et al. (1995).Implementing either of these latter two forms of Fick's law requires amending J [Mass]   diff to include the influence of the temperature gradient on the diffusive flux.This is most easily accomplished with J [Pressure]   diff by combining it with the ideal gas law for water vapor, e v = ρ v RT K /M w , to yield A similar result is produced for J [Fraction]   diff by combining it with the ideal gas law for the dry air component of the soil pore space air, ρ d = P d RT K /M d , with the understanding that P d = P atmos is held constant because P atmos , the external ambient atmospheric pressure during the time of the experiment, is constant and not influenced by changes in vapor pressure within the soil column during the soil heating experiment.
The final expression for J [Fraction]   diff is Comparing Eqs. ( 27) and ( 28) shows that except for the (1 + ρ v /ρ d ) term in the denominator of J [Fraction]   diff these two expressions for J diff are the same.For the purposes of the sensitivity test both Eqs. ( 27) and ( 28) were used with ρ v as the predicted variable (i.e., ρ j +1 vi and ρ j vi are the Crank-Nicolson finite difference variables) and T as the model variable for the linearized coefficients (i.e., T j i ).The resulting finite difference terms were implemented into Eq.( 23) and into the finite difference equivalent of Eq. ( 6).
None of the model solutions resulting from either J [Fraction]   diff or J [Pressure]   diff were useful.They were all either (a) unstable or (b), if stable, physically unrealistic.In other words, no solution associated with either of these other two forms of J diff could be found that was not largely meaningless when compared to the observations.Consequently, the only possible conclusion here is that no improvement to the present model's performance is possible when using either J [Fraction]   diff or J [Pressure]   diff .This in turn supports the notion that for the physical problem considered in the present study J [Mass]   diff appears to be the preferred expression for J diff .

Summary and recommendations
This study has developed and tested a non-equilibrium (liquid to vapor phase change) model for simulating heat and moisture flow in soils during fires, but the model does assume thermal equilibrium.By and large the simulations of soil temperature and moisture are not only credible but also often quite good.In general, all model results showed a significant improvement over all comparable results from the companion equilibrium model of Massman (2012).
The principal reason for the present model's success is the incorporation of a dynamic condensation coefficient, K c (parameterized as a function of temperature and soil water potential), into the non-equilibrium evaporative source term, S v ; both of which are modeled after the Hertz-Knudsen equation.Physically K c suppressed condensation in favor of evaporation at high temperatures and soil water potentials, which in turn insured model stability.Furthermore, the non-equilibrium assumption also seemed to have improved the parameterization (and performance) of the mass transport associated with the advective velocity, u vl , relative to the model's of Massman (2012) and Campbell et al. (1995).The model's performance was further and significantly en-hanced by the inclusion of a hydraulic conductivity function (HCF) for liquid water transport, which was calibrated here by "fitting" the HCF parameters to ensure that the model optimally reproduced the observed temperature and moisture dynamics.This fitting procedure was necessary because no data are (nor will be) available for the soil samples used in the laboratory heating experiments (Campbell et al., 1995).Another important (and novel) feature of the model is the inclusion of a dynamic residual soil moisture θ r , also parameterized as a function of temperature and soil water potential, which is introduced into the model in an attempt to capture the long evaporative tail that seems to require temperatures well beyond 100 • C in order to evaporate at all.Physically θ r is intended to represent the strongly bound soil moisture, which for the present purposes is conceptualized as a monolayer.Including θ r was sometimes, but not always, beneficial to model performance.So it seems worthy of further consideration and possible refinement in any future studies of a similar nature.Finally, the model is also sensitive to the thermal infrared radiation component to the soil's thermal conductivity (λ [2]  s ; Eq. 14), which increases the thermal conductivity within the pore space of the soil as temperature increases.It is recommended that this term also be included and further tested when evaluating any future models describing the heating soils to high temperatures.
In general, the model simulates the observed soil temperatures quite well.It is often slightly less precise for soil moisture and the best simulations were usually a compromise between faithfully representing the observed soil temperatures or the observed soil moisture.Nonetheless, the model does capture reasonably well many observed features of the soil moisture dynamics, viz., it simulates an increase in soil moisture ahead of the drying front (due to the condensation of evaporated soil water at the front) and the hiatus in the soil temperature rise during the strongly evaporative stage of the soil drying.Furthermore, the model also captures the observed rapid evaporation of soil moisture that occurs at relatively low temperatures (50-90 • C), as well as some aspects of the long evaporative tail associated with strongly bound soil moisture.But, the model also displays a tendency to predict a greater depth of the drying front than suggested by the observations.
Sensitivity analyses (SAs) were also performed with different formulations for the water retention curve, soil hydraulic conductivity function, one variant of the present evaporative source term, S v , and different soil types with different initial conditions.The principal conclusion from these SAs is that some care (and testing) must be given to the selection of the WRC and HCF, as not all of them performed equally well.Some further investigations into the modeling benefit of film flow as part of the HCF also seems warranted.The two forms of S v tested here performed about the same, and the model's performance (at least for soil moisture) was poorest compared to the experiment with the highest initial soil moisture content.No obvious explanation for this "underperformance" could be found, so it seems worthwhile to further test the model for high initial saturation conditions.Finally, it is important to test the present model's performance and its associated parameterizations (particularly the WRC and HCF) against laboratory data and field data associated with daily cycles of soil heating and moisture transport.

Figure 1 .
Figure 1.Comparison of measured (symbols) and modeled (lines) soil temperatures during the Quincy Sand heating experiment.Neither simulation includes a dynamic residual soil moisture term, θ r .The solid lines are for a model simulation with R p = 1 mm; the dotted lines corresponds to a simulation for R p = 4 mm.Note as R p

Figure 2 .
Figure 2. of measured (symbols) and modeled (lines) soil moisture contents during the Quincy Sand heating experiment.Neither simulation includes a dynamic residual soil moisture term, θ r .The solid lines are for a model simulation with R p = 1 mm; the dotted lines correspond to a simulation with R p = 4 mm.To compare with the equilibrium model, see Fig.3ofMassman (2012).

Figure 3 .Figure 4 .
Figure 3. Measured soil moisture vs measured soil temperatures for the Quincy Sand heating experiment (see previous two figures).

Figure 5 .
Figure 5.Comparison of the final modeled and measured temperature profiles at the completion of the 90 min Quincy Sand heating experiment.Because the data shown in the measured profile (black) are not precisely coincident in time, the full model results (solid red and blue lines) were sub-sampled in synchrony in time (and coincide in space) with the observations.These time-synchronized model profiles are shown as dashed red and blue lines.To compare with the equilibrium model, see Fig. 6 of Massman (2012).

Figure 6 .
Figure6.Comparison of the final modeled and measured moisture profiles at the completion of the Quincy Sand heating experiment.Because the data shown in the measured profile (black) are not precisely coincident in time, the full model results (solid red and blue lines) were sub-sampled in synchrony in time (and coincide in space) with the observations.These synchronized model profiles are shown as dashed red and blue lines.The observed data (black) suggest that the total water lost during the 90 min experiment was 31 % of the initial amount.The (red) model simulation indicated a 31.4 % loss and the corresponding (red) synchronized model yielded a 33.8 % loss.The (blue) model simulation indicated a 34.6 % loss and the corresponding (blue) synchronized model yielded a 34.2 % loss.Note there is very little recondensing soil moisture ahead of the drying front (at about 40-50 mm depth), in agreement with Figs.2 and 4above and in contrast with the equilibrium model, Fig.7ofMassman (2012), where there was significant recondensation.

Figure 7 .
Figure 7. Final modeled profiles of vapor density [ρ v ], equilibrium vapor density [ρ ve ], and the condensation coefficient (K c ) modified vapor density term [K c ρ v ] used with the non-equilibrium model source term, S v , at the completion of the 90 min model simulation.The three solid lines are for a model simulation with R p = 1 mm; the single dotted line corresponds to a simulation with R p = 4 mm.The maximum vapor density for these two simulations is between about 1.3 and 1.5 times the density of the standard atmosphere (i.e., 1.292 kg m −3 ) and is located near the position of the maximum in the vapor source term, S v .This figure can be compared with the equilibrium model result: Fig.8ofMassman (2012).

Figure 8 .
Figure8.Final modeled profile of vapor pressure at the end of the 90 min model simulation.The solid line is the model simulation with R p = 1 mm and the dotted line corresponds to the simulation with R p = 4 mm.In both cases the maximum vapor pressure occurs at the soil surface and/or near the level of the maximum S v .For these two scenarios the maximum vapor pressure is about 3.2 times the pressure of 1 standard atmosphere (i.e., P ST = 101.325kPa).

Figure 9 .
Figure 9. Example of the hydraulic conductivity, K H , and magnitude of the hydraulic diffusivity, |K n |, as functions of soil moisture, θ , for the Quincy Sand R p = 1 mm scenario.K H corresponds to the HCF for capillary flow from Assouline (2001), Eq. (17), and the model for film flow fromZhang (2011), Eq. (18).Numerically

Figure 10 .Figure 11 .
Figure10compares the measured (symbols) and modeled (lines) of soil temperatures during the Wet Palousse B heating experiment with the model run that includes the residual soil moisture, θ r (dashed), and that which does not (solid).Figure11compares the measured (symbols) and modeled (lines) of soil moisture during the same experiment.Figure12shows the two models' solution trajectories, where the dashed line is the model run that includes θ r and the solid line does not.These results suggest that the inclusion of the θ r in the model does not influence temperatures very

Figure 12 .
Figure 12.Modeled soil moisture vs modeled soil temperatures for the Palousse B Wet heating experiment (see Figs. 10 and 11 above).The solid lines are for a model simulation that does not include the dynamic residual soil moisture, θ r ; the dotted lines correspond to a simulation that includes θ r .This is the solution space representation of the model's solutions.

Figure 13 .Figure 14 .
Figure 13.Observed and modeled soil moisture vs soil temperatures (trajectories) for the Wet Palousse B heating experiment.Solid lines are observed data and the dash-dot lines are from the model that includes the dynamic residual soil moisture, θ r .

Figure 15 .
Figure 15.Comparison of the final modeled and measured moisture profiles at the completion of the 70 min Palousse B Wet heating experiment.Because the data shown in the measured profile (black) are not precisely coincident in time, the full model results (solid red and blue lines) were sub-sampled in synchrony in time (and coincide in space) with the observations.These time-synchronized model profiles are shown as dashed red and blue lines.The observed data (black) suggest that the total water lost during the 70 min experiment was 28.8 % of the initial amount.The (red) model simulation indicated a 14.7 % loss and the corresponding (red) synchronized model yielded a 15.8 % loss.The (blue) model simulation indicated a 27.8 % loss and the corresponding (blue) synchronized model yielded a 29.4 % loss.Note there is very little recondensing soil moisture ahead of the drying front (at about 35 mm depth), in agreement with Figs.11 and 12.
is just a rescaling of K H (see Sect. 2.2.5 for further details) and |K surf n | is derived from the model for V θ,surf from Gawin et al. (1999) (again see Sect.2.2.5).