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

. The net storage heat ﬂux ( (cid:49)Q S ) is important in the urban surface energy balance (SEB) but its determination remains a signiﬁcant challenge. The hysteresis pattern of the diurnal relation between the (cid:49)Q S and net all-wave radiation ( Q ∗ ) has been captured in the Objective Hysteresis Model (OHM) parameterization of (cid:49)Q S . Although successfully used in urban areas, the limited availability of coefﬁ-cients for OHM hampers its application. To facilitate use, and enhance physical interpretations of the OHM coefﬁcients, 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 hy-drometeorological forcing is presented using a stochastic approach (subset simulation). The sensitivity test suggests that the albedo, Bowen ratio and bulk transfer coefﬁcient, solar radiation and wind speed are most critical. AnOHM, driven by local meteorological conditions at ﬁve sites with different land use, is shown to simulate the (cid:49)Q S ﬂux well (RMSE values of ∼ 30 W m − 2 ) . The intra-annual dynamics of OHM coefﬁcients are explored. AnOHM offers signiﬁcant 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 landatmospheric interactions are strongly modulated by the partitioning of solar energy at the land surface (Chen and Dudhia, 2001;McCumber and Pielke, 1981; which can be considered through the surface energy balance (SEB) equation (Oke, 1988): where Q * , Q 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 * , Q H , and Q E ) are more readily measured using standard techniques (e.g. eddycovariance instruments, radiometry) than Q S (Offerle et al., 2005;Pauwels and Daly, 2016;Roberts et al., 2006;Wang, 2012). For Q 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 Published by Copernicus Publications on behalf of the European Geosciences Union.
d. Parameterization as a function of Q * : 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 Q S in urban areas, result in the SEB residual (i.e. Q * +Q F −(Q H + Q E )) 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 * 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): where the a 1 , a 2 , and a 3 coefficients are for individual facets determined by least-square regression between Q S and Q * 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 Q S by Q * compared with other approaches. Despite the shortage of OHM coefficients for the wide range of facet types found in cities, OHM captures the urban Q S overall generally well (Grimmond and Oke, 1999;Järvi et al., 2011Järvi et al., , 2014Karsisto 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., 2011Järvi et al., , 2014Ward 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. (2003Gao et al. ( , 2008 was used to explore the physical relation of OHM coefficients a 1 and a 2 to the phase lag between Q S and Q * . However, insight into a 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(Gao et al., , 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).
where T is the temperature at a reference depth z (positive downward), t is time, λ is the thermal diffusivity, and W = ∂λ/∂z − C W /C g wϕ is the soil water flux density (Ren et al., 2000), with C W the volumetric heat capacity of water, C g the volumetric heat capacity of soil, w the pore water velocity, and ϕ the volumetric soil water content. The steady-periodic solution of Eq. (3) corresponding to the principal Earth rotation frequency (ω = 2π 24, in rad h −1 ), with boundary condition is given by (Gao et al., 2003(Gao et al., , 2010) where M = 2λ +W , N = ω , and with T S , A T S , and γ 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 where δ = arctan M N = arctan 2λω ( +W ) and k is the thermal conductivity. In particular, at the surface z = 0, the ground heat flux G 0 is given by and a simple written form of Q S (if only one surface) can be given as where η = δ − γ and c η = kA T S √ M 2 +N 2 MN . Although the above derivation only considers the land surface made of a single material type, the derived Q S (Eq. 8) can be adapted for surfaces made of composite materials or volumes given appropriate bulk/ensemble properties.
2.2 Parameterization of net all-wave radiation Q * for a land surface Given the parameterizations of incoming longwave radiation L ↓ , outgoing longwave radiation L ↑ , sensible heat flux Q H , latent heat flux Q E , and storage heat flux Q S as follows: The boundary condition imposed by the SEB relation can be rewritten as where the turbulent fluxes Q H and Q E are parameterized as functions of temperature gradient T S − T a with albedo α, bulk transfer coefficient C h , wind speed U , and Bowen ratio (β = Q H /Q E ). Theoretically, the second part of Eq. (10) (i.e. (1 − ε s ) L ↓ ) should be accounted for in the estimation of L ↑ (Oke, 1987); however, given that it is usually less than ∼ 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 ↓ and air temperature T a follow sinusoidal forms through a day as function of the mean value for the day (e.g. K ↓ ) (Sun et al., 2013), and introducing the solar radiation scale, and longwave radiation scale (assuming ε a ≈ ε s ≈ ε as a first-order estimate (as AnOHM is insensitive to this parameter; see Sect. 3.2); see clear sky of ∼ 0.85 (Staley and Jurica, 1972) and urban surfaces of ∼ 0.95 (Kotthaus et al., 2014): where τ denotes phase differences between T a and K ↓ , the f = f L + f T consists of the longwave energy redistribution factor: f L = 4εσ T tained: where ζ = arctan (N * /M * ), γ = ζ + arctan . The net all-wave radiation Q * is parameterized as where ϕ = arctan (χγ sin(γ )−sin(τ )) (f A * K )/(fLA * T )−(χγ cos(γ )−cos(τ )) and

Derivation of AnOHM coefficients
Based on the above parameterizations of Q * (Eq. 21) and Q S (Eq. 8), together with OHM for a specific surface: the coefficients can be readily derived from the parameterization in Sect. 2.2, as 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: 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 where a 3F (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 1 characterizes the ratio of Q S and Q * and depends on the energy scales (i.e. c η and c ϕ ) and their phase difference (i.e. η − ϕ). The energy scales, representing daily amplitudes of Q S and Q * , determine the overall magnitude, while the phase difference moderates the ratio value.
b. a 2 accounts for the temporal changes in Q S and Q * by including the principal Earth rotation frequency ω, in addition to the same determinants of a 1 (i.e. c η , c ϕ , and η − ϕ). The complementary sinusoidal functions, with phase difference (i.e. sin (η − ϕ) and cos (η − ϕ)), in the formulations of a 1 and a 2 are inversely related with a stronger lag effect from a 2 , and less contribution to Q S by Q * (i.e. smaller a 1 ).
c. a 3 (or a 3F ) indicates the baseline Q S determined by energy redistribution factors (i.e. f T and f ) and energy inputs (i.e. K ↓ , and Q F if anthropogenic heat is considered) as well as a 1 . It can be inferred from Eq. (2) that the nocturnal Q S is largely determined by a 3 when the absolute values and variability of Q * are small at night. A larger daytime energy input (i.e. K ↓ , 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 Geosci. Model Dev., 10, 2875-2890, 2017 www.geosci-model-dev.net/10/2875/2017/ 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 1 at which P (Y > y 1 ). Then a Markov chain Monte Carlo (MCMC) procedure is used to generate samples of a given conditional probability p 0 , leading to the exceedance of y 1 in the earlier simulations. This procedure is repeated, for exceedance events F i at which P (Y > y i ) = p i 0 , i = 1, 2, 3, . . ., 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 level = 3) and a conditional probability of p 0 = 0.1i.e. at each level the highest 10 % of the outputs are considered to exceed the intermediate threshold. As such, the threelevel simulation can effectively capture a rare event with the target exceedance probability of 10 −4 (i.e. the probability of occurrence is less than 1 in 10 000) and generate appropriate samples of different conditional probabilities.
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 where i = 1, 2, . . ., N level is the index of conditional sampling level, E [X] is the expectation that the unconditional distribution of a specific uncertainty parameter X, while E X|Y > y i 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 ε are similar to the original proposal distributions (Fig. 1), implying weak dependence of a 1 , a 2 , and a 3 on these properties. However, for albedo (α) both a 2 and a 3 are sensitive, but a 1 is not; changes in inverse Bowen ratio (β −1 ) impact all three coefficients; and the bulk transfer coefficient C h impacts a 1 and a 2 , but has little effect on a 3 .
Using S (Eq. 28) to quantify this, it is found that the surface properties (k, C p , and ε) 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 α sensitivity analysis (Fig. 2), it is apparent that an increase in α will increase a 1 while decreasing a 2 and a 3 , whereas the reverse occurs for β −1 and C h (i.e. their decreases leads to larger a 2 and a 3 values but smaller a 1 ).
From this, the links between the key surface parameters and the storage heat flux can be considered. With an increase in α, there is reduced solar energy in the SEB. This reduces the temporal change in Q S (smaller a 2 ) and decreases the baseline value of Q S (smaller a 3 ); larger β −1 indicates that more available energy is dissipated by Q E than by Q H , leading to decreased T s and Q S (smaller a 1 ); a smaller portion of Q * will be dissipated by Q S (smaller a 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 , T 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 , K ↓ , and U ) and the phase lag τ between K ↓ and T a have large impacts. In terms of the greatest impact on the coefficients (a 1 , a 2 , and a 3 ): A K and U influences a 1 , τ impacts a 2 , and a 3 responds more to A K and K ↓ than the other variables.
Variables that strongly modulate the interactions between Q S and Q * can be informed by the S results (Fig. 4). For instance, a greater range in K ↓ (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 Q S (smaller a 1 ) but higher baseline Q S (larger a 3 ). This is consistent with a reduction inK ↓ having a decrease in a 3 . The temporal change in Q S is highly correlated with the change in τ , an increase in which implies a   Luo et al. (2007) slower response of the surface to solar radiation and an overall decrease in Q S (smaller a 1 , a 2 , and a 3 ). The greater sensitivity to τ of a 2 is a key part of the original hysteresis nature of the heating/cooling of a surface. The sensitivity responses of a 1 , a 2 , and a 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 Q 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.

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;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 Q S is thus estimated as the modified residual of urban energy balance as Q S = Q * + 0.75Q F − 1.2 (Q H + Q E ) (Kotthaus and Grimmond, 2014a, b), which is then used in this evaluation. A similar approach for estimating Q S (i.e. residual of surface energy balance, Q S = Q * + Q F − (Q H + Q E )) is applied at the other (non-urban) sites but with Q F = 0.
AnOHM is first calibrated with observations under sunny conditions, when the assumptions of AnOHM are best satisfied (i.e. diurnal cycles of K ↓ and T a follow sinusoidal forms), to obtain surface properties required by AnOHM (Table 3). As the Bowen ratio β varies daily and monthly (Kotthaus and Grimmond, 2014a, b), β is either determined as the daily value if available, or based on the observation-based monthly climatology (Table 3). The seasonality in albedo α is accounted for also by using its monthly climatology (Table 3). AnOHM is driven by atmospheric forcing (i.e. K ↓ , T a , and U ) and/or their derived scales (A K , K ↓ , A T , T a , and τ ) to generate the OHM coefficients (i.e. a 1 , a 2 , and a 3 , see Geosci. Model Dev., 10, 2875-2890, 2017 www.geosci-model-dev.net/10/2875/2017/ Figure 1. Histograms of conditional samples at different conditional levels for surface property parameters (rows from top: thermal conductivity k in W m −1 K −1 , heat capacity C p in MJ m −3 K −1 , albedo α, emissivity ε, inverse Bowen ratio β −1 , and bulk transfer coefficient C h in J m −3 K −1 ) with AnOHM coefficients as the model output (columns from left: a 1 , a 2 and a 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.  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). Intraannual 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 1 (Fig. 5a) and a 2 (Fig. 5b) but general similarity in a 3 (Fig. 5c). Compared to a 1 and a 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 3 (see Eqs. 25 and 27), such seasonality in a 3 is predicted by AnOHM, and evident  Figure 3. Histograms of conditional samples at different conditional levels for ambient forcing parameters (rows from top: incoming solar radiation amplitude A K in W m −2 and its daytime mean K ↓ in W m −2 , air temperature amplitude A T in • C and its daily mean T a in • C, the phase lag τ in rad between K ↓ and T a , wind speed U in m s −1 , and water flux density W in m s −1 ) with AnOHM coefficients as the model output (columns from left: a 1 , a 2 and a 3 ). As Fig. 1. in the observations (Fig. 5c, also Ward et al., 2013). Larger K ↓ in warm seasons (May-September) will lead to smaller a 3 (Eqs. 25, 27) and vice versa. The AnOHM simulated and observed Q S agree well at the five different land cover sites, with RMSE values of ∼ 30 W m −2 . For comparison purposes, it is noted that the urban land surface model comparison (Best and Grimmond, 2015;Grimmond et al., 2011) found Q S to be the most poorly represented among all the SEB components with the best RMSE values of 53 W m −2 (Lipson et al., 2017). Although the much smaller Q S RMSE obtained by AnOHM uses a prescribed Bowen ratio in the offline evaluation, such improvement indicates the ability of AnOHM to simulate Geosci. Model Dev., 10, 2875-2890, 2017 www.geosci-model-dev.net/10/2875/2017/ Table 3. Surface properties used in AnOHM simulation for the study sites based on calibration. The values of α and β are monthly climatology from January to December and are used when observations are not available (see Table 1 for notation definition).   (Anandakumar, 1999) are shown. The LOESS (Cleveland and Devlin, 2012) fitting is a locally weighted polynomial regression approach. a more consistent Q S with observations. Compared with OHM predictions (orange lines in Fig. 6), AnOHM (blue lines in Fig. 6) better reproduces the seasonality in Q 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 Q S by AnOHM. Overall, the evaluation demonstrates good performance of AnOHM in predicting the long-term Q 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 Q S and Q * . The sensitivity of AnOHM to surface properties and hydrometeorological conditions is analysed through Monte Carlo-based subset simulations (Au and   Table 2 for site information). Statistics include average bias and RMSE (W m −2 ). The OHM coefficients a 1 , a 2 , and a 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). 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 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 ↓ and air temperature T a diurnal cycles are sinusoidal, optimal performance of AnOHM occurs under clear-sky conditions. The current parameterizations of K ↓ and T a within AnOHM only consider the harmonics of principal frequencies for formulation simplicity. More frequencies may potentially resolve more realistic diurnal variations in K ↓ and T a . As the reflected part of L ↓ (i.e. (1 − ε s ) L ↓ ) is assumed negligible, and similar emissivity values are assumed for sky and land surface (i.e. ε s ≈ ε a ≈ ε), 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 ε 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 landatmosphere 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 Q S . Compared to OHM, AnOHM has the benefit of allowing Q 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 ↓ , T a and U ) to other land surface models (LSMs), it can potentially be used within in LSMs to estimate Q 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.
T. Sun et al.: The Analytical Objective Hysteresis Model (AnOHM v1.0) Appendix A: Rationale for a simplified formulation of outgoing longwave radiation In the formulation of outgoing longwave radiation L ↑ , a simplified form (i.e. ε s σ T 4 s ) is used for AnOHM by ignoring part 2 of Eq. (10) (i.e. (1 − ε s ) L ↓ ). The rationale for such simplification is that given ε s is usually larger than 0.9, (1 − ε s ) L ↓ 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 ↓ = ε a σ T 4 a ≈ ε s σ T 4 a ), we conduct a sensitivity analysis of the ratio between the ignored part (i.e. (1 − ε s ) L ↓ ) and total outgoing longwave radiation (i.e. ε s σ T 4 s + (1 − ε s ) L ↓ ) at a constant air temperature of 20 • C and find this ratio is generally less than 5 % given ε s ranges between 0.90 and 0.99 (Fig. A1).
Moreover, if (1 − ε s ) L ↓ is included in the net longwave radiation, the induced effect can be incorporated into a modified sky emissivity ε a = ε s ε a as follows: Then by assuming ε ≈ ε a ≈ ε s , the derivation following Eq. (18) still holds. The sensitivity analysis suggests that the derived coefficients are insensitive to ε (see S for ε in Fig. 2). As such, we deem the omission of (1 − ε s ) L ↓ will not qualitatively change the results of this work.  (Grimmond and Oke, 1999), whereas the green lines show median values derived from results by observation regression at corresponding sites.