A more accurate scheme for calculating Earth’s skin temperature

Abstract

The theoretical framework of the vertical discretization of a ground column for calculating Earth’s skin temperature is presented. The suggested discretization is derived from the evenly heat-content discretization with the optimal effective thickness for layer-temperature simulation. For the same level number, the suggested discretization is more accurate in skin temperature as well as surface ground heat flux simulations than those used in some state-of-the-art models. A proposed scheme (“op(3,2,0)”) can reduce the normalized root–mean–square error (or RMSE/STD ratio) of the calculated surface ground heat flux of a cropland site significantly to 2% (or 0.9 W m⁻²), from 11% (or 5 W m⁻²) by a 5-layer scheme used in ECMWF, from 19% (or 8 W m⁻²) by a 5-layer scheme used in ECHAM, and from 74% (or 32 W m⁻²) by a single-layer scheme used in the UCLA GCM. Better accuracy can be achieved by including more layers to the vertical discretization. Similar improvements are expected for other locations with different land types since the numerical error is inherited into the models for all the land types. The proposed scheme can be easily implemented into state-of-the-art climate models for the temperature simulation of snow, ice and soil.

1 Introduction

Land/ocean skin temperature (ST) is an important parameter for quantifying the energy and water vapor exchanges between land/ocean and the atmosphere. In addition, accurate predictions of the time of temperatures reaching melting point for snow surface and ice surface are important factors for determining the times of snow and ice to melt and water to freeze (e.g., Ek et al. 2003). Despite this importance, many state-of-the-art models fail to accurately predict those times of phase changes. For example, ECHAM (Arpe et al. 1994) and SNOBAL (Link and Marks 1999) lag-predict the snow melting time. Furthermore, NCEP (Kanamitsu et al. 2002) and ECMWF (Simmons and Gibson 2000) are unable to predict the ice periods for the ocean grids in the 40–70° latitude zone (Tsuang et al. 2008). Since 1992, the Project for the Intercomparison of Land Surface Parameterization Schemes (PILPS) (Henderson-Seller 1995) has been systematically intercomparing the performance of land surface schemes (LSSs) (Schlosser 2000; Slater et al. 2001; Luo et al. 2003; van den Hurk and Viterbo 2003; Bowling et al. 2003). By contrast, this study tries to conduct a theoretical analysis of the discretization of the skin layer for various LSSs.

This study proposes the following numerical discretizations for calculating ST (denoted as T ₀), ground temperature at level k (T _k) and the bottom-layer numerical ground temperature (T _m) (Fig. 1) as:

Fig. 1

Grid structure for a ground column of the “op” scheme proposed by this study, where z _k is vertical coordinate of ground temperature T _k . The variables h _e0, h _ek and h _em are the effective thicknesses of the skin layer, of the layer k, and the bottom layer, respectively. The variables D _k , G _k and T _sk are heat diffusivity, ground heat flux and upper boundary ground temperature of level k at ζ = 0, respectively. The origin (ζ = 0) of the vertical coordinate ζ is located at the location of ground heat flux G _k . Of the ζ coordinate, T _k is at ζ = −h _tk and 0 ≤ h _tk  ≤ h _k . The optimal effective thicknesses h _e0, h _ek and h _em are calculated according to Eq. 20, derived by minimizing the error for T ₀, T _k and T _m simulations, respectively. It can be seen that the effective thickness h _e0 of the skin layer is determined under the conditions that the level index k = 0, h _t0 = 0 and h ₀ = −z ₁/2, and the effective thickness h _ek of layer k is determined under the conditions that h _tk  = (z _k−1–z _k )/2 and h _k  = (z _k−1–z _k+1)/2

$$ \left\{ {\begin{array}{*{20}c} {\frac{{\partial T_{0,n} \left( {z = 0,\;t} \right)}}{\partial t} = - \frac{{G_{0,n} - G_{1,n} \left( {z = - h_{0} ,\;t} \right)}}{{\rho_{g} c_{g} h_{e0} }}} \\ \vdots \\ {\frac{{\partial T_{k,n} \left( {z = z_{k} ,\;t} \right)}}{\partial t} = - \frac{{G_{k,n} - G_{k + 1,n} }}{{\rho_{g} c_{g} h_{ek} }}} \\ \vdots \\ {\frac{{\partial T_{m,n} \left( {z = z_{m} ,\;t} \right)}}{\partial t} = - \frac{{G_{m,n} - G_{m + 1,n} \left( {z = - \infty ,\;t} \right)}}{{\rho_{g} c_{g} h_{em} }} = - \frac{{G_{m,n} }}{{\rho_{g} c_{g} h_{em} }}} \\ \end{array} } \right. $$

(1)

where T _0,n , T _k,n and T _m,n are calculated skin temperature, ground temperature at level k and bottom-layer numerical ground temperature (K), respectively; t is time (s); G _k,n are calculated ground heat flux at level k (W m⁻²) (positive upward); ρ _g and c _g are density (kg m⁻³) and specific heat of the surface (J kg⁻¹ K⁻¹), respectively. The subscript “0” denotes a property at the surface, the subscript “n” denotes a property solved by a numerical method, and the subscript “m” denotes a property at the bottom level. The variable h ₀ is the physical thickness of the skin layer, h _e0 is the effective thickness of the skin layer, h _ek is the effective thickness of the level-k numerical layer, and h _em is the effective thickness of the bottom layer. Note that the physical thickness of the bottom layer is infinity, and the ground heat flux approaches zero at the infinite depth, i.e., G _m+1,n (z = −∞, t) = 0 (Carslaw and Jaeger 1959). Setting G _m+1,n  = 0 is also adopted by a few state-of-the-art climate models, such as by ECMWF (Viterbo and Beljaars 1995) and ECHAM (Roeckner et al. 2003). It is a good boundary condition since surface energy budget can be closed, when compared to a force-restored scheme used in some other climate models (Tsuang 2005; Tsuang et al. 2008). Nonetheless, how to determine h _em is not well documented. Note that ρ _g c _g h _e0 (J m⁻² K⁻¹) is also called surface area heat capacity (Grell et al. 1995) or the heat capacity of the surface layer (Roeckner et al. 2003). Similarly, we name ρ _g c _g h _ek as level-k area heat capacity (J m⁻² K⁻¹).

This study tries to derive optimal value for the effective thickness h _ek by minimizing the error for the temperature simulation of T _k . The conventional finite-difference scheme (denoted as “cv”) assumes that the ST is equal to the mean-layer temperature of the uppermost finite-difference layer (e.g., Tsuang and Dracup 1990; Gaspar et al. 1990; Blondin 1991; Hostetler et al. 1993; Sellers et al. 1996; Chia and Wu 1998; Link and Marks 1999), which can be calculated by the first equation in Eq. 1 by setting h _e0 = h ₀. Nonetheless, it is a well known fact that the calculated amplitude of the ST is dependent on vertical resolution, and the amplitude usually decreases with h ₀ (e.g., Gaspar et al. 1990) unless the thickness of the uppermost numerical layer is discretized infinitely thin (e.g., Mote and O’Neill 2000). Consequently, the diurnal cycle of the surface temperature is poorly simulated. Hence, a time lag for determining the timing of snow melt may be produced. However, an infinitely thin thickness is numerically impossible; and the thinner the thickness, the more the number of numerical layers required and the higher the computational cost. On the other hand, Viterbo and Beljaars (1995) set h _e0 at 0 by assuming that the skin layer has no heat capacity (denoted as “nh” scheme); Oleson et al. (2004) set h _e0 at a value <h ₀. Various studies have been conducted to find a suitable discretization to compromise between the accuracy and the computational cost as summarized in Table 1 (Viterbo and Beljaars 1995; Chen and Dudhia 2001; Roeckner et al. 2003; Oleson et al. 2004).

Table 1 Comparison of skin-layer parameterizations in the literature, where k = heat diffusivity, ω = diurnal angular velocity

Minimizing the error for calculating the ST of the dominant frequency component leads to an optimal effective thickness while calculating the ST. The optimal effective thickness of a single-layer discretization (h ₀ → ∞) for a soil column has been introduced by Arakawa and Mintz (1974) (denoted as A74) and Deardorff (1978) (denoted as D78). A74 has been adopted in many land surface parameterizations (Sellers et al. 1986; Tsuang and Yuan 1994; Tsuang and Tu 2002). Tsuang (2003), based on it, derived an analytical solution set for describing both land skin temperature and Planetary Boundary Layer air temperature. Tsuang (2005) reformulated the equation to determine surface ground heat flux from land skin temperature retrieved from satellites.

Since A74 and D78 use only a single layer for representing a soil column, they are very economical. Nonetheless, their accuracies have limitations. Recently, many state-of-the-art land surface parameterizations have started to simulate temperatures at multiple ground layers such as the force-restored scheme (also proposed by Deardorff 1978) used in SiB2 (Sellers et al. 1996), ECHAM (Roeckner et al. 2003), ECMWF (Blondin 1991; Viterbo and Beljaars 1995), MM5 (Grell et al. 1995), NOAH (Chen and Dudhia 2001) and CLM (Oleson et al. 2004; Dickinson et al. 2006). This study extends A74 and D78’s work from single-layer discretization to multiple-layer discretization to further increase accuracy.

Earlier versions of the proposed scheme have been applied for the determination of the phase changes for snow, ice and water by implementing it into a climate model (Tsuang et al. 2001) and a turbulent kinetic energy ocean model (Tu and Tsuang 2005). Chen and Dudhia (2001) and Oleson et al. (2004) also implemented similar concept into their models. Nonetheless, so far these numerical schemes have not been well documented and analyzed.

2 Optimal effective thickness

This section tries to determine the optimal value of the effective thickness (denoted as h ^p _ek ) of layer k under the conditions that the true ground heat fluxes (G _k , G _k+1) at ζ = 0 and ζ = −h _k are known. Nonetheless, it should be noted that due to numerical error, the calculated ground heat fluxes (G _k,n , G _k+1,n ) can be departed from their true fluxes. In this section, the vertical coordinate system ζ is used, where the origin (ζ = 0) of ζ is at G _k (Fig. 1). Then, the ground temperature T _k can be determined as:

$$ \frac{{\partial T_{k,n} \left( {\zeta = - h_{t} ,t} \right)}}{\partial t} = - \frac{{G_{k,n} \left( {\zeta = 0,t} \right) - G_{k + 1,n} \left( {\zeta = - h_{k} ,t} \right)}}{{\rho_{g} c_{g} h_{ek} }} $$

(2)

where h _ek is the effective thickness of layer k, T _k is the temperature at depth h _t below G _k . From Fig. 1 (or Eq. 1), it can be seen that T ₀ is determined under the conditions k = 0, h _t  = 0 and h ₀ = −z ₁/2, of which the effective thickness of the skin layer is h _e0; T _k is determined under the conditions h _t  = (z _k−1–z _k )/2 and h _k  = (z _k−1–z _k+1)/2, of which the effective thickness is h _ek ; and, T _m is determined under the conditions k = m, h _t  = (z _m–1–z _m )/2 and h _m  = ∞, of which the effective thickness of the bottom layer is h _em .

It should be noted that the solutions of the true ground temperature profile and the true ground heat flux profile of a ground column can be determined analytically under the conditions that for an ideal surface, in which the heat diffusion coefficient D is constant and neglecting the horizontal heat transport, the heat transfer in the ground can be assumed to obey the Fourier law of diffusion, and the upper boundary temperature can be expressed in the Fourier series (e.g., Tsuang 2003). Please refer to Appendix 1 for further details. Under the above conditions, the true profiles of ground temperature and heat flux within a ground column can be described analytically as shown in Eqs. 36 and 38, respectively, of which the upper boundary temperature T _sk is written as:

$$ T_{sk} = \overline{{T_{sk} }} + \sum\limits_{j = 1}^{j = \infty } {\Updelta T_{skj} \cos \left( {\omega_{j} \left( {t - t_{kmj} } \right)} \right)} $$

(3)

where the overbar “^–” is the average, the subscript “j” denotes a property of frequency component j, ΔT _skj is the amplitude of the T _sk of the particular frequency component, ω _j is the angular velocity of the frequency component, t _kmj is the time when the highest T _sk of the particular frequency component occurs, and t is local time. Moreover, the surface ground heat flux G ₀ can also be determined from the energy budget of the land surface, where outgoing terrestrial radiation, surface sensible heat flux and surface latent heat flux are functions of ST (e.g., Brutsaert 1982; Garratt 1992). That is, G ₀ is a function of T ₀. Similarly the subsurface ground heat fluxes G _k (k > 0) and G _k+1 are also a function of T _k , according to the Fourier law of diffusion. Hence, the error in calculated T _k causes an error in simulated G _k and G _k+1, which in turn changes the simulated T _k . It is assumed that the optimization can be done layer by layer, and only the dependencies of the heat fluxes G _k and G _k+1 with respect to T _k are considered. The temperatures of the level above and below are ignored. Then, the relation between calculated G _k,n and G _k+1,n can be related to their respective true fluxes G _k and G _k+1 (Eq. 38) at ζ = 0 and –h _k , respectively, by Taylor’s series expansion on T _k as:

$$ \left\{ {\begin{array}{*{20}c} \begin{aligned} G_{k,n} \left( {T_{k,n} } \right) & = G_{k} \left( {T_{k} } \right) + \frac{{\partial G_{k} }}{{\partial T_{k} }}\left( {T_{k,n} - T_{k} } \right) + {\text{H}} . {\text{O}} . {\text{T}} .\\ & = G\left( {\zeta = 0} \right) + x_{k} \left( {T_{k,n} - T\left( {\zeta = - h_{t} } \right)} \right) + {\text{H}} . {\text{O}} . {\text{T}}. \\ \end{aligned} \\ \begin{aligned} G_{k + 1,n} \left( {T_{k,n} } \right) & = G_{k + 1} \left( {T_{k} } \right) + \frac{{\partial G_{k + 1} }}{{\partial T_{k} }}\left( {T_{k,n} - T_{k} } \right) + {\text{H}} . {\text{O}} . {\text{T}} .\\ & = G\left( {\zeta = - h_{k} } \right) - x_{k + 1} \left( {T_{k,n} - T\left( {\zeta = - h_{t} } \right)} \right) + {\text{H}} . {\text{O}} . {\text{T}}. \\ \end{aligned} \\ \end{array} } \right. $$

(4)

where T _k,n is the simulated T _k (K). Note that G _k  = G(ζ = 0), G _k+1 = G(ζ = −h _k ) and T _k  = T(ζ = −h _t ). The second term shows the difference between calculated and true ground heat flux, and x _k is defined (Refer to Appendix 2 for the derivation) as:

$$ x_{0} = \frac{{\partial G_{0} \left( {T_{0} } \right)}}{{\partial T_{0} }} = 4\varepsilon \sigma T_{0}^{3} + \frac{{\rho_{a} c_{a} }}{{r_{a} }} + \frac{{\rho_{a} L_{v} }}{{r_{a} + r_{c} }}\frac{{\partial q^{*} \left( {T_{0} } \right)}}{{\partial T_{0} }} $$

(5a)

$$ x_{k} = - \frac{{\partial G_{k} }}{{\partial T_{k - 1} }} = \frac{{\partial G_{k} }}{{\partial T_{k} }} = \frac{{\rho_{g} c_{g} D}}{{z_{k - 1} - z_{k} }},\quad k = 1,\;m $$

(5b)

$$ x_{m + 1} = - \frac{{\partial G_{m + 1} }}{{\partial T_{m} }} = \frac{{\rho_{g} c_{g} D}}{{z_{m} - z_{m + 1} }} = \frac{{\rho_{g} c_{g} D}}{{z_{m} + \infty }} = 0 $$

(5c)

where ε is emissivity of surface; σ is Stefan-Boltzman constant (~5.67 × 10⁻⁸ W m⁻² K⁻⁴); ρ _a and c _a are density (~1.16 kg m⁻³) and constant pressure heat capacity (~1,005 J kg⁻¹ K⁻¹) of air, respectively; q ^*(T) is saturated specific humidity at temperature T (=0.622 e ^*(T)/P), where P is atmospheric pressure (Pa) and e ^* is saturated vapor pressure (Pa) (Richards 1971); r _a is aerodynamic resistance (s m⁻¹); r _c is canopy resistance (s m⁻¹) for evapotranspiration; and L _v is latent heat of evaporation (~2.5 × 10⁶ J/kg). Note that the third equality of Eq. 5c is derived due to that z _m+1 → −∞. Figure 2 shows ∂G ₀/∂T ₀ as a function of ST for various aerodynamic and canopy resistances with ε = 0.97 and P = 1,013 hPa. It can be seen that ∂G ₀/∂T ₀ increases with ST, but decreases with aerodynamic and canopy resistances. In general, during the daylight hours, aerodynamic and canopy resistances are much lower than those during the nights (Blondin 1991; Tsai et al. 2007). Typical values of daytime aerodynamic and canopy resistances observed in various land covers in the summer can be found in Wilson et al. (2002). Therefore, from the Fig. 2, we can infer that the values of ∂G ₀/∂T ₀ during summer days are much higher than those during winter nights.

Fig. 2

∂G ₀/∂T ₀ as a function of skin temperature (T ₀) for various aerodynamic resistances (r _a) and canopy resistances (r _c) with ε = 0.97 and P = 1,013 hPa

Substituting Eq. 4 into Eq. 2 and neglecting the high order terms (H.O.T.), the governing Eq. 2 can be rewritten as:

$$ \begin{aligned} \frac{{\partial T_{k,n} }}{\partial t} & = - \frac{{G_{k,n} - G_{k + 1,n} }}{{\rho_{g} c_{g} h_{ek} }} \\ & = \frac{ - 1}{{\rho_{g} c_{g} h_{ek} }}\left[ {G_{k} - G_{k + 1} + s_{k} \left( {T_{k,n} - T_{k} } \right) + {\text{H}} . {\text{O}} . {\text{T}}.} \right] \\ & = \frac{ - 1}{{\rho_{g} c_{g} h_{ek} }}\left[ {G\left( {\zeta = 0,t} \right) - G\left( {\zeta = - h_{k} ,t} \right) + s_{k} \left( {T_{k,n} - T\left( {\zeta = - h_{t} ,t} \right)} \right) + {\text{H}} . {\text{O}} . {\text{T}}.} \right] \\ \end{aligned} $$

(6)

where the elasticity s _k is defined as:

$$ s_{k} \equiv \left( {\frac{{\partial G_{k} }}{{\partial T_{k} }} - \frac{{\partial G_{k + 1} }}{{\partial T_{k} }}} \right) = x_{k} + x_{k + 1} $$

(7)

Hereafter, we denote the condition s _k  → ∞ as the “stiff” condition, and the condition s _k  → 0 as the “non-stiff” condition, borrowing the expressions from mass-spring system. Under the stiff condition, the restoring forcing (s _k ) is large, which forces T _k,n  → T _k according to Eq. 6. Under the non-stiff condition, the restoring forcing (s _k ) is weak and T _k,n can be very different from T _k .

The simulated T _k can also be written in the frequency domain as:

$$ T_{k,n} = \overline{{T_{k,n} }} + \sum\limits_{j = 1}^{j = \infty } {\Updelta T_{kj,n} \cos \left( {\omega_{j} \left( {t - t_{kmj} - t_{klj} } \right)} \right)} = \sum\limits_{j = 0}^{j = \infty } {\tau_{kj,n} } $$

(8)

where $\overline{{T_{k,n} }}$ and t _lj are the mean value and the time lag of the simulated T _k , respectively. From the second equality, it can be observed that for j = 0, $\tau_{k0,n} = \overline{{T_{k,n} }} ,$ but j ≥ 1, τ _kj,n  = ΔT _kj,n  cos (ω _j (t − t _kmj  − t _klj )). Substituting Eqs. 36, 38 and 8 into Eq. 6 and neglecting the H.O.T., the governing Eq. 6 can be rewritten in the frequency domain as:

$$ \begin{gathered} \sum\limits_{j = 0}^{j = \infty } {\left[ {\frac{{\partial \tau_{kj,n} }}{\partial t} + \frac{{s_{k} }}{{\rho_{g} c_{g} h_{ek} }}\tau_{kj,n} } \right]} \hfill \\ \quad = - \frac{1}{{\rho_{g} c_{g} h_{ek} }}\sum\limits_{j = 0}^{j = \infty } {\left[ {g_{j} \left( {\zeta = 0,t} \right) - g_{j} \left( {\zeta = - h_{k} ,t} \right) - s_{k} \tau_{j} \left( {\zeta = - h_{t} ,t} \right)} \right]} \hfill \\ \end{gathered} $$

(9)

where g _j (ζ, t) and τ _j (ζ, t) are true ground heat flux and true ground temperature of frequency component j, respectively. Please refer to Eqs. 39 and 37 for details. Due to the orthogonal property between different frequency components (Kreyszig 2006), the above governing equation for determining T _k by a numerical method can be decomposed for each frequency component into:

$$ \frac{{\partial \tau_{kj,n} }}{\partial t} + \frac{{s_{k} }}{{\rho_{g} c_{g} h_{ek} }}\tau_{kj,n} = - \frac{1}{{\rho_{g} c_{g} h_{ek} }}\left[ {g_{j} \left( {\zeta = 0,t} \right) - g_{j} \left( {\zeta = - h_{k} ,t} \right) - s_{k} \tau_{j} \left( {\zeta = - h_{t} ,t} \right)} \right] $$

(10)

Note that from the above equation, it can be solved that for j = 0, $\tau_{k0,n} = \overline{{T_{k,n} }} = \tau_{0} \left( {\zeta = - h_{t} ,t} \right) = \overline{{T_{k} }}$ since ${{\partial \tau_{k0,n} } \mathord{\left/ {\vphantom {{\partial \tau_{k0,n} } {\partial t}}} \right. \kern-\nulldelimiterspace} {\partial t}} = g_{0} \left( {\zeta = 0,t} \right) = g_{0} \left( {\zeta = - h_{k} ,t} \right) = 0.$ Note from Eq. 36, it can be derived that $\overline{{T_{k} }} = \overline{{T_{sk} }} .$ For j ≥ 1, the solution of τ _kj,n is more complicated and has been solved in Appendix 3.

The overall RMSE of the calculated T _k (denoted as e(T _k )), according to its definition, can be written as:

$$ e\left( {T_{k} } \right)^{2} \equiv \overline{{\left( {T_{k,n} - T_{k} } \right)^{2} }} = \sum\limits_{j} {\overline{{\left( {\tau_{kj,n} - \tau_{j} \left( {\zeta = - h_{t} ,\;t} \right)} \right)^{2} }} } = \sum\limits_{j} {e\left( {\tau_{kj} } \right)^{2} } $$

(11)

where the second equality is derived due to the orthogonal property between different frequency components (Kreyszig 2006). According to Appendix 3, the normalized root–mean–square error (NRMSE) (or RMSE/STD ratio) e(τ ^*) of the calculated T _k,n of the frequency component compared to the true T _k can be determined from Eq. 66 as:

$$ e\left( {\tau^{*} } \right) = \sqrt {\frac{{h_{a}^{*2} - 2h_{a}^{*} h_{e}^{*} \exp \left( { - h_{t}^{*} } \right)\cos \left( {t_{a}^{*} - h_{t}^{*} } \right) + \exp \left( { - h_{t}^{*} } \right)^{2} h_{e}^{*2} }}{{h_{e}^{*2} + 0.5s^{*2} }}} $$

(12)

where the asterisks denote the nondimensional variables of the frequency component j at the level k. The variables are made nondimensional, by multiplying time by the angular velocity ω _j of the frequency component, dividing the length by $\sqrt {2D/\omega_{j} } ,$ dividing energy flux by the standard deviation (STD) of the ground heat flux component g _kj at the upper boundary (ζ = 0), i.e., g _kjSTD, and dividing temperature component τ _kj by its STD at the upper boundary (ζ = 0), i.e., τ _kjSTD. That is:

$$ \tau^{*} \equiv \frac{{\tau_{kj} }}{{\tau_{{kj{\text{STD}}}} }} = \frac{{\sqrt 2 \tau_{kj} }}{{\Updelta T_{skj} }} $$

(13)

$$ h_{e}^{*} \equiv \frac{{h_{ek} }}{{\sqrt {2D/\omega_{j} } }} $$

(14)

$$ h_{a}^{*} \left( {h^{*} } \right) = \frac{1}{\sqrt 2 }\sqrt {1 - 2\cos \left( {h^{*} } \right)\exp \left( { - h^{*} } \right) + \exp \left( { - h^{*} } \right)^{2} } $$

(15)

$$ h_{t}^{*} \equiv \frac{{h_{t} }}{{\sqrt {2D/\omega_{j} } }} $$

(16)

$$ s^{*} \equiv \frac{{\tau_{{kj{\text{STD}}}} }}{{g_{{kj{\text{STD}}}} }}s_{k} = \frac{1}{{\rho_{g} c_{g} \sqrt {D\omega_{j} } }}s_{k} $$

(17)

$$ t_{a}^{*} \left( {h^{*} } \right) \equiv \omega_{j} t_{a} \left( {h^{*} } \right) = \frac{\pi }{4} - \tan^{ - 1} \left( {\frac{{\exp \left( { - h^{*} } \right)\sin \left( {h^{*} } \right)}}{{1 - \exp \left( { - h^{*} } \right)\cos \left( {h^{*} } \right)}}} \right) $$

(18)

where $h^{*} \equiv {{h_{k} } \mathord{\left/ {\vphantom {{h_{k} } {\sqrt {2D/\omega_{j} } }}} \right. \kern-\nulldelimiterspace} {\sqrt {2D/\omega_{j} } }}.$ Note that $\tau_{{kj{\text{STD}}}} = {{\Updelta T_{skj} } \mathord{\left/ {\vphantom {{\Updelta T_{skj} } {\sqrt 2 }}} \right. \kern-\nulldelimiterspace} {\sqrt 2 }}$ and $g_{{kj{\text{STD}}}} = \rho_{g} c_{g} \sqrt {D\omega_{j} } {{\Updelta T_{skj} } \mathord{\left/ {\vphantom {{\Updelta T_{skj} } {\sqrt 2 }}} \right. \kern-\nulldelimiterspace} {\sqrt 2 }}$ according to Eqs. 41 and 42, respectively. It shows that e(τ ^*) is a function of h ^*, h ^* _e , h ^* _t and s ^* only. Note that h ^* _a and t ^* _a are functions of h ^* according to Eqs. 15 and 18, respectively. This study tries to determine a better value for h _ek , of which the root–mean–square error (RMSE) of the calculated ground temperature T _k,n of the dominant frequency component j is minimal. That is,

$$ \mathop {\min }\limits_{{h_{ek} }} e\left( {\tau_{kj} } \right) $$

(19)

Solving the above equation, the optimal value of the effective thickness (denoted as h ^p _ek ) can be determined. Please refer to Appendix 3 for the derivation. And its dimensionless form is written (denoted as h ^p* _e ) as:

$$ \begin{aligned} h_{e}^{p*} & \equiv \frac{{h_{ek}^{p} }}{{\sqrt {2D/\omega_{j} } }} \\ & = \frac{{2h_{a}^{*2} - \exp \left( { - h_{t}^{*} } \right)^{2} s^{*2} + \sqrt {4h_{a}^{*4} + 4\cos 2\left( {t_{a}^{*} - h_{t}^{*} } \right)\exp \left( { - h_{t}^{*} } \right)^{2} h_{a}^{*2} s^{*2} + \exp \left( { - h_{t}^{*} } \right)^{4} s^{*4} } }}{{4\cos \left( {t_{a}^{*} - h_{t}^{*} } \right)\exp \left( { - h_{t}^{*} } \right)h_{a}^{*} }} \\ \end{aligned} $$

(20)

We name the numerical scheme, which sets h ^* _e  = h ^p* _e , as the optimal scheme, proposed by this study (denoted as “op”) to calculate ground temperature at each numerical layer. For the skin layer, the optimal effective thickness h _e0 can be determined by putting level index k = 0 and h _t = 0 into Eq. 20 (Fig. 1b) as:

$$ h_{e0}^{p*} = \frac{{2h_{a0}^{*2} - s_{0}^{*2} + \sqrt {4h_{a0}^{*4} + 4\cos \left( {2t_{a0}^{*} } \right)h_{a0}^{*2} s_{0}^{*2} + s_{0}^{*4} } }}{{4\cos \left( {t_{a0}^{*} } \right)h_{a0}^{*} }} $$

(21)

where the subscript “0” denotes a property of the skin layer. For the bottom-numerical layer, the effective thickness h _em can be determined by putting level index k = m and h _m  = ∞ into Eq. 20 (Fig. 1c) as:

$$ h_{em}^{p*} = \frac{{1 - \exp \left( { - h_{tm}^{*} } \right)^{2} s_{m}^{*2} + \sqrt {1 + 2\sin \left( {2h_{tm}^{*} } \right)\exp \left( { - h_{tm}^{*} } \right)^{2} s_{m}^{*2} + \exp \left( { - h_{tm}^{*} } \right)^{4} s_{m}^{*4} } }}{{2\sqrt 2 \cos \left( {{\pi \mathord{\left/ {\vphantom {\pi 4}} \right. \kern-\nulldelimiterspace} 4} - h_{tm}^{*} } \right)\exp \left( { - h_{tm}^{*} } \right)}} $$

(22)

where the subscript “m” denotes a property of the bottom numerical layer, $h_{am}^{*} \left( \infty \right) \to {1 \mathord{\left/ {\vphantom {1 {\sqrt 2 }}} \right. \kern-\nulldelimiterspace} {\sqrt 2 }},$ and t ^* _am (∞) = π/4.

Moreover, under the condition that ∂G ₀/∂T ₀ is constant, the RMSE of the calculated surface ground heat flux (denoted as e(G ₀)) can be derived from Eq. 5a by neglecting the H.O.T. as:

$$ e\left( {G_{0} } \right)^{2} \equiv \overline{{\left( {G_{0,n} - G_{0} } \right)^{2} }} = \sum\limits_{j} {\overline{{\left( {g_{0j,n} - g_{0j} } \right)^{2} }} } = \sum\limits_{j} {e\left( {g_{0j} } \right)^{2} } \approx \left( {\frac{{\partial G_{0} }}{{\partial T_{0} }}} \right)^{2} \sum\limits_{j} {e\left( {\tau_{0j} } \right)^{2} } $$

(23)

It implies that while the error of the simulated ST of the dominant frequency component is minimal, the error of the simulated G ₀ of the frequency component is also minimal. Neglecting the H.O.T. and under the conditions that the true ground heat fluxes (G ₀, G ₁) at z = 0 and z = −h ₀ are known, the NRMSE of the calculated surface ground heat flux e(g ^*₀ ) of a frequency component compared to the true flux can be derived from Eq. 23, by substituting Eqs. 12 and 21 into it, as:

$$ \begin{aligned} e\left( {g_{0}^{*} } \right) \approx & \sqrt {\overline{{\left[ {\frac{{\partial G_{0}^{*} }}{{\partial T_{0}^{*} }}\left( {\tau_{0,n}^{*} - \tau_{0}^{*} } \right)} \right]^{2} }} } = \frac{{\partial G_{0}^{*} }}{{\partial T_{0}^{*} }}\sqrt {\overline{{\left[ {\left( {\tau_{0,n}^{*} - \tau_{0}^{*} } \right)} \right]^{2} }} } = \frac{{\partial G_{0}^{*} }}{{\partial T_{0}^{*} }}e\left( {\tau_{0}^{*} } \right) \\ = & \frac{{\partial G_{0}^{*} }}{{\partial T_{0}^{*} }}\sqrt {\frac{{h_{a0}^{*2} - 2h_{a0}^{*} h_{e0}^{*} \cos \left( {t_{a0}^{*} } \right) + h_{e0}^{*2} }}{{h_{e0}^{*2} + 0.5s_{0}^{*2} }}} \\ \end{aligned} $$

(24)

where h ^* _t  = 0 for the skin layer. It shows that e(g ^*₀ ) is a function of h ^*₀ , h ^* _e0 , k ^*₀ and ∂G ^*₀ /∂T ^*₀ only. Note that h ^* _a0 is a function of h ^*₀ .

Figure 3 shows the accuracies e(g ^*₀ ) of the “A74”, “D78”, “cv”, “nh” and “op” schemes as functions of h ^*₀ under ∂G ^*₀ /∂T ^*₀  = 4 and 0.5. The e(g ^*₀ ) of the “cv”, “nh” and “op” schemes are calculated from Eq. 24 by setting h ^* _e0 at h ^*₀ , 0 and h ^p* _e0 , respectively. The e(g ^*₀ ) of the “A74” and “D78” schemes are calculated from Eq. 24 by setting h ^* _e0 at ${1 \mathord{\left/ {\vphantom {1 {\sqrt 2 }}} \right. \kern-\nulldelimiterspace} {\sqrt 2 }}$ and 1, respectively, under the condition h ^*₀  → ∞. This figure shows that for the same h ^*₀ , the NRMSE simulated by “op” is always the lowest among all the schemes. This is expected since we derive “op” by minimizing the NRMSE. In addition, it shows that the NRMSEs of “cv”, “nh” and “op” decrease with thinner h ^*₀ . Therefore, for achieving a higher accuracy of “op”, a finer thickness of the skin layer is needed. In contrast, the accuracy of the single layer schemes “A74” and “D78” are about as good as “op”, and better than “cv” and “nh” when h ^*₀ is thick (>2). Surprisedly, from the figure, when h ^*₀  < 1, the accuracy of “nh” is worse than “cv”. It implies that using a thinner top layer or inclusion of an infinite thin skin layer like Viterbo and Beljaars (1995) showed a worse behavior than the standard “cv” method. This implies that when h ^*₀  < 1, it is rather calculating the skin temperature by the standard “cv” method than the “nh” scheme.

Fig. 3

Normalized root–mean–square errors (or RMSE/STD ratios) e(g ^*₀ ) of the surface ground heat flux of a particular frequency, as functions of h ^*₀ , under the conditions that the true ground heat fluxes (G ₀, G ₁) at surface (z = 0) and at z = −h ₀ are known, for ∂G ^*₀ /∂T ^*₀  = 4 and 0.5 as determined by the optimal scheme (denoted as “op”), the optimal stiff scheme (denoted as “os”), the non-stiff equal-amplitude scheme (denoted as “ne”), the optimal non-stiff scheme (denoted as “on”), the conventional finite-difference scheme (“cv”), the no-heat capacity scheme (“nh”), the single-layer scheme of Arakawa and Mintz (1974) (“A74”) and the single-layer scheme of Deardorff (1978) (“D78”)

3 Vertical discretization for single frequency component

In the above section, the optimal effective thickness h ^p* _ek has been determined for each numerical layer k under the conditions that z ^* _k−1 , z ^* _k , z ^* _k+1 are known. Note that h ^* _k  = 0.5(z ^* _k−1  − z ^* _k+1 ) and h ^* _tk  = 0.5(z ^* _k−1  − z ^* _k ). In addition, the accuracy e(g ^*₀ ) (Eq. 24) is derived under the conditions that the true ground heat fluxes at z = 0 and z = −h ₀ are known. In fact, we usually do not know the true fluxes. Therefore, the overall RMSE (denoted as e(G _0j )) of the calculated surface ground heat flux of a particular frequency component j can be larger than e(g ^*₀ ). For determining vertical discretization z _k , after various tests, it is suggested to use the evenly heat-content discretized grid (discretization the vertical profile to have the same heat content within z _k and z _k+1) to have the least overall error among our tests for surface ground heat flux simulation.

From Eq. 36, it can be seen that the heat content of a particular frequency component decays exponentially with depth. Of the frequency component, the ratio p of the heat content stored from the surface to the depth d to the total heat content within the ground column can be determined as:

$$ p\left( {d^{*} } \right) = \frac{{\int_{{ - d^{*} }}^{0} {\left\| {\tau^{*} \left( {z^{*} } \right)} \right\|{\text{d}}z^{*} } }}{{\int_{ - \infty }^{0} {\left\| {\tau^{*} \left( {z^{*} } \right)} \right\|{\text{d}}z^{*} } }} = \frac{{\int_{{ - d^{*} }}^{0} {\exp \left( {z^{*} } \right){\text{d}}z^{*} } }}{{\int_{ - \infty }^{0} {\exp \left( {z^{*} } \right){\text{d}}z^{*} } }} = \frac{{\left. {\exp \left( {z^{*} } \right)} \right|_{{ - d^{*} }}^{0} }}{{\left. {\exp \left( {z^{*} } \right)} \right|_{ - \infty }^{0} }} = 1 - \exp \left( { - d^{*} } \right) $$

(25)

Reorganizing the above equation, the function d(p) can be determined as:

$$ d^{ *} \left( p \right) = - \ln \left( {1 - p} \right) $$

(26)

Based on the above criterion, the vertical discretization z _k of a ground column can be determined (denoted as evenly heat-content discretization).

Table 2 lists the values of z ^* _k , h ^* _k , h ^p* _ek and normalized overall error e(G ^*_0j ) (=e(G _0j )/g_0jSTD) of the evenly heat-content discretization from single layer (m = 0) up to 6 layers (m = 5) of a particular frequency component. It is uneven grid with finer discretization near the surface, and infinitely coarse grid for the bottommost layer. The subsurface fluxes G _k,n (for k ≥ 1) can be determined as a function of ground temperatures according to the finite difference Eq. 49. Then, the system of (m + 1) ordinary differential equations for determining the skin temperature (T ₀) as well as ground temperatures (T ₁,…,T _m) can be rewritten from Eq. 1 as:

Table 2 Vertical discretizations of evenly heat-content discretized grids and their normalized overall errors e(G ^*_0j ) (= e(G _0j )/g _0jSTD) of surface ground heat flux of a particular frequency component j when ∂G ^*₀ /∂T ^*₀  = 2.64

$$ \left\{ {\begin{array}{*{20}c} {\frac{{\partial T_{0,n} }}{\partial t} = - \frac{{G_{0,n} }}{{\rho_{g} c_{g} h_{e0} }} - \frac{D}{{h_{e0} }}\left( {\frac{{T_{0,n} - T_{1,n} }}{{z_{0} - z_{1} }}} \right)} \\ \vdots \\ {\frac{{\partial T_{k,n} }}{\partial t} = \frac{D}{{h_{ek} }}\left( {\frac{{T_{k - 1,n} - T_{k,n} }}{{z_{k - 1} - z_{k} }} - \frac{{T_{k,n} - T_{k + 1,n} }}{{z_{k} - z_{k + 1} }}} \right)} \\ \vdots \\ {\frac{{\partial T_{m,n} }}{\partial t} = - \frac{{G_{m,n} }}{{\rho_{g} c_{g} h_{em} }} = \frac{D}{{h_{em} }}\left( {\frac{{T_{m - 1,n} - T_{m,n} }}{{z_{m - 1} - z_{m} }}} \right)} \\ \end{array} } \right. $$

(27)

where the “op” scheme, proposed by this study, sets h _ek  = h ^p _ek according to Eq. 20, and the “cv” scheme sets h _ek  = h _k . The normalized overall error e(G ^*_0j ) in the table is determined by solving the above ordinary-differential-equation system by the Runge–Kutta method using the software Mathematica (http://www.wolfram.com) by setting ∂G ^*₀ /∂T ^*₀ at 2.64. It can be seen that e(G ^*_0j ) reduces from 68% of the single-layer discretization (m = 0) to 2% of the 4-layer discretization (m = 3). Thus, the more the layers, the lower the overall error is. The effective thickness in each layer is thinner than its respective physical thickness, especially for the bottommost-numerical layer. Comparing to the normalized error e(g ^*₀ ) of the skin layer, the normalized overall error e(G ^*_0j ) for m in the range of 1–5 is a factor of 3 larger than that of e(g ^*₀ ) for the same h ₀. This is reasonable since the error in subsurface flux (z = −h ₀) is not included in e(g ^*₀ ). A tri-diagonal matrix for solving the temperatures by the “op” scheme is illustrated in the Appendix (See Appendix 4 for details), which can be solved directly without any iteration by numerical solvers (e.g., Press et al. 1992).

Figure 4 shows the optimal effective thicknesses h ^p* _e0 of various ∂G ^*₀ /∂T ^*₀ (0, 1, 2, 4, ∞), h ^*₀ and h ^* _a0 as functions of h ^*₀ . It shows that when h ^*₀  > 0.6, the value of h ^p* _e0 is scattered and decreased with ∂G ^*₀ /∂T ^*₀ . But when h ^*₀  < 0.6, the value of h ^p* _e0 is almost independent to ∂G ^*₀ /∂T ^*₀ . For example, when h ^*₀  = d ^*(25%) = 0.144, h ^p* _e0 is within [0.133530, 0.133533], changing very little with ∂G ^*₀ /∂T ^*₀ (bottom panel of the Fig. 4). From Table 2, it can be seen h ^*₀  < 0.35 for the evenly heat-content discretization with layers ≥2 (or m ≥ 1). Therefore, it is suggested to determine h ^p* _e0 using some typical values of ∂G ^*₀ /∂T ^*₀ , such as the annual median value of 2.64 of the Bondville site (shown in the later section) when the value of ∂G ^*₀ /∂T ^*₀ is not immediately available.

Fig. 4

Top panel optimal effective thicknesses h ^p* _e0 of various ∂G ^*₀ /∂T ^*₀ and h ^* _a0 as functions of the skin-layer thickness h ^*₀ , where the numbers in the legend denote the values of ∂G ^*₀ /∂T ^*₀ ; bottom panel h ^p* _e0 as a function of ∂G ^*₀ /∂T ^*₀ for h ^*₀  = d ^*(25%) = 0.14381. All the variables are expressed in dimensionless units. Please see Eq. 52 for the definitions

4 Vertical discretization for multiple frequency components

Table 2 shows the vertical discretizations of evenly heat-content discretized grids for a particular frequency component. Nonetheless, there are uncertainties with the frequency chosen. In fact, diurnal, annual and many other frequency components are present in ST. In the followings, we use the subscripts “d”, “y” and “s” denoting the properties of “diurnal”, “annual” and “sun-spot” frequency components, respectively. Besides, it is well known that in the deep ground, the dominant frequency component usually is much slower than in the surface (Huang et al. 2000). Hence, the dominant frequency component of the bottommost-numerical layer can be slower than the diurnal frequency. Therefore, the locations z _k of a ground column should be chosen to be able to record the temperature evolutions of multiple dominant frequency components. This section uses the skin temperature data at a cropland site, Bondville, IL, USA (40.01°N, 88.29°W) to design a better discretization for multiple frequency components.

At the Bondville site, the soil is silt loam. The heat diffusivity D for the soil with water content within wilting point and field capacity is 6.2 × 10⁻⁷ m² s⁻¹ and its volume heat capacity ρ _g c _g is about 2.4 × 10⁶ J m⁻³ K⁻¹ (de Vries 1975). A more detailed description of the site can be found in Meyers and Hollinger (2004) and Tsuang (2005). Figure 5 shows the time series of the observed ST at the cropland site from 1997 to 2000, and its frequency spectrum. It can be seen that its major frequency components include 1/y, 1/d, 1/12h, 1/8h and 1/6h. It can be fitted by cosine functions. The result after truncating terms with ΔT ₀ ≤ 1 K or Δg ₀ ≤ 1 W m⁻² is as:

Fig. 5

Top skin temperature (T ₀) observed at a cropland site, Bondville, IL, USA (40.01°N, 88.29°W) from 1997–2000, middle the amplitudes ΔT of skin temperature, and bottom the amplitudes Δg of surface ground heat flux of various frequency components

$$ \begin{gathered} T_{0} \left( t \right) = 285.15 - 1.14\cos \left( {\frac{2\pi }{{4{\text{y}}}}t} \right) + 11.88\cos \left( {\frac{2\pi }{{1{\text{y}}}}\left( {t - 199d} \right)} \right) + 3.44\cos \left( {\frac{2\pi }{1d}\left( {t - 14h} \right)} \right) \\ + 0.94\cos \left( {\frac{4\pi }{1d}\left( {t - 1h} \right)} \right) + 0.25\cos \left( {\frac{6\pi }{1d}\left( {t - 5h} \right)} \right) + 0.10\cos \left( {\frac{8\pi }{1d}\left( {t - 3h} \right)} \right) \\ \end{gathered} $$

(28)

where t is time starting from 0:00 LT, 1 January 1997 and the unit of the ST is in K. Three dominant frequency components are identified. Their frequencies are 1/d, 1/y and 1/4y. The analytical solutions of their corresponding temperature profiles and ground heat fluxes of the ST (Eq. 28) observed at the site can be determined from Eqs. 36 and 38, directly. A comparison between observed monthly ground heat flux with the analytical solution can be found in Tsuang (2005). It shows that the correlation coefficient is 0.8 with RMSE at 5.6 W m⁻². Table 3 shows the STDs of ST and G ₀ of each frequency component. It can be seen that the STD of G ₀ of the diurnal frequency component is the largest among all the frequency components although it’s STD of ST is less than the annual frequency component. Therefore, calculating the optimal h _e0 based on the diurnal frequency component is a proper choice.

Table 3 Characteristics of the dominant frequency components of the skin temperature observed at the Bondville site (except for j = 7) during 1997–2000

Solar radiation is the dominant surface energy component for heating Earth’s surface (Tsuang 2003) and solar radiation has strong diurnal and annual cycles. In addition, the sun-spot cycle of 11 years might be of interest for climate study. The variation caused by the sunspot cycle to solar output is on the order of 0.1% of the solar constant (a peak-to-trough range of 1.3 W m⁻² compared to 1,366 W m⁻² for the average solar constant) (Ramaswamy et al. 2001). This range is slightly smaller than the change in radiative forcing caused by the increase in atmospheric CO₂ since the eighteenth century.

At the Bondville cropland site, from Table 3, the STD of G ₀ of the diurnal component was 39 W m⁻², and that of the annual frequency component was 7 W m⁻². To have an accuracy at 1 W m⁻², the level index “m” should be equal to or larger than 3 for the diurnal frequency component, and 2 for the annual frequency component according to Table 2. In respect to the sun-spot frequency component there is no need to reserve layers for the component since its STD of G ₀ was 0.46 W m⁻², which is only half of 1 W m⁻². If we choose m = 3 for the diurnal component, m = 2 for the annual component and m = 0 for the sun-spot component (denoted as the “op(3,2,0)” scheme), then, the overall error can be roughly estimated as:

$$ \begin{aligned} e\left( {G_{0} } \right) = & \sqrt {\sum\limits_{j} {\left( {e\left( {G_{0j}^{*} } \right)g_{{0j{\text{STD}}}} } \right)^{2} } } \\ \approx & \sqrt {\left( {e\left( {G_{0}^{p*} \left( {m = 3} \right)} \right)g_{{0d{\text{STD}}}} } \right)^{2} + \left( {e\left( {G_{0}^{p*} \left( {m = 2} \right)} \right)g_{{0y{\text{STD}}}} } \right)^{2} + \left( {e\left( {G_{0}^{n*} \left( {m = 0} \right)} \right)g_{{0s{\text{STD}}}} } \right)^{2} } \\ = & \sqrt {\left( {0.020 \times 39.2} \right)^{2} + \left( {0.048 \times 7.08} \right)^{2} + \left( {1 \times 0.46} \right)^{2} } \\ = & 1.09\;{\text{Wm}}^{ - 2} \\ \end{aligned} $$

(29)

where e(G ^p*₀ ) and e(G ^n*₀ ) denote the normalized overall errors of surface ground heat flux simulation of the “op” scheme and the no-heat-capacity scheme (“nh”), respectively. Note that the value of e(G ^n*_0j ) is 1 for the single-layer (m = 0, h ^*₀  → ∞) “nh” scheme, of which e(G ^*_0j ) can be determined by setting h ^* _e0  = 0 and h ^*₀  → ∞ into Eq. 24. Note that for the single layer scheme e(G ^*_0j ) = e(g ^*₀ ) since the subsurface ground heat flux at z = −h ₀ (→ −∞) is known (=0). For the “op(3,2,0)” scheme, since there is no middle layer recording the temperature evolution of the sun-spot frequency component, h ^*_0s of the sun-spot frequency component approaches infinity. Therefore, the effective thickness h ^* _e0s of the sun-spot cycle is close to 0 comparing to h ^*_0s . ($h_{e0s}^{*} = {{h_{e0d}^{p*} \sqrt {{{2D} \mathord{\left/ {\vphantom {{2D} {\omega_{d} }}} \right. \kern-\nulldelimiterspace} {\omega_{d} }}} } \mathord{\left/ {\vphantom {{h_{e0d}^{p*} \sqrt {{{2D} \mathord{\left/ {\vphantom {{2D} {\omega_{d} }}} \right. \kern-\nulldelimiterspace} {\omega_{d} }}} } {\sqrt {{{2D} \mathord{\left/ {\vphantom {{2D} {\omega_{s} }}} \right. \kern-\nulldelimiterspace} {\omega_{s} }}} }}} \right. \kern-\nulldelimiterspace} {\sqrt {{{2D} \mathord{\left/ {\vphantom {{2D} {\omega_{s} }}} \right. \kern-\nulldelimiterspace} {\omega_{s} }}} }} = 0.000 2 7 5).$ Therefore, the accuracy of the sun-spot frequency component is close to that of “nh” for h ^*₀  → ∞. Hence, e(G ₀)of “op(3,2,0)” can be approximated as the above Eq. 29.

Table 4 lists the grid structures and effective thicknesses of “op(3,2,0)” and others (such as “A74”, “D78”, “ECHAM”). The vertical coordinates (z _k ) of the 5 layers of “op(3,2,0)” are at 25, 50 and 75% of the heat storage of the diurnal component, and 33 and 67% of the annual component. That is, they are at $- d^{ *} \left( { 0. 2 5} \right)\sqrt {{{ 2 {\text{D}}} \mathord{\left/ {\vphantom {{ 2 {\text{D}}} {\omega_{\text{d}} }}} \right. \kern-\nulldelimiterspace} {\omega_{\text{d}} }}} ,\; - d^{ *} \left( { 0. 5} \right)\sqrt {{{ 2 {\text{D}}} \mathord{\left/ {\vphantom {{ 2 {\text{D}}} {\omega_{\text{d}} }}} \right. \kern-\nulldelimiterspace} {\omega_{\text{d}} }}} ,\; - d^{ *} \left( { 0. 7 5} \right)\sqrt {{{ 2 {\text{D}}} \mathord{\left/ {\vphantom {{ 2 {\text{D}}} {\omega_{\text{d}} }}} \right. \kern-\nulldelimiterspace} {\omega_{\text{d}} }}} ,\; - d^{ *} \left( { 0. 3 3} \right) \sqrt {{{ 2 {\text{D}}} \mathord{\left/ {\vphantom {{ 2 {\text{D}}} {\omega_{\text{y}} }}} \right. \kern-\nulldelimiterspace} {\omega_{\text{y}} }}} \;{\text{and}}\; - d^{ *} \left( { 0. 6 7} \right)\sqrt {{{ 2 {\text{D}}} \mathord{\left/ {\vphantom {{ 2 {\text{D}}} {\omega_{\text{y}} }}} \right. \kern-\nulldelimiterspace} {\omega_{\text{y}} }}} ,$ respectively; or −0.038, −0.090, −0.181, −1.012, and −2.742 m, respectively, for the heat diffusivity D of 6.2 × 10⁻⁷ m² s⁻¹ (silt loam). Then, h _e0 is determined to be 0.017 m by minimizing the error of the diurnal frequency component when ∂G ₀/∂T ₀ = 42 W m⁻² K⁻¹ (the annual median value of a cropland site in Bondville, USA, shown in the later section), and h _em to be 2.579 m by minimizing the error of the annual frequency component.

Table 4 Vertical discretization of various schemes for heat diffusivity at 6.2 × 10⁻⁷ m² s⁻¹, when ∂G ₀/∂T ₀ is fixed at 42 W m⁻² K⁻¹ (the annual median value of a cropland site in Bondville, USA)

5 Case study—a cropland site

The above section estimates that “op(3,2,0)” can simulate surface ground heat flux with an accuracy at about 1 W m⁻². This section applies the above discretization for determining skin temperature at a cropland site, Bondville, IL, USA (40.01°N, 88.29°W), for a case study. Seven cases of cascade frequency components based on the frequency component data observed at this study site are used to construct the cases (Table 5).

Table 5 Analytical skin temperature used for case studies

Figure 6 shows the time series of observed ∂G ₀/∂T ₀ at a cropland site, Bondville, IL, USA (40.01°N, 88.29°W) in 2000. In addition, ∂G ₀ ^*/∂T ^*₀ of the diurnal frequency component is also shown in the figure. At the site, the daytime aerodynamic and canopy resistances in the summer are about 37 and 60 s m⁻¹, respectively (Wilson et al. 2002). Following the method of Wilson et al. (2002), r _a is derived from the sum of the resistance to momentum transport and an excess resistance for scalar fluxes (Verma et al. 1986) and r _c is determined according to the Penman–Monteith approximation to the big leaf equations (Jarvis and McNaughton 1986; Shuttleworth et al. 1984). It can be seen that at this site, typical value for winter nighttime ∂G ^*₀ /∂T ^*₀ (∂G ₀/∂T ₀) was about 1.8 (29.04 W m⁻² K⁻¹), and that of summer day was about 4.4 (70.93 W m⁻² K⁻¹), of which the annual mean was 2.86 (46 W m⁻² K⁻¹) with the median value at 2.64 (42 W m⁻² K⁻¹), when T ₀ during wintertime (November–January) was usually <273 K.

Fig. 6

Time series of observed ∂G ₀/∂T ₀ and hourly composite values at a cropland site, Bondville, IL, USA (40.01°N, 88.29°W) in 2000, where ∂G ^*₀ /∂T ^*₀ is making nondimensional by the diurnal frequency component. The midday period is from 10 to 14 LT and the midnight period is from 22 to 02 LT. Summer period is defined from June to August and winter period is defined from November to January

Table 6 shows the bias, RMSE and normalized RMSE of skin temperature, and the bias, RMSE and normalized RMSE for surface ground heat of Case 7 simulated by “A74”, “D78”, “op(0,0,0)”, “op(1,0,0)”, “op(1,1,0)”, “op(2,1,0)”, “echam”, “ecmwf”, “ecmwf + op”, “cv(3,2,0)”, “op(3,1,0)”, “echam + op”, “cv(3,2,0) + nh”, “cv(3,2,0) + op”, “op(3,2,0)” and “op(3,2,1)” by solving Eq. 27 by the Runge–Kutta method. Figure 7 illustrates the time series and the spectrums of calculated skin temperature T _0,n and surface ground heat flux G _0,n among “A74”, “echam”, “cv(3,2,0)” and “op(3,2,0)” for Case 7. Here, the parentheses (d,y,s) denotes the number “m” (Table 2) of evenly heat-content discretized layers used for recording the diurnal temperature profile (d), the annual profile (y) and the sun-spot frequency profile (11 y) (s), respectively. The same initial condition (T(z _k ,0)) and the same boundary conditions (G(0,t), G(−∞,t)) are used for all the calculations. The conditions are prescribed as those of their corresponding analytical forms (Eqs. 36 and 37). The numerical results are compared with the analytical T ₀ of Eq. 28 and the analytical G ₀ (Eq. 39). The “cv + nh” scheme denotes using the “cv” scheme for the middle layers, but using the “nh” scheme for the skin layer; similarly, the “cv + op” scheme denotes using the “cv” scheme for the middle layers, but using the “op” scheme for the skin layer.

Table 6 Bias (bias(T)), RMSE (e(T ₀)) and normalized RMSE (e(T ^*₀ )) of skin temperature simulation, and bias (bias(G)), RMSE (e(G ₀)) and normalized RMSE (e(G ^*₀ )) for surface ground heat flux simulation of various schemes (as listed in Tables 2 and 4) for Case 7 (as listed in Table 6) during 6-day integration starting from t = 0. During the period, ∂G ₀/∂T ₀ is fixed at 42 W m⁻² K⁻¹

Fig. 7

Time series of calculated skin temperature T _0,n and surface ground heat flux G _0,n among “A74”, “echam”, “cv(3,2,0)” and “op(3,2,0)” for Case 7, with frequency components observed at the Bondville cropland site. The vertical discretizations of the schemes are listed in Table 4. Same initial (T(z,0)) and boundary (G(0,t)) conditions with zero-heat flux at the bottom boundary layer are used for all the calculations, where T(z,0) and G(0,t) are prescribed as those of their corresponding analytical forms (Eqs. 36 and 39). The results are compared with the analytical T ₀ of Eq. 28 and its corresponding analytical G ₀ (Eq. 39)

For Case 7, it can be seen that the normalized overall error of the simulated G ₀ by “op(3,2,0)” is 0.89 W m⁻², which is close to the value 1.09 W m⁻² roughly estimated from Eq. 29. This shows the error can be approximated as Eq. 29. For Case 7, it can be seen that the normalized overall errors of the simulated G ₀ by “op(3,2,0)” and “op(3,2,1)” are the lowest group. They are at 2%, or at 0.9 W m⁻². Those by “A74”, “D78” and “op(0,0,0)” are the highest group, varying from 70 to 90%, or within 30–39 W m⁻²; it is expected since there is only one layer (m = 0) in these schemes. In addition, it can be seen that for the same vertical structure (3,2,0), “op” is better than “cv + op”, “cv + op” is better than “cv + nh”, and “cv + nh” is better than “cv”. For Case 7, e(G ^*₀ ) of “op(3,2,0)” is at 2%, that of “cv(3,2,0) + op” 7%, that of “cv(3,2,0) + nh” 12%, and that of “cv(3,2,0)” 28%. This shows that using the optimal effective thickness (determined by Eq. 20) for each numerical layer does improve the model system for surface ground heat flux simulation, by more than an order comparing to the conventional finite difference scheme, and by about 6 times comparing to the no-heat capacity scheme used in some state-of-the-art models.

In addition, the difference between simulated results and the analytical solutions are expressed in the frequency spectrum in Fig. 8. It can be seen that the errors simulated by “op(3,2,0)” have been reduced for the component with frequency of the diurnal cycle or higher comparing to “A74”, “echam” and “cv(3,2,0)”.

Fig. 8

Same as Fig. 7, but these differences between numerical solution and analytical solution are expressed in the frequency spectrum

6 Discussion

6.1 Multiple frequency components

In respect to frequency components, Table 7 shows the normalized overall errors e(G ^*₀ ) of the 7 cases simulated by “op(0,0,0)”, “op(1,0,0)”, “op(2,0,0)”, “op(3,0,0)”, “op(4,0,0)”, “op(5,0,0)”, “op(1,1,0)”, “op(2,1,0)”, “op(3,1,0)”, “op(3,2,0)” and “op(3,2,1)”. Cases 1–4 have diurnal and higher frequency components, but no annual frequency component. Cases 5–7 contain diurnal as well as annual frequency components. In the table, for Cases 1–4, the sequence of the values of e(G ^*₀ ) from low to high is “op(5,0,0)” < “op(4,0,0)” < “op(3,0,0)” < “op(2,0,0)” < “op(1,0,0)” < “op(0,0,0)”. It shows that the more layers of z _k for keeping tracks the diurnal temperature profile are, its accuracy for calculating G ₀ of the frequency component becomes higher. Nonetheless, although “op(5,0,0)” can calculate G ₀ of the diurnal frequency component accurately, it is not able to handle G ₀ of the annual (or lower) frequency component. The normalized overall errors e(G ^*₀ ) of Cases 1–4 by “op(5,0,0)” are at about 0.8%, but the errors increase to about 16% for Cases 5–7. In contrast, those for Cases 5–7 by “op(3,2,0)” are only 2%. Note that “op(3,2,0)” consists two layers for recording the annual profile of ground temperature. Therefore, it is necessary to allocate a few layers of z _k to simulate the annual frequency component properly. Nonetheless, in respect to the sun-spot frequency component, which is in Case 7, it shows that the normalized overall error e(G ^*₀ ) of op(3,2,1) for Case 7 is only slightly better than that of “op(3,2,0)” by 0.01% (or 0.004 W m⁻²). Therefore, allocation of an additional layer for recording sun-spot frequency is not crucial for having accuracy at 1 W m⁻² for G ₀ simulation.

Table 7 Same as Table 6 but for the overall normalized RMSE error (e(G ^*₀ )) (%) for surface ground heat flux simulation of various schemes for 7 cases (as listed in Table 5)

6.2 Lag-predict problem of snow melting time

The lag-predict of snow melting time has been found in many state-of-the-art models (e.g., Arpe et al. 1994; Link and Marks 1999). One likely reason is caused by the “cv” scheme used for their snow ST simulation. Figure 9 shows the STs simulated by the “cv” schemes of various discretizations (“cv(1,1,0)”, “cv(2,1,0)”, “cv(3,2,0)”) and by the “op(3,2,0)” scheme, suggested by this study. It can be seen that all the “cv” schemes underestimate the range of ST. They systematically underestimate the ST during the day and overestimate it during the night. Since snow usually melts during the day, therefore, the “cv” scheme can produce a time lag for snowmelt simulation. Besides, the underestimation increases with the numerical thickness h ₀ of the skin layer. The “cv(1,1,0)” scheme, of which h ₀ is 0.91 m, underestimates the ST by 4 K at noon, while “cv(3,2,0)”, of which h ₀ is 0.064 m, underestimates it by 0.5 K. In contrast, the “op(3,2,0)” scheme accurately predicts the ST with RMSE at 0.02 K. Therefore, the error of the time-lag problem in snowmelt simulation can be reduced.

Fig. 9

Same as Fig. 7, but for the last day of the simulation (day 6) among “cv(1,1,0)”, “cv(2,1,0)”, “cv(3,2,0)” and “op(3,2,0)”, of which the physics thicknesses h ₀ of the skin layers are 0.91, 0.098, 0.064 and 0.019 m, respectively

6.3 Other effective thickness

The above section shows that setting effective thicknesses of the skin layer at the optimal effective thicknesses is the most accurate for ST simulation. Nonetheless, there are uncertainties with the value for ∂G ₀/∂T ₀ although it is suggested using a typical value of 42 W m⁻² K⁻¹ (the annual median value of a cropland site in Bondville, USA) for the purpose. Here, various effective thicknesses are derived.

For example, under two extreme elasticities (s ^* → 0, s ^* → ∞), i.e., (the non-stiff condition and the stiff condition), the effective thickness h ^* _e0 of the skin layer can be derived from Eq. 21 as:

$$ h_{e0}^{n*} = \mathop {\lim }\limits_{{s_{0}^{*} \to 0}} h_{e0}^{p*} = \frac{{h_{a0}^{*} }}{{\cos \left( {t_{a0}^{*} } \right)}} $$

(30)

$$ h_{e0}^{s*} = \mathop {\lim }\limits_{{s_{0}^{*} \to \infty }} h_{e0}^{p*} = \cos \left( {t_{a0}^{*} } \right)h_{a0}^{*} $$

(31)

where the first superscripts of h ^n* _e0 and h ^s* _e0   denote the optimal thickness derived under the “non-stiff” and “stiff” conditions, respectively. We name the effective thickness h ^n* _e0 as the optimal non-stiff thickness, and denote the numerical scheme, which sets h ^* _e0  = h ^n* _e0 , as “on”. And, we name the thickness h ^s* _e0   as the optimal stiff thickness, and denote the numerical scheme, which sets h ^* _e0  = h ^s* _e0  , as “os”. In addition, From the above two equations, it can be seen that h ^n* _e0  > h ^* _a0  > h  ^s* _e0 . Since h ^* _a0 is also not a function of ∂G ^*₀ /∂T ^*₀ , therefore, setting h ^* _e0 at h ^* _a0 , is another choice for the effective thickness. We name the numerical scheme, which sets h ^* _e0  = h ^* _a0 , as the non-stiff equal-amplitude scheme (denoted as “ne”) since its simulated amplitude of ST of the concerned frequency component j is equal to its true amplitude $\left( {\Updelta T^{*} = \sqrt 2 } \right)$ under the non-stiff condition according to Eq. 61.

Figure 3 also shows the accuracies of the “ne”, “on” and “os” schemes as functions of h ^*₀ under ∂G ^*₀ /∂T ^*₀  = 4 and 0.5. Table 8 lists the equations for determining the effective thickness of various schemes (“cv”, “nh”, “ne”, “on”, “os” and “op”) and their corresponding accuracies. The e(g ^*₀ ) of the “ne”, “on” and “os” schemes are calculated from Eq. 24 by setting h ^* _e0 at h ^n* _e0 and h ^s* _e0 , respectively. This figure shows that calculating STs of a particular frequency component using the schemes of “op” “ne”, “on” and “os” are always better than “cv” and “nh”. In addition, for high ∂G ₀ ^*/∂T ^*₀ (= 4 in this case study), the accuracy of “os” is the closest to “op”; but for low ∂G ^*₀ /∂T ^*₀ (= 0.5 in this case study), the accuracy of “on” is the closest to “op”. Moreover, the figure shows that for both high and low ∂G ^*₀ /∂T ^*₀ , the accuracy simulated by “ne” is only slightly worse than “op”. Therefore, the “ne” scheme is a good scheme for implementing to models for which ∂G ^*₀ /∂T ^*₀ is not immediately available. It can be shown that “A74” is a special case of “ne” for h ^*₀  → ∞, and “D78” is a special case of “on” for h ^*₀  → ∞.

Table 8 Characteristics of various effective thickness parameterizations for skin layer

6.4 Adding the skin layer onto the conventional finite difference scheme “cv + op”

Table 6 also shows that “ecmwf + op” is better than “ecmwf”, and “echam + op” is better than “echam”. For Case 7, e(G ^*₀ ) of “ecmwf + op” is at 10% while that of “ecmwf” is at 11%, and e(G ^*₀ ) of “echam + op” 15% while that of “echam” is at 19%. This shows that using the optimal effective thickness for the skin layer also improves the model system for surface ground heat flux simulation for both the state-of-the-art models “ECMWF” and “ECHAM”. This “cv + op” type scheme can be important for a land type such as ocean which is not suitable for using the effective thicknesses for all the ground layers, or for a model of which the model structure does not want to be changed significantly. In addition to being more accurate than “cv”, the benefits of using “cv + op” or “cv + ne” include: (1) keeping track of the energy budget of a land column in the layer-mean temperatures of a land column, and (2) retaining the same memory allocation system of most climate models. Note that most of the models store ST as well as the layer-mean temperatures of a land column [e.g., ECMWF (Viterbo and Beljaars 1995), ECHAM (Roeckner et al. 2003), and NCEP (Kanamitsu et al. 2002)]. The scheme “cv + ne” proposed in this study has been used for a turbulent kinetic energy ocean model (Tu and Tsuang 2005).

7 Conclusion

This study introduces a differential equation for calculating skin temperature, and derives the optimal effective thickness analytically by minimizing the error for the temperature simulation at each numerical layer. The optimal effective thickness of each numeral layer can be determined from Eq. 20. It shows that the effective thickness is always thinner than its physical thickness. The value of the effective thickness h ^p _e0 of the skin layer is a function of its physical thickness h ₀ and the temperature derivate ∂G ₀/∂T ₀ of surface ground heat flux (Fig. 4). The value of h ^p _e0 approaches to h ₀ when h ₀ → 0, and h ^p _e0 is fixed at a range within $\left( {0.5\sqrt {{{2D} \mathord{\left/ {\vphantom {{2D} {\omega_{d} }}} \right. \kern-\nulldelimiterspace} {\omega_{d} }}} ,\;\sqrt {{{2D} \mathord{\left/ {\vphantom {{2D} {\omega_{d} }}} \right. \kern-\nulldelimiterspace} {\omega_{d} }}} } \right),$ varying with ∂G ₀/∂T ₀, when h ₀ → ∞. Therefore, the assumptions of low or no heat capacity (i.e., h _e0 ≪ h ₀) for the skin layer when h ₀ → 0 are not good approximations for a homogeneous ground column. Nonetheless, it may be valid where the land surface is covered by vegetation or it consists of organic soil, of which the heat diffusivity is usually much lower than that of the underneath soil (Best 1998). The characteristics such as the thickness and the accuracies of the various effective thicknesses (“cv”, “nh”, “ne”, “on”, “os” and “op”), especially under h ^*₀  → 0 and h ^*₀  → ∞, are listed in Table 8.

The most beneficial scheme is “op”. Table 4 lists the vertical discretizations of the “op” scheme of 1–6 numerical levels for calculating Earth’s skin temperature. The suggested discretization is derived from the evenly heat-content discretization with the optimal effective thickness for layer-temperature simulation. For the same level number, the suggested discretization is more accurate in skin temperature as well as surface ground heat flux simulations than those used in some state-of-the-art models. The proposed scheme (“op(3,2,0)”) has shown to be more accurate than the schemes used in state-of-art climate models including ECMWF (Viterbo and Beljaars 1995), ECHAM (Roeckner et al. 2003) and the UCLA GCM (Arakawa and Mintz 1974). The profiles of diurnal and annual ground temperatures are recoded in the middle layers of “op(3,2,0)”. This type of arrangement is important for reducing the error of the corresponding frequency component. In addition, it is found that “cv” systematically underestimates ST during the day. The underestimation can be as high as 4 K for h ₀ at 0.91 m. Since snow usually melts during the day, the “cv” scheme can cause a time lag for snowmelt. In contrast, the “op(3,2,0)” scheme, proposed by this study, accurately predicts the ST with RMSE at 0.02 K. Therefore, the error in the time lag can be reduced. Nonetheless, it should be noted that we are not able to prove the evenly heat-content discretization to be the optimal vertical discretization. A better vertical discretization may exist. A tri-diagonal matrix for solving the temperatures by the “op” scheme is illustrated in the Appendix 4.

The introduction of an additional differential equation for calculating skin temperature is also found beneficial for the temperature simulation. For the same vertical structure (3,2,0), “cv + op” is better than “cv + nh”, and “cv + nh” is better than “cv”. In addition, “ecmwf + op” is better than “ecmwf”, and “echam + op” is better than “echam” (Table 6). Although the proposed “op” scheme can be easily implemented into state-of-the-art climate models for the temperature simulation of snow, ice sheet and soil, nonetheless, it should be reminded that the effective thickness is derived based on the assumption that heat source is from the surface only. If horizontal advected heat flux is important, such as for ocean water temperature simulation, “cv + op” can be a better option than “op”. In addition, the introduction of effective thickness for each numerical layer, which is different from the real layer thickness, causes the energy conservation equation to be different from the conventional form. From Eq. 1, it can be easily proved that the heat content H (J) of an entire ground column from infinite depth (z = −∞) to the surface (z = 0) can be determined using a modified form as:

$$ \frac{\partial H}{\partial t} = \frac{\partial }{\partial t}\left( {\int\limits_{ - \infty }^{0} {\rho_{g} c_{g} T{\text{d}}z} } \right) = \frac{\partial }{\partial t}\left( {\sum\limits_{k = 0}^{m} {\rho_{g} c_{g} T_{k} h_{ek} } } \right) $$

(32)

Appendices

Appendix 1: Analytical ground temperature and heat flux equations

Considering an ideal surface, in which the heat diffusion coefficient is constant and neglecting the horizontal heat transport, the heat transfer in the ground can be assumed to obey the Fourier law of diffusion (e.g., Hillel 1982) as:

$$ \frac{\partial T}{\partial t} = \frac{ - 1}{{\rho_{g} c_{g} }}\frac{\partial G}{\partial z} = \frac{\partial }{\partial z}D\frac{\partial T}{\partial z} = D\frac{{\partial^{2} T}}{{\partial z^{2} }} $$

(33)

and

$$ G = - \rho_{g} c_{g} D\frac{\partial T}{\partial z} $$

(34)

where T is ground temperature (K), t is time (s), G is ground heat flux (W m⁻²) (positive upward), D is heat diffusivity of ground surface (m² s⁻¹), and z is the vertical coordinate system (m) (positive upward).

Considering the upper boundary temperature (denoted as T _s) at z = 0, i.e., the upper boundary condition of Eq. 33, can be described in the frequency domain (e.g., Tsuang 2003) as:

$$ T_{s} = \overline{{T_{s} }} + \sum\limits_{j = 1}^{j = \infty } {\Updelta T_{sj} \cos \left( {\omega_{j} \left( {t - t_{mj} } \right)} \right)} $$

(35)

where the overbar “^–” is the average, the subscript “j” denotes a property of frequency component j, ΔT _sj is the amplitude of T _s of the particular frequency component, ω _j is the angular velocity of the frequency component, t _mj is the time when the highest T _s of the particular frequency component occurs, and t is local time.

The analytical solution of the temperature profile of Eq. 33 can be determined (after Carslaw and Jaeger 1959) as:

$$ \begin{gathered} T\left( {z,t} \right) = \overline{{T_{s} }} + \sum\limits_{j = 1}^{j = \infty } {\Updelta T_{sj} \exp \left( {\sqrt {\frac{{\omega_{j} }}{2D}} z} \right)\cos \left( {\omega_{j} \left( {t - t_{mj} } \right) + \sqrt {\frac{{\omega_{j} }}{2D}} z} \right)} \\ = \sum\limits_{j = 0}^{j = \infty } {\tau_{j} \left( {z,t} \right)} \\ \end{gathered} $$

(36)

where τ _j is the ground temperature of frequency component j. For j = 0, $\tau_{0} = \overline{{T_{s} }} .$ For j ≥ 1, τ _j can be determined by observing the second equality of the above equation as

$$ \tau_{j} \left( {z,t} \right) = \Updelta T_{sj} \exp \left( {\sqrt {\frac{{\omega_{j} }}{2D}} z} \right)\cos \left( {\omega_{j} \left( {t - t_{mj} } \right) + \sqrt {\frac{{\omega_{j} }}{2D}} z} \right) $$

(37)

Note that Carslaw and Jaeger (1959) only proved the above solution (A4) with a single frequency component. Nonetheless, the above form (A4) is also valid for multiple frequency components. It can be proven by substituting the above equation to Eq. 33 to check the equality between the right hand side (RHS) and the left hand side (LHS) of Eq. 33.

Substituting the above equation into Eq. 34, the analytical solution of ground heat flux can be determined as:

$$ \begin{aligned} G\left( {z,t} \right) = & \sum\limits_{j = 1}^{j = \infty } { - \rho_{g} c_{g} \Updelta T_{sj} \sqrt {D\omega_{j} } \exp \left( {\sqrt {\frac{{\omega_{j} }}{2D}} z} \right)\cos \left( {\omega_{j} \left( {t - t_{mj} } \right) + \sqrt {\frac{{\omega_{j} }}{2D}} z + \frac{\pi }{4}} \right)} \\ = & \sum\limits_{j = 0}^{j = \infty } {g_{j} \left( {z,t} \right)} \\ \end{aligned} $$

(38)

where g _j is the ground heat flux of frequency component j. For j = 0, g ₀ = 0. For j ≥ 1, g _j can be determined by observing the second equality of the above equation as

$$ g_{j} \left( {z,t} \right) = - \Updelta g_{j} \exp \left( {\sqrt {\frac{{\omega_{j} }}{2D}} z} \right)\cos \left( {\omega_{j} \left( {t - t_{mj} } \right) + \sqrt {\frac{{\omega_{j} }}{2D}} z + \frac{\pi }{4}} \right) $$

(39)

and

$$ \Updelta g_{j} \equiv \rho_{g} c_{g} \Updelta T_{sj} \sqrt {D\omega_{j} } $$

(40)

From the above equation, the dominant frequency component can be identified by choosing the maximum value among Δg _j of various frequency components j.

The STD of T _s of frequency component j (denoted as τ _jSTD) can be determined from Eq. 36 as:

$$ \tau_{{0j{\text{STD}}}} \equiv \sqrt {\overline{{\left( {\tau_{j} \left( {0,t} \right) - \overline{{\tau_{j} \left( {0,t} \right)}} } \right)^{2} }} } = \frac{1}{\sqrt 2 }\Updelta T_{sj} $$

(41)

And the STD of surface ground heat flux, where z = 0, of the frequency component (denoted as g _jSTD) can be determined from Eq. 39 as:

$$ g_{{0j{\text{STD}}}} \equiv \sqrt {\overline{{\left( {g_{j} \left( {0,t} \right) - \overline{{g_{j} \left( {0,t} \right)}} } \right)^{2} }} } = \sqrt {\frac{1}{2}\Updelta g_{j}^{2} } = \frac{{\rho_{g} c_{g} \sqrt {D\omega_{j} } \Updelta T_{sj} }}{\sqrt 2 } = \rho_{g} c_{g} \sqrt {D\omega_{j} } \tau_{{j{\text{STD}}}} $$

(42)

Appendix 2: Temperature derivative of ground heat flux

This appendix shows the form of the temperature derivative of ground heat flux ∂G/∂T. G ₀ can be determined from the energy budget of the land surface (positive upward) (e.g., Brutsaert 1982) as:

$$ G_{0} = - \left[ {\left( {1 - \alpha } \right)R_{s} + R_{\text{ld}} - R_{\text{lu}} - H - {\rm LE}} \right] $$

(43)

where α is albedo, R _s is incoming solar radiation (positive downward), R _ld is incoming atmospheric radiation (positive downward), R _lu is outgoing terrestrial radiation (positive upward), H is surface sensible heat flux (positive upward) and LE is surface latent heat flux (positive upward). R _lu, H and LE are functions of ST. They can be determined (e.g., Brutsaert 1982; Garratt 1992) as:

$$ R_{\text{lu}} = \varepsilon \sigma T_{0}^{4} $$

(44)

$$ H = \frac{{\rho_{a} c_{a} }}{{r_{a} }}\left( {T_{0} - T_{a} } \right) $$

(45)

$$ {\text{LE}} = \frac{{\rho_{a} L_{v} }}{{r_{a} + r_{c} }}\left( {q^{*} \left( {T_{0} } \right) - q_{a} } \right) $$

(46)

where ε is emissivity of surface; σ is Stefan–Boltzman constant (~5.67 × 10⁻⁸ W m⁻² K⁻⁴); ρ _a, c _a, T _a and q _a are density (~1.16 kg m⁻³), constant pressure heat capacity (~1,005 J kg⁻¹ K⁻¹), temperature (K) and specific humidity (kg kg⁻¹) of air, respectively; q ^*(T) is saturated specific humidity at temperature T (= 0.622 e ^*(T)/P), where P is atmospheric pressure (Pa) and e ^* is saturated vapor pressure (Pa) (Richards 1971); r _a is aerodynamic resistance (s m⁻¹); r _c is canopy resistance (s m⁻¹) for evapotranspiration; and L _v is latent heat of evaporation (~2.5 × 10⁶ J/kg). Substituting the above three equations into Eq. 43, G ₀ can be written as a function of ST as:

$$ G_{0} \left( {T_{0} } \right) = - \left[ {\left( {1 - \alpha } \right)R_{s} + R_{\text{ld}} - \varepsilon \sigma T_{0}^{4} - \frac{{\rho_{a} c_{a} }}{{r_{a} }}\left( {T_{0} - T_{a} } \right) - \frac{{\rho_{a} L_{v} }}{{r_{a} + r_{c} }}\left( {q^{*} \left( {T_{0} } \right) - q_{a} } \right)} \right] $$

(47)

The above equation shows surface ground heat flux is a function of ST. Therefore, ∂G ₀/∂T ₀ can be derived by taking the derivative of the above equation on ST as:

$$ \frac{{\partial G_{0} \left( {T_{0} } \right)}}{{\partial T_{0} }} = 4\varepsilon \sigma T_{0}^{3} + \frac{{\rho_{a} c_{a} }}{{r_{a} }} + \frac{{\rho_{a} L_{v} }}{{r_{a} + r_{c} }}\frac{{\partial q^{*} \left( {T_{0} } \right)}}{{\partial T_{0} }} $$

(48)

In addition, from Eq. 34, the ground heat flux at level k can be determined conventionally by the finite difference scheme (Fig. 1) as:

$$ G_{k} = - \rho_{g} c_{g} D\frac{{T_{k - 1} - T_{k} }}{{z_{k - 1} - z_{k} }} + {\text{H}} . {\text{O}} . {\text{T}}. $$

(49)

Therefore, ∂G _k /∂T _k−1 and ∂G _k /∂T _k can be derived by taking the derivative of the above equation on T _k−1 and T _k as:

$$ \frac{{\partial G_{k} }}{{\partial T_{k - 1} }} = - \frac{{\rho_{g} c_{g} D}}{{z_{k - 1} - z_{k} }} = - s_{k} $$

(50)

$$ \frac{{\partial G_{k} }}{{\partial T_{k} }} = \frac{{\rho_{g} c_{g} D}}{{z_{k - 1} - z_{k} }} $$

(51)

Appendix 3: Nondimensionalization and optimal effective thickness

To derive a general form for a particular frequency component j ≥ 1, equations of the frequency component are made nondimensional, by multiplying time by the angular velocity ω _j of the frequency component, dividing the length by $\sqrt {2D/\omega_{j} } ,$ dividing energy flux by the STD of the ground heat flux component g _kj at the upper boundary (ζ = 0), i.e., g _kjSTD, and dividing temperature component τ _kj by its STD at the upper boundary (ζ = 0), i.e., τ _kjSTD. In addition denoting nondimensional quantities of the frequency component with asterisks, but without the level index k and the frequency index j for simplicity (although occasionally we will recover these indices k and j in the text for clarity), we obtain:

$$ \tau^{*} \equiv \frac{{\tau_{kj} }}{{\tau_{{kj{\text{STD}}}} }} = \frac{{\sqrt 2 \tau_{kj} }}{{\Updelta T_{skj} }},\quad t^{*} \equiv \omega_{j} t,\quad g^{*} \equiv \frac{{g_{kj} }}{{g_{{kj{\text{STD}}}} }} = \frac{{\sqrt 2 g_{kj} }}{{\rho_{g} c_{g} \sqrt {D\omega_{j} } \Updelta T_{skj} }},\quad \zeta^{*} \equiv \frac{\zeta }{{\sqrt {2D/\omega_{j} } }} $$

(52)

where $\tau_{{kj{\text{STD}}}} = {{\Updelta T_{skj} } \mathord{\left/ {\vphantom {{\Updelta T_{skj} } {\sqrt 2 }}} \right. \kern-\nulldelimiterspace} {\sqrt 2 }}$ and $g_{{kj{\text{STD}}}} = \rho_{g} c_{g} \sqrt {D\omega_{j} } {{\Updelta T_{skj} } \mathord{\left/ {\vphantom {{\Updelta T_{skj} } {\sqrt 2 }}} \right. \kern-\nulldelimiterspace} {\sqrt 2 }}$ according to Eqs. 41 and 42, respectively. As a result, the dimensionless forms of the analytical solution (Eq. 37) of the ground temperature of the frequency component and its analytical solution (Eq. 39) of the ground heat flux profile of the component can be rewritten using dimensionless variables as:

$$ \tau^{*} \left( {\zeta^{*} ,t^{*} } \right) = \sqrt 2 \exp \left( {\zeta^{*} } \right)\cos \left( {t^{*} - t_{m}^{*} + \zeta^{*} } \right) $$

(53)

$$ g^{*} \left( {\zeta^{*} ,t^{*} } \right) = - \sqrt 2 \exp \left( {\zeta^{*} } \right)\cos \left( {t^{*} - t_{m}^{*} + \zeta^{*} + \frac{\pi }{4}} \right) $$

(54)

Note that from Eq. 53, it can be seen that the amplitude of the true dimensionless skin temperature component τ ^* is $\sqrt 2 .$ And the dimensionless forms of ∂G ₀/∂T ₀ and the elasticity s _k can be written as:

$$ x_{0}^{*} = \frac{{\partial G_{0}^{*} }}{{\partial T_{0}^{*} }} = \frac{{\tau_{{0j{\text{STD}}}} }}{{g_{{0j{\text{STD}}}} }}\frac{{\partial G_{0} }}{{\partial T_{0} }} = \frac{1}{{\rho_{g} c_{g} \sqrt {D\omega_{j} } }}\frac{{\partial G_{0} }}{{\partial T_{0} }} $$

(55a)

$$ x_{k}^{*} = - \frac{{\partial G_{k}^{*} }}{{\partial T_{k - 1}^{*} }} = \frac{{\partial G_{k}^{*} }}{{\partial T_{k}^{*} }} = \frac{{\rho_{g} c_{g} D}}{{\rho_{g} c_{g} \sqrt {D\omega_{j} } \left( {z_{k - 1} - z_{k} } \right)}} = \frac{1}{{\sqrt 2 \left( {z_{k - 1}^{*} - z_{k}^{*} } \right)}},\quad {{k}} = 1,\;{\text{m}} $$

(55b)

$$ x_{m + 1}^{*} = 0 $$

(55c)

and

$$ s^{*} \equiv x_{k}^{*} + x_{k + 1}^{*} $$

(56)

where the second equalities are derived by substituting Eqs. 41 and 42 into 5. The above Eq. 55a implies that if ∂G ₀/∂T ₀ remains constant, the lower the frequency, the higher the value of ∂G ^*₀ /∂T ^*₀ will be, i.e., approaching the stiff condition.

In addition, the dimensionless form of the governing Eq. 10 for determining T _k of the frequency component by a numerical method can be derived by multiplying it by ${{\sqrt 2 } \mathord{\left/ {\vphantom {{\sqrt 2 } {\omega_{j} \Updelta T_{skj} }}} \right. \kern-\nulldelimiterspace} {\omega_{j} \Updelta T_{skj} }}.$ The form, then, can be written as:

$$ \begin{gathered} \frac{{\partial \tau_{n}^{*} }}{{\partial t^{*} }} + \frac{{s^{*} }}{{\sqrt 2 h_{e}^{*} }}\tau_{n}^{*} = - \frac{1}{{\sqrt 2 h_{e}^{*} }}\left[ {g^{*} \left( {0,t^{*} } \right) - g^{*} \left( { - h^{*} ,t^{*} } \right) - s^{*} \tau^{*} \left( { - h_{t}^{*} ,t^{*} } \right)} \right] \hfill \\ \quad = \frac{1}{{\sqrt 2 h_{e}^{*} }}\left[ \begin{gathered} \sqrt 2 \cos \left( {t^{*} - t_{m}^{*} + \frac{\pi }{4}} \right) - \sqrt 2 \exp \left( { - h^{*} } \right)\cos \left( {t^{*} - t_{m}^{*} - h^{*} + \frac{\pi }{4}} \right) \hfill \\ + \sqrt 2 s^{*} \exp \left( { - h_{T}^{*} } \right)\cos \left( {t^{*} - t_{m}^{*} - h_{t}^{*} } \right) \hfill \\ \end{gathered} \right] \hfill \\ \quad = - \sqrt 2 \frac{{h_{a}^{*} }}{{h_{e}^{*} }}\sin \left( {t^{*} - t_{m}^{*} - t_{a}^{*} } \right) + \frac{{s^{*} }}{{h_{e}^{*} }}\exp \left( { - h_{T}^{*} } \right)\cos \left( {t^{*} - t_{m}^{*} - h_{t}^{*} } \right) \hfill \\ \end{gathered} $$

(57)

where $h^{*} \equiv {{h_{k} } \mathord{\left/ {\vphantom {{h_{k} } {\sqrt {2D/\omega_{j} } }}} \right. \kern-\nulldelimiterspace} {\sqrt {2D/\omega_{j} } }}.$ Note that the second equality is derived by substituting Eqs. 53 and 54 into it, which assumes that the ground heat fluxes at the surface and at depth h _k are accurately provided as inputs. And,

$$ h_{a}^{*} \left( {h^{*} } \right) = \frac{1}{\sqrt 2 }\sqrt {1 - 2\cos \left( {h^{*} } \right)\exp \left( { - h^{*} } \right) + \exp \left( { - h^{*} } \right)^{2} } $$

(58)

$$ t_{a}^{*} \left( {h^{*} } \right) = \frac{\pi }{4} - \tan^{ - 1} \left( {\frac{{\exp \left( { - h^{*} } \right)\sin \left( {h^{*} } \right)}}{{1 - \exp \left( { - h^{*} } \right)\cos \left( {h^{*} } \right)}}} \right) $$

(59)

Under the condition that s ^* is constant, Eq. 57 is a first-order linear ordinary differential equation with constant coefficients (e.g., Kreyszig 2006), which can be solved analytically as:

$$ \begin{aligned} \tau_{n}^{*} & = \frac{1}{{2h_{e}^{*2} + s^{*2} }}\left[ {2\sqrt 2 h_{a}^{*} h_{e}^{*} \cos \left( {t^{*} - t_{m}^{*} - t_{a}^{*} } \right) - 2s^{*} h_{a}^{*} \sin \left( {t^{*} - t_{m}^{*} - t_{a}^{*} } \right)} \right. \\ & \quad \left. { + \sqrt 2 s^{*2} \exp \left( { - h_{t}^{*} } \right)\cos \left( {t^{*} - t_{m}^{*} - h_{t}^{*} } \right) + 2s^{*} h_{e}^{*} \exp \left( { - h_{t}^{*} } \right)\sin \left( {t^{*} - t_{m}^{*} - h_{t}^{*} } \right)} \right] \\ & = \Updelta T_{n}^{*} \cos \left( {t^{*} - t_{m}^{*} - t_{l}^{*} } \right) \\ \end{aligned} $$

(60)

where ΔT _n ^* and t ^* _l are the amplitude and the phase lag of the simulated T _k,n of the frequency component. They are

$$ \Updelta T_{n}^{*} = \sqrt 2 \sqrt {\frac{{2h_{a}^{*2} + 2\sqrt 2 \sin \left( {t_{a}^{*} - h_{t}^{*} } \right)h_{a}^{*} s^{*} \exp \left( { - h_{t}^{*} } \right) + s^{*2} \exp \left( { - h_{t}^{*} } \right)^{2} }}{{2h_{e}^{*2} + s^{*2} }}} $$

(61)

$$ t_{l}^{*} = \tan^{ - 1} \left[ {\frac{\begin{gathered} 2\sqrt 2 h_{a}^{*} h_{e}^{*} \sin \left( {t_{a}^{*} } \right) - 2s^{*} h_{a}^{*} \cos \left( {t_{a}^{*} } \right) \hfill \\ - \sqrt 2 s^{*2} \exp \left( { - h_{t}^{*} } \right)\sin \left( { - h_{t}^{*} } \right) + 2s^{*} h_{e}^{*} \exp \left( { - h_{t}^{*} } \right)\cos \left( { - h_{t}^{*} } \right) \hfill \\ \end{gathered} }{\begin{gathered} 2\sqrt 2 h_{a}^{*} h_{e}^{*} \cos \left( {t_{a}^{*} } \right) + 2s^{*} h_{a}^{*} \sin \left( {t_{a}^{*} } \right) \hfill \\ + \sqrt 2 s^{*2} \exp \left( { - h_{t}^{*} } \right)\cos \left( { - h_{t}^{*} } \right) + 2s^{*} h_{e}^{*} \exp \left( { - h_{t}^{*} } \right)\sin \left( { - h_{t}^{*} } \right) \hfill \\ \end{gathered} }} \right] $$

(62)

As a result, the difference (denoted as ɛ(τ ^*)) between calculated T _k,n and the true T _k of the frequency component can be written as a cosine function from Eqs. 60 and 53 as:

$$ \begin{aligned} \varepsilon \left( {\tau^{*} } \right) \equiv & \tau_{n}^{*} - \tau^{*} \left( {h_{t}^{*} ,t^{*} } \right) \\ = & \; \frac{1}{{2h_{e}^{*2} + s^{*2} }}\left[ \begin{gathered} 2\sqrt 2 h_{a}^{*} h_{e}^{*} \cos \left( {t^{*} - t_{m}^{*} - t_{a}^{*} } \right) - 2s^{*} h_{a}^{*} \sin \left( {t^{*} - t_{m}^{*} - t_{a}^{*} } \right) \hfill \\ + \sqrt 2 s^{*2} \exp \left( { - h_{t}^{*} } \right)\cos \left( {t^{*} - t_{m}^{*} - h_{t}^{*} } \right) + 2s^{*} h_{e}^{*} \exp \left( { - h_{t}^{*} } \right)\sin \left( {t^{*} - t_{m}^{*} - h_{t}^{*} } \right) \hfill \\ \end{gathered} \right] \\ & - \sqrt 2 \exp \left( { - h_{t}^{*} } \right)\cos \left( {t^{*} - t_{m}^{*} - h_{t}^{*} } \right) \\ = & \; \Updelta T_{\varepsilon }^{*} \cos \left( {t^{*} - t_{m}^{*} - t_{\varepsilon }^{*} } \right) \\ \end{aligned} $$

(63)

where the second equality is derived by substituting Eqs. 60 and 53 into it, and ΔT _ɛ ^* and t ^* _ɛ are determined as:

$$ \Updelta T_{\varepsilon }^{*} = 2\sqrt {\frac{{h_{a}^{*2} - 2h_{a}^{*} h_{e}^{*} \exp \left( { - h_{t}^{*} } \right)\cos \left( {t_{a}^{*} - h_{t}^{*} } \right) + \exp \left( { - h_{t}^{*} } \right)^{2} h_{e}^{*2} }}{{2h_{e}^{*2} + s^{*2} }}} $$

(64)

$$ t_{\varepsilon }^{*} = \tan^{ - 1} \left[ {\frac{\begin{gathered} 2\sqrt 2 h_{a}^{*} h_{e}^{*} \sin \left( {t_{a}^{*} } \right) - 2s^{*} h_{a}^{*} \cos \left( {t_{a}^{*} } \right) \hfill \\ + \sqrt 2 \exp \left( { - h_{t}^{*} } \right)2h_{e}^{*2} \sin \left( { - h_{t}^{*} } \right) + 2s^{*} h_{e}^{*} \exp \left( { - h_{t}^{*} } \right)\cos \left( { - h_{t}^{*} } \right) \hfill \\ \end{gathered} }{\begin{gathered} 2\sqrt 2 h_{a}^{*} h_{e}^{*} \cos \left( {t_{a}^{*} } \right) + 2s^{*} h_{a}^{*} \sin \left( {t_{a}^{*} } \right) \hfill \\ - \sqrt 2 \exp \left( { - h_{t}^{*} } \right)2h_{e}^{*2} \cos \left( { - h_{t}^{*} } \right) + 2s^{*} h_{e}^{*} \exp \left( { - h_{t}^{*} } \right)\sin \left( { - h_{t}^{*} } \right) \hfill \\ \end{gathered} }} \right] $$

(65)

Therefore, the normalized root–mean–square error (NRMSE) (or RMSE/STD ratio) e(τ ^*) of the calculated T _k,n of the frequency component compared to the true T _k can be determined from Eq. 63 as:

$$ \begin{gathered} e\left( {\tau^{*} } \right) \equiv \sqrt {\overline{{\left( {\tau_{n}^{*} - \tau^{*} \left( {h_{t}^{*} ,t^{*} } \right)} \right)^{2} }} } = \sqrt {\overline{{\varepsilon \left( {\tau^{*} \left( {h_{t}^{*} ,t^{*} } \right)} \right)^{2} }} } = \sqrt {\overline{{\left[ {\Updelta T_{\varepsilon }^{*} \cos \left( {t^{*} - t_{m}^{*} - t_{\varepsilon }^{*} } \right)} \right]^{2} }} } = \sqrt {0.5\Updelta T_{\varepsilon }^{*2} } \\ = \sqrt {\frac{{h_{a}^{*2} - 2h_{a}^{*} h_{e}^{*} \exp \left( { - h_{t}^{*} } \right)\cos \left( {t_{a}^{*} - h_{t}^{*} } \right) + \exp \left( { - h_{t}^{*} } \right)^{2} h_{e}^{*2} }}{{h_{e}^{*2} + 0.5s^{*2} }}} \\ \end{gathered} $$

(66)

where the fourth equality is derived, since $\overline{{\cos \left( {t^{*} - t_{m}^{*} - t_{\varepsilon }^{*} } \right)^{2} }} = 0. 5.$ It shows that e(τ ^*) is a function of h ^*, h ^* _e , h ^* _t and s ^* only. Note that h ^* _a and t ^* _a are functions of h ^* according to Eqs. 58 and 59.

In order to minimize the NRMSE τ ^*, the optimal effective thickness h ^* _e can be determined by taking the derivative of e(τ ^*) from Eq. 66 and letting:

$$ \frac{{\partial e\left( {\tau^{*} } \right)}}{{\partial h_{e}^{*} }} = 0 $$

(67)

Solving the above equation, the optimal value of h ^* _e (denoted as h ^p* _e ) is determined as:

$$ h_{e}^{p*} = \frac{{2h_{a}^{*2} - \exp \left( { - h_{t}^{*} } \right)^{2} s^{*2} + \sqrt {4h_{a}^{*4} + 4\cos 2\left( {t_{a}^{*} - h_{t}^{*} } \right)\exp \left( { - h_{t}^{*} } \right)^{2} h_{a}^{*2} s^{*2} + \exp \left( { - h_{t}^{*} } \right)^{4} s^{*4} } }}{{4\cos \left( {t_{a}^{*} - h_{t}^{*} } \right)\exp \left( { - h_{t}^{*} } \right)h_{a}^{*} }} $$

(68)

where the negative value is discarded.

Appendix 4: Numerical matrix for the “optimal” scheme

This appendix describes the numerical parameterization in this study. A numerical grid is chosen to discretize the land column as shown in Fig. 1. That is, T is determined at the center of the grid, and flux G _k is determined at the boundaries of the grid. The ST and the ground temperature of each layer are parameterized according to Eq. 27. The resulting temperature matrix is:

$$ \begin{gathered} \left[ {\begin{array}{*{20}c} {1 + \beta y_{0} } & { - \beta y_{0} } & {} & {} & {} & {} & {} & {} & {} \\ {} & { - \beta x_{1} } & {1 + \beta x_{1} + \beta y_{1} } & { - \beta y_{1} } & {} & {} & {} & {} & {} \\ {} & {} & {} & {} & {} & {} & {} & {} & {} \\ {} & {} & {} & {} & {} & {} & {} & {} & {} \\ {} & {} & {} & { - \beta x_{k} } & {1 + \beta x_{k} + \beta y_{k} } & { - \beta y_{k} } & {} & {} & {} \\ {} & {} & {} & {} & {} & {} & {} & {} & {} \\ {} & {} & {} & {} & {} & {} & {} & {} & {} \\ {} & {} & {} & {} & {} & {} & {} & {} & {} \\ {} & {} & {} & {} & {} & {} & {} & { - \beta x_{m} } & {1 + \beta x_{m} } \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} {T_{0}^{j + 1} } \\ {T_{1}^{j + 1} } \\ {} \\ {T_{k - 1}^{j + 1} } \\ {T_{k}^{j + 1} } \\ {T_{k + 1}^{j + 1} } \\ {} \\ {T_{m - 1}^{j + 1} } \\ {T_{m}^{j + 1} } \\ \end{array} } \right] = \hfill \\ \left[ {\begin{array}{*{20}c} {T_{0}^{j} - \frac{{\Updelta t \times G_{0} }}{{\rho_{g} c_{g} h_{e0} }} - \left( {1 - \beta } \right)y_{0} \left( {T_{0}^{j} - T_{1}^{j} } \right)} \\ {T_{1}^{j} + \left( {1 - \beta } \right)\left[ {x_{1} \left( {T_{0}^{j} - T_{1}^{j} } \right) - y_{1} \left( {T_{1}^{j} - T_{2}^{j} } \right)} \right]} \\ {} \\ {} \\ {T_{k}^{j} + \left( {1 - \beta } \right)\left[ {x_{k} \left( {T_{k - 1}^{j} - T_{k}^{j} } \right) - y_{k} \left( {T_{k}^{j} - T_{k + 1}^{j} } \right)} \right]} \\ {} \\ {} \\ {} \\ {T_{m}^{j} + \left( {1 - \beta } \right)x_{m} \left( {T_{m - 1}^{j} - T_{m}^{j} } \right)} \\ \end{array} } \right] \hfill \\ \end{gathered} $$

(69)

where the variables in Eq. 69 are described as follows:

$$ y_{0} = \frac{\Updelta t}{{h_{e0} }}\frac{D}{{z_{0} - z_{1} }} $$

(70)

$$ x_{k} = \frac{\Updelta t}{{h_{ek} }}\left( {\frac{D}{{z_{k - 1} - z_{k} }}} \right),\quad {\text{for}}\;k = 1,\;m - 1$$

(71)

$$ y_{k} = \frac{\Updelta t}{{h_{ek} }}\left( {\frac{D}{{z_{k} - z_{k + 1} }}} \right),\quad {\text{for}}\;k = 1,\;m - 1$$

(72)

$$ x_{m} = \frac{\Updelta tD}{{h_{em} \left( {z_{m - 1} - z_{m} } \right)}} $$

(73)

where the coordinates z _k are listed in Table 4, varying with discretization chosen, and the effective thickness of h _ek is parameterized according to Eq. 20. The values of h _k and h _tk is determined as:

$$ h_{0} = 0.5\left( {z_{0} - z_{1} } \right) , $$

(74)

$$ h_{k} = 0.5\left( {z_{k - 1} - z_{k + 1} } \right),\quad {\text{for}}\;k = 1,\;m - 1 $$

(75)

$$ h_{t0} = 0 $$

(76)

$$ h_{tk} = 0.5\left( {z_{k - 1} - z_{k} } \right),\quad {\text{for}}\;k = 1,\;m $$

(77)

Note that β = 0 for the forward scheme, β = 0.5 for the Crank-Nicolson scheme, and β = 1 for the backward scheme. The backward scheme is found to be desirable since it is numerical unconditionally stable, and the backward scheme has been implemented into a climate model (Tsuang et al. 2001). The matrix is a tri-diagonal matrix and can be easily and efficiently solved by the LU method (e.g., Press et al. 1992).