The challenges of an in situ validation of a nonequilibrium model of soil heat and moisture dynamics during fires

With the increasing frequency and severity of fire, there is an increasing desire to better manage fuels and minimize, as much as possible, the impacts of fire on soils and other natural resources. Piling and/or burning slash is one method of managing fuels and reducing the risk and consequences of wildfire, but the repercussions to the soil, although very localized, can be significant and often irreversible. In an effort to provide a tool to better understand the impact of fire on soils, this study outlines the improvements to and the in situ validation of a nonequilibrium model for simulating the coupled interactions and transport of heat, moisture and water vapor during fires. Improvements to the model eliminate the following two important (but heretofore universally overlooked) inconsistencies: one that describes the relationship between evaporation and condensation in the parameterization of the nonequilibrium vapor source term, and the other that is the incorrect use of the apparent thermal conductivity in the soil heat flow equation. The first of these made a small enhancement in the stability and performance of the model. The second is an important improvement in the physics underpinning the model but had less of an impact on the model's performance and stability than the first. This study also (a) develops a general heating function that describes the energy input to the soil surface by the fire and (b) discusses the complexities and difficulties of formulating the upper boundary condition from a surface energy balance approach. The model validation uses (in situ temperature, soil moisture and heat flux) data obtained in a 2004 experimental slash pile burn. Important temperature-dependent corrections to the instruments used for measuring soil heat flux and moisture are also discussed and assessed. Despite any possible ambiguities in the calibration of the sensors or the simplicity of the parameterization of the surface heating function, the difficulties and complexities of formulating the upper boundary condition and the obvious complexities of the dynamic response of the soil's temperature and heat flux, the model produced at least a very credible, if not surprisingly good, simulation of the observed data. This study then continues with a discussion and sensitivity analysis of some important feedbacks (some of which are well known and others that are more hypothetical) that are not included in the present (or any extant) model, but that undoubtedly are dynamically influencing the physical properties of the soil in situ during the fire and, thereby, modulating the behavior of the soil temperature and moisture. This paper concludes with a list of possible future observational and modeling studies and how they would advance the research and findings discussed here.

1 Introduction

Fire has been a largely beneficial part of the landscape in most areas of the world for millennia (Harrison et al., 2010). But, over the past few decades, fire has increased significantly in frequency, extent and severity – to the point that it now poses substantial risks to most of the world's wildlands and forested ecosystems and the goods and services they provide (e.g., Kasischke and Turetsky, 2006; Mortiz et al., 2012; Abatzogloua and Williams, 2016; Stambaugh et al., 2018; San-Miguel-Ayanz et al., 2019). Consequently, there is an increased desire to reduce wildfire risk through better fuels management, mitigate the consequences of fire by improving management methods, and promote the recovery of soils and vegetation after fire (e.g., Millar et al., 2007; McCaffrey et al., 2015; Vallejo and Alloza, 2015; Schoennagel et al., 2017; Dey and Schweitzer, 2018). Therefore, to achieve any of these desired outcomes for managing the ecological effects of fire, it is necessary to improve our understanding of the impacts that fire and extreme heating can have on soils. An important and extremely useful tool in this effort is having better models of soil heating during fires. The benefit of a model lies in its ability to vary the amount and duration of heating (e.g., Steward et al., 1990) that characterize different fire types and to better judge how fire impacts different soil types. And, just as there are a variety of different soil types with differing thermophysical properties, there are also a variety fire types characterized by duration, intensity and, for the present purposes, the ability to heat the soil's surface, i.e., the soil heat flux (W m⁻²). For example, the prescribed burning of understory vegetation, which is an example of spreading or dynamic fire used to reduce understory fuels, is usually low intensity and will heat the soil for only a few minutes. But as the fuel loading increases, so does the fire's duration and intensity (e.g., Massman et al., 2003). On the other hand, prescribed slash pile burns, like the one studied here, are stationary fires that can that can sustain a soil heat flux of 1–10 kW m⁻² for tens of hours (e.g., Massman and Frank, 2004; Massman et al., 2008). Wildfires, and crown fires in particular, are dynamic and often fast moving and relatively brief by comparison to slash pile burns, but they can also be extraordinarily intense (10–100 kW m⁻²). Yet, despite their often extreme intensity, dynamic fires can produce a counterintuitive spatial pattern of soil surface heating, i.e., negligible heating near the region of the most intense fire to significant heating in areas of much lower fire intensity (Stoof et al., 2013). It is also possible that a fast-moving crown fire can cause burning material, e.g., tree boles or other woody material, to come into direct contact with the soil, producing a stationary fire that can continue to burn or smolder long after the fire front itself has passed. With the present modeling approach, it makes it relatively easy (just by changing the model's soil surface boundary forcing) to estimate the depth of penetration of critical temperature thresholds (e.g., Massman et al., 2010a, Fig. 1) for differing soil types. Here I summarize the changes made to the HMV (Heat–Moisture–Vapor) model (Massman, 2015) and assess the improvements they made to the model's performance by comparing modeled and observed (in situ) soil temperatures, heat fluxes and changes in soil moisture during a 2004 slash pile burn (Massman et al., 2008). Note that this model has been tested on several other instrumented slash pile burns and a few dynamic wildfires and controlled laboratory fires. Although the results will not be presented here (principally because the large majority of these fires were not instrumented with in situ soil moisture probes), these other fires did provide insights into soil heating dynamics during the fires and additional tests for building confidence in the model's performance.

The HMV model is a 1D (soil depth) model with three time-dependent predictive variables, namely temperature, T_K (K) or T (^∘C), soil water potential, ψ (J kg⁻¹) (where ψ<0 and is relatable to volumetric soil moisture θ (m³ m⁻³) through a water retention curve), and soil vapor density, ρ_v (kg m⁻³). At any specific depth, the model assumes thermal equilibrium between the soil matrix, the soil vapor and the soil moisture. However, it is termed a nonequilibrium model because it does not assume a priori that the soil moisture and soil vapor are in equilibrium, contrary to the equilibrium approach that has been the basis of virtually all models of coupled heat and moisture flow in soils since Philip and de Vries (1957) and de Vries (1958). Although the equilibrium assumption has led to many insights into the nature of soil heat and moisture transport processes in the last six decades, it must fail at some point as the soil dries out for the simple reason that it is difficult to maintain vapor in equilibrium with soil moisture when there is little to no soil moisture (Novak, 2012; Massman, 2015). In the case of rapid soil heating and drying during fires, Massman (2015) further indicates that at the drying front, where local soil evaporation rates are highest, θ and ρ_v are forced out of equilibrium as soil moisture rapidly decreases and the soil vapor rapidly increases. Novak (2019) also demonstrates (under less extreme conditions than during fires, i.e., $\mathrm{0}<T\le \mathrm{60}$ ^∘C) that the greatest departure from equilibrium occurs at the drying front. The equilibrium model cannot capture this evaporative disequilibrium, which may explain why soil evaporation is better modeled with a nonequilibrium approach (Smits et al., 2011; Ouedraogo et al., 2013; Massman, 2015; Borujerdi et al., 2019). In fact, the most important change/improvement in the HMV model (detailed in the next section) is in the parameterization of the vapor source term, S_v (kg m⁻³ s⁻¹), which is the essence of the nonequilibrium approach and its ability to capture the evaporative disequilibrium. As with its predecessor (Massman, 2015), the present vapor source term is formulated on the basis of the Hertz–Knudsen equation, which Trautz et al. (2015) have suggested describes evaporation better than other nonequilibrium models of S_v. Nonetheless, all extant models of S_v have overlooked (and 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 in 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 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 soil physical properties. The fifth section discusses possible future directions for modeling and observational studies. The final section summarizes this study.

2 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 paper also uses the same notation and the 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 will be updated 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.

2.1 Conservation of mass and energy

The HMV model is composed of the three conservation equations, namely 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 as follows:

$$\begin{array}{}\text{(1)}& C_{\mathrm{s}}{\displaystyle \frac{\partial T}{\partial t}}-{\displaystyle \frac{\partial }{\partial z}}\left[{\mathit{\lambda }}_{\mathrm{s}}{\displaystyle \frac{\partial T}{\partial z}}\right]=-L_{\mathrm{v}}{S}_{\mathrm{v}}+W{S}_{\mathrm{w}}\equiv -L_{\mathrm{v}}^{*}{S}_{\mathrm{v}},\end{array}$$

where C_s (J m⁻³ K⁻¹) is the volumetric specific heat of the soil, such that $C_{\mathrm{s}}=C_{\mathrm{s}}(T,\mathit{\theta })$ is a function of both temperature and volumetric soil moisture. t (s) is time, z (m) is soil depth, and λ_s (W m⁻¹ K⁻¹) is soil thermal conductivity, such that ${\mathit{\lambda }}_{\mathrm{s}}={\mathit{\lambda }}_{\mathrm{s}}(T,\mathit{\theta },{\mathit{\rho }}_{\mathrm{v}})$; L_v=L_v(T_K) (J kg⁻¹) is the enthalpy of vaporization and −L_vS_v represents the change in enthalpy associated with evaporation/condensation. $S_{\mathrm{v}}=S_{\mathrm{v}}(T_{\mathrm{K}},\mathit{\theta },\mathit{\psi },{\mathit{\rho }}_{\mathrm{v}})$ is the source term for water vapor and is discussed in more detail in the following section. WS_w is the change in enthalpy associated with the heat of wetting (also termed the heat of immersion), where W (J kg⁻¹) is the heat of wetting and S_w (kg m⁻³ s⁻¹) is the source term for water liquid, or equivalently the sink term for water vapor, i.e., $S_{\mathrm{w}}\equiv -S_{\mathrm{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=-\mathit{\psi }$ and 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 $-\mathit{\psi }\to \mathrm{\infty }$). Finally, with this identification for W and the above identity between S_w and S_v, it follows that $L_{\mathrm{v}}^{*}\equiv L_{\mathrm{v}}-\mathit{\psi }$.

The conservation of liquid water is as follows:

$$\begin{array}{}\text{(2)}& {\mathit{\rho }}_{\mathrm{w}}{\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}-{\mathit{\rho }}_{\mathrm{w}}{\displaystyle \frac{\partial }{\partial z}}\left[K_{\mathrm{n}}{\displaystyle \frac{\partial {\mathit{\psi }}_{\mathrm{n}}}{\partial z}}+K_H-V_{\mathit{\theta },\mathrm{surf}}\right]=-S_{\mathrm{v}},\end{array}$$

where ρ_w=ρ_w(T_K) (kg m⁻³) is the density of liquid water, and ψ_n (dimensionless) is the nondimensional form of ψ, i.e., ${\mathit{\psi }}_{\mathrm{n}}=\mathit{\psi }/{\mathit{\psi }}_{*}$, where ${\mathit{\psi }}_{*}=-{\mathrm{10}}^{\mathrm{6}}$ J kg⁻¹ is the nominal soil water potential of ovendried soil (Campbell et al., 1995). (Note that ψ_n is used interchangeably with ψ throughout this paper.) $K_{\mathrm{n}}=K_{\mathrm{n}}(T_{\mathrm{K}},{\mathit{\psi }}_{\mathrm{n}},\mathit{\theta })$ (m² s⁻¹) is the hydraulic diffusivity, $K_{\mathrm{H}}=K_{\mathrm{H}}(T_{\mathrm{K}},{\mathit{\psi }}_{\mathrm{n}},\mathit{\theta })$ (m s⁻¹) is the hydraulic conductivity, and $V_{\mathit{\theta },\mathrm{surf}}=V_{\mathit{\theta },\mathrm{surf}}(T_{\mathrm{K}},\mathit{\theta })$ (m s⁻¹) is the velocity of liquid water associated with surface diffusion of water. The hydraulic conductivity functions, $K_{\mathrm{n}}(T_{\mathrm{K}},{\mathit{\psi }}_{\mathrm{n}},\mathit{\theta })$ and $K_{\mathrm{H}}(T_{\mathrm{K}},{\mathit{\psi }}_{\mathrm{n}},\mathit{\theta })$, are given as follows:

$$\begin{array}{}\text{(3)}& K_{\mathrm{n}}={\displaystyle \frac{K_{\mathrm{I}}{K}_{\mathrm{R}}{\mathit{\rho }}_{\mathrm{w}}}{{\mathit{\mu }}_{\mathrm{w}}}}{\mathit{\psi }}_{*}\mathrm{and}{K}_{\mathrm{H}}={\displaystyle \frac{K_{\mathrm{I}}{K}_{\mathrm{R}}{\mathit{\rho }}_{\mathrm{w}}}{{\mathit{\mu }}_{\mathrm{w}}}}g,\end{array}$$

where μ_w=μ_w(T_K) (Pa s) (Huber et al., 2009) is the viscosity of water, and g=9.81 m s⁻² is the acceleration due to gravity. K_I (m²) 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). $K_{\mathrm{R}}=K_{\mathrm{R}}(\mathit{\theta },{\mathit{\psi }}_{\mathrm{n}},T_{\mathrm{K}})$ (dimensionless) is the relative hydraulic conductivity (used to describe capillary flow in soils). The model for intrinsic permeability is taken from Bear (1972) and is $K_{\mathrm{I}}=(\mathrm{6.17}\times {\mathrm{10}}^{-\mathrm{4}})d_{\mathrm{g}}^{\mathrm{2}}$, 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, as in Massman (2015), film flow did not really impact the model's performance.

The conservation of water vapor is as follows:

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\partial (\mathit{\eta }-\mathit{\theta }){\mathit{\rho }}_{\mathrm{v}}}{\partial t}}-{\displaystyle \frac{\partial }{\partial z}}\left[D_{\mathrm{ve}}{\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{v}}}{\partial z}}-(\mathit{\eta }-\mathit{\theta })u_{\mathrm{vl}}{\mathit{\rho }}_{\mathrm{v}}\right]=S_{\mathrm{v}},\end{array}$$

where η (m³ m⁻³) is the total soil porosity, assumed to be temporally constant and spatially uniform (note that η is obtained from the soil's bulk and particle densities, which are the actual model input variables, and because both of these variables are assumed constant and uniform, so also is η), and (η−θ) is the soil's air filled porosity. $D_{\mathrm{ve}}=D_{\mathrm{ve}}(T_{\mathrm{K}},\mathit{\psi },{\mathit{\rho }}_{\mathrm{v}})$ (m² s⁻¹) 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 et al. (1995) and detailed in Massman (2012). u_vl (m s⁻¹) is the advective velocity induced by the change in volume associated with the rapid volatilization of soil moisture (e.g., Ki et al., 2005), which is given as follows:

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\partial u_{\mathrm{vl}}}{\partial z}}={\displaystyle \frac{S_{\mathrm{v}}}{(\mathit{\eta }-\mathit{\theta }){\mathit{\rho }}_{\mathrm{v}}}}.\end{array}$$

The final model equations result from preserving Eq. (4) and eliminating S_v from Eqs. (1) and (2), such that Eq. (1) is replaced with the following:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{C}_{\mathrm{s}}& {\displaystyle \frac{\partial T}{\partial t}}-{\displaystyle \frac{\partial }{\partial z}}\left[{\mathit{\lambda }}_{\mathrm{s}}{\displaystyle \frac{\partial T}{\partial z}}\right]\\ & -L_{\mathrm{v}}^{*}{\mathit{\rho }}_{\mathrm{w}}\left({\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}-{\displaystyle \frac{\partial }{\partial z}}\left[K_{\mathrm{n}}{\displaystyle \frac{\partial {\mathit{\psi }}_{\mathrm{n}}}{\partial z}}+K_{\mathrm{H}}-V_{\mathit{\theta },\mathrm{surf}}\right]\right)=\mathrm{0},\end{array}\end{array}$$

and Eq. (2) is replaced with the following:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\mathit{\rho }}_{\mathrm{w}}& {\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}-{\mathit{\rho }}_{\mathrm{w}}{\displaystyle \frac{\partial }{\partial z}}\left[K_{\mathrm{n}}{\displaystyle \frac{\partial {\mathit{\psi }}_{\mathrm{n}}}{\partial z}}+K_{\mathrm{H}}-V_{\mathit{\theta },\mathrm{surf}}\right]\\ & +{\displaystyle \frac{\partial (\mathit{\eta }-\mathit{\theta }){\mathit{\rho }}_{\mathrm{v}}}{\partial t}}+{\displaystyle \frac{\partial }{\partial z}}\left[D_{\mathrm{ve}}{\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{v}}}{\partial z}}-(\mathit{\eta }-\mathit{\theta })u_{\mathrm{vl}}{\mathit{\rho }}_{\mathrm{v}}\right]=\mathrm{0}.\end{array}\end{array}$$

2.2 Improvements in nonequilibrium vapor source term

Massman (2015, Equation (10)) adapted the Hertz–Knudsen equation to develop the following formulation for the vapor source term, S_v:

$$\begin{array}{}\text{(8)}& S_{\mathrm{v}}=S_{*}{A}_{\mathrm{wa}}\sqrt{{\displaystyle \frac{R{T}_{\mathrm{K}}}{M_{\mathrm{w}}}}}\left({\mathcal{K}}_{\mathrm{e}}{\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}-{\mathcal{K}}_{\mathrm{c}}{\mathit{\rho }}_{\mathrm{v}}\right),\end{array}$$

where S_* is an empirical dimensionless parameter, which is tuned as necessary to ensure model stability. A_wa (m² m⁻³ or m⁻¹) is the volume-normalized soil water–air interfacial surface area, which Massman (2015) parameterized as A_wa=A_wa(θ). R=8.314 J mol⁻¹ K⁻¹ is the universal gas constant, M_w=0.01802 kg mol⁻¹ is the molar mass of water vapor, and 𝒦_e (dimensionless) is the mass accommodation (or evaporation) coefficient, which Massman (2015) sets ≡1. ${\mathcal{K}}_{\mathrm{c}}={\mathcal{K}}_{\mathrm{c}}(T_{\mathrm{K}},{\mathit{\psi }}_{\mathrm{n}})$ (dimensionless) is the thermal accommodation (or condensation) coefficient, which Massman (2015) parameterizes as a physicochemical (Arrhenius) function as follows: ${\mathcal{K}}_{\mathrm{c}}(T_{\mathrm{K}},{\mathit{\psi }}_{\mathrm{n}})=e^{\frac{E_{\mathrm{av}}-M_{\mathrm{w}}\mathit{\psi }}{R}\left(\frac{\mathrm{1}}{T_{\mathrm{K}}}-\frac{\mathrm{1}}{T_{\mathrm{K},\mathrm{in}}}\right)}$, where E_av−M_wψ (J mol⁻¹) is an empirical surface condensation/evaporation activation energy for which E_av≈30–40 kJ mol⁻¹ was determined empirically, and T_K,in is the initial soil temperature. ρ_v,eq (kg m⁻³) is the equilibrium vapor density, defined as ${\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}=a_{\mathrm{w}}{\mathit{\rho }}_{\mathrm{v},\mathrm{sat}}\left(T_{\mathrm{K}}\right)$, where $a_{\mathrm{w}}=e^{\frac{M_{\mathrm{w}}{\mathit{\psi }}_{*}}{R{T}_{\mathrm{K}}}{\mathit{\psi }}_{\mathrm{n}}}$ is the dimensionless water activity, modeled here with the Kelvin equation, and ρ_v,sat(T_K) (kg m⁻³) is the saturated vapor density, which is a function only of T_K.

But the model of S_v, embodied by Eq. (8), assumes that the interfacial surfaces appropriate for evaporation and condensation are the same, i.e., that A_wa is the same for both evaporation $\left(\sqrt{R{T}_{\mathrm{K}}/M_{\mathrm{w}}}{\mathcal{K}}_{\mathrm{e}}{\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}\right)$ and condensation $\left(\sqrt{R{T}_{\mathrm{K}}/M_{\mathrm{w}}}{\mathcal{K}}_{\mathrm{c}}{\mathit{\rho }}_{\mathrm{v}}\right)$. 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 monolayer of liquid water (e.g., Novak, 2019). But it is physically 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 nonequilibrium models of heat and moisture flow in dry soils.

The new version of the nonequilibrium model parameterizes S_v as follows:

$$\begin{array}{}\text{(9)}& S_{\mathrm{v}}=S_{*}\sqrt{{\displaystyle \frac{R{T}_{\mathrm{K}}}{M_{\mathrm{w}}}}}\left(A_{\mathrm{wa}}\left(\mathit{\theta }\right){\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}-A_{\mathrm{wa},\mathrm{dry}}{\mathcal{K}}_{\mathrm{c}}{\mathit{\rho }}_{\mathrm{v}}\right),\end{array}$$

where 𝒦_e≡1 has been retained – as has the original formulation for A_wa (Massman, 2015), as follows:

$$\begin{array}{}\text{(10)}& A_{\mathrm{wa}}\left(\mathit{\theta }\right)={\mathcal{S}}_{\mathrm{w}}{\left(\mathrm{1}-{\mathcal{S}}_{\mathrm{w}}\right)}^{a_{\mathrm{1}}}+a_{\mathrm{2}}{\left[{\mathcal{S}}_{\mathrm{w}}\left(\mathrm{1}-{\mathcal{S}}_{\mathrm{w}}\right)\right]}^{a_{\mathrm{3}}},\end{array}$$

where ${\mathcal{S}}_{\mathrm{w}}=\mathit{\theta }/\mathit{\eta }$ is the soil water saturation and a₁=50 (rather than the original value of 40), a₂=0.003 and $a_{\mathrm{3}}=\mathrm{1}/\mathrm{8}$. This particular value for the parameter a₁ was chosen so that the maximum value of A_wa occurs at 𝒮_w≈0.02 ($=\mathrm{1}/a_{\mathrm{1}}$) and is assumed to be where the soil surfaces are covered by a monolayer of water (Brusseau et al., 2006). $A_{\mathrm{wa},\mathrm{dry}}\equiv A_{\mathrm{wa}}\left(\mathit{\theta }\right)$, as long as ${\mathcal{S}}_{\mathrm{w}}>\mathrm{1}/a_{\mathrm{1}}$, and $A_{\mathrm{wa},\mathrm{dry}}\equiv \max \left(A_{\mathrm{wa}}\right)$ whenever ${\mathcal{S}}_{\mathrm{w}}\le \mathrm{1}/a_{\mathrm{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 monolayer of water. Empirical tuning of S_* and E_av, after implementing the other changes, yielded $\mathrm{0.01}\le S_{*}\le \mathrm{0.1}$ and E_av=10 kJ mol⁻¹. Together, these changes in S_v improved the model's stability and robustness and its fidelity to the observed soil moisture during the 2004 burn (detailed later).

2.3 Corrections and improvements in soil thermal conductivity

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; 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, in the following:

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& {\mathit{\lambda }}_{\mathrm{s}}=\\ & {\displaystyle \frac{k_{\mathrm{w}}\mathit{\theta }{\mathit{\lambda }}_{\mathrm{w}}\left(T_{\mathrm{K}},{\mathit{\rho }}_{\mathrm{w}}\right)+k_{\mathrm{a}}[\mathit{\eta }-\mathit{\theta }]{\mathit{\lambda }}_{\mathrm{a}}^{*}\left(\mathit{\theta },T_{\mathrm{K}},{\mathit{\rho }}_{\mathrm{v}}\right)+k_{\mathrm{m}}[\mathrm{1}-\mathit{\eta }]{\mathit{\lambda }}_{\mathrm{m}}}{k_{\mathrm{w}}\mathit{\theta }+k_a[\mathit{\eta }-\mathit{\theta }]+k_{\mathrm{m}}[\mathrm{1}-\mathit{\eta }]}}\\ & +\mathrm{3.8}\mathit{\sigma }{N}^{\mathrm{2}}{R}_{\mathrm{p}}{T}_{\mathrm{K}}^{\mathrm{3}},\end{array}\end{array}$$

where k_w, k_a and k_m (dimensionless) are the Campbell et al. (1994) generalized formulations of the de Vries (1963) weighting factors. λ_w(T_K,ρ_w) (W m⁻¹ K⁻¹) is the thermal conductivity of liquid water, and ${\mathit{\lambda }}_{\mathrm{a}}^{*}(\mathit{\theta },T_{\mathrm{K}},{\mathit{\rho }}_{\mathrm{v}})$ (W m⁻¹ K⁻¹) is the apparent thermal conductivity of moist air, which is the sum of the true thermal conductivity of moist air, ${\mathit{\lambda }}_{\mathrm{a}}(\mathit{\theta },T_{\mathrm{K}},{\mathit{\rho }}_{\mathrm{v}})$ (W m⁻¹ K⁻¹) and the term, ${\mathit{\lambda }}_{\mathrm{v}}^{*}(\mathit{\theta },T_{\mathrm{K}},{\mathit{\rho }}_{\mathrm{v}})$ (W m⁻¹ K⁻¹), 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 (W m⁻¹ K⁻¹) is the thermal conductivity of the mineral components of the soil, σ (W m⁻² K⁻⁴) is the Stefan–Boltzmann constant, $N=N\left(\mathit{\theta }\right)=\mathrm{1}+\mathit{\theta }/\left(\mathrm{3}\mathit{\eta }\right)$ (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 Eq. (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), as follows:

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{k}_{*}=& k_{*}(\mathit{\theta },T)={\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\\ & \left[{\displaystyle \frac{\mathrm{2}}{\mathrm{1}+\left(\frac{{\mathit{\lambda }}_{*}}{{\mathit{\lambda }}_{\mathrm{f}}}-\mathrm{1}\right)g_{\mathrm{a}}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\left(\frac{{\mathit{\lambda }}_{*}}{{\mathit{\lambda }}_{\mathrm{f}}}-\mathrm{1}\right)\left(\mathrm{1}-\mathrm{2}{g}_{\mathrm{a}}\right)}}\right],\end{array}\end{array}$$

where the subscript ^* refers to water (w), air (a) or mineral (m). g_a is the de Vries (1963) shape factor, an empirically determined model parameter. In general, g_a≈0.1 (Campbell et al., 1994). ${\mathit{\lambda }}_{\mathrm{f}}={\mathit{\lambda }}_{\mathrm{f}}(\mathit{\theta },T)$ is a weighted mixture of the thermal conductivities of air and water as follows:

$$\begin{array}{}\text{(13)}& {\mathit{\lambda }}_{\mathrm{f}}(\mathit{\theta },T)={\mathit{\lambda }}_{\mathrm{a}}\left(\mathit{\theta },T,{\mathit{\rho }}_{\mathrm{v}}\right)+f_{\mathrm{w}}(\mathit{\theta },T)\left[{\mathit{\lambda }}_{\mathrm{w}}\left(T\right)-{\mathit{\lambda }}_{\mathrm{a}}\left(\mathit{\theta },T,{\mathit{\rho }}_{\mathrm{v}}\right)\right],\end{array}$$

and

$$\begin{array}{}\text{(14)}& f_{\mathrm{w}}(\mathit{\theta },T)={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+{\left(\frac{\mathit{\theta }}{{\mathit{\theta }}_{\mathrm{w}}}\right)}^{-q_{\mathrm{w}}\left(T\right)}}},\end{array}$$

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_{\mathrm{w}}\left(T_{\mathrm{K}}\right)=q_{\mathrm{w}\mathrm{0}}(T_{\mathrm{K}}/\mathrm{303}{)}^{\mathrm{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.

A total of three changes have been made to this original formulation. First, ${\mathit{\lambda }}_{\mathrm{v}}^{*}$ is no longer included because to do so is to double 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 nonequilibrium models of soil heating include the vapor distillation term, either explicitly through the conservation of mass of water vapor (nonequilibrium models) or implicitly through the conservation of mass of liquid water (equilibrium models). Consequently, it is unnecessary and redundant to include ${\mathit{\lambda }}_{\mathrm{v}}^{*}$ in Eq. (11).

Second, λ_m now includes an explicit temperature dependency of the form ${\mathit{\lambda }}_{\mathrm{m}}={\mathit{\lambda }}_{\mathrm{m}}\left(T_{\mathrm{K}}\right)={\mathit{\lambda }}_{\mathrm{m}\mathrm{0}}\left(\mathrm{8}\exp \right(-\mathrm{0.008}(T_{\mathrm{K}}-\mathrm{300}))+\mathrm{3})/\mathrm{11}$, where λ_m0 (W m⁻¹ K⁻¹) is basically an adjustable parameter. The results for α quartz (α referring to a specific crystalline structure of quartz) from (Yoon et al., 2004, Fig. 4) suggests that it is reasonable to expect that λ_m0≤15 W m⁻¹ K⁻¹. The temperature function $\left(\mathrm{8}\exp \right(-\mathrm{0.008}(T_{\mathrm{K}}-\mathrm{300}))+\mathrm{3})/\mathrm{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 and 600 K, as shown by Kanamori et al. (1968) and Yoon et al. (2004). At the Manitou Experimental Forest (MEF) burn site, and more generally throughout the region surrounding MEF, sand quartz is likely to be the dominant soil mineral, but other crystalline mineral forms, such as granite and feldspar, are also common (e.g., Retzer, 1949; Mathews, 1900; Smith et al., 1999). Although the thermal conductivities associated with these other mineral forms also decrease with increasing temperature (Heuze, 1983; Mottaghy et al., 2008; Miao et al., 2014), overall they decrease with about half the slope ($\partial {\mathit{\lambda }}_{\mathrm{m}}/\partial T$) of the sand quartz parameterization used here. They also tend to show a much lower λ_m0 than quartz sand. Here λ_m0=4.42 W m⁻¹ K⁻¹, which should help accommodate these other nonsand quartz mineral forms.

Third, R_p is now estimated as $R_{\mathrm{p}}=\mathrm{0.408}{d}_{\mathrm{g}}\sqrt{({\mathit{\rho }}_{\mathrm{p}}/{\mathit{\rho }}_{\mathrm{b}})-\mathrm{1}}$ from Arya et al. (1999), where d_g (m) is the mean or effective soil particle diameter, ρ_p (Mg m⁻³) is the particle density (which can usually be assumed to be about 2.65 Mg m⁻³) and ρ_b (Mg m⁻³) 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 soils, even during fires. Massman (2015) reached a similar conclusion.

2.4 Evaluation of changes in 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 in S_v and λ_s to those shown in Massman (2015). The results indicated that (a) the model's ability to faithfully reproduce the data was very similar to those in Massman (2015), and (b) the model's stability was noticeably 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 in λ_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.

2.5 Complexities and challenges

2.5.1 Surface heating, surface energy balance and upper boundary condition

The forcing function is the energy that is input to the soil at its surface, denoted here by Q_F(t) (W m⁻²). 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 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 burning slash pile, surface heating may be combination of radiation and conduction, it may change over time as the pile burns and as the ash accumulates and, at 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 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 Campbell et al. (1995), which more or less dictated the following forcing function: $Q_{\mathrm{F}}\left(t\right)=Q_{\mathrm{Fmax}}(\mathrm{1}-\exp (-t/\mathit{\tau }\left)\right)$, where τ (s) and Q_Fmax (W m⁻²) are adjustable parameters that describe the rate of rise of the forcing function (τ) and its steady state asymptotic (maximal) value Q_Fmax. This study uses the following modified form of the “BFD curve” (Barnett, 2002) as the forcing function as follows:

$$\begin{array}{}\text{(15)}& Q_{\mathrm{F}}\left(t\right)=Q_{\mathrm{Fin}}+\left(Q_{\mathrm{Fmax}}-Q_{\mathrm{Fin}}\right)e^{-\mathit{\alpha }{\left(\ln \left(t/t_{\mathrm{m}}\right)\right)}^{\mathrm{2}}},\end{array}$$

where $e^{-\mathit{\alpha }\left(\ln \right(t/t_{\mathrm{m}}){)}^{\mathrm{2}}}$ is the modified BFD curve, t_m (s) is the time at which the maximum forcing occurs, α (dimensionless) = $\mathrm{2}\ln \left(\mathrm{10}\right)/\left({\sinh }^{-\mathrm{1}}\right[\mathrm{0.5}{t}_{\mathrm{d}}/t_{\mathrm{m}}]{)}^{\mathrm{2}}$, 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_Fmax, t_m and t_d are adjustable input parameters to the model that define the forcing function, much the same way that Q_Fmax and τ are for the Campbell et al. (1995) forcing function. Q_Fin (W m⁻²), on the other hand, is not an adjustable input parameter. It is determined from other considerations of the soil surface energy balance and is defined later. The boundary conditions for temperature and vapor pressure, e_v (Pa), have similar functional forms as Eq. (15) as follows: $V_{\mathrm{a}}=V_{\mathrm{in}}+V_{\max }{e}^{-\mathit{\alpha }\left(\ln \right(t/t_{\mathrm{m}}){)}^{\mathrm{2}}}$, 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 that is tuned so that the model matches (as much as possible) the soil observational data and (if necessary) helps ensure model stability. The vapor density boundary condition is derived from these latter two boundary conditions using the ideal gas law. The upper boundary condition on soil moisture is $(\partial \mathit{\theta }/\partial z{)}_{\mathrm{0}}=\mathrm{0}$ and is discussed further in Massman (2015).

The energy balance at the soil surface used with the present study formulates the net infrared radiation loss at the surface as a balance between the outgoing and incoming infrared radiation. This is different from both Massman (2012) and Massman (2015), which did not include the possibility of incoming environmental infrared radiation being absorbed by the soil's surface. Here the surface energy balance is expressed as follows:

$$\begin{array}{}\text{(16)}& \begin{array}{rl}{\mathit{ϵ}}_{\mathrm{0}}\left({\mathit{\theta }}_{\mathrm{0}}\right)Q_{\mathrm{F}}\left(t\right)& ={\mathit{ϵ}}_{\mathrm{0}}\left({\mathit{\theta }}_{\mathrm{0}}\right)\mathit{\sigma }\left[T_{\mathrm{K}\mathrm{0}}^{\mathrm{4}}-{\mathit{ϵ}}_{\mathrm{a}}\left({\mathit{\rho }}_{\mathrm{va}}\right)T_{\mathrm{Ka}}^{\mathrm{4}}\right]\\ & +{\mathit{\rho }}_{\mathrm{a}}{c}_{\mathrm{pa}}{C}_{\mathrm{H}}\left[T_{\mathrm{0}}-T_{\mathrm{a}}\right]+L_{\mathrm{v}\mathrm{0}}^{*}{E}_{\mathrm{0}}+G_{\mathrm{0}},\end{array}\end{array}$$

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 the absorptivity and emissivity of the soil are the same), and the first term on the right-hand side is the net infrared heat loss (where the term $\propto {\mathit{ϵ}}_{\mathrm{a}}\left({\mathit{\rho }}_{\mathrm{va}}\right)T_{\mathrm{Ka}}^{\mathrm{4}}$ was not included in Massman (2012) or Massman (2015)). The second term is the convective heat loss, the third is the rate of evaporation and the last is the soil conductive heat. ϵ₀(θ₀) is the soil emissivity and is a function of the soil moisture, θ₀. $\mathit{\sigma }=\mathrm{5.67}\times {\mathrm{10}}^{-\mathrm{8}}$ W m⁻² K⁻⁴ 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 ambient vapor density ρ_va, following the clear sky parameterization of Brutsaert (1984, Eq. (6.18)). ${\mathit{\rho }}_{\mathrm{a}}={\mathit{\rho }}_{\mathrm{a}}\left(T_{\mathrm{K}\mathrm{0}}\right)=\mathrm{1.29}(P_{\mathrm{a}}/P_{\mathrm{ST}})(T_{\mathrm{ST}}/T_{\mathrm{K}\mathrm{0}})$ (kg m⁻³) 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=101 325 Pa and T_ST=273.15 K are the standard atmospheric pressure and temperature. C_H=0.032 (m s⁻¹) 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) (K), is the ambient temperature somewhere above the soil surface (upper boundary condition), $L_{\mathrm{v}\mathrm{0}}^{*}{E}_{\mathrm{0}}$ (W m⁻²) is the rate of soil water evaporation, E₀ (kg m⁻² s⁻¹) is the evaporative mass flux at the surface, and G₀ (W m⁻²) is the soil conductive heat flux and the upper boundary condition for the modeled soil temperatures (Eq. 6). E₀ is parameterized as the sum of a diffusional component and an advective component (Massman, 2012) as follows:

$$\begin{array}{}\text{(17)}& E_{\mathrm{0}}=C_{\mathrm{E}}{h}_{\mathrm{s}\mathrm{0}}\left[{\mathit{\rho }}_{\mathrm{v}\mathrm{0}}-{\mathit{\rho }}_{\mathrm{va}}\left(t\right)\right]+C_{\mathrm{U}}{u}_{\mathrm{vl}\mathrm{0}}{\mathit{\rho }}_{\mathrm{v}\mathrm{0}},\end{array}$$

where C_E (m s⁻¹) and C_U (dimensionless) are adjustable model transfer coefficients which were determined empirically to maximize E₀ without destabilizing the model. C_U=0.125 is associated with the jet of volatilized air emanating from the soil with velocity u_vl0 and $C_{\mathrm{E}}={\mathrm{10}}^{-\mathrm{4}}$ m s⁻¹. The surface humidity, h_s0 (dimensionless), is assumed to be the same as a_w0, which is the water activity at the surface.

Although Eq. (16) completely accounts for the energy exchange between the atmosphere and the soil, there are two key issues that need to be addressed when initializing the model for applications to fires. The first was mentioned earlier; that is, 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. Further complicating a proper formulation of the soil surface infrared (IR) and convective heat flow terms are the structure and porosity of the pile, which will influence the IR impinging on the soil surface, the airflow ventilating the pile interior, and the size and efficiencies of any convective eddies that may act to transfer heat toward or away from the soil surface. In this study, the upper boundary condition for temperature, T_a(t), is chosen to ensure that the convective heat flow is away from the soil surface. 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 these later two terms from Eq. (16) yields the following simplified version of the soil surface energy balance:

$$\begin{array}{}\text{(18)}& {\mathit{ϵ}}_{\mathrm{0}}\left({\mathit{\theta }}_{\mathrm{0}}\right)Q_{\mathrm{F}}\left(t\right)=L_{\mathrm{v}\mathrm{0}}^{*}{E}_{\mathrm{0}}+G_{\mathrm{0}}.\end{array}$$

The second issue is basically a feature of all forcing functions. In the case of the BFD curve, t_m≥0.5 h 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_Fin. The choice of Q_Fin depends on whether Eq. (16) or Eq. (18) is used. In the case of Eq. (18), Q_Fin=0. In this case, the soil conductive heat flux becomes $G_{\mathrm{0}}=-L_{\mathrm{v}\mathrm{0}}^{*}{E}_{\mathrm{0}}$, and because the source term and the evaporation rates are very nearly at equilibrium (i.e., S_v≈0 and E₀≈0) during this ramp-up, it follows that G₀≈0 and that $\partial T/\partial t\approx \mathrm{0}$. (Note that the soil is initialized to be isothermal at the temperature obtained at the soil surface just before initiating the burn.) But for Eq. (16), the full surface energy balance equation, this equilibrium condition does not occur during the ramp-up if Q_Fin=0 because the net IR term in Eq. (16) is not initially in equilibrium. This is true despite initializing T_a=T₀. Without Q_Fin, Eq. (16) reduces to $G_{\mathrm{0}}\approx -{\mathit{ϵ}}_{\mathrm{0}}\left({\mathit{\theta }}_{\mathrm{0}}\right)\mathit{\sigma }{T}_{\mathrm{K}\mathrm{0}}^{\mathrm{4}}(\mathrm{1}-{\mathit{ϵ}}_{\mathrm{a}}({\mathit{\rho }}_{\mathrm{va}}\left)\right)<\mathrm{0}$, which induces a transient in the solution that causes the soil temperature to drop slightly (i.e., $\partial T/\partial t<\mathrm{0}$). But assigning $Q_{\mathrm{Fin}}\approx \mathit{\sigma }{T}_{\mathrm{K}\mathrm{0}}^{\mathrm{4}}(\mathrm{1}-{\mathit{ϵ}}_{\mathrm{a}}({\mathit{\rho }}_{\mathrm{va}}\left)\right)$ eliminates this transient and ensures that G₀≈0 and $\partial T/\partial t\approx \mathrm{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 the best way to prevent this initial transient. Note that this initial period of disequilibrium 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) and with the Campbell forcing function used in Massman (2012) and Massman (2015).

2.5.2 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 (2012) and Massman (2015), i.e., the second derivative, ${\partial }^{\mathrm{2}}/\partial z^{\mathrm{2}}$, of the three 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 model predictions within the upper few centimeters of the soil too much. The advective velocity, u_vl, requires only one boundary condition – see Eq. (5) – and is u_vl=0 at the bottom boundary.

Although this combination of the pass-through lower boundary condition and placement depth (0.6 m) is fairly general, it is still possible that they may influence the model's solution within the region of greatest interest (the top 20 cm or so). To test for this possibility, a sensitivity test was performed by increasing the domain depth to 1 m. The resulting change in the model's solution was negligible throughout the upper 0.6 m. A total of two further simulations were done each with a zero-flux lower boundary condition for soil moisture (i.e., $\partial \mathit{\theta }/\partial z{|}_{z_{\mathrm{bottom}}}=\mathrm{0}$ or, more specifically, $\partial \mathit{\psi }/\partial z{|}_{z_{\mathrm{bottom}}}=\mathrm{0}$), where one had a domain depth of 0.6 m and the other had 1 m depth. The impact of these changes on the model simulation was, for the present purposes, also negligible. This zero-flux lower boundary condition for θ would seem to be a credible, but apparently unnecessary, and alternative to the present pass-through boundary condition, which is preferred here because it allows all model variables and their associated fluxes to be dynamically modeled. θ was chosen for this sensitivity analysis because it can be plausibly argued that there is enough known about a site's soil structure to make a credible guess at the depth at which the zero-flux boundary condition could be assumed. Soil temperature is better dynamically modeled at the lower boundary because the heat pulse for the fire can extend far deeper than 0.6 m, regardless of the soil structure (e.g., Massman and Frank, 2004). Finally, the pass-through formulation is also the logical default for vapor density because (a) ρ_v is more strongly controlled locally by the balance between evaporation and condensation than by the lower boundary condition, and (2) simply decreasing the effective vapor diffusivity, D_ve , is sufficient for mimicking regions of highly impermeable soil without the need to explicitly incorporate a specific non-pass-through boundary condition.

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 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, which is discussed next.

2.6 The water retention curve

Unlike previous studies (Massman, 2012, 2015) which employed soils for which the water retention curve was unavailable, the present study employs a water retention curve appropriate to the soils at the MEF experimental burn site. Figure 1 shows the water retention curves (WRCs) for preburn (red) and postburn (black) soil. Note that both pre- and postburn soils are included in this study because they are part of the model sensitivity analysis. The data shown in this figure were provided by Gregory Butters (personal communication, 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 fitted with the Fredlund and Xing (1994) model as follows:

$$\begin{array}{}\text{(19)}& \mathit{\theta }/\mathit{\eta }={\mathcal{S}}_{\mathrm{w}}=\left(\mathrm{1}-{\displaystyle \frac{\ln \left(\mathrm{1}+a{\mathit{\psi }}_{\mathrm{n}}\right)}{\ln (\mathrm{1}+a)}}\right){\left(\ln \left(e+{\left(b{\mathit{\psi }}_{\mathrm{n}}\right)}^n\right)\right)}^{-m},\end{array}$$

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: ${\mathit{\rho }}_{\mathrm{b},\mathrm{pre}}=\mathrm{1.30}$ Mg m⁻³ and ${\mathit{\rho }}_{\mathrm{b},\mathrm{post}}=\mathrm{1.46}$ Mg m⁻³. 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 fitted to the observations.

Figure 1 Pre- and postburn water retention curves at the 2004 Manitou Experimental Forest (MEF) burn site. Observations were fitted with the 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 preburn (red), and it is about 0.45 for the postburn (black).

For the HMV model, the choice of WRC is important. The Fredlund–Xing model did provide the best fit to the data, and its impact on model performance was judged to be 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, Eq. (19) and Campbell–Shiozawa model and (the class of models represented here by) the van Genuchten and the Groenevelt–Grant models is the WRC's description and mathematical behavior when the soil is extremely or completely dry. Under these conditions, Eq. (19) when θ≈0, the soil water potential remains bounded, i.e., $\mathit{\psi }\approx -{\mathrm{10}}^{\mathrm{6}}$ J kg⁻¹. But, with the van Genuchten and the Groenevelt–Grant models, a completely dry soil is impossible to achieve because θ=0 can only occur when $\mathit{\psi }=-\mathrm{\infty }$. In the HMV model an unbounded function for the WRC does result in the degradation in model performance, i.e., some loss in the physical realism in the simulation of soil moisture, and can introduce model instabilities.

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, 2015). In the case of Massman (2012), when θ=0, $\mathit{\psi }=-{\mathrm{10}}^{\mathrm{6}}$ J kg⁻¹ and $\partial \mathit{\psi }/\partial t=\mathrm{0}$. But as long as the temperature keeps rising, $\partial T/\partial t>\mathrm{0}$, it follows that $\partial {\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}/\partial t>\mathrm{0}$. For the equilibrium model (${\mathit{\rho }}_{\mathrm{v}}\equiv {\mathit{\rho }}_{\mathrm{v},\mathrm{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), 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 led to the improved parameterization of S_v embodied in Eq. (9) of the present study.

Additionally, temperature can also have significant effects on the WRC (e.g. Olivella and Gens, 2000; Salager et al., 2010), which is a consequence of the influence temperature has on the surface tension of water (σ_w=σ_w(T_K) N m⁻¹ (Vargaftik et al., 1983). Massman (2012) found that these effects were relatively small, even at very high temperatures, but his 2012 model did not include any possibility of soil water movement (i.e., K_R=0). Milly (1984) includes the temperature effects on both the WRC (σ_w) and the hydraulic conductivity functions (μ_w) and, likewise, found little effect, but his study was restricted to relatively normal soil temperatures (i.e., T≤60 ^∘C). Here the impact that this temperature dependency (σ_w) has on the WRC and the model solution was investigated by multiplying the model variable ψ_n of Eq. (19) by the function ${\mathit{\varphi }}_{\mathrm{n}}\left(T_{\mathrm{K}}\right)={\mathit{\sigma }}_{\mathrm{w}}\left(T_{\mathrm{K},\mathrm{in}}\right)/{\mathit{\sigma }}_{\mathrm{w}}\left(T_{\mathrm{K}}\right)$, which ensures that $\partial \mathit{\theta }/\partial T_{\mathrm{K}}<\mathrm{0}$ is in accordance with observations (e.g. Olivella and Gens, 2000; Salager et al., 2010), and that ϕ_n is consistent with theoretical considerations and similar to other model parameterizations (e.g. Milly, 1984; Zhou et al., 2014). It also ensures that ϕ_n=1 at the beginning of the model simulation. The results suggested that temperature has only a minor influence on the WRC and the model solution, and because the model's sensitivity to ϕ_n is small compared to the model's sensitivity to changes in other parameters and parameterizations, this issue will not be considered any further in this study.

2.7 The hydraulic conductivity functions

The hydraulic conductivity functions, K_n and K_H, are both basically determined by K_R, as discussed above (Eq. 3). Like Massman (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, in the following:

$$\begin{array}{}\text{(20)}& K_{\mathrm{R}}\left(\mathit{\theta }\right)={\left(\mathrm{1}-{\left[\mathrm{1}-{\left({\displaystyle \frac{\mathit{\theta }}{\mathit{\eta }}}\right)}^{\frac{\mathrm{1}}{m_k}}\right]}^{m_k}\right)}^{n_k},\end{array}$$

where m_k and n_k are parameters ($\mathrm{0}<m_k<\mathrm{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 no data available to determine K_R for the soil at the burn site.

3 Manitou Experimental Forest – burn site and instrumentation

Much of the general description of climate and physical characteristics of the MEF 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.

3.1 General site and soil description

The burn experiment is located within the MEF (39^∘04^′ N, 105^∘04^′ W), a dry montane ponderosa pine (Pinus ponderosa) forest in the central Rocky Mountains, about 45 km west of Colorado Springs, CO, USA. MEF has a mean elevation of about 2400 m above sea level (a.s.l.) and an annual mean temperature of about 5 ^∘C. The annual precipitation is about 400 mm. Soils within MEF tend to have a low available water-holding capacity and moderately high permeability. The dominant parent materials of the soils within MEF are primarily Pikes Peak granite and secondarily weathered red arkosic sandstone.

The area surrounding the burn site is a large grassy opening that had been created in the surrounding ponderosa pine forest in 2001 when several trees were cut in an effort to reduce the amount of mistletoe in the area. The vegetation within this opening is predominantly grasses, forbs and shrubs (including some nonnative invasives). At the time of the experiment (fall 2003–spring 2004), this opening was covered primarily by senescent bunchgrasses. The soil within the general area of the burn pile is a deep (>1 m), fine loamy, mixed, frigid, Pachic Argiustoll (or a Luvic Phaeozem (Pachic); FAO UN, 2014) and is typical of soils throughout this experimental area. 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 Mg m⁻³. Soil organic matter comprises about 1 %–2 % of the soil by mass and is more or less uniform through at least the top 10 cm of soil (Jiménez Esquilín et al., 2008). Previous grazing and mechanical harvesting throughout the area has resulted in a moderately disturbed soil.

Figure 2 Layout of the data system and slash pile of the MEF 26 April 2004 experimental burn about 25 min after ignition. Image credit: John M. Frank.

Figure 3 The slash pile of the MEF 26 April 2004 experimental burn about 75 min after ignition. Image credit: John M. Frank.

Figure 4 The slash pile of the MEF 26 April 2004 experimental burn about 195 min after ignition. Image credit: John M. Frank.

3.2 The slash pile burn – description and instrumentation

The burn site instrumentation was installed in August 2003 at two control plots and two slash plots (under the center and the edge of the slash pile). Data include soil temperatures, soil moisture, soil heat flux and soil CO₂ at several soil depths. All these in situ sensors and their associated connectors and cables were buried several centimeters deep and connected to data loggers (CR23X data logger; Campbell Scientific, Inc., Logan, UT, USA) via a 27 m trench that had been backfilled to protect the data communications from the heat of the fire. The slash pile (located at 3.11439^∘ N, 105.10284^∘ W) was mechanically constructed in March 2004 and measured about 42 m in circumference, about 6 m in height and covered an elliptically shaped area of about 130 m². The fuel loading was estimated to be between 450 and 600 kg m⁻². The burn was initiated a few minutes after 10:00 Mountain daylight time (MDT) on 26 April 2004. Figures 2–4 are a sequence of photographs taken during the first 3 h or so of the burn. Figure 2 shows the burning slash pile about 25 min after the fire was initiated and shows the deployment of the data loggers, CO₂ pumps and analyzers and the supporting infrastructure. Figure 3 shows the burning slash pile about 75 min after initiation, and Fig. 4 was taken about 195 min after initiation. Note that although the soil CO₂ data 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.

3.2.1 Soil temperature

Soil temperatures were measured with thermocouples (Omega Engineering Inc., Norwalk, CT, USA) and sampled every 2 min during the fire and for the week following the day the fire was initiated. To insure electrical isolation, all thermocouple junctions were coated with epoxy (Omegabond 101) prior to insertion into the soil. Thermocouples were placed at the soil surface (0.00 m) and at 0.02, 0.05, 0.10, 0.15, 0.20 and 0.50 m depths. The four uppermost sensors, including the soil surface, were K-type (rated to 704 ^∘C), J-type (rated to 260 ^∘C), which was used at 0.15 m, and the bottom two depths were T-type (rated to 100 ^∘C) sensors.

3.2.2 Soil heat flux

The soil heat fluxes were measured at three depths (0.02, 0.10 and 0.20 m) under the center of the pile with a high-temperature probe (HFT) of AlumaCore with an exterior ceramic glaze (Thermonetics Corporation, La Jolla, CA, USA). They were also sampled at the same rate and time as the soil thermocouples. These HFTs are rated to 775 ^∘C and have a nominal sensitivity between 1250 and 1750 W m⁻² mV⁻¹. These HFTs were attached to a data logger by a chromel extension wire (Omega Engineering Inc.; TFCH-020; rated to 260 ^∘C).

Because these 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. But the details concerning the thermal conductivity of these HFTs are very relevant here and need highlighting. The thermal conductivity of these HFTs, λ_p (W m⁻¹ K⁻¹), is ${\mathit{\lambda }}_{\mathrm{p}}=\mathrm{0.7}+\mathrm{0.003}T$ (where T is degrees Celsius). Knowledge of λ_p is important to correct for the discrepancy between the true soil heat flux, G_soil (W m⁻²), 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 the Philip correction, is given as follows:

$$\begin{array}{}\text{(21)}& G_{\mathrm{soil}}=\left[\mathrm{1}-\mathit{\beta }r\left(\mathrm{1}-{\mathit{ϵ}}^{-\mathrm{1}}\right)\right]G_{\mathrm{m}},\end{array}$$

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 ${\mathit{ϵ}}^{-\mathrm{1}}={\mathit{\lambda }}_{\mathrm{s}}/{\mathit{\lambda }}_{\mathrm{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 likely to occur during a fire. But this correction also requires in situ knowledge of λ_s, which was not (and probably could never have been) measured during this or any fire. So the model's predicted λ_s is used in Eq. (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.

3.2.3 Soil moisture and CO₂

All soil moisture and CO₂ data were measured at 0.05 and 0.15 m depths at both the center and edge position under the slash pile. They were both sampled every 30 min for the week during and after the burn. Soil moisture was measured using a specially designed high-temperature time domain reflectometer (TDR; Zostrich Geotechnical, Oroville, 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 to continue operating and providing reliable data at temperatures exceeding 250 ^∘C. To ensure this, external portions of the coaxial cables that were likely to be exposed to such high temperatures were wrapped in silicon tape. The calibration factor for this TDR is temperature dependent and is discussed in more detail in the Appendix of Massman et al. (2010a). As 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₂ was measured by drawing a continuous sample for approximately 0.5 min through $\mathrm{3}/\mathrm{8}$ in. (9.53 mm internal diameter) Dekabon tubing into a LI-820 (LI-COR Inc., Lincoln, NE, USA) that was housed several meters from the slash pile, as shown in Fig. 2.

4 Model performance and sensitivity analysis

4.1 Surface energy balance

Figures 5 and 6 show the details of the surface energy balance for Eq. (16), the model's upper boundary condition that includes the sensible heat and the net infrared terms. Included in Fig. 5 is the forcing function, Q_F(t), which shows the shape of the BFD curve. Figure 7 shows the surface energy balance for Eq. (18), which excludes the sensible heat and the net infrared terms. Comparing Figs. 6 and 7 suggests that, once tuned appropriately, either formulation of the surface energy balance will give very similar simulations of the evaporative flux, $L_{\mathrm{v}\mathrm{0}}^{*}{E}_{\mathrm{0}}$. Likewise, during the period of soil heating (i.e., when G₀>0), both formulations give similar and reasonable simulations for G₀. But the simplified model of the surface energy balance, Eq. (18), does not capture (and in fact cannot capture) G₀ during the period when the soil is cooling (i.e., when G₀<0, which starts at about 22 h in Figs. 5 and 6). Without the possibility of radiative and convective cooling of the surface, Eq. (18) does not reproduce G₀<0.

Figure 5 Soil surface energy balance components for the 26 April 2004 MEF burn resulting from the full surface energy balance model (Eq. 16). The parameters for the forcing function, Q_F(t) or Eq. (15), are Q_in=120 W m⁻², Q_Fmax=18 kW m⁻², t_m=13.5 h and t_d=35 h. The net forcing is defined as the difference between the radiative terms, i.e., ${\mathit{ϵ}}_{\mathrm{0}}\left({\mathit{\theta }}_{\mathrm{0}}\right)Q_{\mathrm{F}}\left(t\right)-{\mathit{ϵ}}_{\mathrm{0}}\left({\mathit{\theta }}_{\mathrm{0}}\right)\mathit{\sigma }[T_{K\mathrm{0}}^{\mathrm{4}}-{\mathit{ϵ}}_{\mathrm{a}}({\mathit{\rho }}_{\mathrm{va}}\left)T_{\mathrm{Ka}}^{\mathrm{4}}\right]$, which, from Eq. (16), is equal to the sum of the three nonradiative terms, i.e., ${\mathit{\rho }}_{\mathrm{a}}{c}_{\mathrm{pa}}{C}_H[T_{\mathrm{0}}-T_{\mathrm{a}}]+L_{\mathrm{v}\mathrm{0}}^{*}{E}_{\mathrm{0}}+G_{\mathrm{0}}$. Figure 6 is an magnified version of the net forcing and these three energy components.

Figure 6 Replotted from Fig. 5 to facilitate comparison with the surface energy terms as simulated with the simplified surface energy balance, Eq. (18), and shown in Fig. 7.

Figure 7 Soil surface energy balance components for the 26 April 2004 MEF burn resulting from the simplified surface energy balance model (Eq. 18). The parameters for the forcing function, Q_F(t) or Eq. (15), are Q_in=0 W m⁻², Q_Fmax=2.7 kW m⁻², t_m=9.5 h and t_d=27.5 h. Here the net forcing is the sum $L_{\mathrm{v}\mathrm{0}}^{*}{E}_{\mathrm{0}}+G_{\mathrm{0}}$ because the surface sensible heat flux is not included as part the simplified surface energy balance.

Another important aspect about the surface energy balance that Eq. (18) does not capture as well as Eq. (16) is the time lag between the maximum in soil heat flux and the maximum in soil temperature. At 0.02 m depth, the observed time lag is about 5 h. Predictions of the time lag with the full surface energy balance, Eq. (16), agree almost exactly with the observed time lag. But the simplified model predicts a time lag of about 4.6 h. This difference in the time lags is not the result of the different choices of t_m, t_d or Q_Fmax between the two models of the surface energy balance. Instead, Eq. (18) inherently constrains the soil heat flux more strongly than Eq. (16).

These apparent limitations of the simplified surface energy balance, Eq. (18), do not lessen the argument that it may be appropriate for some slash pile burns or any time that burning material is in direct contact with the soil, and the concomitant soil heating is overwhelmingly conductive. Instead, these present comparisons suggest that a hybrid of Eqs. (18) and (16) may be more appropriate. Such a hybrid would employ Eq. (18) in the early part of the burn (before significant loss of mass from the slash pile due to combustion) and Eq. (16) later after the fire intensity has peaked (i.e., sometime after t_d), when the soil surface or possibly an ash surface is more exposed to the ambient environment. The time course of the experimental burn, shown in Figs. 2–4, is consistent with a hybrid formulation of the boundary forcing. Nonetheless, it is beyond the intent of the present paper to explore this concept any further.

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.

4.2 Soil temperature

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 Fig. 8, 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 Fig. 8 shows some slight discrepancies in the model's performance. (1) The model does not capture the early temperature rise beginning at about 4 h, (2) nor does it capture the secondary maximum temperature at about 20 h. (3) It appears to overestimate the maximum 2 cm soil temperatures (at about 21 h), and (4) the soil appears to cool off faster than observed (most obvious after about 28 h). Discrepancies (1) and (2) are not unexpected because it is impossible for a simple forcing function like the BFD curve and Eq. (15) to produce an exact or even a near-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, and it could not be easily generalized from one fire to the next. The (relatively significant) overestimation by nearly 60 ^∘C of the 2 cm soil temperatures is not fully understood, but it may be a consequence of a mismeasurement of the installation depth of the soil temperature sensors and/or poor contact with the soil, meaning that the soil air in contact with the temperature sensor would be acting as an insulator relative to the heat conducting soil particles in contact with the sensor. It is also possible that λ_s possess greater vertical structure than is included in the model. This could easily be the case for MEF soils because the present model of λ_s neither includes the observed vertical structure in soil bulk density nor its relationship to the vertical structure of the soil's thermal conductivity (Massman et al., 2008). Finally, the issue involving the difference between the observed and modeled rates of cooling, discrepancy (4), seems to be characteristic of all model simulations and not just the simulation shown in Fig. 8. Some of this is undoubtedly related to the simple shape of the BFD curve (forcing function) and the constraints it imposes on the upper boundary conditions. Another likely contributing factor to discrepancy (4), and the other three 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.

Figure 8 Soil temperatures during the 26 April 2004 MEF burn. Observed data are denoted with symbols, and modeling results (solid lines) are simulated with the full version of the soil surface energy balance (Eq. 16).

4.3 Soil heat flux

Comparisons between measured and modeled heat fluxes, shown in Fig. 9, are an independent check of the mathematical structure and tuning parameters of the surface forcing function, Q_F(t). For this study, Q_F(t) was tuned primarily for soil temperatures and secondarily for soil heat flux. Figure 9 compares the observed heat fluxes (blue color-filled area and symbols) to the modeled heat fluxes (solid lines). The upper boundary of the (2 cm) blue color-filled region is the heat flux measured by the heat flux plate without the Philip correction. The lower boundary of this region is the measured heat flux with the Philip correction, indicating that, for this experiment, the Philip correction reduced the amplitude of the uncorrected flux. The Philip correction is not shown for the two lower heat flux measurements (symbols) because it made virtually no change in the uncorrected fluxes. In general, the model does appear to capture many features of the observations. Of particular interest is the fact that the amplitude of the modeled heat flux more closely resembles 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 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.

Figure 9 Soil heat flux during the 26 April 2004 MEF burn, with modeled (solid lines) and observed (blue color-filled area and red and black symbols) measurements. The upper boundary of the (2 cm) blue color-filled region is the heat flux measured by the heat flux plate without the Philip correction (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.

4.4 Soil moisture

The modeled soil moisture is shown as a function of time in both Figs. 10 and 11 for the same depths (and color coding) as shown for temperature in Fig. 8. Included in each figure is the soil moisture measured at 5 cm (red) and 15 cm (magenta) depths, but Fig. 10 includes the temperature-corrected calibration of the soil moisture probe, whereas Fig. 11 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, 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 of the uncertainty inherent in these soil moisture measurements.

Figure 10 Volumetric soil moisture, after temperature correction to the TDR probe (Massman et al., 2010a, Appendix A), during the 26 April 2004 MEF burn, with modeled (solid lines) and observed (symbols) measurements. The color coding for depth is the same as that in Fig. 8, but observations of soil moisture are taken only at two depths, i.e., 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.

Figure 11 Volumetric soil moisture, without temperature correction to the TDR (Massman et al., 2010a, Appendix A) during the 26 April 2004 MEF burn, with modeled (solid lines) and observed (symbols) measurements. The color coding for depth is the same as that in Fig. 8, but observations of soil moisture are taken only at two depths, i.e., 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.

Figure 12 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 of than 5 cm, the soil moisture begins to evaporate (decrease) between 50 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 performed to see if the model could be made to better emulate the observed θ−T trajectory, including different formulations and parameter values for K_R; different parameter values for the forcing function, Q_F(t), and the simplified surface energy boundary condition, Eq. (18); and different WRCs, variations in the parameters of the source term, S_v, and removing the enhancement factor from the vapor diffusivity, D_ve and Eq. (4). Only changes in S_v and D_ve made any significant positive difference to the model's θ−T trajectory. Changes in either of these model factors can strongly influence the amount of vapor within the soil pores, which in turn influences 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 one), the resulting decrease in vapor diffusivity causes the vapor within the soil pores to increase such that it feeds back to soil moisture via the condensation term of S_v (proportional to ρ_v; Eq. 9), 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).

Figure 12 Modeled and observed θ−T trajectory. Observed data are denoted with symbols and modeling results are solid lines. The observed volumetric soil moisture has been temperature corrected.

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 K_R or changing the heating rates or 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 with respect 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, i.e., fidelity to the soil's thermal response to heating by fire and the soil's moisture and vapor response.

4.5 Dynamic feedbacks

Fire changes soil, and the more intense the fire and the greater amount of soil heating, the greater the changes in the soil will be. Although this statement should be self-evident, it is probably more precise to say that changes in the soil result from the integrated effects of the surface heat flux (delivered to the soil), the area of soil in contact with or influenced by the fire and the duration of the fire. That is, the probability of the change in soil due to fire is $\propto \underset{\mathrm{0}}{\overset{t_{\mathrm{d}}}{\int }}\stackrel{\mathrm{Area}}{\int }G(z=\mathrm{0},\mathrm{Area},t)\mathrm{d}\left(\mathrm{Area}\right)\mathrm{d}t$. This more precise approach clearly suggests that changes in the soil are more likely to occur during spatially stationary fires (e.g., slash pile burns or relatively low-intensity but long-duration biomass burning) than during a fast-moving, dynamic fire (such as a grassland fire or a low shrubland fire or even an extremely high-intensity crown fire).

Examples of some of the changes that are relevant to the present study include (a) changes in the soil bulk density (Gregory Butters, personal communication, 2009; Kojima et al., 2018), which will alter the WRC (Fig. 1) and hydraulic properties (Tian et al., 2018), (b) changes in the thermal conductivity of the soil (Massman et al., 2008; Kojima et al., 2018) and the volumetric specific heat (Gregory Butters, personal communication, 2009), (c) changes in the heat and vapor fluxes resulting from (a) and (b) (Kojima et al., 2018), and (d) changes in 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 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 the soil heat and moisture transport, this section details the results of a model simulation based on observed or inferred changes in 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 to 1.46 Mg m⁻³ (a 12 % increase, as per Fig. 1); simultaneously, the thermal conductivity of the mineral fraction, λ_m0, increases from 4.42 to 8 Wm K⁻¹; the de Vries shape factor, g_a, decreased from 0.123 to 0.06; the Campbell et al. (1994) parameter q_w0 (which determines when the 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 Gregory Butters, personal communication, 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 soil thermal conductivity, the WRC, the hydraulic conductivity function and the source term, S_v. On its own this decrease in η yields an increase in S_v. But under extreme drying it is more reasonable to expect A_wa (and S_v) to decrease (Chen et al., 2018); so, to compensate, S_* is reduced when η is reduced. 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 in 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, (i.e., $\partial {\mathit{\lambda }}_{\mathrm{s}}/\partial \mathit{\theta }$ decreased as a result of the fire) and increased λ_s when the soil was 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 fire-affected soils. The conductivity of the mineral fraction, λ_m, was increased to capture the creation of thin, but highly conductive, coatings of mineral oxides (MnO₂ 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., O'Brian et al., 2013).

The base case model simulation (Figs. 8–12) and the feedback simulation are compared in Figs. 13–15. The upper boundaries of the color-filled areas in Fig. 13 (temperature) and Fig. 14 (heat flux) correspond to the feedback simulation, whereas the lower boundary in Fig. 15 (soil volumetric moisture) is associated with the feedback simulation. These comparisons demonstrate that the feedbacks overwhelmingly act to increase the soil temperatures and heat fluxes throughout the 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 feeds back (nonlinearly) upon itself, creating conditions that allow the heat to penetrate more deeply into the soil and allowing the soil temperatures to exceed what they would have been without the feedbacks.

Figure 13 Potential impacts of dynamic feedbacks on the soil temperatures during the 26 April 2004 MEF burn. Observed data are denoted with symbols and modeling results are solid filled areas. The lower boundary of the solid filled areas corresponds to the base case simulation and the upper boundary to the feedback simulation.

Figure 14 Potential impacts of dynamic feedbacks on the soil heat flux during the 26 April 2004 MEF burn. The observed data are shown as the (2 cm) darker blue color-filled area and red and black symbols. The upper boundary of the observed heat flux includes the Philip correction, whereas 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. The upper boundary of the light blue, red and black color-filled areas corresponds to the feedback simulation and the lower boundary to the base case simulation.

Figure 15 Potential impacts of dynamic feedbacks on the volumetric soil moisture during the 26 April 2004 MEF burn. Observed data, which has been temperature corrected, are denoted with symbols, and modeling results are solid filled areas. The upper boundary of the solid filled areas corresponds to the base case simulation. The lower boundary corresponds to the feedback simulation.

5 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 these should be useful to consider for future studies of the soil's 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 is needed. Also, there needs to be 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 than just including ρ_b because the observed increase in λ_s for a dry 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 volatilization of the soil moisture, i.e., u_vl and Eq. (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. 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. The results of Zeng et al. (2011) for less extreme conditions support this notion. 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 Eq. (1). This energy transport term has often been included when modeling the daily cycle of energy flow through soils. Because some fires, especially slash pile burns, can continue for a week or more, it seems appropriate to investigate the influence that this type of energy transport could have on model solutions.

• c.

  The 2D and 3D effects – there are (at least) two physical processes that cannot be fully represented in a 1D model, i.e., (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, b; Nobles et al., 2010). (i) In an earlier section, mention was made of the second installation of soil sensors at the edge of the pile. With 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 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 it further underscores the nature of the challenges and empiricism inherent in the 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 2D or 3D in nature and it appears to carry combustion products into the soil (Massman et al., 2010b, Fig. 4). Such an advective current is hypothesized to have caused the extremely rapid 20-fold increase in soil CO₂ which was observed to have occurred within the first 30 min 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 overwhelm 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 1D model like the HMV model, but such an adjustment should probably be guided by observational studies and experiments designed to establish the existence, nature and dynamics of u_adv.

6 Concluding summary

This study describes the continuing development, improvement and validation of the HMV model, which is a nonequilibrium model of the coupled transport of heat, moisture and water vapor in soils during surface fires. The purpose of the research supporting 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. In fact, a version of the HMV model is currently available to land managers through the USDA Forest Service's First Order Fire Effects Model (FOFEM) version 6.7 (USDA Forest Service, 2020). Key improvements to the model, which noticeably improved its stability and performance, include more physically realistic parameterizations of the nonequilibrium 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 discussions of the complexities and difficulties regarding formulating the surface (or upper) boundary condition. The model is validated using in situ measurements of soil temperatures, heat flux and soil moisture obtained during a 2004 experimental burn carried out at MEF (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 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 in 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 to 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 2D or 3D advective flows in soils induced by fires and the possibility of 2D or 3D heat flow.

Appendix A: Soil apparent thermal conductivity

A1 Mass gradient diffusional flux ($-D_{\mathrm{ve}}\partial {\mathit{\rho }}_{\mathrm{v}}/\partial z$)

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 the modeling of soil heat and 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 Eqs. (6), (7) and (4). But there is a valid alternative expression for Eq. (6), which results from the combination of the conservation of enthalpy, Eq. (1), with the conservation of water vapor, Eq. (4), and the equilibrium assumption (${\mathit{\rho }}_{\mathrm{v}}={\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}=a_{\mathrm{w}}{\mathit{\rho }}_{\mathrm{v},\mathrm{sat}}$), and is expressed as follows:

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}{C}_{\mathrm{s}}{\displaystyle \frac{\partial T}{\partial t}}& -{\displaystyle \frac{\partial }{\partial z}}\left[{\mathit{\lambda }}_{\mathrm{s}}{\displaystyle \frac{\partial T}{\partial z}}\right]+L_{\mathrm{v}}^{*}{\displaystyle \frac{\partial (\mathit{\eta }-\mathit{\theta }){\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}}{\partial t}}\\ & -L_{\mathrm{v}}^{*}{\displaystyle \frac{\partial }{\partial z}}\left[D_{\mathrm{ve}}{\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}}{\partial z}}\right]=\mathrm{0},\end{array}\end{array}$$

which, after some simple mathematical manipulation, can also be written as follows:

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}{C}_{\mathrm{s}}{\displaystyle \frac{\partial T}{\partial t}}& -L_{\mathrm{v}}^{*}{\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}{\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}+L_{\mathrm{v}}^{*}(\mathit{\eta }-\mathit{\theta }){\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}}{\partial t}}-{\displaystyle \frac{\partial }{\partial z}}\left[{\mathit{\lambda }}_{\mathrm{s}}{\displaystyle \frac{\partial T}{\partial z}}\right.\\ & \left.+L_{\mathrm{v}}^{*}{D}_{\mathrm{ve}}{\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}}{\partial z}}\right]+D_{\mathrm{ve}}{\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}}{\partial z}}{\displaystyle \frac{\partial L_{\mathrm{v}}^{*}}{\partial z}}=\mathrm{0},\end{array}\end{array}$$

where the last term on the right-hand side of this equation accounts for moving $L_{\mathrm{v}}^{*}$ inside the gradient operator $\partial /\partial z$. As discussed in the main text, $a_{\mathrm{w}}=a_{\mathrm{w}}(T_{\mathrm{K}},\mathit{\psi })$ is modeled by the Kelvin equation, $a_{\mathrm{w}}=e^{\frac{M_{\mathrm{w}}{\mathit{\psi }}_{*}}{R{T}_{\mathrm{K}}}{\mathit{\psi }}_n}$ and ${\mathit{\rho }}_{\mathrm{v},\mathrm{sat}}={\mathit{\rho }}_{\mathrm{v},\mathrm{sat}}\left(T_{\mathrm{K}}\right)$, so that $\partial {\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}/\partial t$ can be expanded in terms of $\partial T/\partial t$, $\partial {\mathit{\psi }}_n/\partial t$ and $\partial {\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}/\partial z$ in terms of $\partial T/\partial z$ and $\partial {\mathit{\psi }}_{\mathrm{n}}/\partial z$. The following shows $\partial {\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}/\partial z$. (Note that $\partial {\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}/\partial t$ can be found by substituting t for z).

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}{\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}}{\partial z}}& =\left[a_{\mathrm{w}}{\mathrm{\Delta }}_{\mathrm{sat}}-a_{\mathrm{w}}\left({\displaystyle \frac{M_{\mathrm{w}}{\mathit{\psi }}_{*}{\mathit{\psi }}_n}{R{T}_{\mathrm{K}}}}\right){\displaystyle \frac{{\mathit{\rho }}_{\mathrm{v},\mathrm{sat}}}{T_{\mathrm{K}}}}\right]{\displaystyle \frac{\partial T}{\partial z}}\\ & +\left[a_{\mathrm{w}}\left({\displaystyle \frac{M_{\mathrm{w}}{\mathit{\psi }}_{*}}{R{T}_{\mathrm{K}}}}\right){\mathit{\rho }}_{\mathrm{v},\mathrm{sat}}\right]{\displaystyle \frac{\partial {\mathit{\psi }}_n}{\partial z}},\end{array}\end{array}$$

where Δ_sat (kg m⁻³ K⁻¹) is the slope of the saturation vapor curve, $\mathrm{d}{\mathit{\rho }}_{\mathrm{v},\mathrm{sat}}/\mathrm{d}T$, and $\mathrm{7}\times {\mathrm{10}}^{-\mathrm{4}}<{\mathrm{\Delta }}_{\mathrm{sat}}<\mathrm{70}\times {\mathrm{10}}^{-\mathrm{4}}$ is for temperatures between about 5 and 60 ^∘C. Furthermore, unless the soil is extremely dry (${\mathit{\psi }}_n\approx \mathrm{1})-[\left(M_{\mathrm{w}}{\mathit{\psi }}_{*}{\mathit{\psi }}_n\right)/\left(R{T}_{\mathrm{K}}\right)\left]\right[{\mathit{\rho }}_{\mathrm{v},\mathrm{sat}}/T_{\mathrm{K}}]\ll {\mathrm{\Delta }}_{\mathrm{sat}}$, implying that the second term that multiplies the temperature gradient of Eq. (A3) can be ignored relative to the first term for most applications. Otherwise, for an extremely dry soil, $-\left[\right(M_{\mathrm{w}}{\mathit{\psi }}_{*}\mathit{\psi }n)/(R{T}_{\mathrm{K}}\left)\right][{\mathit{\rho }}_{\mathrm{v},\mathrm{sat}}/T_{\mathrm{K}}]\approx \mathrm{0.5}{\mathrm{\Delta }}_{\mathrm{sat}}$. Henceforth, ${\mathrm{\Delta }}_{\mathrm{sat}}^{*}\equiv {\mathrm{\Delta }}_{\mathrm{sat}}-\left[\right(M_{\mathrm{w}}{\mathit{\psi }}_{*}{\mathit{\psi }}_n)/(R{T}_{\mathrm{K}}\left)\right][{\mathit{\rho }}_{\mathrm{v},\mathrm{sat}}/T_{\mathrm{K}}]$. In addition, it is always true that $({\mathit{\psi }}_n/T_{\mathrm{K}})\partial T/\partial z\gg \partial {\mathit{\psi }}_n/\partial z$, which means that the rightmost term in Eq. (A3) can also be ignored. Therefore, Eq. (A3) is now as follows:

$$\begin{array}{}\text{(A4)}& {\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}}{\partial z}}=a_{\mathrm{w}}{\mathrm{\Delta }}_{\mathrm{sat}}^{*}{\displaystyle \frac{\partial T}{\partial z}}.\end{array}$$

Substituting Eq. (A4) into Eq. (A2) yields the following:

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}{C}_{\mathrm{s}}{\displaystyle \frac{\partial T}{\partial t}}& -L_{\mathrm{v}}^{*}{\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}{\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}-{\displaystyle \frac{\partial }{\partial z}}\left[\left({\mathit{\lambda }}_{\mathrm{s}}+a_{\mathrm{w}}{L}_{\mathrm{v}}^{*}{D}_{\mathrm{ve}}{\mathrm{\Delta }}_{\mathrm{sat}}^{*}\right){\displaystyle \frac{\partial T}{\partial z}}\right]\\ & +D_{\mathrm{ve}}{\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}}{\partial z}}{\displaystyle \frac{\partial L_{\mathrm{v}}^{*}}{\partial z}}=\mathrm{0},\end{array}\end{array}$$

where the term $L_{\mathrm{v}}^{*}(\mathit{\eta }-\mathit{\theta })\partial {\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}/\partial t$ has been dropped from Eq. (A2) because $\left|L_{\mathrm{v}}^{*}\right(\mathit{\eta }-\mathit{\theta })\partial {\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}/\partial t|\ll |C_{\mathrm{s}}\partial T/\partial t|$ (although it is a bit tedious to show). Equation (A5) yields the following identification for λ_s,app:

$$\begin{array}{}\text{(A6)}& {\mathit{\lambda }}_{\mathrm{s},\mathrm{app}}={\mathit{\lambda }}_{\mathrm{s}}+a_{\mathrm{w}}{L}_{\mathrm{v}}^{*}{D}_{\mathrm{ve}}{\mathrm{\Delta }}_{\mathrm{sat}}^{*}={\mathit{\lambda }}_{\mathrm{s}}+{\mathit{\lambda }}_{\mathrm{s},\mathrm{dis}},\end{array}$$

where the term describing the effect of the vapor distillation due to temperature gradients on the soil's thermal conductivity is $a_{\mathrm{w}}{L}_{\mathrm{v}}^{*}{D}_{\mathrm{ve}}{\mathrm{\Delta }}_{\mathrm{sat}}^{*}$ (identified in Eq. (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 the conservation of enthalpy in soils (e.g., Hillel, 2004, Eq. (12.24), p. 226; 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, Eq. (A4), into Eq. (7). This yields the following:

$$\begin{array}{}\text{(A7)}& \begin{array}{rl}\left({\mathit{\rho }}_{\mathrm{w}}\right.& \left.-{\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}\right){\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}+(\mathit{\eta }-\mathit{\theta }){\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}}{\partial t}}\\ & +{\displaystyle \frac{\partial }{\partial z}}\left[a_{\mathrm{w}}{D}_{\mathrm{ve}}{\mathrm{\Delta }}_{\mathrm{sat}}^{*}{\displaystyle \frac{\partial T}{\partial z}}\right]=\mathrm{0},\end{array}\end{array}$$

or to a very good approximation (which follows from Eq. (A7) because ${\mathit{\rho }}_{\mathrm{w}}\gg {\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}$), as follows:

$$\begin{array}{}\text{(A8)}& {\mathit{\rho }}_{\mathrm{w}}{\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}+(\mathit{\eta }-\mathit{\theta }){\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}}{\partial t}}+{\displaystyle \frac{\partial }{\partial z}}\left[a_{\mathrm{w}}{D}_{\mathrm{ve}}{\mathrm{\Delta }}_{\mathrm{sat}}^{*}{\displaystyle \frac{\partial T}{\partial z}}\right]=\mathrm{0}.\end{array}$$

Equations (A5) and (A8) now form an equilibrium model with two independent predictive variables (T and θ or ψ) that describe 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, be 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 and vapor; Eq. (A8) above) with the simplified version of the conservation of energy equation, Eq. (6), as discussed in the development of the nonequilibrium model. This simplified equation is as follows:

$$\begin{array}{}\text{(A9)}& C_{\mathrm{s}}{\displaystyle \frac{\partial T}{\partial t}}-{\displaystyle \frac{\partial }{\partial z}}\left[{\mathit{\lambda }}_{\mathrm{s}}{\displaystyle \frac{\partial T}{\partial z}}\right]-L_{\mathrm{v}}^{*}{\mathit{\rho }}_{\mathrm{w}}{\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}=\mathrm{0}.\end{array}$$

It is important to reiterate that the paired Eqs. (A8) and (A5) are as valid a basis for describing coupled heat and moisture flow in soils, as are the paired Eqs. (A8) and (A9), although Eqs. (A8) and (A9) are much more extensively used. But the conservation of energy equation, as expressed with Eq. (A9), does not explicitly include the effect of the vapor distillation due to temperature gradients (λ_s,dis), whereas Eq. (A5) does. Instead, the effects of the distillation of water vapor are implicit in Eq. (A9) because they are included in $\partial \mathit{\theta }/\partial t$ (via Eq. A8). Therefore, to substitute λ_s,app for λ_s in Eq. (A9), or any θ-based variant of Eq. (A9) (as opposed to a ρ_v-based variant like Eq. 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 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 Eq. (A5), which already explicitly includes λ_s,dis.

The nonequilibrium form of λ_s,app follows from the same general methodology as the equilibrium form from Eq. (A5) does. First, combining Eqs. (1) and (4) and then simplifying yields the following:

$$\begin{array}{}\text{(A10)}& \begin{array}{rl}{C}_{\mathrm{s}}{\displaystyle \frac{\partial T}{\partial t}}& -L_{\mathrm{v}}^{*}{\mathit{\rho }}_{\mathrm{v}}{\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}-{\displaystyle \frac{\partial }{\partial z}}\left[{\mathit{\lambda }}_{\mathrm{s}}{\displaystyle \frac{\partial T}{\partial z}}-L_{\mathrm{v}}^{*}{D}_{\mathrm{ve}}{\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{v}}}{\partial z}}\right]\\ & +D_{\mathrm{ve}}{\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{v}}}{\partial z}}{\displaystyle \frac{\partial L_{\mathrm{v}}^{*}}{\partial z}}=\mathrm{0}.\end{array}\end{array}$$

Next, employing the ideal gas law for ρ_v, i.e., $e_{\mathrm{v}}={\mathit{\rho }}_{\mathrm{v}}R{T}_{\mathrm{K}}/M_{\mathrm{w}}$, yields the following:

$$\begin{array}{}\text{(A11)}& {\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{v}}}{\partial z}}=-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{v}}}{T_{\mathrm{K}}}}{\displaystyle \frac{\partial T}{\partial z}}+{\displaystyle \frac{M_{\mathrm{w}}}{R{T}_{\mathrm{K}}}}{\displaystyle \frac{\partial e_{\mathrm{v}}}{\partial z}}.\end{array}$$

Finally, introducing Eq. (A11) into Eq. (A10) yields the following:

$$\begin{array}{}\text{(A12)}& \begin{array}{rl}{C}_{\mathrm{s}}{\displaystyle \frac{\partial T}{\partial t}}& -L_{\mathrm{v}}^{*}{\mathit{\rho }}_{\mathrm{v}}{\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}-{\displaystyle \frac{\partial }{\partial z}}\left[\left({\mathit{\lambda }}_{\mathrm{s}}+L_{\mathrm{v}}^{*}{D}_{\mathrm{ve}}{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{v}}}{T_{\mathrm{K}}}}\right){\displaystyle \frac{\partial T}{\partial z}}\right]\\ & +{\displaystyle \frac{\partial }{\partial z}}\left[\left({\displaystyle \frac{L_{\mathrm{v}}^{*}{D}_{\mathrm{ve}}{M}_{\mathrm{w}}}{R{T}_{\mathrm{K}}}}\right){\displaystyle \frac{\partial e_{\mathrm{v}}}{\partial z}}\right]+D_{\mathrm{ve}}{\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{v}}}{\partial z}}{\displaystyle \frac{\partial L_{\mathrm{v}}^{*}}{\partial z}}=\mathrm{0},\end{array}\end{array}$$

where the nonequilibrium form of λ_s,dis is now $L_{\mathrm{v}}^{*}{D}_{\mathrm{ve}}{\mathit{\rho }}_{\mathrm{v}}/T_{\mathrm{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 nonequilibrium 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.

To summarize and conclude, the conservation of enthalpy (Eq. 1, or one of its many derivatives, e.g., Eqs. 6, A5 and 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_{\mathrm{v}}^{*}{S}_{\mathrm{v}}$. Therefore, all effects associated with the distillation of water are a result of S_v. This is true whether discussing either an equilibrium or nonequilibrium model. Imposing or assuming liquid/vapor equilibrium or nonequilibrium has mathematical consequences with regard 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, among which there are soil volumetric water content and soil vapor density, i.e., ${\mathit{\lambda }}_{\mathrm{s}}={\mathit{\lambda }}_{\mathrm{s}}(\mathit{\theta },{\mathit{\rho }}_{\mathrm{v}},\mathrm{\cdots })$ (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 (Eq. 1). To do so violates the conservation of enthalpy by effectively double counting λ_s,dis.

A2 Mass mixing ratio diffusional flux ($-{\mathit{\rho }}_{\mathrm{a}}{D}_{\mathrm{ve}}\partial {\mathit{\chi }}_{\mathrm{v}}/\partial z$)

The discussion in this section of the appendix complements the discussion of λ_app and λ_s,dis in Sect. A1 above. It does not change the final outcome or conclusions reached in Sect. A1. Section 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 Sect. A1 above develops λ_s,app more or less along the traditional lines, using the mass gradient form of the diffusional flux, $-D_{\mathrm{ve}}\partial {\mathit{\rho }}_{\mathrm{v}}/\partial z$. But the diffusional flux is also represented (and sometimes more appropriately represented) in terms of the gradient of the mass mixing ratio, i.e., $-{\mathit{\rho }}_{\mathrm{a}}{D}_{\mathrm{ve}}\partial {\mathit{\chi }}_{\mathrm{v}}/\partial z$, where ρ_a (kg m⁻³) is the total (dry air and vapor) soil gas density and χ_v (kg kg⁻¹) = ${\mathit{\rho }}_{\mathrm{v}}/{\mathit{\rho }}_{\mathrm{a}}$. This section of the Appendix shows that, for a relatively moist soil (i.e., a_w≈1 and ${\mathit{\chi }}_{\mathrm{v},\mathrm{eq}}\approx {\mathit{\chi }}_{\mathrm{v},\mathrm{sat}}$), λ_s,dis developed in Sect. A1 above differs between +2.5 % (for T≈5 ^∘C) and −4.5 % (for temperatures T≈60 ^∘C) from that given in Eq. (A6). All that this requires is a demonstration of how ${\mathit{\rho }}_{\mathrm{a}}\partial {\mathit{\chi }}_{\mathrm{v},\mathrm{eq}}/\partial z$ generalizes Eq. (A4).

From the identity ${\mathit{\rho }}_{\mathrm{a}}={\mathit{\rho }}_{\mathrm{d}}+{\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}$, where ρ_d (kg m⁻³) is the soil dry air density, it follows that:

$$\begin{array}{}\text{(A13)}& \begin{array}{rl}{\mathit{\rho }}_{\mathrm{a}}{\displaystyle \frac{\partial {\mathit{\chi }}_{\mathrm{v},\mathrm{eq}}}{\partial z}}& =\left({\mathit{\rho }}_{\mathrm{d}}+{\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}\right){\displaystyle \frac{\partial }{\partial z}}\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}}{{\mathit{\rho }}_{\mathrm{d}}+{\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}}}\right)\\ & ={\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}}{\partial z}}-{\mathit{\chi }}_{\mathrm{v},\mathrm{eq}}{\displaystyle \frac{\partial \left({\mathit{\rho }}_{\mathrm{d}}+{\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}\right)}{\partial z}}.\end{array}\end{array}$$

Next, assuming that ρ_d is an ideal gas, it follows that:

$$\begin{array}{}\text{(A14)}& {\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{d}}}{\partial z}}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{d}}}{p_{\mathrm{d}}}}{\displaystyle \frac{\partial p_{\mathrm{d}}}{\partial z}}-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{d}}}{T_{\mathrm{K}}}}{\displaystyle \frac{\partial T}{\partial z}}.\end{array}$$

Because the pressure gradient term, $({\mathit{\rho }}_{\mathrm{d}}/p_{\mathrm{d}})\partial p_{\mathrm{d}}/\partial z$, is not relevant for the present purposes (and it is unlikely to contribute much to $\partial {\mathit{\rho }}_{\mathrm{d}}/\partial z$ anyway), it can be dropped from Eq. (A14). This yields the following:

$$\begin{array}{}\text{(A15)}& {\mathit{\rho }}_{\mathrm{a}}{\displaystyle \frac{\partial {\mathit{\chi }}_{\mathrm{v},\mathrm{eq}}}{\partial z}}=\left(\mathrm{1}-{\mathit{\chi }}_{\mathrm{v},\mathrm{eq}}\right){\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}}{\partial z}}+{\mathit{\chi }}_{\mathrm{v},\mathrm{eq}}{\displaystyle \frac{\left({\mathit{\rho }}_{\mathrm{a}}-{\mathit{\rho }}_{\mathrm{v},\mathrm{eq}}\right)}{T_{\mathrm{K}}}}{\displaystyle \frac{\partial T}{\partial z}},\end{array}$$

which, after combining with Eq. (A4) and some further simplification, yields the following:

$$\begin{array}{}\text{(A16)}& {\mathit{\rho }}_{\mathrm{a}}{\displaystyle \frac{\partial {\mathit{\chi }}_{\mathrm{v},\mathrm{eq}}}{\partial z}}=\left(\mathrm{1}-{\mathit{\chi }}_{\mathrm{v},\mathrm{eq}}\right)a_{\mathrm{w}}\left({\mathrm{\Delta }}_{\mathrm{sat}}^{*}+{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{v},\mathrm{sat}}}{T_{\mathrm{K}}}}\right){\displaystyle \frac{\partial T}{\partial z}},\end{array}$$

from which it follows that:

$$\begin{array}{}\text{(A17)}& {\mathit{\lambda }}_{\mathrm{s},\mathrm{dis}}=L_{\mathrm{v}}^{*}{D}_{\mathrm{ve}}{a}_{\mathrm{w}}\left(\mathrm{1}-{\mathit{\chi }}_{\mathrm{v},\mathrm{eq}}\right)\left({\mathrm{\Delta }}_{\mathrm{sat}}^{*}+{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{v},\mathrm{sat}}}{T_{\mathrm{K}}}}\right).\end{array}$$

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 ($\mathrm{1}-{\mathit{\chi }}_{\mathrm{v},\mathrm{eq}}$) instead of 1 and ${\mathrm{\Delta }}_{\mathrm{sat}}^{*}+{\mathit{\rho }}_{\mathrm{v},\mathrm{sat}}/T_{\mathrm{K}}$ instead of ${\mathrm{\Delta }}_{\mathrm{sat}}^{*}$. At about 5 ^∘C, χ_v≈0.005 and ${\mathit{\rho }}_{\mathrm{v},\mathrm{sat}}/T_{\mathrm{K}}\approx \mathrm{0.03}{\mathrm{\Delta }}_{\mathrm{sat}}^{*}$, whereas, at 60 ^∘C, χ_v≈0.10 and ${\mathit{\rho }}_{\mathrm{v},\mathrm{sat}}/T_{\mathrm{K}}\approx \mathrm{0.06}{\mathrm{\Delta }}_{\mathrm{sat}}^{*}$. These results imply that, at 5 ^∘C, Eq. (A17) will yield a value for λ_s,dis that is about 2.5 % higher than ${\mathit{\lambda }}_{\mathrm{s},\mathrm{dis}}=L_{\mathrm{v}}^{*}{D}_{\mathrm{ve}}{a}_{\mathrm{w}}{\mathrm{\Delta }}_{\mathrm{sat}}^{*}$ and that, near 60 ^∘C, Eq. (A17) yields a value for λ_s,dis that is about 4.5 % lower. Unless a correction between +2 % and −4.5 % is important for estimating λ_s,dis, it seems that the mass mixing ratio formulation for the diffusional mass flux (Eq. A17) adds very little value to the original formulation (Eq. A6).