The Analytical Objective Hysteresis Model (AnOHM v1.0): methodology to determine bulk storage heat flux coefficients

Abstract

The net storage heat flux ($\mathrm{\Delta }{Q}_{\mathrm{S}})$ is important in the urban surface energy balance (SEB) but its determination remains a significant challenge. The hysteresis pattern of the diurnal relation between the $\mathrm{\Delta }{Q}_{\mathrm{S}}$ and net all-wave radiation ($Q^{\ast })$ has been captured in the Objective Hysteresis Model (OHM) parameterization of $\mathrm{\Delta }{Q}_{\mathrm{S}}$. Although successfully used in urban areas, the limited availability of coefficients for OHM hampers its application. To facilitate use, and enhance physical interpretations of the OHM coefficients, an analytical solution of the one-dimensional advection–diffusion equation of coupled heat and liquid water transport in conjunction with the SEB is conducted, allowing development of AnOHM (Analytical Objective Hysteresis Model). A sensitivity test of AnOHM to surface properties and hydrometeorological forcing is presented using a stochastic approach (subset simulation). The sensitivity test suggests that the albedo, Bowen ratio and bulk transfer coefficient, solar radiation and wind speed are most critical. AnOHM, driven by local meteorological conditions at five sites with different land use, is shown to simulate the $\mathrm{\Delta }{Q}_{\mathrm{S}}$ flux well (RMSE values of $\sim$ 30 W m${}^{-\mathrm{2}})$. The intra-annual dynamics of OHM coefficients are explored. AnOHM offers significant potential to enhance modelling of the surface energy balance over a wider range of conditions and land covers.

Introduction

The essential role of an integrated land surface model is to physically predict the land–atmosphere interactions by resolving the transfer of energy, water, and trace gases (Katul et al., 2012; Liang et al., 1994; Sellers et al., 1997). Such land–atmospheric interactions are strongly modulated by the partitioning of solar energy at the land surface (Chen and Dudhia, 2001; McCumber and Pielke, 1981; Yang and Wang, 2014) which can be considered through the surface energy balance (SEB) equation (Oke, 1988): $${\displaystyle }{Q}^{\ast }-\mathrm{\Delta }{Q}_{\mathrm{S}}=Q_H+Q_E,$$ where $Q^{\ast }$, $\mathrm{\Delta }{Q}_{\mathrm{S}}$, $Q_H$, and $Q_E$ are the net all-wave radiation, net storage, turbulent sensible, and latent heat fluxes, respectively. Equation (1) distinguishes the available energy at the land surface (left-hand side) from the heat transfer through turbulent transport (right-hand side).

The turbulent and radiative fluxes ($Q^{\ast }$, $Q_H$, and $Q_E)$ are more readily measured using standard techniques (e.g. eddy-covariance instruments, radiometry) than $\mathrm{\Delta }{Q}_{\mathrm{S}}$ (Offerle et al., 2005; Pauwels and Daly, 2016; Roberts et al., 2006; Wang, 2012). For $\mathrm{\Delta }{Q}_{\mathrm{S}}$ the net energy stored or released by changes in sensible heat within the canopy air layer, roughness elements (e.g. vegetation, buildings in an urban environment), and the ground all have to be considered. The volume of interest extends from the top of the roughness sub-layer to the depth in the ground where the daily averaged vertical net heat conduction is zero (see Fig. 2 in Masson et al., 2002); this presents very significant challenges of spatial sampling.

Knowledge of $\mathrm{\Delta }{Q}_{\mathrm{S}}$ is crucial to a wide range of processes and applications: from modelling turbulent heat transfer and boundary layer development to predicting soil thermal fields. In rural sites, or simple bare soil sites, the flux may be a small fraction of the net all-wave radiation (Oliphant et al., 2004). However, in areas where there is more mass, such as cities, the term becomes much more significant (Kotthaus and Grimmond, 2014a) and a key element of the SEB and well-known effects such as the urban heat island.

In urban systems a wide range of techniques have been used to estimate $\mathrm{\Delta }{Q}_{\mathrm{S}}$ (Grimmond et al., 1991; Roberts et al., 2006). These include the following:

• a.

  Heat conduction approach: the weighted average of heat flows through all urban materials and surfaces by solving heat conduction equations – e.g. buildings, streets, vegetated lands (Offerle et al., 2005; Wang et al., 2012; Yang et al., 2014).

• b.

  Thermal mass scheme: the storage heat is inferred from the changes in thermal mass of all components of the urban system (Kerschgens and Kraus, 1990).

• c.

  Heat flux plates: combined measurements from grass and paved surfaces (Kerschgens and Drauschke, 1986; Kerschgens and Hacker, 1985).

• d.

  Parameterization as a function of $Q^{\ast }$: either as a linear function (Oke et al., 1981), hyperbolic (cotangent, secant) function (Doll et al., 1985), or hysteresis relation (Camuffo and Bernardi, 1982). The last of these is used in the Objective Hysteresis Model (OHM) (Grimmond et al., 1991).

• e.

  Residual: practical difficulties of direct measurement of $\mathrm{\Delta }{Q}_{\mathrm{S}}$ in urban areas, result in the SEB residual (i.e. $Q^{\ast }+Q_F-\left(Q_H,+,Q_E\right)$) frequently being the “preferred” observations (Ao et al., 2016; Ching et al., 1983; Doll et al., 1985; Li et al., 2015; Oke and Cleugh, 1987) (where $Q_F$ is the anthropogenic heat flux).

The focus here is on the OHM approach, which is forced by $Q^{\ast }$ and accounts for the diversity of the surface materials (sub-facets $i)$ in the measurement source area of interest with weightings ($f)$ for their two- or three-dimensional extent (Grimmond et al., 1991): $${\displaystyle }\mathrm{\Delta }{Q}_{\mathrm{S}}=\sum _i{f}_i\left(a_{\mathrm{1},i},Q_i^{\ast },+,a_{\mathrm{2},i},{\displaystyle \frac{\partial Q_i^{\ast }}{\partial t}},+,a_{\mathrm{3},i}\right),$$ where the $a_{\mathrm{1}}$, $a_{\mathrm{2}}$, and $a_{\mathrm{3}}$ coefficients are for individual facets determined by least-square regression between $\mathrm{\Delta }{Q}_{\mathrm{S}}$ and $Q^{\ast }$ using results from observations (e.g. asphalt road (Anandakumar, 1999), wetlands (Souch et al., 1998), forests (Oliphant et al., 2004), or numerical modelling – e.g. urban canyons (Arnfield and Grimmond, 1998) and roofs (Meyn and Oke, 2009). These coefficients capture the net behaviour of a facet type in a typical setting, rather than being required to identify the component materials within a facet (e.g. multiple materials making up a roof, wall, with varying thermal connectivity and individual properties). As such, OHM is one of the less demanding parameterizations, yet does capture a more realistic understanding of the relation between $\mathrm{\Delta }{Q}_{\mathrm{S}}$ by $Q^{\ast }$ compared with other approaches. Despite the shortage of OHM coefficients for the wide range of facet types found in cities, OHM captures the urban $\mathrm{\Delta }{Q}_{\mathrm{S}}$ overall generally well (Grimmond and Oke, 1999; Järvi et al., 2011, 2014; Karsisto et al., 2015; Roth and Oke, 1995).

OHM is a cornerstone in the urban land surface models SUEWS (Surface Urban Energy And Water Balance Scheme; Järvi et al., 2011, 2014; Ward et al., 2016) and LUMPS (Local-Scale Urban Meteorological Parameterization Scheme; Grimmond and Oke, 2002), and plays an essential role in determining the initial energy partitioning at each time step of the models' simulations. Previous modelling studies (Arnfield and Grimmond, 1998; Meyn and Oke, 2009) have led to better understanding of the OHM coefficients. Solution of the one-dimensional advection–diffusion equation of coupled heat and liquid water transport by Gao et al. (2003, 2008) was used to explore the physical relation of OHM coefficients $a_{\mathrm{1}}$ and $a_{\mathrm{2}}$ to the phase lag between $\mathrm{\Delta }{Q}_{\mathrm{S}}$ and $Q^{\ast }$. However, insight into $a_{\mathrm{3}}$ remain unclear (Sun et al., 2013).

In this paper, the solutions of the one-dimensional advection–diffusion equation of coupled heat and liquid water transport (Gao et al., 2003, 2008) are employed with the SEB (Eq. 1) to investigate more fully the three OHM coefficients, the outcomes of which lead to development of the Analytical Objective Hysteresis Model (AnOHM) (Sect. 2). The Monte Carlo-based subset simulation (Au and Beck, 2001) approach is then used to undertake a sensitivity analysis of AnOHM to surface properties and hydrometeorological conditions (Sect. 3). An offline evaluation of AnOHM's performance for five sites with different land covers (Sect. 4) provides evidence that this is an alternative approach to obtain OHM coefficients. Given that this allows applications across a much wider range of environments and meteorological conditions, we conclude that AnOHM has important implications for land surface modelling (urban and non-urban).

Model development

Parameterization of storage heat flux $\mathrm{\Delta }{Q}_{\mathrm{S}}$ for a land surface

For a given land surface (e.g. bare soil), the governing heat conduction–advection equation can be written (Gao et al., 2003, 2010) as $${\displaystyle \frac{\partial T}{\partial t}}=\mathit{\lambda }{\displaystyle \frac{{\partial }^{\mathrm{2}}T}{\mathrm{d}{z}^{\mathrm{2}}}}+W{\displaystyle \frac{\partial T}{\partial z}},$$ where $T$ is the temperature at a reference depth $z$ (positive downward), $t$ is time, $\mathit{\lambda }$ is the thermal diffusivity, and $W=\partial \mathit{\lambda }/\partial z-\left(C_{\mathrm{W}},/,C_g\right)w\mathit{\phi }$ is the soil water flux density (Ren et al., 2000), with $C_{\mathrm{W}}$ the volumetric heat capacity of water, $C_g$ the volumetric heat capacity of soil, $w$ the pore water velocity, and $\mathit{\phi }$ the volumetric soil water content.

The steady-periodic solution of Eq. (3) corresponding to the principal Earth rotation frequency ($\mathit{\omega }=\mathrm{2}\mathit{\pi }/\phantom{\mathrm{2}\mathit{\pi }\mathrm{24}}\mathrm{24}$, in $\mathrm{rad}{\mathrm{h}}^{-\mathrm{1}})$, with boundary condition $${\displaystyle }{T}_{\mathrm{S}}=A_{T_{\mathrm{S}}}\sin \left(\mathit{\omega },t,-,\mathit{\gamma }\right)+{\stackrel{\mathrm{‾}}{T}}_{\mathrm{S}}$$ is given by (Gao et al., 2003, 2010) $${\displaystyle }T\left(z,,,t\right)=A_{T_{\mathrm{S}}}\exp \left(-,z,/,M\right)\sin \left(\mathit{\omega },t,-,z,/,N,-,\mathit{\gamma }\right)+{\stackrel{\mathrm{‾}}{T}}_{\mathrm{S}},$$ where $M=\frac{\mathrm{2}\mathit{\lambda }}{\mathrm{\Delta }+W}$, $N=\frac{\mathrm{\Delta }}{\mathit{\omega }}$, and $\mathrm{\Delta }=\sqrt{\frac{W^{\mathrm{2}}+\sqrt{W^{\mathrm{4}}+\mathrm{16}{\mathit{\lambda }}^{\mathrm{2}}{\mathit{\omega }}^{\mathrm{2}}}}{\mathrm{2}}}$; with ${\stackrel{\mathrm{‾}}{T}}_{\mathrm{S}}$, $A_{T_{\mathrm{S}}}$, and $\mathit{\gamma }$ denoting the daily mean value, amplitude, and initial phase of surface temperature, respectively, which need to be determined by the boundary conditions imposed by the SEB.

From Fourier's law, the soil heat flux is then given by $$\begin{array}{cc}{\displaystyle }G\left(z,,,t\right)\equiv -{\displaystyle \frac{k\partial T}{\partial z}}& {\displaystyle }=k{A}_{T_{\mathrm{S}}}{\displaystyle \frac{\sqrt{M^{\mathrm{2}}+N^{\mathrm{2}}}}{MN}}\exp \left(-,{\displaystyle \frac{z}{M}}\right)\\ & {\displaystyle }& {\displaystyle }\sin \left(\mathit{\omega },t,-,{\displaystyle \frac{z}{N}},-,\mathit{\gamma },+,\mathit{\delta }\right),\end{array}$$ where $\mathit{\delta }=\arctan \left(\frac{M}{N}\right)=\arctan \left[\frac{\mathrm{2}\mathit{\lambda }\mathit{\omega }}{(\mathrm{\Delta }+W)\mathrm{\Delta }}\right]$ and $k$ is the thermal conductivity. In particular, at the surface $z=\mathrm{0}$, the ground heat flux $G_{\mathrm{0}}$ is given by $${\displaystyle }{G}_{\mathrm{0}}\left(t\right)=k{A}_{T_{\mathrm{S}}}{\displaystyle \frac{\sqrt{M^{\mathrm{2}}+N^{\mathrm{2}}}}{MN}}\sin \left(\mathit{\omega },t,-,\mathit{\gamma },+,\mathit{\delta }\right),$$ and a simple written form of $\mathrm{\Delta }{Q}_{\mathrm{S}}$ (if only one surface) can be given as $${\displaystyle }\mathrm{\Delta }{Q}_{\mathrm{S}}=G_{\mathrm{0}}=c_{\mathit{\eta }}\sin \left(\mathit{\omega },t,+,\mathit{\eta }\right),$$ where $\mathit{\eta }=\mathit{\delta }-\mathit{\gamma }$ and $c_{\mathit{\eta }}=k{A}_{T_{\mathrm{S}}}\frac{\sqrt{M^{\mathrm{2}}+N^{\mathrm{2}}}}{MN}$ .

Although the above derivation only considers the land surface made of a single material type, the derived $\mathrm{\Delta }{Q}_{\mathrm{S}}$ (Eq. 8) can be adapted for surfaces made of composite materials or volumes given appropriate bulk/ensemble properties.

Parameterization of net all-wave radiation $Q^{\ast }$ for a land surface

Given the parameterizations of incoming longwave radiation $L_{\downarrow }$, outgoing longwave radiation $L_{\uparrow }$, sensible heat flux $Q_H$, latent heat flux $Q_E$, and storage heat flux $\mathrm{\Delta }{Q}_{\mathrm{S}}$ as follows: $$\begin{array}{}& {\displaystyle }& {\displaystyle }{L}_{\downarrow }={\mathit{\epsilon }}_{\mathrm{a}}\mathit{\sigma }{T}_{\mathrm{a}}^{\mathrm{4}},& {\displaystyle }& {\displaystyle }{L}_{\uparrow }={\underbrace{{\mathit{\epsilon }}_{\mathrm{s}}\mathit{\sigma }{T}_{\mathrm{S}}^{\mathrm{4}}}}_{\mathrm{1}}+{\underbrace{(\mathrm{1}-{\mathit{\epsilon }}_{\mathrm{s}})L_{\downarrow }}}_{\mathrm{2}}& {\displaystyle }& {\displaystyle }{Q}_H=C_{h}U\left(T_{\mathrm{S}},-,T_{\mathrm{a}}\right),& {\displaystyle }& {\displaystyle }{Q}_E=Q_H/\mathit{\beta },& {\displaystyle }& {\displaystyle }\mathrm{\Delta }{Q}_{\mathrm{S}}=G_{\mathrm{0}}.\end{array}$$ The boundary condition imposed by the SEB relation can be rewritten as $$\begin{array}{cc}{\displaystyle }\left(\mathrm{1},-,\mathit{\alpha }\right)K_{\downarrow }& {\displaystyle }+{\mathit{\epsilon }}_{\mathrm{a}}\mathit{\sigma }{T}_{\mathrm{a}}^{\mathrm{4}}-{\mathit{\epsilon }}_{\mathrm{s}}\mathit{\sigma }{T}_{\mathrm{S}}^{\mathrm{4}}\\ & {\displaystyle }& {\displaystyle }=C_{h}U\left(\mathrm{1},+,{\mathit{\beta }}^{-\mathrm{1}}\right)\left(T_{\mathrm{S}},-,T_{\mathrm{a}}\right)+G_{\mathrm{0}},\end{array}$$ where the turbulent fluxes $Q_H$ and $Q_E$ are parameterized as functions of temperature gradient $T_{\mathrm{S}}-T_{\mathrm{a}}$ with albedo $\mathit{\alpha }$, bulk transfer coefficient $C_h$, wind speed $U$, and Bowen ratio ($\mathit{\beta }=Q_H/Q_E)$. Theoretically, the second part of Eq. (10) (i.e. $\left(\mathrm{1},-,{\mathit{\epsilon }}_{\mathrm{s}}\right)L_{\downarrow })$ should be accounted for in the estimation of $L_{\uparrow }$ (Oke, 1987); however, given that it is usually less than $\sim \mathrm{5}$ % of the first part of the equation (see full discussion in Appendix A) for most land covers (Oke, 1987), here it is omitted from consideration and in the development of AnOHM.

By assuming that the incoming solar radiation $K_{\downarrow }$ and air temperature $T_{\mathrm{a}}$ follow sinusoidal forms through a day as function of the mean value for the day (e.g. ${\stackrel{\mathrm{‾}}{K}}_{\downarrow })$ (Sun et al., 2013), $$\begin{array}{}& {\displaystyle }& {\displaystyle }{K}_{\downarrow }=A_K\sin \left(\mathit{\omega },t\right)+{\stackrel{\mathrm{‾}}{K}}_{\downarrow },& {\displaystyle }& {\displaystyle }{T}_{\mathrm{a}}=A_T\sin \left(\mathit{\omega },t,-,\mathit{\tau }\right)+{\stackrel{\mathrm{‾}}{T}}_{\mathrm{a}},\end{array}$$ and introducing the solar radiation scale, $${\displaystyle }{A}_K^{\ast }=\left(\mathrm{1},-,\mathit{\alpha }\right)A_K,$$ and longwave radiation scale (assuming ${\mathit{\epsilon }}_{\mathrm{a}}\approx {\mathit{\epsilon }}_{\mathrm{s}}\approx \mathit{\epsilon }$ as a first-order estimate (as AnOHM is insensitive to this parameter; see Sect. 3.2); see clear sky of $\sim$ 0.85 (Staley and Jurica, 1972) and urban surfaces of $\sim$ 0.95 (Kotthaus et al., 2014): $$\begin{array}{cc}{\displaystyle }{A}_T^{\ast }& {\displaystyle }=\left(\mathrm{4},\mathit{\epsilon },\mathit{\sigma },{\stackrel{\mathrm{‾}}{T}}_{\mathrm{a}}^{\mathrm{3}},+,\left(\mathrm{1},+,{\mathit{\beta }}^{-\mathrm{1}}\right),C_h,U\right)A_T\\ & {\displaystyle }& {\displaystyle }=(f_L+f_T)A_T=f{A}_T,\end{array}$$ where $\mathit{\tau }$ denotes phase differences between $T_{\mathrm{a}}$ and $K_{\downarrow }$, the $f=f_L+f_T$ consists of the longwave energy redistribution factor: $f_L=\mathrm{4}\mathit{\epsilon }\mathit{\sigma }{\stackrel{\mathrm{‾}}{T}}_{\mathrm{a}}^{\mathrm{3}}$ and a turbulent energy redistribution factor: $f_T=\left(\mathrm{1},+,{\mathit{\beta }}^{-\mathrm{1}}\right)C_{h}U$. Linearizing the fourth-order longwave expressions of temperature at mean daily air temperature $\stackrel{\mathrm{‾}}{T_{\mathrm{a}}}$ (Sun et al., 2013), the values of $\stackrel{\mathrm{‾}}{T_{\mathrm{S}}}$ and $A_{T_{\mathrm{S}}}$ are obtained: $$\begin{array}{}& {\displaystyle }& {\displaystyle }{\overline{T}}_{\mathrm{S}}={\displaystyle \frac{\mathrm{1}-\mathit{\alpha }}{f}}{\stackrel{\mathrm{‾}}{K}}_{\downarrow }+{\stackrel{\mathrm{‾}}{T}}_{\mathrm{a}}& {\displaystyle }& {\displaystyle }\begin{array}{l}{A}_{T_{\mathrm{S}}}={\displaystyle \frac{fMN\sin \left(\mathit{\tau }\right)}{N\left(f,M,+,k\right)\sin \left(\mathit{\gamma }\right)-kM\cos \left(\mathit{\gamma }\right)}}{A}_T\\ ={\displaystyle \frac{\mathrm{1}}{\sqrt{M_{\ast }^{\mathrm{2}}+N_{\ast }^{\mathrm{2}}}}}{\displaystyle \frac{\sin \left(\mathit{\tau }\right)}{\sin \left(\mathit{\gamma },-,\mathit{\zeta }\right)}}{A}_T\\ ={\mathit{\chi }}_{\mathit{\gamma }}{A}_T,\end{array}\end{array}$$ where $\mathit{\zeta }=\arctan \left(N_{\ast },/,M_{\ast }\right)$, $\mathit{\gamma }=\mathit{\zeta }+\arctan \left(\frac{\sin \left(\mathit{\tau }\right)}{\cos \left(\mathit{\tau }\right)+A_K^{\ast }/A_T^{\ast }}\right)$, $M_{\ast }=\mathrm{1}+k/\left(fM\right)$, $N_{\ast }=k/\left(fN\right)$ and ${\mathit{\chi }}_{\mathit{\gamma }}=\frac{\mathrm{1}}{\sqrt{M_{\ast }^{\mathrm{2}}+N_{\ast }^{\mathrm{2}}}}\frac{\sin \left(\mathit{\tau }\right)}{\sin \left(\mathit{\gamma },-,\mathit{\zeta }\right)}$.

The net all-wave radiation $Q^{\ast }$ is parameterized as $${\displaystyle }\begin{array}{ll}{Q}^{\ast }& =\left(\mathrm{1},-,\mathit{\alpha }\right)K_{\downarrow }+\mathit{\epsilon }\mathit{\sigma }{T}_{\mathrm{a}}^{\mathrm{4}}-\mathit{\epsilon }\mathit{\sigma }{T}_{\mathrm{S}}^{\mathrm{4}}\\ & =\left(\mathrm{1},-,\mathit{\alpha }\right)\left(A_K,\sin ,\left(\mathit{\omega },t\right),+,{\stackrel{\mathrm{‾}}{K}}_{\downarrow }\right)+f_L\left(T_{\mathrm{a}},-,T_{\mathrm{S}}\right)\\ & =c_{\mathit{\phi }}\sin \left(\mathit{\omega },t,+,\mathit{\phi }\right)+{\displaystyle \frac{f_L}{f}}\left(\mathrm{1},-,\mathit{\alpha }\right){\stackrel{\mathrm{‾}}{K}}_{\downarrow },\end{array}$$ where $\mathit{\phi }=\arctan \left[\frac{\left({\mathit{\chi }}_{\mathit{\gamma }},\sin ,\left(\mathit{\gamma }\right),-,\sin ,\left(\mathit{\tau }\right)\right)}{\left(f,A_K^{\ast }\right)/\left(f_L,A_T^{\ast }\right)-\left({\mathit{\chi }}_{\mathit{\gamma }},\cos ,\left(\mathit{\gamma }\right),-,\cos ,\left(\mathit{\tau }\right)\right)}\right]$ and $${\displaystyle }{c}_{\mathit{\phi }}=\sqrt{{\left[{\displaystyle \frac{{\left(f,A_K^{\ast }\right)}^{\mathrm{2}}}{{\left(f_L,A_T^{\ast }\right)}^{\mathrm{2}}}},-,\left({\mathit{\chi }}_{\mathit{\gamma }},\cos ,\left(\mathit{\gamma }\right),-,\cos ,\left(\mathit{\tau }\right)\right)\right]}^{\mathrm{2}}+{\left[{\mathit{\beta }}_{\mathit{\gamma }},\sin ,\left(\mathit{\gamma }\right),-,\sin ,\left(\mathit{\tau }\right)\right]}^{\mathrm{2}}}.$$

Derivation of AnOHM coefficients

Based on the above parameterizations of $Q^{\ast }$ (Eq. 21) and $\mathrm{\Delta }{Q}_{\mathrm{S}}$ (Eq. 8), together with OHM for a specific surface: $${\displaystyle }\mathrm{\Delta }{Q}_{\mathrm{S}}=a_{\mathrm{1}}{Q}^{\ast }+a_{\mathrm{2}}{\displaystyle \frac{\partial Q^{\ast }}{\partial t}}+a_{\mathrm{3}},$$ the coefficients can be readily derived from the parameterization in Sect. 2.2, as $$\begin{array}{}& {\displaystyle }& {\displaystyle }{a}_{\mathrm{1}}={\displaystyle \frac{c_{\mathit{\eta }}}{c_{\mathit{\phi }}}}\cos \left(\mathit{\eta },-,\mathit{\phi }\right),& {\displaystyle }& {\displaystyle }{a}_{\mathrm{2}}={\displaystyle \frac{c_{\mathit{\eta }}}{\mathit{\omega }{c}_{\mathit{\phi }}}}\sin \left(\mathit{\eta },-,\mathit{\phi }\right),& {\displaystyle }& {\displaystyle }\begin{array}{l}{a}_{\mathrm{3}}=-{\displaystyle \frac{c_{\mathit{\eta }}}{c_{\mathit{\phi }}}}\cos \left(\mathit{\eta },-,\mathit{\phi }\right)\cdot {\displaystyle \frac{f_T}{f}}\left(\mathrm{1},-,\mathit{\alpha }\right){\stackrel{\mathrm{‾}}{K}}_{\downarrow }\\ =-a_{\mathrm{1}}\cdot {\displaystyle \frac{f_T}{f}}\left(\mathrm{1},-,\mathit{\alpha }\right){\stackrel{\mathrm{‾}}{K}}_{\downarrow }.\end{array}\end{array}$$ In the densest parts of cities, the anthropogenic heat ($Q_F)$ often has a large influence on the SEB and it needs to be accounted for (Allen et al., 2011; Chow et al., 2014; Nie et al., 2014; Sailor, 2011). This requires the governing SEB relation (Eq. 14) to be rewritten: $$\begin{array}{cc}{\displaystyle }\left(\mathrm{1},-,\mathit{\alpha }\right)K_{\downarrow }& {\displaystyle }+\mathit{\epsilon }\mathit{\sigma }{T}_{\mathrm{a}}^{\mathrm{4}}-\mathit{\epsilon }\mathit{\sigma }{T}_{\mathrm{S}}^{\mathrm{4}}+Q_F\\ & {\displaystyle }& {\displaystyle }=C_{h}U\left(\mathrm{1},+,{\mathit{\beta }}^{-\mathrm{1}}\right)\left(T_{\mathrm{S}},-,T_{\mathrm{a}}\right)+G_{\mathrm{0}}.\end{array}$$

Assuming $Q_F$ is diurnally invariant (as a first-order estimate – e.g. Best and Grimmond, 2016), the derivation (Sect. 2.2) can be extended to include a first-order estimate of $Q_F$ to obtain $$\begin{array}{cc}{\displaystyle }{a}_{\mathrm{3}F}& {\displaystyle }=-{\displaystyle \frac{c_{\mathit{\eta }}}{c_{\mathit{\phi }}}}\cos \left(\mathit{\eta },-,\mathit{\phi }\right)\cdot {\displaystyle \frac{f_T}{f}}\left(\mathrm{1},-,\mathit{\alpha }\right){\stackrel{\mathrm{‾}}{K}}_{\downarrow }-Q_F\\ & {\displaystyle }& {\displaystyle }=-a_{\mathrm{1}}\cdot {\displaystyle \frac{f_T}{f}}\left(\mathrm{1},-,\mathit{\alpha }\right){\stackrel{\mathrm{‾}}{K}}_{\downarrow }-Q_F,\end{array}$$ where $a_{\mathrm{3}F}$ (subscript “$F$” indicates the inclusion of $Q_F)$. The other two coefficients remain unchanged.

Physical interpretations of AnOHM coefficients

Based on the parameterizations of AnOHM coefficients (Eqs. 23, 24, 25/27), physical interpretations can be more fully described compared with OHM:

• a.

  $a_{\mathrm{1}}$ characterizes the ratio of $\mathrm{\Delta }{Q}_{\mathrm{S}}$ and $Q^{\ast }$ and depends on the energy scales (i.e. $c_{\mathit{\eta }}$ and $c_{\mathit{\phi }})$ and their phase difference (i.e. $\mathit{\eta }-\mathit{\phi })$. The energy scales, representing daily amplitudes of $\mathrm{\Delta }{Q}_{\mathrm{S}}$ and $Q^{\ast }$, determine the overall magnitude, while the phase difference moderates the ratio value.

• b.

  $a_{\mathrm{2}}$ accounts for the temporal changes in $\mathrm{\Delta }{Q}_{\mathrm{S}}$ and $Q^{\ast }$ by including the principal Earth rotation frequency $\mathit{\omega }$, in addition to the same determinants of $a_{\mathrm{1}}$ (i.e. $c_{\mathit{\eta }}$, $c_{\mathit{\phi }}$, and $\mathit{\eta }-\mathit{\phi })$. The complementary sinusoidal functions, with phase difference (i.e. $\sin \left(\mathit{\eta },-,\mathit{\phi }\right)$ and $\cos \left(\mathit{\eta },-,\mathit{\phi }\right)$), in the formulations of $a_{\mathrm{1}}$ and $a_{\mathrm{2}}$ are inversely related with a stronger lag effect from $a_{\mathrm{2}}$, and less contribution to $\mathrm{\Delta }{Q}_{\mathrm{S}}$ by $Q^{\ast }$ (i.e. smaller $a_{\mathrm{1}})$.

• c.

  $a_{\mathrm{3}}$ (or $a_{\mathrm{3}F})$ indicates the baseline $\mathrm{\Delta }{Q}_{\mathrm{S}}$ determined by energy redistribution factors (i.e. $f_T$ and $f)$ and energy inputs (i.e. ${\stackrel{\mathrm{‾}}{K}}_{\downarrow }$, and $Q_F$ if anthropogenic heat is considered) as well as $a_{\mathrm{1}}$. It can be inferred from Eq. (2) that the nocturnal $\mathrm{\Delta }{Q}_{\mathrm{S}}$ is largely determined by $a_{\mathrm{3}}$ when the absolute values and variability of $Q^{\ast }$ are small at night. A larger daytime energy input (i.e. ${\stackrel{\mathrm{‾}}{K}}_{\downarrow }$, and $Q_F$ if anthropogenic heat is considered) suggests more heat released at night.

Sensitivity analysis

Given the complex dependence of AnOHM coefficients on surface properties and meteorological forcing (Sect. 2.3), the impacts of these coefficients are assessed further by a sensitivity analysis.

Subset simulation

To improve the computational efficiency of undertaking Monte Carlo sensitivity analyses, subset simulation is used (Au and Beck, 2001). This is an adaptive stochastic simulation procedure with particular efficiency in analysing the short-tail of a distribution probability (while also adaptable to long-tail scenarios) (Wang et al., 2011).

If the probability that a critical response $Y$ exceeds a threshold $y,P(Y>y$), a range of exceedance regions can be specified and sampled using Markov chains. Initially a direct Monte Carlo method is used to choose possible values for the parameter of interest in the anticipated range with a specified distribution (or probability distribution function, PDF) of the uncertainty. From this (level 0), the first exceedance level probability is determined, $F_{\mathrm{1}}$ at which $P\left(Y,>,y_{\mathrm{1}}\right)$. Then a Markov chain Monte Carlo (MCMC) procedure is used to generate samples of a given conditional probability $p_{\mathrm{0}}$, leading to the exceedance of $y_{\mathrm{1}}$ in the earlier simulations. This procedure is repeated, for exceedance events $F_i$ at which $P\left(Y,>,y_i\right)=p_{\mathrm{0}}^i$, $i=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\dots }$, until simulations reach a target exceedance probability, e.g. associated with rare events or risk analysis. Further details of this subset simulation process are provided in Wang et al. (2011).

Subset simulation efficiently generates conditional samples with Metropolis algorithms (Hastings, 1970; Metropolis et al., 1953). This is the basis of MCMC. To generate samples that successively approach a certain conditional probability, a specific Markov chain is designed with the target PDF as its limiting stationary distribution trend as its length increases. The selection of a distribution is key as this controls the next sample generated from the current one. Ideally, the distribution selection would be automatic but this has an efficiency cost relative to the robustness benefit. For the surface parameters (Table 1a) and hydrometeorological forcing (Table 1b) analyses a normal distribution PDF is used (Au and Beck, 2003; Au et al., 2007), with three conditional levels ($N_{\mathrm{level}}=\mathrm{3}$) and a conditional probability of $p_{\mathrm{0}}=\mathrm{0.1}$ – i.e. at each level the highest 10 % of the outputs are considered to exceed the intermediate threshold. As such, the three-level simulation can effectively capture a rare event with the target exceedance probability of ${\mathrm{10}}^{-\mathrm{4}}$ (i.e. the probability of occurrence is less than 1 in 10 000) and generate appropriate samples of different conditional probabilities.

Range of values used as basis for the sensitivity analysis: (a) surface parameters and (b) hydrometeorological variables. All are assumed to have normal PDF. Values of surface parameters are based on values reported in Stull (1988).

Parameter/variable		Unit	Min	Max	Mean	Standard deviation
(a) Surface
Thermal conductivity	$k$	W m${}^{-\mathrm{1}}$ K${}^{-\mathrm{1}}$	0	3	1.2	0.1
Bulk material heat capacity	$C_p$	MJ m${}^{-\mathrm{3}}$ K${}^{-\mathrm{1}}$	0	4	2.0	0.04
Albedo	$\mathit{\alpha }$	–	0	1	0.27	0.07
Emissivity	$\mathit{\epsilon }$	–	0.8	1.0	0.93	0.025
Midday* mean Bowen ratio (inverse)	${\mathit{\beta }}^{-\mathrm{1}}$	–	0	20	0.05	0.05
Bulk transfer coefficient	$C_h$	J m${}^{-\mathrm{3}}$ K${}^{-\mathrm{1}}$	0	8	4	0.5
(b) Hydrometeorological
Amplitude or range of the daily incoming shortwave radiation	$A_K$	W m${}^{-\mathrm{2}}$	0	1200	800	200
Mean daytime incoming shortwave radiation	${\stackrel{\mathrm{‾}}{K}}_{\downarrow }$	W m${}^{-\mathrm{2}}$	0	500	200	50
Amplitude or range of the daily air temperature	$A_T$	${}^{\circ }$C	0	15	8	2
Mean daily air temperature	${\stackrel{\mathrm{‾}}{T}}_{\mathrm{a}}$	${}^{\circ }$C	0	40	30	7.5
Phase lag between radiation and air temperature	$\mathit{\tau }$	rad	0	$\mathit{\pi }/\mathrm{2}$	$\mathit{\pi }/\mathrm{4}$	$\mathit{\pi }/\mathrm{10}$
Mean daytime wind speed	$U$	m s${}^{-\mathrm{1}}$	0	4	2	0.5
Mean daily water flux density	$W$	10${}^{-\mathrm{7}}$ m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$ m${}^{-\mathrm{2}}$	0	100	10	5

* midday period: 1000–1400 local standard time.

Characteristics of the flux towers at the study sites.

Site	UK-Ldn	US-Wlr	CA-NS5	US-SRM	US-SO4
Location	51.50${}^{\circ }$ N, 0.12${}^{\circ }$ W	37.52${}^{\circ }$ N, 96.86${}^{\circ }$ W	55.86${}^{\circ }$ N, 98.49${}^{\circ }$ W	31.82${}^{\circ }$ N, 110.87${}^{\circ }$ W	33.38${}^{\circ }$ N, 116.64${}^{\circ }$ W
Land cover classification	Urban/built-up	Grassland	Evergreen needleleaf forest	Woody savannas	Closed shrublands
Land cover code	URB	GRA	ENF	WSA	CSH
Study year	2011	2003	2004	2004	2005
Reference	Kotthaus and Grimmond (2014a, b)	Klazura et al. (2006), Coulter et al. (2006)	Goulden et al. (2006)	Scott et al. (2009)	Luo et al. (2007)

The metric $S$ (in %), used to indicate the sensitivity of the model output $Y$ to a specific uncertainty parameter $X$ (Wang et al., 2011), is $${\displaystyle }S=\left[{\displaystyle \frac{\mathrm{1}}{N_{\mathrm{level}}}},\sum _{i=\mathrm{1}}^{N_{\mathrm{level}}},{\displaystyle \frac{E\left[X,\mathrm{|},Y,>,y_i\right]-E\left[X\right]}{E\left[X\right]}}\right]\times \mathrm{100},$$ where $i=\mathrm{1},\mathrm{2},\mathrm{\dots },N_{\mathrm{level}}$ is the index of conditional sampling level, $E\left[X\right]$ is the expectation that the unconditional distribution of a specific uncertainty parameter $X$, while $E\left[X,\mathrm{|},Y,>,y_i\right]$ is the expectation of $X$ at conditional level $i$. A positive (negative) $S$ indicates an increase will lead to increase (decrease) in simulated value. Hence the sign of $S$ indicates the impact of a change in parameter uncertainty. The absolute magnitude of $S$ indicates the sensitivity.

This assessment does not consider if the simulated values have low probability. Later analyses (Sect. 4) consider the simulation results relative to observed fluxes.

Impacts of surface properties

Following the sensitivity analysis of AnOHM coefficients to the surface properties, the distributions of conditional samples for thermal conductivity $k$, bulk heat capacity $C_p$, and emissivity $\mathit{\epsilon }$ are similar to the original proposal distributions (Fig. 1), implying weak dependence of $a_{\mathrm{1}}$, $a_{\mathrm{2}}$, and $a_{\mathrm{3}}$ on these properties. However, for albedo ($\mathit{\alpha })$ both $a_{\mathrm{2}}$ and $a_{\mathrm{3}}$ are sensitive, but $a_{\mathrm{1}}$is not; changes in inverse Bowen ratio (${\mathit{\beta }}^{-\mathrm{1}})$ impact all three coefficients; and the bulk transfer coefficient $C_h$ impacts $a_{\mathrm{1}}$ and $a_{\mathrm{2}}$, but has little effect on $a_{\mathrm{3}}$.

Histograms of conditional samples at different conditional levels for surface property parameters (rows from top: thermal conductivity $k$ in W m${}^{-\mathrm{1}}$ K${}^{-\mathrm{1}}$, heat capacity $C_p$ in MJ m${}^{-\mathrm{3}}$ K${}^{-\mathrm{1}}$, albedo $\mathit{\alpha }$, emissivity $\mathit{\epsilon }$, inverse Bowen ratio ${\mathit{\beta }}^{-\mathrm{1}}$, and bulk transfer coefficient $C_h$ in J m${}^{-\mathrm{3}}$ K${}^{-\mathrm{1}})$ with AnOHM coefficients as the model output (columns from left: $a_{\mathrm{1}}$, $a_{\mathrm{2}}$ and $a_{\mathrm{3}})$. Each subplot $x$ axis is the parameter value and $y$ axis is the PDF value. The original proposal distribution (dashed line) and simulation levels (different colours) are shown.

[Figure omitted. See PDF]

Using $S$ (Eq. 28) to quantify this, it is found that the surface properties ($k$, $C_p$, and $\mathit{\epsilon })$ have less sensitivity, with less skewed conditional samples between levels, so $S$ values close to 0 (Fig. 2). The $S$ of $k$ is the largest of the three. From the $S$ results for the $\mathit{\alpha }$ sensitivity analysis (Fig. 2), it is apparent that an increase in $\mathit{\alpha }$ will increase $a_{\mathrm{1}}$ while decreasing $a_{\mathrm{2}}$ and $a_{\mathrm{3}}$, whereas the reverse occurs for ${\mathit{\beta }}^{-\mathrm{1}}$ and $C_h$ (i.e. their decreases leads to larger $a_{\mathrm{2}}$ and $a_{\mathrm{3}}$ values but smaller $a_{\mathrm{1}})$.

Relative variation in sensitivity ($S$, %, Eq. 28) to surface parameters. See Fig. 1 for further details.

[Figure omitted. See PDF]

From this, the links between the key surface parameters and the storage heat flux can be considered. With an increase in $\mathit{\alpha }$, there is reduced solar energy in the SEB. This reduces the temporal change in $\mathrm{\Delta }{Q}_{\mathrm{S}}$ (smaller $a_{\mathrm{2}})$ and decreases the baseline value of $\mathrm{\Delta }{Q}_{\mathrm{S}}$ (smaller $a_{\mathrm{3}})$; larger ${\mathit{\beta }}^{-\mathrm{1}}$ indicates that more available energy is dissipated by $Q_E$ than by $Q_H$, leading to decreased $T_{\mathrm{s}}$ and $\mathrm{\Delta }{Q}_{\mathrm{S}}$ (smaller $a_{\mathrm{1}})$; a smaller portion of $Q^{\ast }$ will be dissipated by $\mathrm{\Delta }{Q}_{\mathrm{S}}$ (smaller $a_{\mathrm{1}})$ as the increased $C_h$ can facilitate the turbulent convection and thus increase the total turbulent fluxes.

Impacts of hydrometeorological conditions

Similarly, the sensitivity of AnOHM to hydrometeorological variables is explored (Fig. 3). The air temperature (range, mean) and water flux density related variables (i.e. $A_T$, ${\stackrel{\mathrm{‾}}{T}}_{\mathrm{a}}$, and $W)$ have minimal influence on the skewness of the conditional samples. In contrast, the incoming shortwave (solar) radiation (range, mean) and wind-related variables (i.e. $A_K$, ${\stackrel{\mathrm{‾}}{K}}_{\downarrow }$, and $U)$ and the phase lag $\mathit{\tau }$ between $K_{\downarrow }$ and $T_{\mathrm{a}}$ have large impacts. In terms of the greatest impact on the coefficients ($a_{\mathrm{1}}$, $a_{\mathrm{2}}$, and $a_{\mathrm{3}})$: $A_K$ and $U$ influences $a_{\mathrm{1}}$, $\mathit{\tau }$ impacts $a_{\mathrm{2}}$, and $a_{\mathrm{3}}$ responds more to $A_K$ and ${\stackrel{\mathrm{‾}}{K}}_{\downarrow }$ than the other variables.

Histograms of conditional samples at different conditional levels for ambient forcing parameters (rows from top: incoming solar radiation amplitude $A_K$ in W m${}^{-\mathrm{2}}$ and its daytime mean ${\stackrel{\mathrm{‾}}{K}}_{\downarrow }$ in W m${}^{-\mathrm{2}}$, air temperature amplitude $A_T$ in ${}^{\circ }$C and its daily mean ${\stackrel{\mathrm{‾}}{T}}_{\mathrm{a}}$ in ${}^{\circ }$C, the phase lag $\mathit{\tau }$ in rad between $K_{\downarrow }$ and $T_{\mathrm{a}}$, wind speed $U$ in m s${}^{-\mathrm{1}}$, and water flux density $W$ in m s${}^{-\mathrm{1}})$ with AnOHM coefficients as the model output (columns from left: $a_{\mathrm{1}}$, $a_{\mathrm{2}}$ and $a_{\mathrm{3}})$. As Fig. 1.

[Figure omitted. See PDF]

Variables that strongly modulate the interactions between $\mathrm{\Delta }{Q}_{\mathrm{S}}$ and $Q^{\ast }$ can be informed by the $S$ results (Fig. 4). For instance, a greater range in $K_{\downarrow }$ (i.e. larger $A_K)$ will occur with larger energy input from solar radiation, leading to stronger heating of the near-surface atmosphere and a smaller portion to $\mathrm{\Delta }{Q}_{\mathrm{S}}$ (smaller $a_{\mathrm{1}})$ but higher baseline $\mathrm{\Delta }{Q}_{\mathrm{S}}$ (larger $a_{\mathrm{3}})$. This is consistent with a reduction in $\overline{K_{\downarrow }}$ having a decrease in $a_{\mathrm{3}}.$ The temporal change in $\mathrm{\Delta }{Q}_{\mathrm{S}}$ is highly correlated with the change in $\mathit{\tau }$, an increase in which implies a slower response of the surface to solar radiation and an overall decrease in $\mathrm{\Delta }{Q}_{\mathrm{S}}$ (smaller $a_{\mathrm{1}}$, $a_{\mathrm{2}}$, and $a_{\mathrm{3}})$. The greater sensitivity to $\mathit{\tau }$ of $a_{\mathrm{2}}$ is a key part of the original hysteresis nature of the heating/cooling of a surface. The sensitivity responses of $a_{\mathrm{1}}$, $a_{\mathrm{2}}$, and $a_{\mathrm{3}}$ to $U$ are very consistent with those to $C_h$, suggesting the similar pathway that turbulent fluxes (i.e. $Q_H$ and $Q_E)$ modulate $\mathrm{\Delta }{Q}_{\mathrm{S}}$. As $W$ mostly influences the heat conduction–diffusion in the underlying surface as thermal properties (i.e. $C_p$ and $k)$, less dependence is observed on it. This is similar with $C_p$ and $k$.

Relative variation in sensitivity ($S$, %, Eq. 28) to forcing parameters. See Fig. 3 for further details.

[Figure omitted. See PDF]

Model evaluation

In this section, the actual ability of AnOHM to determine the storage heat flux relative to observations is evaluated using 30 min observations from five sites of different land use/covers (Table 2). The measurements include turbulent sensible and latent fluxes, along with incoming and outgoing shortwave and longwave radiation and basic meteorological variables (see Kotthaus and Grimmond, 2014a, b; Klazura et al., 2006; Coulter et al., 2006; Goulden et al., 2006; Scott et al., 2009; Luo et al., 2007, for details). Anthropogenic heat flux $Q_F$ at the urban site (i.e. UK-Ldn) is estimated using the GreaterQF model (Iamarino et al., 2011); the heat storage flux $\mathrm{\Delta }{Q}_{\mathrm{S}}$ is thus estimated as the modified residual of urban energy balance as $\mathrm{\Delta }{Q}_{\mathrm{S}}=Q^{\ast }+\mathrm{0.75}{Q}_F-\mathrm{1.2}\left(Q_H,+,Q_E\right)$ (Kotthaus and Grimmond, 2014a, b), which is then used in this evaluation. A similar approach for estimating $\mathrm{\Delta }{Q}_{\mathrm{S}}$ (i.e. residual of surface energy balance, $\mathrm{\Delta }{Q}_{\mathrm{S}}=Q^{\ast }+Q_F-\left(Q_H,+,Q_E\right)$) is applied at the other (non-urban) sites but with $Q_F=\mathrm{0}$.

Surface properties used in AnOHM simulation for the study sites based on calibration. The values of $\mathit{\alpha }$ and $\mathit{\beta }$ are monthly climatology from January to December and are used when observations are not available (see Table 1 for notation definition).

Parameter	Unit	Site
		UK-Ldn	US-Wlr	CA-NS5	US-SRM	US-SO4
$k$	W m${}^{-\mathrm{1}}$ K${}^{-\mathrm{1}}$	2.8	0.43	0.51	0.41	0.56
$C_p$	MJ m${}^{-\mathrm{3}}$ K${}^{-\mathrm{1}}$	2.4	0.31	0.36	0.56	0.27
$\mathit{\alpha }$	–	0.24, 0.24, 0.22, 0.20, 0.14, 0.13, 0.12, 0.14, 0.18, 0.24, 0.24, 0.18	0.29, 0.29, 0.17, 0.18, 0.18, 0.12, 0.11, 0.10, 0.19, 0.13, 0.24, 0.35	0.30, 0.29, 0.22, 0.15, 0.10, 0.10, 0.10, 0.11, 0.22, 0.24, 0.28, 0.30	0.13, 0.17, 0.16, 0.14, 0.13, 0.12, 0.13, 0.15, 0.14, 0.19, 0.13, 0.18	0.22, 0.11, 0.11, 0.10, 0.11, 0.10, 0.10, 0.10, 0.11, 0.10, 0.17, 0.24
$\mathit{\epsilon }$	–	0.92	0.93	0.95	0.95	0.92
$\mathit{\beta }$	–	6.1, 5.1, 8.3, 7.9, 5.4, 3.9, 5.3, 4.2, 5.2, 4.3, 4.8, 3.2	2.9, 0.8, 7.6, 2.7, 0.3, 0.3, 0.3, 0.8, 0.5, 0.7, 2.3, 2.3	6.1, 6.0, 8.7, 8.0, 1.9, 1.6, 0.7, 0.7, 1.3, 1.4, 3.1, 8.0	1.9, 5.5, 3.3, 2.0, 10.1, 9.7, 2.0, 0.9, 3.0, 4.3, 10.0, 3.3	1.5, 1.4, 1.9, 3.0, 1.4, 1.4, 2.1, 1.2, 2.8, 1.9, 2.1, 4.1
$C_h$	J m${}^{-\mathrm{3}}$ K${}^{-\mathrm{1}}$	4.3	1.9	5.1	3.6	3.9

AnOHM is first calibrated with observations under sunny conditions, when the assumptions of AnOHM are best satisfied (i.e. diurnal cycles of $K_{\downarrow }$ and $T_{\mathrm{a}}$ follow sinusoidal forms), to obtain surface properties required by AnOHM (Table 3). As the Bowen ratio $\mathit{\beta }$ varies daily and monthly (Kotthaus and Grimmond, 2014a, b), $\mathit{\beta }$ is either determined as the daily value if available, or based on the observation-based monthly climatology (Table 3). The seasonality in albedo $\mathit{\alpha }$ is accounted for also by using its monthly climatology (Table 3). AnOHM is driven by atmospheric forcing (i.e. $K_{\downarrow }$, $T_{\mathrm{a}}$, and $U)$ and/or their derived scales ($A_K$, ${\stackrel{\mathrm{‾}}{K}}_{\downarrow }$, $A_T$, ${\stackrel{\mathrm{‾}}{T}}_{\mathrm{a}}$, and $\mathit{\tau }$) to generate the OHM coefficients (i.e. $a_{\mathrm{1}}$, $a_{\mathrm{2}}$, and $a_{\mathrm{3}}$, see Fig. 5), from which the net heat storage flux $\mathrm{\Delta }{Q}_{\mathrm{S}}$ can be predicted (Fig. 6) using the observed $Q^{\ast }$ with Eq. (2).

Intra-annual variations of OHM coefficients: (a) $a_{\mathrm{1}}$, (b) $a_{\mathrm{2}}$, and (c) $a_{\mathrm{3}}$. LOESS fits (solid lines) through the daily values predicted by AnOHM and daily values (squares) measured at an asphalt road site (Anandakumar, 1999) are shown. The LOESS (Cleveland and Devlin, 2012) fitting is a locally weighted polynomial regression approach.

[Figure omitted. See PDF]

Monthly median (line) diurnal cycles and interquartile range (shaded) values of $\mathrm{\Delta }{Q}_{\mathrm{S}}$ for AnOHM predictions (blue), OHM predictions (orange) and observations (green) at (a) UK-Ldn (URB), (b) US-Wlr (GRA), (c) CA-NS5 (ENF), (d) US-SRM (WSA), and (e) US-SO4 (CSH) (see Table 2 for site information). Statistics include average bias and RMSE (W m${}^{-\mathrm{2}})$. The OHM coefficients $a_{\mathrm{1}}$, $a_{\mathrm{2}}$, and $a_{\mathrm{3}}$ used for different land covers are: 0.553, 0.303, and $-$37.6 at the urban site (UK-Ldn) (Ward et al., 2016), 0.32, 0.54, and $-$27.4 at the grass-covered sites (US-Wlr and US-SRM) (Grimmond and Oke, 1999), and 0.11, 0.11, and $-$12.3 at the forest-covered sites (CA-NS5 and US-SO4) (Grimmond and Oke, 1999).

[Figure omitted. See PDF]

To examine the seasonality of the OHM coefficients, rather than the daily variations in hydrometeorological forcing, LOESS (LOcally wEighted Scatter-plot Smoother; Cleveland and Devlin, 1988) curves are obtained to filter out day-to-day variations in the OHM coefficients (see Appendix B for a direct comparison of these coefficients by different modelling and observational regression approaches). Intra-annual variations are found in all the three OHM coefficients (Fig. 5), indicating the strong impact of seasonality of meteorological conditions. These controls, as indicated by Eqs. (23)–(25/27), are complex and will vary with local conditions. For instance, comparison of OHM coefficients between the AnOHM predictions (LOESS fitted solid lines in Fig. 5) and observations at an asphalt road site in Alland, Austria, reported in Anandakumar (1999) (empty squares in Fig. 5) demonstrates differences in $a_{\mathrm{1}}$ (Fig. 5a) and $a_{\mathrm{2}}$ (Fig. 5b) but general similarity in $a_{\mathrm{3}}$ (Fig. 5c). Compared to $a_{\mathrm{1}}$ and $a_{\mathrm{2}}$, it is noteworthy that, in addition to the $S$ results (see Fig. 4) given the more explicit mechanism by which the atmospheric conditions moderate $a_{\mathrm{3}}$ (see Eqs. 25 and 27), such seasonality in $a_{\mathrm{3}}$ is predicted by AnOHM, and evident in the observations (Fig. 5c, also Ward et al., 2013). Larger ${\stackrel{\mathrm{‾}}{K}}_{\downarrow }$ in warm seasons (May–September) will lead to smaller $a_{\mathrm{3}}$ (Eqs. 25, 27) and vice versa.

The AnOHM simulated and observed $\mathrm{\Delta }{Q}_{\mathrm{S}}$ agree well at the five different land cover sites, with RMSE values of $\sim$ 30 W m${}^{-\mathrm{2}}$. For comparison purposes, it is noted that the urban land surface model comparison (Best and Grimmond, 2015; Grimmond et al., 2011) found $\mathrm{\Delta }{Q}_{\mathrm{S}}$ to be the most poorly represented among all the SEB components with the best RMSE values of 53 W m${}^{-\mathrm{2}}$ (Lipson et al., 2017). Although the much smaller $\mathrm{\Delta }{Q}_{\mathrm{S}}$ RMSE obtained by AnOHM uses a prescribed Bowen ratio in the offline evaluation, such improvement indicates the ability of AnOHM to simulate a more consistent $\mathrm{\Delta }{Q}_{\mathrm{S}}$ with observations. Compared with OHM predictions (orange lines in Fig. 6), AnOHM (blue lines in Fig. 6) better reproduces the seasonality in $\mathrm{\Delta }{Q}_{\mathrm{S}}$ but gives larger bias at two sites with natural land covers (i.e. US-SRM and US-SO4). This can be attributed to the overestimates of nocturnal $\mathrm{\Delta }{Q}_{\mathrm{S}}$ by AnOHM. Overall, the evaluation demonstrates good performance of AnOHM in predicting the long-term $\mathrm{\Delta }{Q}_{\mathrm{S}}$ with clear seasonality reproduced across a wide range of surface types.

Discussion and concluding remarks

In this study, the Analytical Objective Hysteresis Model (AnOHM) is developed to obtain OHM coefficients across a wide range of surface and meteorological conditions and to improve physical understanding of the interactions between $\mathrm{\Delta }{Q}_{\mathrm{S}}$ and $Q^{\ast }$. The sensitivity of AnOHM to surface properties and hydrometeorological conditions is analysed through Monte Carlo-based subset simulations (Au and Beck, 2001). The results highlight the importance of the albedo, the Bowen ratio, and the bulk transfer coefficient, and the importance of solar radiation and wind speed in regulating the heat storage. The importance of albedo in modulating the heat storage was also found by Wang et al. (2011), who also used the same subset simulation approach with the single-layer urban canopy model (SLUCM; for details see Kusaka et al., 2001). This demonstrates the consistency in heat storage modelling between AnOHM and SLUCM. From the sensitivity results, variations in OHM coefficients of a similar size may arise from either surface property parameters or hydrometeorological forcing that are associated with the same physical processes (see bulk transfer coefficient $C_h$ in Fig. 2 and wind speed $U$ in Fig. 4). This supports the ability of AnOHM in representing physical processes. An offline evaluation of AnOHM using flux observations from five sites with different land covers demonstrates its ability to predict the intra-annual dynamics of OHM coefficients and shows good agreement between simulated and observed storage heat fluxes. In particular, the seasonality in the OHM coefficient $a_{\mathrm{3}}$ observed in a previous study (Anandakumar, 1999) is well predicted by AnOHM.

The limitations of AnOHM are important to consider. First, given the assumption that the incoming solar radiation $K_{\downarrow }$ and air temperature $T_{\mathrm{a}}$ diurnal cycles are sinusoidal, optimal performance of AnOHM occurs under clear-sky conditions. The current parameterizations of $K_{\downarrow }$ and $T_{\mathrm{a}}$ within AnOHM only consider the harmonics of principal frequencies for formulation simplicity. More frequencies may potentially resolve more realistic diurnal variations in $K_{\downarrow }$ and $T_{\mathrm{a}}$. As the reflected part of $L_{\downarrow }$ (i.e. $\left(\mathrm{1},-,{\mathit{\epsilon }}_{\mathrm{s}}\right)L_{\downarrow }$) is assumed negligible, and similar emissivity values are assumed for sky and land surface (i.e. ${\mathit{\epsilon }}_{\mathrm{s}}\approx {\mathit{\epsilon }}_{\mathrm{a}}\approx \mathit{\epsilon })$, the outgoing longwave radiation is underestimated. These simplifications greatly facilitate the AnOHM formulation without qualitatively changing the final results as the sensitivity analyses (see the minimal $S$ values for $\mathit{\epsilon }$ in Fig. 2) demonstrate. The inclusion of water flux density $W$ equips AnOHM with an ability to investigate the hydrological impacts of the underlying surface on land–atmosphere interactions. However, estimation of $W$ remains challenging (Wang, 2014) and the resulting uncertainty in the final results warrants caution in conducting simulations over land covers with strong soil moisture dynamics (e.g. grassland with high soil moisture under clear-sky condition). Despite these limitations, AnOHM does permit improved modelling of the surface energy balance through its physically based parameterization scheme for storage heat flux $\mathrm{\Delta }{Q}_{\mathrm{S}}$. Compared to OHM, AnOHM has the benefit of allowing $\mathrm{\Delta }{Q}_{\mathrm{S}}$ to be simulated for land covers for which coefficients are not available and to allow for seasonal variability to be accounted for. As AnOHM shares similar hydrometeorological forcing inputs (i.e. $K_{\downarrow }$, $T_{\mathrm{a}}$ and $U)$ to other land surface models (LSMs), it can potentially be used within in LSMs to estimate $\mathrm{\Delta }{Q}_{\mathrm{S}}$, or if turbulent fluxes are included to be a complete LSM. The overall improvements from adopting AnOHM in modelling land surface processes will be presented in forthcoming work in the SUEWS–AnOHM framework.

The Fortran source code for AnOHM can be obtained from the corresponding authors upon request.

Rationale for a simplified formulation of outgoing longwave radiation

In the formulation of outgoing longwave radiation $L_{\uparrow }$, a simplified form (i.e. ${\mathit{\epsilon }}_{\mathrm{s}}\mathit{\sigma }{T}_{\mathrm{s}}^{\mathrm{4}})$ is used for AnOHM by ignoring part 2 of Eq. (10) (i.e. $\left(\mathrm{1},-,{\mathit{\epsilon }}_{\mathrm{s}}\right)L_{\downarrow })$. The rationale for such simplification is that given ${\mathit{\epsilon }}_{\mathrm{s}}$ is usually larger than 0.9, $\left(\mathrm{1},-,{\mathit{\epsilon }}_{\mathrm{s}}\right)L_{\downarrow }$ contributes a relatively small portion to the total longwave component (Oke, 1987) and omission of this part is well accepted in the parameterization of outgoing longwave radiation for land surface modelling across various land covers (Bateni and Entekhabi, 2012; Lee et al., 2011; Stensrud, 2007).

Using the parameterization of incoming longwave radiation in the AnOHM framework (i.e. $L_{\downarrow }={\mathit{\epsilon }}_{\mathrm{a}}\mathit{\sigma }{T}_{\mathrm{a}}^{\mathrm{4}}\approx {\mathit{\epsilon }}_{\mathrm{s}}\mathit{\sigma }{T}_{\mathrm{a}}^{\mathrm{4}})$, we conduct a sensitivity analysis of the ratio between the ignored part (i.e. $\left(\mathrm{1},-,{\mathit{\epsilon }}_{\mathrm{s}}\right)L_{\downarrow })$ and total outgoing longwave radiation (i.e. ${\mathit{\epsilon }}_{\mathrm{s}}\mathit{\sigma }{T}_{\mathrm{s}}^{\mathrm{4}}+\left(\mathrm{1},-,{\mathit{\epsilon }}_{\mathrm{s}}\right)L_{\downarrow })$ at a constant air temperature of 20 ${}^{\circ }$C and find this ratio is generally less than 5 % given ${\mathit{\epsilon }}_{\mathrm{s}}$ ranges between 0.90 and 0.99 (Fig. A1).

Moreover, if $\left(\mathrm{1},-,{\mathit{\epsilon }}_{\mathrm{s}}\right)L_{\downarrow }$ is included in the net longwave radiation, the induced effect can be incorporated into a modified sky emissivity ${\mathit{\epsilon }}_{\mathrm{a}}^{\prime }={\mathit{\epsilon }}_{\mathrm{s}}{\mathit{\epsilon }}_{\mathrm{a}}$ as follows: $${\displaystyle }\begin{array}{ll}{L}_{\mathrm{net}}& =L_{\downarrow }-L_{\uparrow }\\ & =L_{\downarrow }-\left({\mathit{\epsilon }}_{\mathrm{s}},\mathit{\sigma },T_{\mathrm{s}}^{\mathrm{4}},+,\left(\mathrm{1},-,{\mathit{\epsilon }}_{\mathrm{s}}\right),L_{\downarrow }\right)\\ & ={\mathit{\epsilon }}_{\mathrm{s}}{L}_{\downarrow }-{\mathit{\epsilon }}_{\mathrm{s}}\mathit{\sigma }{T}_{\mathrm{s}}^{\mathrm{4}}\\ & ={\mathit{\epsilon }}_{\mathrm{s}}{\mathit{\epsilon }}_{\mathrm{a}}\mathit{\sigma }{T}_{\mathrm{a}}^{\mathrm{4}}-{\mathit{\epsilon }}_{\mathrm{s}}\mathit{\sigma }{T}_{\mathrm{s}}^{\mathrm{4}}\\ & ={\mathit{\epsilon }}_{\mathrm{a}}^{\prime }\mathit{\sigma }{T}_{\mathrm{a}}^{\mathrm{4}}-{\mathit{\epsilon }}_{\mathrm{s}}\mathit{\sigma }{T}_{\mathrm{s}}^{\mathrm{4}}\end{array}.$$ Then by assuming $\mathit{\epsilon }\approx {\mathit{\epsilon }}_{\mathrm{a}}^{\prime }\approx {\mathit{\epsilon }}_{\mathrm{s}}$, the derivation following Eq. (18) still holds. The sensitivity analysis suggests that the derived coefficients are insensitive to $\mathit{\epsilon }$ (see $S$ for $\mathit{\epsilon }$ in Fig. 2).

As such, we deem the omission of $\left(\mathrm{1},-,{\mathit{\epsilon }}_{\mathrm{s}}\right)L_{\downarrow }$ will not qualitatively change the results of this work.

Ratio between the second part of Eq. (10) (i.e. $\left(\mathrm{1},-,{\mathit{\epsilon }}_{\mathrm{s}}\right)L_{\downarrow })$ and total outgoing longwave radiation (i.e. ${\mathit{\epsilon }}_{\mathrm{s}}\mathit{\sigma }{T}_{\mathrm{s}}^{\mathrm{4}}+\left(\mathrm{1},-,{\mathit{\epsilon }}_{\mathrm{s}}\right)L_{\downarrow })$ at a constant air temperature of 20 ${}^{\circ }$C.

[Figure omitted. See PDF]

Comparison in OHM coefficients between different modelling approaches and observation regression

The comparison in OHM coefficients by different modelling and observational regression approaches (Fig. B1) indicate AnOHM generally follows the results by observation regression, whereas the typical coefficient values adopted by OHM do not.

Comparison of OHM coefficients (left, central and right columns for $a_{\mathrm{1}}$, $a_{\mathrm{2}}$ and $a_{\mathrm{3}}$, respectively) between different modelling approaches and observation regression at five sites: UK-Ldn (a, b, c), US-Wlr (d, e, f), CA-NS5 (g, h, i), US-SRM (j, k, l) and US-SO4 (m, n, o). The blue dots denote the paired values between AnOHM and observation regression. The orange lines represent the reference value used in OHM simulations for land covers of grass and tree (Grimmond and Oke, 1999), whereas the green lines show median values derived from results by observation regression at corresponding sites.

[Figure omitted. See PDF]

The authors declare that they have no conflict of interest.

Acknowledgements

Funding is acknowledged from Met Office/Newton Fund CSSP- China (SG), National Science Foundation of China (51679119, TS), and U.S. National Science Foundation (CBET-1435881, ZHW). The authors thank Ivan Au (University of Liverpool) for providing the Subset Simulation package. The authors acknowledge the large number of people who have contributed to the data collection, the agencies that have provided sites and the agencies that funded the research at the individual sites. The US Department of Energy's Office of Science funded AmeriFlux data (ameriflux-data.lbl.gov) are from US-Wlr (PIs: David Cook and Richard L. Coulter), CA-NS5 (PI: Mike Goulden), US-SRM (PI: Russell Scott) and US-SO4 (PI: Walt Oechel, funded by San Diego State University and SDSU Field Stations Program). The London data are supported by NERC ClearfLo (NE/H003231/1), NERC/Belmont TRUC (NE/L008971/1), EUf7 BRIDGE (211345), H2020 UrbanFluxes (637519), King's College London and University of Reading. In particular, the authors thank Simone Kotthaus (University of Reading) for her detailed preparation of the UK-Ldn site data. For access to the UK-Ldn site data, please contact [email protected] by: Chiel van Heerwaarden Reviewed by: three anonymous referees