The challenges of an in situ validation of a non-equilibrium model of soil heat and moisture dynamics during ﬁres

. With the increasing frequency and severity of ﬁre there is an increasing desire to better manage fuels and minimize, as much as possible, the impacts of ﬁre on soils and other natural resources. Piling and/or burning slash is one method of managing fuels and reducing the risk and consequences of wildﬁre, but the repercussions to the soil, although very localized, can be signiﬁcant and often irreversible. In an effort to provide a tool to better understand the impact of ﬁre on soils, this study outlines the improvements to and the in-situ validation of a non-equilibrium model for simulating the coupled interactions and transport 5 of heat, moisture and water vapor during ﬁres. Improvements to the model eliminate two important (but heretofore universally overlooked) inconsistencies: one that describes the relationship between evaporation and condensation in the parameterization of the non-equilibrium vapor source term and the other, is the incorrect use of the apparent thermal conductivity in the soil heat ﬂow equation. The ﬁrst of these enhanced the stability and performance of the model. The second is an important improvement in the model’s physical realism, but had less of an impact on the model’s performance and stability than the ﬁrst. The model 10 validation uses (in-situ temperature, soil moisture, and heat ﬂux) data obtained in a 2004 experimental slash pile burn. Important temperature dependent corrections to the instruments used for measuring soil heat ﬂux and moisture are also discussed and assessed. Despite any possible ambiguities in the calibration the sensors or the simplicity of the parameterization of the surface heating function, the difﬁculties and complexities of formulating the upper boundary condition, and the obvious complexities of the dynamic response of the soil’s temperature and heat ﬂux, the model produced at least a very credible, if not surprisingly 15 good, simulation of the observed data. This study then continues with a discussion and sensitivity analysis of some important feedbacks (some

therefore include) an implicit and incorrect assumption about soil evaporation that is addressed and corrected in the present study.
The following section also discusses other changes to the HMV-model, including: (a) eliminating the use of the apparent soil thermal conductivity in the soil heat flow equation (also further discussed and justified in the Appendix) and (b) improving the 60 parameterization of the surface energy balance and the upper (soil surface) boundary condition (including the development of a generic soil heating function for use with prescribed burns or wildfires). In addition, this study also discusses the subtleties and difficulties of formulating a universal surface energy balance for soil heating by fire. The third section reviews the site, soils and data of the experimental slash pile burn and the fourth section compares the observations with the model simulations and explores some of the consequences (to the simulations) of some dynamic feedbacks and interactions between the fire and 65 soil physical properties. The fifth section discusses possible future directions for modeling and observational studies. The final section summarizes this study.

Model Description
Similar to version 1 of the HMV-model (Massman , 2015), the present version employs a linearized Crank-Nicolson finite difference scheme and was coded and run using MatLab version 2017b. This manuscript also uses the same notation and the 70 same functional parameterizations for the supporting thermodynamic and physical variables as Massman (2015). For this study details concerning these functional parameterizations will be summarized as necessary for clarity and to update with new information or data. Otherwise many details covered by Massman (2015) will not be repeated here. The remainder of this section discusses the physical fundamentals of the changes made to the HMV-model.

75
The HMV-Model is composed of the three conservation equations: the conservation of energy (or maybe more properly the conservation of enthalpy), the conservation of soil liquid water and the conservation of soil water vapor. The conservation of energy is where C s (Jm −3 K −1 ) is the volumetric specific heat of the soil, such that C s = C s (T, θ) is a function of both temperature 80 and volumetric soil moisture; t (s) is time; z (m) is soil depth; λ s (Wm −1 K −1 ) is soil thermal conductivity, such that λ s = λ s (T, θ, ρ v ); L v = L v (T K ) (Jkg −1 ) is the enthalpy of vaporization and −L v S v represents the change in enthalpy associated with evaporation/condensation; S v = S v (T K , θ, ψ, ρ v ) is the source term for water vapor and is discussed in more detail in the following section; and W S w is the change in enthalpy associated with the heat of wetting (also termed the heat of immersion), where W (Jkg −1 ) is the heat of wetting and S w (kgm −3 s −1 ) is the source term for water liquid or equivalently the sink term 85 for water vapor, i.e., S w ≡ −S v . W is discussed by de Vries (1958) and for the present purposes W can be interpreted as that additional enthalpy of vaporization that is required to break the electrostatic bonds between molecular water and the soil mineral surfaces. In general the wetting reaction is exothermic, i.e., W > 0 and a function of temperature (Grant (2003); Prunty and Bell (2005)). Massman (2012) investigated the effects of temperature on ψ and W , but found that it had little impact on the modeling results. For this study the HMV-model follows Campbell et al. (1995) and assumes that W = −ψ and 90 ignores any temperature dependencies of ψ and W . Note that W is only significant at high temperatures (as L v → 0) and for extremely dry soil (as −ψ → ∞). Finally with this identification for W and the above identity between S w and S v , it follows that L * v ≡ L v − ψ. The conservation of liquid water is where ρ w = ρ w (T K ) (kgm −3 ) is the density of liquid water; ψ n (dimensionless) is the non-dimensional form of ψ, i.e., ψ n = ψ/ψ * , where ψ * = −10 6 Jkg −1 is the nominal soil water potential of oven-dried soil ( (Campbell et al. , 1995)). (Note ψ n is used interchangeably with ψ throughout this manuscript.) K n = K n (T K , ψ n , θ) (m 2 s −1 ) is the hydraulic diffusivity; K H = K H (T K , ψ n , θ) (ms −1 ) is the hydraulic conductivity; and V θ,surf = V θ,surf (T K , θ) (ms −1 ) is the velocity of liquid water associated with surface diffusion of water. The hydraulic conductivity functions, K n (T K , ψ n , θ) and K H (T K , ψ n , θ), are 100 given as follows: where µ w = µ w (T K ) (Pas) is the viscosity of water; g = 9.81 ms −2 is the acceleration due to gravity; K I (m 2 ) is the intrinsic permeability of the soil here assumed to be constant and uniform throughout the soil profile, but does, in fact, vary with the concentration and type of solutes in soil water, (e.g., Lutz and Kemper , 1959); and K R = K R (θ, ψ n , T K ) (dimensionless) is 105 the relative hydraulic conductivity (used to describe capillary flow in soils). The model for intrinsic permeability is taken from Bear (1972) and is K I = (6.17 × 10 −4 )d 2 g ; where d g (m) is the mean or 'effective' soil particle diameter. Note that switching variables from ψ < 0, to ψ n produces ψ n > 0 and K n < 0. The present model considers only capillary flow and will ignore film flow, because like Massman (2015) film flow did not really impact the model's performance.

The conservation of water vapor is
where η (m 3 m −3 ) is the total soil porosity, assumed to be temporally constant and spatially uniform, and (η − θ) is the soil's air filled porosity; and D ve = D ve (T K , ψ, ρ v ) (m 2 s −1 ) is the (equivalent) molecular diffusivity associated with the diffusive transport of water vapor in the soil's air-filled pore space, where D ve includes the enhancement factor developed by Campbell associated with the rapid volitalization of soil moisture, which is given as follows: The final model equations result by preserving Equation (4) and eliminating S v from Equations (1) and (2), such that Equation (1) is replaced with 120 and Equation (2) is replaced with

Improvements in Non-equilibrium Vapor Source Term
Massman (2015, Equation (10)) adapted the Hertz-Knudsen Equation to develop the following formulation for the vapor source where S * is an empirical dimensionless parameter, which is "tuned" as necessary to ensure model stability; A wa (m 2 m −3 or m −1 ) is the volume-normalized soil water-air interfacial surface area, which Massman (2015) parameterized as A wa = A wa (θ); R = 8.314 Jmol −1 K −1 is the universal gas constant; M w = 0.01802 kgmol −1 is the molar mass of water vapor; K e (dimensionless) is the mass accommodation (or evaporation) coefficient, which Massman (2015) sets ≡ 1; K c = K c (T K , ψ n )

130
(dimensionless) is the thermal accommodation (or condensation) coefficient, which Massman (2015) parameterizes as a physicochemical (Arrhenius) function: is an empirical surface condensation/evaporation activation energy for which E av ≈ 30 − 40 kJmol −1 was determined empirically and T K,in is the initial soil temperature; and ρ v,eq (kgm −3 ) is the equilibrium vapor density, defined as ρ v,eq = a w ρ v,sat (T K ) where a w = e Mw ψ * RT K ψn is the dimensionless water activity, modeled here with the Kelvin Equation, and ρ v,sat (T K ) (kgm −3 ) is the 135 saturated vapor density, which is a function only of T K .
But the model of S v embodied by Equation (8) assumes that the interfacial surfaces appropriate to evaporation and condensation are the same, i.e., that A wa is the same for both evaporation RT K /M w K e ρ v,eq and condensation In general this is not a priori the case unless one assumes that soil moisture (θ) never drops below the point at which the soil's interfacial surface area is completely covered by a thin film or mono-layer of liquid water (e.g., Novak , 2019). But it is physically 140 more realistic, at least for very dry soils (which are likely to occur during fires), to assume that condensation can occur even in the absence of liquid water. Otherwise models of S v would impose a physically unrealistic constraint on non-equilibrium models of heat and moisture flow in dry soils.
The new version of the non-equilibrium model parameterizes S v as follows: where K e ≡ 1 has been retained as has the original formulation for A wa (Massman , 2015): where S w = θ/η is the soil water saturation and a 1 = 50 (rather than the original value of 40), a 2 = 0.003, and a 3 = 1/8.
This particular value for the parameter a 1 was chosen so that the maximum value of A wa occurs at S w ≈ 0.02 (= 1/a 1 ), and is assumed to be where the soil surfaces are covered by a mono-layer of water (Brusseau et al, , 2006). A wa,dry ≡ A wa (θ) as long 150 as S w > 1/a 1 , and A wa,dry ≡ max (A wa ) whenever S w ≤ 1/a 1 . In other words, A wa,dry differs from A wa whenever the soil moisture is so low that the soil particle surfaces are covered by, at most, a mono-layer of water. Empirical "tuning" of S * and E av after implementing the other changes yielded S * = 0.1 and E av = 10 kJmol −1 . Together these changes to S v improved the model's stability and robustness, as well as, it's fidelity to the observed soil moisture during the 2004 burn (detailed later).

155
The present model of λ s retains (a) the general structure of the original Campbell-de Vries model for thermal conductivity (see Campbell et al. (1994), de Vries (1963 and Massman (2012)) and (b) the additional Bauer term associated with the high-temperature thermal (infrared) radiant energy transfer within the soil pore space (Bauer , 1993). That is where k w , k a , and k m (dimensionless) are the Campbell et al. (1994) generalized formulations of the de Vries (1963) is the apparent thermal conductivity of moist air, which is the sum of the true thermal conductivity of moist air, λ a (θ, T K , ρ v ) (Wm −1 K −1 ) and the term, λ * v (θ, T K , ρ v ) (Wm −1 K −1 ), which embodies the effects of latent heat transfer or "the effect of the vapor distillation due to temperature gradients" (de Vries (1958); Appendix A of the present study); λ m (Wm −1 K −1 ) is the thermal conductivity of the mineral components of the soil; σ (Wm −2 K −4 ) is the Stefan-Boltzmann constant; N = N (θ) = 165 1 + θ/(3η) (dimensionless) is the pore gas index of refraction; and R p (m) is the soil's pore space volumetric radius. The first term on the right hand side of Equation (11) is the Campbell-de Vries model and the second is the Bauer term.
The weighting factors, k w , k a , and k m all have the same general form (Campbell et al. , 1994): where the subscript * refers to water (w), air (a), or mineral (m); g a is the de Vries (1963) shape factor, an empirically 170 determined model parameter. In general g a ≈ 0.1 (Campbell et al. , 1994). λ f = λ f (θ, T ) is a weighted mixture of the thermal conductivities of air and water: and 175 where q w (T K ) is a dimensionless parameter that describes the water content at which water begins to influence λ s . It is defined as q w (T K ) = q w0 (T K /303) 2 with q w0 another empirically determined parameter. g a and q w0 are important for the present study because they are an important part of the sensitivity analysis (discussed in a later section) that explores the interactions between the fire and the soil physical properties.
Three changes have been made to this original formulation. First, λ * v is no longer included because to do so is to double 180 count the vapor "distillation" (de Vries , 1958) term that accounts for the influence that evaporation, transport and condensation of water vapor can have on the apparent thermal conductivity (Appendix A). In other words as shown in Appendix A, both equilibrium and non-equilibrium models of soil heating include the vapor distillation term, either explicitly through the conservation of mass of water vapor (non-equilibrium models) or implicitly through the conservation of mass of liquid water (equilibrium models). Consequently, it is unnecessary and redundant to include λ * v in Equation (11). Second, λ m now in-185 cludes an explicit temperature dependency of the form λ m = λ m (T K ) = λ m0 (8 exp(−0.008(T K − 300)) + 3)/11; where λ m0 (Wm −1 K −1 ) is basically an adjustable parameter. The results for α-quartz (α referring to a specific crystalline structure of quartz) from (Yoon et al. , 2004, Figure 4) suggests that it is reasonable to expect that λ m0 ≤ 15 Wm −1 K −1 . The temperature function (8 exp(−0.008(T K − 300)) + 3)/11 was chosen to emulate the approximately 70% decrease in the thermal conductivity of α-quartz (here used as a substitute for sand quartz) between about 300 K and 600 K shown by Kanamori et al. (1968) 190 and Yoon et al. (2004). Third, R p is now estimated as R p = 0.408d g (ρ p /ρ b ) − 1 from Arya et al. (1999); where d g (m) is the mean or "effective" soil particle diameter, ρ p (Mgm −3 ) is the particle density (which can usually be assumed to be about 2.65 Mgm −3 and ρ b (Mgm −3 ) is the soil bulk density. This formulation for R p yields more physically realistic estimates of R p (i.e., 5 µm ≤ R p ≤ 200 µm for the soils tested in this study) than the default value of 1000 µm used in the previous version of the HMV-model. It also suggests that the thermal infrared contribution to the soil thermal conductivity is negligible in most 195 soils, even during fires.

Evaluation of changes to S v and λ s
Assessing the consequences of these alterations to the source term and soil thermal conductivity to the model's performance was done using the laboratory data of Campbell et al. (1995) and by comparing the simulations with the changes to S v and λ s to those shown in Massman (2015). The results indicated that (a) the model's ability to faithfully reproduce the data were 200 very similar to those in Massman (2015) and (b) the model's stability was significantly improved. The reason for (b) is (almost exclusively attributable to) the change in S v and is very much a positive benefit to the model. The changes to λ s , did require some adjustments to the values of some of the other parameters included in the model for λ s , but these were minimal and not particularly significant. For the sake of brevity, none of these comparisons are included in this study. The forcing function is the energy that is input to the soil at its surface, denoted here by Q F (t) (Wm −2 ). How that energy is divided between net infrared heat loss, convective heat loss, evaporation, and soil conductive heating is expressed by the surface energy balance and the upper boundary conditions. Although relatively simple in concept, in practice, for application to fires the forcing function, the surface energy balance and the upper boundary all are at best difficult to formulate precisely 210 and at worst potentially a fiction. For example, in the case of a wildfire or soil heating within a few meters outside the physical perimeter of a slash pile surface forcing is primarily radiant energy. See Fig.2 of Massman et al. (2010a) for an example of this second case. On the other hand, when burning material makes direct contact with the soil, it is reasonable to assume that the forcing at the soil surface is more likely to be conduction rather than radiant energy. Beneath a slash pile burn surface forcing may be combination of radiation and conduction, it may change over time as the pile burns, as the ash accumulates and, at 215 later stages of the burn, as the pile collapses. In the case of a moving fire front, the forcing can be highly variable. Radiant energy is clearly a major driver, but in addition, thermal instabilities drive circulations ahead and behind the fire that input energy into the soil when these circulations force hot air into contact with the soil, which in turn causes direct ignition of soil biomass ahead of the flame front (Finney et al. , 2015;Pearce et al. , 2019;Linn , 2019). As the fire front passes the forcing is likely to be a combination of conduction and radiation and possibly convection. Whereas after the fire front conduction is 220 the major forcing in areas covered with burning biomass, and radiant energy and possibly convection in areas free of burning biomass. Finally, in the case of burning duff, forcing is likely to be solely conductive in nature; but complications arise because in this situation duff is a highly porous burning insulator. Parameterizing the forcing in this case is problematic because of the extremely limited (empirical and theoretical) knowledge concerning burning duff. Massman (2015) assessed the performance of the previous version of the HMV-model against the laboratory data of Camp-225 bell et al. (1995), which more or less dictated the following forcing function: Q F (t) = Q F max (1 − exp(−t/τ )); where τ (s) and Q F max (Wm −2 ) are adjustable parameters that describe the rate of rise of the forcing function (τ ) and its stead state asymptotic (maximal) value Q F max . This study uses the following modified form of the BFD curve (Barnett , 2002) as the forcing function

Complexities and Challenges
where e −α(ln(t/tm)) 2 is the modified BFD-curve; t m (s) is the time at which the maximum forcing occurs, and α (dimensionless) = 2 ln(10)/(sinh −1 [0.5t d /t m ]) and t d (s) is the time interval between when the forcing first reaches 1% of its maximum and when it has decayed to 1% of its maximum. Q F max , t m and t d are adjustable input parameters to the model that define the forcing function, much the same way that Q F max and τ are for the Campbell et al. (1995) forcing function. Q F in (Wm −2 ), on the other hand, is not an adjustable input parameter. It is determined from other considerations of the soil surface energy 235 balance and is define later. The boundary conditions for temperature and vapor pressure, e v (Pa), have similar functional forms as Equation (15): where V a refers to the ambient atmospheric temperature (T a ) or vapor pressure (ρ v,a ) at the soil surface and the subscript 'in' refers to the initial value of that variable (taken from observations near the time and location of the fire). V max is an adjustable model input parameter 'tuned' so that the model matches (as much as possible) the soil observational data and (if necessary) to help ensure model stability. The vapor density boundary condition is 240 derived from these latter two boundary conditions using the ideal gas law. The upper boundary condition on soil moisture is (∂θ/∂z) 0 = 0 and is discussed further in Massman (2015).
The energy balance at the soil surface used with the present study is slightly different from either Massman (2012) or Massman (2015). Here it is expressed as where the '0' subscript refers to soil surface and the term on the left hand side of this equation is the energy absorbed by the soil (and assumes that absorptivity and emissivity of the soil are the same) and the first term on the right hand side is the net infrared heat loss, the second is the convective heat loss, the third is rate of evaporation, and the last is the soil conductive heat; 0 (θ 0 ) is the soil emissivity and is a function of the soil moisture, θ 0 ; σ = 5.67 × 10 −8 Wm −2 K −4 is the Stefan-Boltzmann constant; a (ρ va ) is the emissivity of the ambient atmosphere exposed to the soil surface during the fire and is a function of the 250 ambient vapor density ρ va following the 'clear sky' parameterization of Brutsaert (1984, Equation (6.18) is the mass density of the ambient air at the soil surface temperature, T K0 , where P a (Pa) is the ambient pressure (at the time and location of the fire and is a model input variable) and P ST = 101325 Pa and T ST = 273.15 K are the standard atmospheric pressure and temperature; C H = 0.032 (ms −1 ) is the transfer coefficient for convective heat from the surface (see Massman (2012) or Massman (2015)). T a = T a (t) (C), or equivalently T Ka = T Ka (t)

255
(K), is the ambient temperature somewhere above the soil surface (upper boundary condition); L * v0 E 0 (Wm −2 ) is the rate of soil water evaporation; E 0 (kgm −2 s −1 ) is the evaporative mass flux at the surface; and G 0 (Wm −2 ) is the soil conductive heat flux and the upper boundary condition for the modeled soil temperatures, Equation (6). E 0 is parameterized as the sum of a diffusional component and an advective component (Massman , 2012): where C E (ms −1 ) and C U (dimensionless) are adjustable model transfer coefficients, which were determined empirically to maximize E 0 without destabilizing the model. C U = 0.125 is associated with the jet of volitalized air emanating from the soil with velocity u vl0 and C E = 10 −4 ms −1 . The surface humidity, h s0 (dimensionless), is assumed to be the same as a w0 , the water activity at the surface.
Although Equation (16) provides a complete accounting of the energy exchange between the atmosphere and the soil, there 265 are two key issues that need to be addressed when initializing the model for applications to fires. The first was mentioned earlier, i.e., under a burning slash pile or a soil covered (partially or completely) with an ash layer it is not clear what the infrared and convective heat environments really are or what role they may play in the soil energy balance. This issue is relevant here because the data used in this study was obtained under a burning slash pile and it is addressed in the present study in a later section on sensitivity analysis by comparing model simulations with and without the IR and convective heat terms. Removing 270 these later two terms from Equation (16) yields the following simplified version of the soil surface energy balance: The second issue is basically a feature of all forcing functions. In the case of the BFD-curve, t m ≥ 0.5 hr is always true and because a typical model time step is between 1 and 4 s, there will always an initial period between a few minutes to several tens of minutes long where the simulated burning time t t m . During this "ramp-up" period Q F ≈ Q F in . The choice 275 of Q F in depends on whether Equation (16) or Equation (18) is used. In the case of Equation (18), Q F in = 0. In this case the soil conductive heat flux becomes G 0 = −L * v0 E 0 , and because the source term and the evaporation rates are very nearly at equilibrium (i.e., S v ≈ 0 and E 0 ≈ 0) during this ramp-up, it follows that G 0 ≈ 0 and that ∂T /∂t ≈ 0. (Note: the soil is initialized to be isothermal at the temperature obtained at the soil surface just before initiating the burn.) But for Equation (16), the full surface energy balance equation, this equilibrium condition does not occur during the ramp-up if Q F in = 0, because the 280 net IR term in Equation (16) is not initially in equilibrium. This is true despite initializing T a = T 0 . Without Q F in Equation (16) reduces to G 0 ≈ − 0 (θ 0 )σT 4 K0 (1 − a (ρ va )) < 0, which induces a transient in the solution that causes the soil temperature to drop slightly (i.e., ∂T /∂t < 0). But assigning Q F in ≈ σT 4 K0 (1− a (ρ va )) eliminates this transient and ensures that G 0 ≈ 0 and ∂T /∂t ≈ 0 during the ramp-up time. Several other methods were tested for eliminating this unrealistic solution, but the present approach proved to be the least intrusive and best way to prevent this initial transient. Note: this initial period of disequilibrium 285 in the surface energy balance is not unique to the modified BFD curve. It also occurred with other forcing functions that were tested (e.g., Blagojević and Pešić , 2011) as well as with the Campbell forcing function used in Massman (2012) and Massman (2015).

The Lower Boundary Condition and Initial Conditions
The lower boundary condition is the same "pass through" or "extrapolative" boundary condition that was used in Massman 290 (2012) and Massman (2015), i.e., the second derivative, ∂ 2 /∂z 2 , of the 3 model variables = 0. But for the present study (which is devoted to wildfires and slash pile burns, rather than the laboratory experiments of Campbell et al. (1995)), the lower boundary is placed at 0.60 m below the surface (well below 0.20 m used in these previous studies). This pass through boundary condition is used because for field-based applications the lower boundary condition will never be known (or knowable without extraordinary effort before hand) so it must be fairly general and placed at a depth where it will not influence the 295 model predictions within the upper few centimeters of soil too much. The advective velocity, u vl , requires only one boundary condition, see Equation (5), and is u vl = 0 at the bottom boundary.
The initial conditions for soil temperature and moisture are taken from measurements made just before igniting slash pile.
The soil temperature is assumed to be uniform throughout the vertical domain and is taken to be the observed soil surface temperature. The initial soil vapor pressure is estimated to be about 40% of the saturation vapor pressure at the initial soil 300 temperature. The ideal gas law is then used to estimate the initial soil vapor density profile. The initial soil water potential is obtained from the soil moisture data using the water retention curve, discussed next.

The Water Retention Curve
Unlike the previous studies (Massman (2012), Massman (2015)), which employed soils for which the water retention curve was unavailable, the present study does employ a water retention curve appropriate to the soils at the site of the experimental 305 burn (Manitou Experimental Forest: MEF). Figure 1 shows the WRCs for pre-burn (red) and post-burn (black) soils the MEF burn site. Note: both pre-and post-burn soils are included in this study because they are part of the model sensitivity analysis.
The data shown in this figure were provided by Butters (2009) and were obtained from a 2008 burn study performed at a site about 60 m away from the 2004 burn discussed in this study. The data were fit with the Fredlund and Xing (1994) model: where a > 0, b > 0, n > 0, and m > 0 are fitting parameters, e is Euler's Number, and the total porosity η was established beforehand such that η pre ≈ 0.51 and η post ≈ 0.45. These values of η correspond to the following values of soil bulk density: This change in bulk density is revisited in the model sensitivity analysis.
Although the Fredlund-Xing model provided the best fit to the data, other models of the WRC were fit to the observations.
For the HMV-model the choice of WRC is important. The Fredlund-Xing model did provide the best fit to the data and its 315 impact on model performance was judged the best of the models tested. Other models for the WRC that were tested include the Campbell-Shiozawa model (Campbell and Shiozawa , 1992), which was used in both Massman (2012) and Massman (2015), the van Genuchten model (van Genuchten , 1980), and the Groenevelt-Grant model (Groenevelt and Grant , 2004). For the purposes of model performance, the key difference between Fredlund-Xing model, Equation (19) and Campbell-Shiozawa description and mathematical behavior when the soil is extremely or completely dry. Under these conditions Equation (19) when θ ≈ 0 the soil water potential remains bounded, i.e., ψ ≈ −10 6 Jkg −1 . But with the van Genuchten and the Groenevelt-Grant models a completely dry soil is impossible to achieve because θ = 0 can only occur when ψ = −∞. But in the HMVmodel an unbounded function for the WRC does result in the degradation in model performance, some loss of physical realism in the simulation of soil moisture, and can introduce model instabilities.

325
On the other hand, using a WRC that remains bounded produces a logical inconsistency when used with the Kelvin Equation to describe water activity and the equilibrium vapor density. This issue plagued earlier modeling attempts (Massman , 2012(Massman , , 2015. In the case of Massman (2012), when θ = 0, ψ = −10 6 Jkg −1 and ∂ψ/∂t = 0. But as long as the temperature keeps rising, ∂T /∂t > 0, it follows that ∂ρ v,eq /∂t > 0. For the equilibrium model (ρ v ≡ ρ v,eq ) this means that when the soil moisture is completely evaporated and the soil is dry, the model autonomously creates water vapor. In the case of Massman (2015), 330 the same inconsistency produces a source term that does not allow condensation to occur (i.e., S v < 0) on the surface of a soil particle that is completely dry. It is this issue that lead to the improved parameterization of S v embodied in Equation (9) of the present study.

The Hydraulic Functions
The hydraulic functions, K n and K H , are both basically determined by K R , as discussed above (Equation (3)). Like Massman

335
(2015) the present study employs the Assouline model Assouline (2001) for K R , except here K R does not include a residual soil moisture term. Therefore, where m k and n k are parameters (0 < m k < 1, and n k > 1) that are adjusted or tuned to ensure that the present model produces a reasonable and physically realistic simulation of the observed soil moisture dynamics during the fire. Unfortunately, there are 340 no data available to determine K R for the soil at the MEF burn site.

Manitou Experimental Forest: Burn Site and Instrumentation
Much of the general description of climate and physical characteristics of Manitou Experimental Forest and the burn site have been published previously (Massman and Frank , 2004;Massman et al. , 2006). Nonetheless, for the present purposes they do bear repeating. Soils within the burn area are Pendant cobbly loam and range between 60-65% sand, 20-25% silt, and 10-15% clay with bulk densities that usually increase with depth (Massman et al. , 2008) and range between 1.1 and 1.5 Mgm −3 . Soil organic material comprises about 1-2% of the soil by volume. Previous grazing and mechanical harvesting throughout the area has resulted in a moderately disturbed soil.  Figure 2 shows the burning slash pile shortly after the fire was initiated and the deployment of the data loggers, CO 2 pumps and analyzers, and the supporting infrastructure. Note: although the soil CO 2 data 370 and the data from the edge of the pile are somewhat peripheral to this study, they are included here because they offer some important insights into some of the assumptions underlying the present model, as is discussed in a later section. Otherwise the measured soil temperatures, soil moisture and soil heat fluxes are compared to the model's predictions and thereby assess and validate the model's performance. were K-type (rated to 704 C), the J-type (rated to 260 C) was used at 0.15 m, and the bottom two depths were T-type (rated to 380 100 C).

Soil Heat Flux
The soil heat fluxes were measured at three depths (0.02 m, 010 m, 0.20 m) under the center of the pile with a high-temperature probe (HFT) with an aluma core and an exterior ceramic glaze (Thermonetics Corporation, La Jolla, CA). They were also sampled at the same rate and time as the soil thermocouples. These high temperature HFTs are rated to 775 C and have a 385 nominal sensitivity between 1250 and 1750 Wm −2 mV −1 . These HFTs were attached to a data logger by a cromel extension wire (Omega Engineering TFCH-020, rated to 260 C).
Because these high temperature HFTs are exposed to such a wide range of temperatures (potentially anywhere between about -10 and 700 C) it is important to account for the effects of temperature on the sensors' thermal conductivity and calibration factors. Details concerning the calibration factor are discussed by Massman and Frank (2004) and need not be repeated here.

390
But the details concerning the thermal conductivity of these HFTs are very relevant here and need highlighting. The thermal conductivity of these HFTs, λ p (Wm −1 K −1 ), is λ p = 0.7+0.003T (where T is degrees Celsius). Knowledge of λ p is important to correct for the discrepancy between the true soil heat flux, G soil (Wm −2 ), and the measured soil heat flux, G m , that results whenever λ p differs from the soil's thermal conductivity, λ s , (Philip , 1961;Sauer et al. , 2003;Tong et al. , 2019). This relationship, known as Philip's correction, is given as where β is a dimensionless factor related to sensor shape (β = 1.31 for the square HFTs), r is the sensor's aspect ratio, i.e., the ratio of the sensor's thickness to its horizontal length (r = 0.19 for the HFTs), and −1 = λ s /λ p . If a sensor is perfectly matched to its soil environment, then λ p = λ s , G m = G soil , and there would be no need to correct for this discrepancy. But, in general, this is unlikely to occur very often at normal daytime or nighttime soil temperatures and so is, therefore, even less 400 likely during a fire. But this correction also requires in situ knowledge of λ s , which was not (and probably could never have been) measured during the fire. So the model's predicted λ s is used in Equation (21). As a consequence, the model's predicted G soil will be evaluated against the measured heat flux, G m , with and without the Philip correction.

Soil Moisture and CO 2
All soil moisture and CO 2 data were measured at 0.05 and 0.15 m depths at both the center and edge position under the slash 405 pile. They were both sampled every half hour for the week during and after the burn. Soil moisture was measured using a specially designed high-temperature TDR (Zostrich Geotechnical; Pullman, WA, USA). The design of this particular probe is fairly standard, but the material used to house the steel needles and the connectors attaching them to the coaxial (data/signal) cables had a much higher melting temperature than normal. Additionally, those external portions of the coaxial cables that were likely to be exposed to high temperatures were wrapped in silicon tape. The calibration factor for this TDR is temperature 410 dependent and is discussed in more detail in the Appendix of Massman et al. (2010a). Like with the soil heat flux, the model's predictions will be evaluated against the TDR measurements with and without the temperature dependent calibration factor.
Soil CO 2 was measured by drawing a continuous sample for approximately 0.5 minutes through 3/8-inch (id) decabond tubing into a LI-820 (LI-COR Inc.; Lincoln, NE, USA) that was housed several m from the slash pile, as shown in Figure 2.  Figures 3 and 4 show the details of the surface energy balance for Equation (16), the model's upper boundary condition that includes the sensible heat and the net infrared terms. Included in Figure 3 is the forcing function, Q F (t), which shows the shape of the BFD-curve. Figure 5 shows the surface energy balance for Equation (18), which excludes the sensible heat and the net infrared terms. Comparing Figures 4 and 5 suggests that, once tuned appropriately, either formulation of the surface 420 energy balance will give very similar simulations of the evaporative flux, L * v0 E 0 . Likewise during the period of soil heating (i.e., when G 0 > 0), both formulations give similar and reasonable simulations for G 0 . But the simplified model of the surface energy balance, Equation (18), does not capture (and in fact cannot capture) G 0 during the period when the soil is cooling (i.e., when G 0 < 0, which starts at about 22 hrs in Figures 3 and 4). Without the possibility of radiative and convective cooling of the surface, Equation (18) does not reproduce G 0 < 0.

425
Another important aspect about the surface energy balance that Equation (18) does not capture as well as Equation (16)  These apparent limitations of the simplified surface energy balance, Equation (18), do not lessen the argument that it may be appropriate for some slash pile burns. Rather these present comparisons suggest that a hybrid of Equations (18) and (16) may be more appropriate. Such a hybrid would employ Equation (18) in the early part of the burn (before significant loss of mass from the slash pile due to combustion) and Equation (16) later after the fire intensity has peaked (i.e., sometime after t d ) when 435 the soil surface or possibly an ash surface is more exposed to the ambient environment.
Other than these two discrepancies the two forms of the surface energy balance give very similar simulations for soil temperature, moisture and other model variables. Nevertheless, because the full surface energy balance provides the more physically realistic simulation it will be used throughout the remainder of this study.

440
The principal aim of this model is to simulate reasonable and realistic soil temperatures during fires. The comparison between modeled and observed soil temperatures, shown in Figure 6, suggests that the HMV-model with the current set of 'tuned' parameters is reasonably good at this task. Nonetheless, a careful examination of the features of this figure shows some slight discrepancies in the model's performance: (1) the model does not capture the early temperature rise beginning at about 4 hrs, (2) nor does it capture the secondary maximum temperature at about 20 hrs, (3) it appears to overestimate the maximum 2 445 cm soil temperatures (at about 21 hrs), and (4) the soil appears to cool off faster than observed (most obvious after about 28 hrs). Discrepancies (1) and (2) are not unexpected because it is impossible for a simple forcing function like the BFD-curve and Equation (15) to produce an exact or even a nearly exact simulation of the observed temperature dynamics during a fire.
Any real physical surface forcing will always be far more (dynamically) complex than the BFD-curve, nor could it be easily generalized from one fire to the next. The (relatively significant) overestimation of the 2 cm soil temperatures is not fully 450 understood, but it may be a consequence of mis-measurement of the installation depth of the soil temperature sensors. Half of the (approximately) 60 C over-prediction can be accounted for if the sensor was installed at 2.6 cm rather than 2 cm. It is also possible that λ s possess greater vertical structure than is included in the model. Finally the issue involving the difference between the observed and modeled rates of cooling, discrepancy (4), seems to be characteristic of all model simulations, not just the simulation shown in Figure 6. Some of this is undoubtedly related to the simple shape of the BFD-curve (forcing 455 function) and the constraints it imposes on the upper boundary conditions. Another likely contributing factor to (4) (as well as the other 3 issues) is the ash layer that formed during the fire. A month after the fire, the ash layer measured between 0.5 and 8 cm deep within the burn area and was about 2 cm deep over the area where the sensors were buried. Because the ash layer would have insulated the soil surface, it would have acted to slow the rate of cooling as the fire died out. In general, the model does appear to capture many features of the observations. Of particular interest is that the amplitude of the modeled heat flux more closely agrees well with the Philip-corrected heat flux, providing confidence in both the model's performance and the quality of the soil heat flux data. On the other hand, the modeled heat flux peaks several hours before the 470 measured heat fluxes do. Without significantly increasing the complexity of the forcing function and the concomitant tuning effort, it is (at best) unlikely (if possible at all) to improve much on the model's ability to reproduce the observed temperature and heat flux data for this burn.

Soil Moisture
The modeled soil moisture is shown as a function of time in both Figures 8 and 9 for the same depths (and color-coding) as 475 shown for temperature in Figure 6. Included on each figure is the soil moisture measured at 5 cm (red) and 15 cm (magenta) depths; but Figure 8 includes the temperature-corrected calibration of the soil moisture probe, whereas Figure 9 does not. The red (5 cm) stars are interpolated values, which replace values that were flagged by the noise filter as questionable. Nonetheless, confidence in the fidelity of these interpolated values is high. But confidence in the magenta (15 cm) stars (also interpolated) is less. Rather interestingly the modeled 5 cm soil moisture agrees better with temperature-corrected soil moisture, whereas 480 at 15 cm the model resembles the uncorrected 15 cm soil moisture. Consequently, both figures are included here to show the importance of the temperature effects on the high-temperature TDR and to provide an estimate on the uncertainty inherent in these soil moisture measurements. Figure 10 shows the observed (temperature-corrected) soil moisture versus the observed temperature along with the model's solutions of θ vs T (or the model's θ − T trajectory). Here the model predicts that for depths less than 5 cm the soil moisture 485 begins to evaporate (decrease) between 50 C and 90 C and that these initial evaporative temperatures increase deeper into the soil. Nonetheless, the model suggests that the soil does not dry out completely until temperatures have reached about 200 C.
The observations at 5 cm show a similar pattern to the model, except these observations suggest that the initial evaporation stage proceeds more slowly and that overall evaporation occurs at a higher temperature. The same observation was made by Massman (2012) and Massman (2015) regarding the laboratory data of Campbell et al. (1995). Several test simulations were influence the amount of vapor within the soil pores, in turn influencing the soil moisture dynamics. This should not be too surprising in the case of S v , because S v is formulated within the present model as a balance between the physical processes that govern evaporation and condensation to and from soil surfaces, and therefore influences the balance between soil moisture and soil vapor. In the case of reducing the enhancement factor (i.e., assigning it a value of 1), the resulting decrease in vapor diffusivity causes the vapor within the soil pores to increases such that it feedbacks to soil moisture via the condensation term 500 of S v (proportional to ρ v , Equation (9)), and again influencing the balance between θ and ρ v . This exploration of the modeled and observed θ−T trajectory has yielded some insights into what maintains the long evaporative tail (where some soil moisture persists well past the boiling point of water), an issue that was less understood in Massman (2012) and Massman (2015).
Taken in total the model's simulation of soil moisture agrees with the observations fairly well. It is, of course, possible to improve on the model's fidelity to the observations by adjusting the hydraulic function (K R ) or changing in the heating rates or 505 duration and magnitude, Q F (t), and most importantly, by changing the values of the parameters of S v and D ve that influence soil vapor, ρ v ; but usually this comes at some expense to the model's fidelity to the measured soil temperatures. The same issue was noted in Massman (2012) and Massman (2015). Thus the present choice of model parameters is a compromise between two or three somewhat conflicting goals: fidelity to the soil's thermal response to heating by fire and the soil's moisture and vapor response.

Dynamic Feedbacks
Fire changes soil and the more intense the fire and the greater amount of soil heating, the greater will be the changes in the soil. Examples of some of the changes that are relevant to the present study include (a) changes the soil bulk density (Butters , 2009;Kojima et al. , 2018), which will alter the WRC (figure 1) and hydraulic properties (Tian et al. , 2018), (b) changes to the thermal conductivity of the soil (Massman et al. , 2008;Kojima et al. , 2018) and the volumetric specific heat (Butters ,515 2009), (c) changes to the heat and vapor fluxes resulting from (a) and (b) (Kojima et al. , 2018), and (d) changes to the soil's specific surface area, A wa , and particle surface potentials, which result when the fire induces water repellency in soils (Chen et al. , 2018). None of the physical/chemical processes causing these phenomena are included in any extant model of soil heating during fires. Part of the difficulty is that these are dynamic feedbacks occur during fires, consequently they influence soil heating and moisture transport in-situ during the fire. To give some idea of how these dynamic feedbacks may influence 520 the soil heat and moisture transport, this section details the results of a model simulation based on observed or inferred changes to the bulk density and key thermophysical parameters.
The following model parameters were changed for this sensitivity analysis to feedbacks: soil bulk density increases from 1.30 Mgm −3 to 1.46 Mgm −3 (a 12% increase as per figure 1); the thermal conductivity of the mineral fraction, λ m0 , increases from 4.42 WmK −1 to 8 WmK −1 , the de Vries shape factor, g a , decreased from 0.123 to 0.06, the Campbell et al. (1994) 525 parameter q w0 (which determines when water content starts to influence the soil's thermal conductivity) decreased from 0.03 to 0.02, the soil's volumetric specific heat increases by 10% (in accordance with the observations made by Butters (2009)), the overall soil thermal conductivity, λ s , increases by 15%, and finally the source term coefficient, S * , decreases from 0.1 to 0.08 (specifically chosen to be a 20% decrease). This increase in bulk density yields a concomitant decrease in soil porosity, η, which is simply carried over in a purely linear fashion to the WRC, the hydraulic function, and the source term, S v . On its 530 own this decrease in η yields an increase in S v , so S * is reduced to compensate for the decrease in A wa suggested by Chen et al. (2018). It very natural to expect the soil's thermal conductivity to increase with the increase in bulk density, but the present model of thermal conductivity does not explicitly include any dependency on ρ b . So λ s is increased by 15% to accord with the increase predicted by the model developed by Johansen (1975). Changes to g a and q w0 are intended to capture the observation (Massman et al. , 2008) that the fire decreased the sensitivity of the soil's thermal conductivity, λ s , to soil moisture, 535 (i.e., ∂λ s /∂θ decreased as a result of the fire) and increased λ s when the soil is dry (θ < 0.08). Decreasing g a and q w0 was the best way to capture this pair of observations. In addition a decrease in g a was also observed by Smits et al. (2016) for fireaffected soils. The conductivity of the mineral fraction, λ m , was increased to capture the creation of thin, but highly conductive, coatings of mineral oxides (MnO 2 in particular) on the soil particle surfaces during the MEF experimental burn (Massman et al. , 2010b;Nobles et al. , 2010). Such thin coatings are synthetically created for routine application to nanotechnologies (e.g.,

540
O' Brian et al. , 2013). The base case model simulation (Figs. 6,7,8,9,10) and the feedback simulation are compared in Figures 11, 12 and 13. The upper boundaries of the color-filled areas in Figures 11 (temperature) and 12 (heat flux) correspond to the feedback simulation.
Whereas, the lower boundary in Figure 13 (soil volumetric moisture) is associated with the feedback simulation. These com-soil and to significantly increase evaporation and evaporative losses of soil moisture (at least within the upper few centimeters of soil). However, the change in bulk density alone is responsible for about half of the differences shown in these figures.
Overall though, these last three figures indicate that these dynamic feedbacks are potentially quite significant to the magnitude and depth of the soil heating. Or paraphrasing somewhat, during a fire soil temperatures are not a direct linear response to the surface forcing function, rather at any given moment the soil heating feedbacks (non-linearly) upon itself creating conditions 550 that allow the heat to penetrate more deeply into the soil and the soil temperatures to exceed what they would have been without the feedbacks.

Improving physical realism: Future observational and modeling studies
Although the HMV-model gives reasonably realistic simulations of soil temperatures and moisture during fires, there are several enhancements that may further improve its performance and that should be useful to consider for future studies of the soil's 555 response to heating during fires.
(a) Dynamic Feedbacks and Soil Thermal Conductivity. A better physical understanding of the dynamical processes that govern how extreme heating during fires changes ρ b and the development of a model of λ s that captures this dynamic. Also a better formulation of how soil structure influences g a and λ s and how soil structure changes during a fire. This issue of improving model parameterizations of λ s is more complex that just including ρ b because the observed increase in λ s for a dry 560 soil was 200-300% (Massman et al. , 2008), which far exceeds the 15-17% predicted by the model of Johansen (1975).
(b) Mass Transport. The present version of the HMV-model assumes an advective transport of water vapor induced by the rapid volitalization of the soil moisture, u vl and Equation (5). But most formulations of advective transport in soils are based on Darcy's law. For application to fire u vl from Darcy's law would result from the rapidly evolving vapor pressure and temperature gradients. Additionally it would be worthwhile to include the dry air density, ρ d , as a separate model variable.

565
Certainly in any real fire the temperature and pressure of the dry air within the soil pore spaces would respond dynamically to heating. But including ρ d as a dynamic variable should yield a more physically realistic simulation of the diffusional and advective transport of water vapor during the fire. Finally, given the potential for extreme gradients in soil water potential and temperature during fires, it may also be worthwhile to include the heat transported by the movement of water (e.g., Stallman , 1965;Pasquale et al. , 2014) in Equation (1). This energy transport term has often been included when modeling the daily cycle 570 of energy flow through soils. Because some fires, especially slash pile burns, can continue a week or more, it seems appropriate to investigate the influence this type of energy transport could have on model solutions.
(c) 2-and 3-Dimensional Effects. There are (at least) two physical processes that cannot be fully represented in a 1-d model: (i) horizontal heat flux (G hor ) and (ii) possible advective currents (here characterized by an advective velocity u adv ) induced in the soil shortly after the pile is ignited ( (Massman et al. , 2010a;Nobles et al. , 2010;Massman et al. , 2010b)). (i) In an earlier 575 section mention was made of the second installation of soil sensors at the edge of the pile. These data used in conjunction with the soil temperature data at the center of the pile it is possible to obtain a crude estimate of the horizontal temperature gradient.
Then assuming that horizontal and vertical thermal conductivities (at the same depth) are about the same in magnitude, it is possible to show that G hor during the fire is approximately 10% of the vertical soil heat flux. This estimate of G hor was confirmed with data taken during another experimental burn performed in the fall of 2004, in which the 20 or so temperature 580 sensors were placed in a horizontal array that allowed about 20 different horizontal gradients to be compared with about the same number of vertical gradients. This is relevant to the present study because it may help explain some of the divergence between modeled and observed soil temperatures and (vertical) heat fluxes and further underscores the nature of the challenges and empiricism inherent in modeling of soil heating during fires. (ii) Far more important for any further studies is the possibility of a fire-induced u adv because it appears to be inherently 2-or 3-dimensional in nature and to carry combustion products into 585 the soil (Massman et al. , 2010b, figure 4). Such an advective current is hypothesized to have caused the extremely rapid 20-fold increase in soil CO 2 observed to have occurred within the first 30 minutes of the Manitou burn being modeled here (Massman et al. , 2010a). Consequently, any combustion-produced water vapor will also be transported. If this is the case it would likely overwhelms u vl and contradict assumptions made about evaporation at the soil surface and the transport of soil water vapor during (at least) some portion of the burn. It may be possible to just impose u adv in a 1-d model like the HMV-model, but such 590 an adjustment should probably be guided by observational studies and experiments designed to establish the existence, nature and dynamics of u adv .

Concluding Summary
This study describes the continuing development, improvement, and validation of the HMV-model, a non-equilibrium model of the coupled transport of heat, moisture, and water vapor in soils during surface fires. The purpose of the research supporting 595 this study is to provide a tool to aid in the management (and reduction, if possible) of the physical and ecological affects of fire on soils. Key improvements to the model, which noticeably improved its stability and performance, include more physically realistic parameterizations of the non-equilibrium vapor source term, S v , and soil thermal conductivity, λ s . Integral to the validation of this model are the development of a general surface heating (forcing) function and the discussions of the complexities and difficulties regarding formulating the surface (or upper) boundary condition. The model is validated using 600 in-situ measurements of soil temperatures, heat flux, and soil moisture obtained during a 2004 experimental burn carried out at Manitou Experimental Forest (in the central Rocky Mountains of Colorado). Despite any possible ambiguities in the calibration of the soil moisture and heat flux sensors and given the simplicity of the modeled forcing function and the complexities of the true forcing that can be inferred from the dynamics of the soil data, the model's ability to reproduce the observations is at least reasonable, and maybe even surprisingly good. But as with Massman (2012) and Massman (2015), tuning the model 605 parameters requires navigating the somewhat divergent goals of achieving a "best" fit to either the soil temperature observations or soil moisture observations. Absent to the model are important fire-induced feedbacks, in particular in-situ changes to the soil's thermal and hydraulic properties that are inevitable when the fire (as is often observed) causes the soil's bulk density to increase. This important dynamic was investigated with a model sensitivity analysis by making logical and credible changes in the (appropriate) model parameters. The net affect of this feedback is that the heat pulse will propagate much deeper into the soil than it would have otherwise, pointing to the need for further observational and modeling studies of this phenomenon. This study closes by highlighting other areas of research needed to improve our understanding of and ability to model the physical processes that occur in soils during fires, including issues involving the transport of both soil moisture and soil vapor during fires, the potentially very significant 2-or 3-dimensional advective flows in soils induced by fires, and the possibility of 2-or 3-dimensional heat flow.

615
Code and data availability. The computer code used in this study was developed using MatLab version 2017b and is publicly available along with any output data at the Forest Service Research Data Archive https://doi.org/10.2737/RDS-2020-YYYY. Prior to its availability online the code and any output is freely available from the author.
Appendix A: Soil Apparent Thermal Conductivity A1 Mass Gradient Diffusional Flux: −D ve ∂ρ v /∂z 620 de Vries (1958) defines the soil's apparent thermal conductivity, λ s,app , as the medium's thermal conductivity "including the effect of the vapor distillation due to temperature gradients". But in a modeling context, the use of λ s,app has to be treated with some care because it is inappropriate to externally introduce λ s,app into the equation for conservation of energy (employed universally as part of any model of heat and moisture flow in soils). This appendix develops both an equilibrium and a nonequilibrium model for λ s,app to illuminate its thermodynamic origins and to clarify its proper role in modeling of soil heat and 625 moisture flow in soils. Before proceeding, note that all terms relating to advective vapor flow (i.e., u vl ) and liquid water flow (i.e., K n , K H and V θ,surf ) can be ignored because they are superfluous for purposes of this appendix.
The basic conservation equations employed in this study are Equations (6), (7) and (4). But there is a valid alternative expression for Equation (6), which results by combining the conservation of enthalpy, Equation (1), with the conservation of water vapor, Equation (4), and the equilibrium assumption (ρ v = ρ v,eq = a w ρ v,sat ), and is expressed as which after some simple mathematical manipulation can also be written as where the last term on the right hand side of this equation accounts for moving L * v inside the gradient operator ∂/∂z. As discussed in the main text a w = a w (T K , ψ) is modeled by the Kelvin Equation, a w = e Mw ψ * RT K ψn , and ρ v,sat = ρ v,sat (T K ), so that ∂ρ v,eq /∂t can be expanded in terms of ∂T /∂t and ∂ψ n /∂t and ∂ρ v,eq /∂z in terms of ∂T /∂z and ∂ψ n /∂z. The following shows ∂ρ v,eq /∂z. (Note ∂ρ v,eq /∂t can be found by substituting t for z).
where ∆ sat (kgm −3 K −1 ) is the slope of the saturation vapor curve, dρ v,sat /dT , and 7×10 −4 < ∆ sat < 70×10 −4 for temperatures between about 5 and 60C. Furthermore, unless the soil is extremely dry ∆ sat ; implying that the second term multiplying the temperature gradient of Equation (A3) can be ignored relative to the first term for most applications. Otherwise for an extremely dry soil − In addition, it is always true that (ψ n /T K )∂T /∂z ∂ψ n /∂z, which means that the rightmost term in Equation (A3) can also be ignored. Therefore, Equation (A3) is now where the term L * v (η − θ)∂ρ v,eq /∂t has been dropped from Equation (A2) because |L * v (η − θ)∂ρ v,eq /∂t| |C s ∂T /∂t| (although it is a bit tedious to show). Equation (A5) yields the following identification for λ s,app : where the term describing "the effect of the vapor distillation due to temperature gradients" on the soil's thermal conductivity is a w L * v D ve ∆ * sat (identified in Equation (A6) as λ s,dis ). This expression for λ s,app is often used as the justification for including the effects of vapor transfer on λ s , i.e., the substitution of λ s,app for λ s in the equation of conservation of enthalpy in soils ((e.g., Hillel , 2004, Equation (12.24), p. 226) and (Campbell et al. , 1995;Smits et al. , 2011;Massman , 2015)). But as this appendix shows, in a modeling context, this substitution is usually unnecessary (and inappropriate) because these vapor transfer effects are already either directly or indirectly embedded in the equation for conservation of enthalpy.
Next is the introduction of the equilibrium assumption, Equation (A4), into Equation (7). This yields or to a very good approximation (which follows from Equation (A7) because ρ w ρ v,eq ) Equations (A5) and (A8) now form an equilibrium model with two independent predictive variables (T and θ or ψ) that describes the coupled heat and moisture flow in soils. This particular model clearly results in and is completely consistent with the soil's apparent thermal conductivity, λ s,app . On the other hand, there is an equally valid model of this coupled dual-variable model for which λ s,app is not only unnecessary, but it would, in fact, an error to invoke its use with the model. This other coupled heat and moisture flow model combines the conservation of mass for water (liquid + vapor), Equation (A8) above, with the simplified version of the conservation of energy equation, Equation (6), discussed in the development of the non-equilibrium model. This simplified equation is It is important to reiterate that the paired Equations (A8) (A5)) on the assumption that λ s should explicitly include the effects of water distillation, is incorrect because to do so is to double count the 675 effects of λ s,dis on soil thermal energy transport. This double counting is even more obvious (using the same argument that λ s should include the λ s,dis term) when λ s,app is substituted for λ s in Equation (A5), which already explicitly includes λ s,dis .
The non-equilibrium form of λ s,app follows from the same general methodology as the equilibrium form does from Equation (A5). First, combining Equations (1) and (4) and then simplifying yields Next employing the ideal gas law for ρ v , i.e., e v = ρ v RT K /M w yields where now the non-equilibrium form of λ s,dis is L * v D ve ρ v /T K . For practical applications this is obviously not as convenient or as useful as the equilibrium form because λ s,dis is now a function of ρ v , which for most experimental settings in soils is difficult (if not impossible) to measure directly. Nonetheless, for the purposes of this appendix, it suffices to show that a non-equilibrium model of λ s,dis can be defined and that in a modeling context it is just as unnecessary and inappropriate to use as is the better known equilibrium version of λ s,dis .

690
(A5), (A9)), is fundamental to all models of heat and moisture flow in soils. All forms of this conservation law also explicitly include the effects of the phase change of water through the vapor source term, L * v S v . Therefore, all effects associated with the distillation of water are a result of S v . This is true whether discussing an equilibrium or non-equilibrium model. Imposing or assuming liquid/vapor equilibrium or non-equilibrium has mathematical consequences to how ρ v and S v are parameterized and can also have indirect influence on the soil's intrinsic thermal conductivity (λ s ), because λ s is a function of several variables, 695 among which are soil volumetric water content and soil vapor density, i.e., λ s = λ s (θ, ρ v , · · · ) (de Vries (1963); Campbell et al. (1994); Tian et al. (2016)). But within a modeling context there is no justification for substituting the apparent thermal conductivity, λ s,app , for the soil's intrinsic thermal conductivity, λ s , in the equation of conservation of enthalpy, Equation (1).
To do so violates the conservation of enthalpy by effectively double counting λ s,dis .

A2 Mass Mixing Ratio Diffusional
The discussion in this section of the appendix complements the discussion of λ app and λ s,dis in A1 above. It does not change the final outcome or conclusions reached in A1. A2 is included here only to further refine λ s,dis . Furthermore, results of this section are limited to relatively moist soils, for which λ s,dis is of greatest interest and significance.
The discussion in A1 above develops λ s,app more-or-less along the traditional lines, using the mass gradient form of the diffusional flux, −D ve ∂ρ v /∂z. But the diffusional flux is also represented (and sometimes more appropriately represented) 705 in terms of the gradient of the mass mixing ratio, i.e., −ρ a D ve ∂χ v /∂z; where ρ a (kgm −3 ) is the total (dry air + vapor) soil gas density and χ v (kgkg −1 ) = ρ v /ρ a . This section of the appendix shows that for a relatively moist soil (i.e., a w ≈ 1 and χ v,eq ≈ χ v,sat ) λ s,dis developed in A1 above differs between +2.5% (for T ≈ 5 C) and -4.5% (for temperatures T ≈ 60 C) from that given in Equation (A6). All that this requires is to show how ρ a ∂χ v,eq /∂z generalizes Equation (A4).
From the identity ρ a = ρ d + ρ v,eq , where ρ d (kgm −3 ) is the soil dry air density, it follows that Next, assuming that ρ d is an ideal gas it follows that Because the pressure gradient term, (ρ d /p d )∂p d /∂z, is not relevant for the present purposes (and it is unlikely to contribute much to ∂ρ d /∂z anyway) it can be dropped from Equation (A14). This yields 715 ρ a ∂χ v,eq ∂z = (1 − χ v,eq ) ∂ρ v,eq ∂z + χ v,eq (ρ a − ρ v,eq ) T K ∂T ∂z (A15) which after combining with Equation (A4) and some further simplification yields from which it follows that This last expression indicates that there are two (moderately) compensating terms to λ s,dis that do not occur in the diffusional mass flux form for the apparent thermal conductivity. These are (1 − χ v,eq ) instead of 1 and ∆ * sat + ρ v,sat /T K instead of ∆ * sat . At about 5 C χ v ≈ 0.005 and ρ v,sat /T K ≈ 0.03∆ * sat . Whereas at 60 C χ v ≈ 0.10 and ρ v,sat /T K ≈ 0.06∆ * sat . These results imply that at 5 C Equation (A17) will yield a value for λ s,dis that is about 2.5% higher than λ s,dis = L * v D ve a w ∆ * sat and that near 60 C Equation (A17) yields a value for λ s,dis that is about 4.5% lower. Unless a correction between +2% and -4.5% is 725 important for estimating λ s,dis , then it seems that the mass mixing ratio formulation for the diffusional mass flux, Equation (A17), adds very little value to the original formulation, Equation (A6).
Author contributions. W. J. Massman is the sole author of this paper.
Competing interests. The author declares no competing interests.
Acknowledgements. The author would like to thank Dr M. Novak for his continued interest in this research and for many discussions on 730 various aspects of this work. I would also like to thank Dr G. Kuitenberg for sharing his time for discussions on all aspects of soil physics relevant to the research discussed in this manuscript.  Fredlund and Xing (1994) model. The only difference between these two curves is that the total or air-filled porosity is about 0.51 for the pre-burn (red) and it is about 0.45 for the post-burn (black).   (16). The parameters for the forcing function, QF (t) or Equation (15), are Qin = 120 Wm −2 , QF max = 18 kWm −2 , tm = 13.5 hrs and t d = 35 hrs. The Net Forcing is defined as the difference between the radiative terms, i.e., 0(θ0)QF (t) − 0(θ0)σ T 4 K0 − a(ρva)T 4 Ka , which from Equation (16) is equal to the sum of the three non-radiative terms: ρacpaCH [T0 − Ta]+L * v0 E0+ G0. Figure 4 is an expanded version of the Net Forcing and these three energy components.   (15), are Qin = 0 Wm −2 , QF max = 2.7kWm −2 , tm = 9.5 hrs and t d = 27.5 hrs. Here the Net Forcing is the sum L * v0 E0 + G0 because the surface sensible heat flux is not included as part the simplified surface energy balance.   (Philip , 1961). The lower boundary of this region is the measured heat flux after applying the Philip correction. The Philip correction is not shown for the two lower heat flux measurements because it made virtually no difference to these measurements. April 2004 Manitou Experimental Forest burn: modeled (solid lines) and observed (symbols). The color coding for depth is the same as that in Figure 6, but observations of soil moisture are taken only at two depths 5 cm (red) and 15 cm (magenta). The red stars are interpolated values and appear to be fairly trustworthy. The magenta stars are less trustworthy because the performance of the 15 cm TDR during the early part of the burn was not completely satisfactory.