A parameterization of sub-grid topographical effects on solar radiation in the E3SM Land Model (version 1.0): implementation and evaluation over the Tibetan Plateau

Topography exerts significant influences on the incoming solar radiation at the land surface. A few stand-alone regional and global atmospheric models have included parameterizations for sub-grid topographic effects on solar radiation. However, nearly all Earth system models (ESMs) that participated in the Coupled Model Intercomparison Project (CMIP6) use a plane-parallel (PP) radiative transfer scheme that assumes that the terrain is flat. In this study, we incorporated a well-validated sub-grid topographic (TOP) parameterization in the Energy Exascale Earth System Model (E3SM) Land Model (ELM) version 1.0 to quantify the effects of sub-grid topography on solar radiation flux, including the shadow effects and multi-scattering between adjacent terrain. We studied the role of sub-grid topography by performing ELM simulations with the PP and TOP schemes over the Tibetan Plateau (TP). Additional ELM simulations were performed at multiple spatial resolutions to investigate the role of spatial scale on sub-grid topographic effects on solar radiation. The Moderate Resolution Imaging Spectroradiometer (MODIS) data was used to compare with the ELM simulations. The results show that topography has non-negligible effects on surface energy budget, snow cover, snow depth, and surface temperature over the TP. The absolute differences in surface energy fluxes for net solar radiation, latent heat flux, and sensible heat flux between TOP and PP exceed 20, 10, and 5 W m⁻², respectively. The differences in land surface albedo, snow cover fraction, snow depth, and surface temperature between TOP and PP exceed 0.1, 0.1, 10 cm, and 1 K, respectively. The magnitude of the sub-grid topographic effects is dependent on seasons and elevations and is also sensitive to the spatial scales. Although the sub-grid topographic effects on solar radiation are larger with more spatial details at finer spatial scales, they cannot be simply neglected at coarse spatial scales. When compared to MODIS data, incorporating the sub-grid topographic effects overall reduces the biases of ELM in simulating surface energy balance, snow cover, and surface temperature, especially in the high-elevation and snow-covered regions over the TP. The inclusion of sub-grid topographic effects on solar radiation parameterization in ELM will contribute to advancing our understanding of the role of the surface topography on terrestrial processes over complex terrain.

1 Introduction

Earth system models (ESMs), which simulate the interactions between atmosphere, land, ocean, and cryosphere systems, are powerful tools for understanding, reconstructing and projecting the Earth's climate (Bonan and Doney, 2018). Land surface models (LSMs) in ESMs represent the terrestrial water, energy, and carbon cycles (Dickinson et al., 2006). However, most of the state-of-the-art LSMs necessarily adopt some oversimplified and unrealistic schemes to treat the transfer of radiation, heat, water, and carbon. For example, lateral transport of water and energy in the subsurface and sub-grid topographic effects on solar radiation are neglected (Y. Fan et al., 2019). These simplifications could lead to large uncertainties, especially at finer spatial scales (Fisher and Koven, 2020; Prentice et al., 2015).

The Energy Exascale Earth System Model (E3SM), a relatively new fully coupled ESM supported by the U.S. Department of Energy (DOE), aims to tackle the grand challenge of actionable predictions of Earth system variability and change (Leung et al., 2020; Golaz et al., 2019). With the capacity to run at relatively high resolution (Caldwell et al., 2019) and include more realistic human–natural processes (Zhou et al., 2020), E3SM provides a good opportunity to better understand the complex Earth system processes and their interactions. However, improving the representations of the complex, multi-scale processes in the Earth system is important to more fully realize the benefits of high-resolution modeling.

As the horizontal grid spacing of ESMs increases, topography is expected to exert more significant influences on many land surface processes, including surface energy balance, surface hydrology, and snowmelt. The incoming and reflected solar radiations, as well as their direct and diffuse components, depend on surface topography (Dubayah and Rich, 1995; Hao et al., 2019a, b). Topography modifies the direct radiation reaching the Earth surface through self-shadowing or blocking by adjacent topography. Topography also decreases the diffuse radiation from sky by decreasing the portion of the visible sky and increases the reflected radiation from adjacent topography due to the multi-scattering effects (Dubayah, 1992; Proy et al., 1989). The changes in net solar radiation due to topography significantly influence surface energy budget (Gu et al., 2012; Lee et al., 2019; Liou et al., 2007), surface hydrology (Lee et al., 2015; Zhang et al., 2018), snowmelt (Zaramella et al., 2019), precipitation (Gu et al., 2020), and vegetation distribution (Alexander et al., 2016). Incorporating the sub-grid topographic effects on solar radiation into LSMs such as the E3SM Land Model (ELM) is key to enhancing our understanding and modeling of surface processes and land–atmosphere interactions in regions of complex terrain, with potential remote effects through excitation of Rossby waves in the atmosphere (Koster et al., 2016).

However, nearly all ESMs (including E3SM) that participated in the Coupled Model Intercomparison Project Phase 6 (CMIP6) neglect the sub-grid topographic effects on solar radiation. Sub-grid topographic effects have been recognized and parameterized in a few regional weather and climate models (Arthur et al., 2018; Gu et al., 2020) and global climate models (Lee et al., 2015). Most CMIP6-class ESMs adopt simple plane-parallel (PP) radiative transfer schemes based on a two-stream approximation, which assumes that topography is flat (Dai et al., 2004; Dickinson, 1983; Sellers, 1985). Such simplified radiation parameterizations do not account for sub-grid topographic effects and can lead to large systematic biases in simulating land surface processes over complex terrain (X. Fan et al., 2019; Lee et al., 2019; Song et al., 2020). Song et al. (2020) reported that both CLM4.5 and CLM5.0 failed to capture the asymmetric diurnal cycles of solar radiation, surface albedo and carbon fluxes in a mountainous rainforest in Costa Rica. Lee et al. (2019) showed that accounting for the sub-grid topographic effects in the Community Land Model (CLM) 4.0 with a spatial resolution of 0.9^∘ × 1.25^∘ reduced the biases of reflected solar radiation in winter over the Tibetan Plateau (TP). However, the sub-grid topographic effects on solar radiation at a fine spatial resolution (e.g., 0.125^∘); the contributions of different factors; the sensitivity to elevations, seasons, and spatial scales; and the consistencies with high-resolution observation data from remote sensing over the TP still need further investigations.

Sub-grid topographic parameterizations for solar radiation in the LSMs need to account for the effects of sub-grid topography without significantly increasing the computational cost. Sub-grid radiation fluxes can be explicitly calculated using a high-resolution digital elevation model (DEM) and then averaged to derive grid scale radiation fluxes (Dubayah, 1992). However, this approach involves a vast data processing and thus introduces substantial computational burden (Helbig and Löwe, 2012). Parameterizations for sub-grid topography based on the statistical characteristics of sub-grid topography (Dubayah, 1990; Essery and Marks, 2007; Gu et al., 2020; Helbig and Löwe, 2012; Lee et al., 2011; Müller and Scherer, 2005) provide a computationally efficient approach for LSMs. Lee et al. (2011) used 3D Monte Carlo photon tracing simulations to develop a parameterization scheme where a set of multiple linear regression equations associate the sub-grid topographic effects on solar radiation with the domain-averaged topographic factors. The parameterization scheme developed by Lee et al. (2011) is computationally efficient because the domain-averaged topographic factors can be calculated a priori based on high-resolution DEM. This parameterization has been successfully applied in the Weather Research and Forecasting (WRF) model (Gu et al., 2012; Liou et al., 2013), CLM4.0 (Lee et al., 2015, 2019), and Taiwan Earth System Model Version 1 (TaiESM) (Lee et al., 2020) and is promising for incorporation in ELM.

The objective of this study is to update and evaluate the radiative transfer scheme to account for sub-grid topographic effects on solar radiation in ELM. We implemented the computationally efficient and physically realistic sub-grid parameterization scheme for solar radiation of Lee et al. (2011) into ELM. ELM simulations over the TP were performed with and without the sub-grid topographic parameterizations for solar radiation from 2000 to 2010 at multiple spatial resolutions. The sub-grid topographic effects on surface energy balance, snow cover and depth, and surface temperature were investigated based on the ELM simulations. The contributions of different factors to the sub-grid topographic effects and the dependence of the sub-grid topographic effects on seasons, elevations, and spatial scales were also analyzed. A suite of remotely sensed data from the Moderate Resolution Imaging Spectroradiometer (MODIS) were used to compare with the ELM simulations with different parameterizations for solar radiation in different seasons.

2 Materials and methods

2.1 Model overview

ELM (Version 1.0) is based on the Community Land Model Version 4.5 (CLM4.5) (Golaz et al., 2019). ELM calculates canopy radiation flux using the two-stream approximation methods; snow albedo using the Snow, Ice, and Aerosol Radiative Model (SNICAR) model (Flanner et al., 2007); and snow cover fraction based on snow water equivalent (Swenson and Lawrence, 2012). ELM also represents the snow hydrological processes including snowfall accumulation, melting, refreezing, compaction, aging, water transfer across layers, etc. New features in ELM to better represent land surface processes include an updated representation of soil hydrology, improved treatment of ecosystem carbon dynamics, a novel topography-based sub-grid spatial structure, and an irrigation scheme constrained by water management (Bisht et al., 2018; Tang and Riley, 2018; Tesfa and Leung, 2017; Zhou et al., 2020).

2.2 Sub-grid topographic parameterizations for solar radiation

The incoming solar radiation for a flat surface is composed of direct radiation ($F_{\mathrm{dir}}^{\mathrm{PP}}$) from sun, diffuse radiation ($F_{\mathrm{dif}}^{\mathrm{PP}}$) from sky, and coupled radiation ($F_{\mathrm{couple}}^{\mathrm{PP}}$) that represents surface reflected radiation that is further reflected or scattered by atmospheric particles. ELM v1.0 assumes flat surfaces and accounts for $F_{\mathrm{dir}}^{\mathrm{PP}}$ and $F_{\mathrm{dif}}^{\mathrm{PP}}$, while neglecting $F_{\mathrm{couple}}^{\mathrm{PP}}$. The solar radiation scheme of ELM v1.0 uses the two-stream approximations (Oleson et al., 2013). In contrast, the incoming solar radiation parametrization of Lee et al. (2011) over mountainous regions includes five components (illustrated in Fig. 1): (1) direct flux ($F_{\mathrm{dir}}^{\mathrm{TOP}}$) represents photons that are transmitted from the sun to the ground surface without encountering any reflection or scattering, (2) the direct-reflected flux ($F_{\mathrm{rdir}}^{\mathrm{TOP}}$) represents photons that are not scattered photons reflected by surrounding terrain, (3) diffuse flux ($F_{\mathrm{dif}}^{\mathrm{TOP}}$) represents photons that are scattered by atmospheric particles but are not reflected by the ground surface, (4) diffuse-reflected flux ($F_{\mathrm{rdif}}^{\mathrm{TOP}}$) represents scattered photons reflected by surrounding terrain, and (5) coupled flux ($F_{\mathrm{couple}}^{\mathrm{TOP}}$) represents remaining photons that are reflected multiple times or scattered by ground surface and atmospheric particles. $F_{\mathrm{dir}}^{\mathrm{TOP}}$ is different from $F_{\mathrm{dir}}^{\mathrm{PP}}$ because of the adjustment of solar illumination geometry and shadowing effects. $F_{\mathrm{dif}}^{\mathrm{PP}}$ is different from $F_{\mathrm{dif}}^{\mathrm{TOP}}$ because the sky hemisphere is occluded by adjacent terrain. Lee et al. (2011) used the radiation fluxes over flat surfaces (i.e., $F_{\mathrm{dir}}^{\mathrm{PP}}$ and $F_{\mathrm{dif}}^{\mathrm{PP}}$) to calculate the radiation fluxes over mountainous terrain based on sub-grid topographic factors. The relative deviation (f_dir) of direct flux between flat surfaces and mountains under the same atmospheric condition is defined as follows:

$$\begin{array}{}\text{(1)}& f_{\mathrm{dir}}={\displaystyle \frac{F_{\mathrm{dir}}^{\mathrm{TOP}}-F_{\mathrm{dir}}^{\mathrm{PP}}}{F_{\mathrm{dir}}^{\mathrm{PP}}}}.\end{array}$$

The relative deviation (f_rdir) of direct-reflected flux over mountains to direct flux over flat surfaces is defined as follows:

$$\begin{array}{}\text{(2)}& f_{\mathrm{rdir}}={\displaystyle \frac{F_{\mathrm{rdir}}^{\mathrm{TOP}}}{F_{\mathrm{dir}}^{\mathrm{PP}}}}.\end{array}$$

Similarly, the relative deviations (f_dif and f_rdif) of diffuse and diffuse-reflected fluxes are expressed as follows:

$$\begin{array}{}\text{(3)}& {\displaystyle }{f}_{\mathrm{dif}}={\displaystyle \frac{F_{\mathrm{dif}}^{\mathrm{TOP}}-F_{\mathrm{dif}}^{\mathrm{PP}}}{F_{\mathrm{dif}}^{\mathrm{PP}}}},\text{(4)}& {\displaystyle }{f}_{\mathrm{rdif}}={\displaystyle \frac{F_{\mathrm{rdif}}^{\mathrm{TOP}}}{F_{\mathrm{dif}}^{\mathrm{PP}}}}.\end{array}$$

In theory, these four relative deviations (i.e., f_dir, f_rdir, f_dif and f_rdif) depend on solar illumination geometry and sub-grid topographic distribution. Based on a series of 3D Monte Carlo photon tracing simulations, Lee et al. (2011) built a multiple linear regression parameterization to predict these four relative deviations well. The parameterization of Lee et al. (2011) uses four variables that include the standard deviation of elevation (σ_h) within a grid cell, grid-averaged values of cosine of the local solar incident angle ($\stackrel{\mathrm{‾}}{\mathit{\mu }}$), sky view factor ($\stackrel{\mathrm{‾}}{V_{\mathrm{d}}}$), and terrain configuration factor ($\stackrel{\mathrm{‾}}{C_{\mathrm{T}}}$). Lee et al. (2011) parameterization is given as follows:

$$\begin{array}{}\text{(5)}& {\left[f_{\mathrm{dir}}{f}_{\mathrm{rdif}}{f}_{\mathrm{rdir}}{f}_{\mathrm{rdif}}\right]}^T=\mathbf{A}\cdot {\left[\stackrel{\mathrm{‾}}{\mathit{\mu }}{\mathit{\sigma }}_{\mathrm{h}}\stackrel{\mathrm{‾}}{V_{\mathrm{d}}}\stackrel{\mathrm{‾}}{C_{\mathrm{T}}}\mathrm{1}\right]}^T,\end{array}$$

where A represents the fitted parameter matrix, which was obtained based on the data generated by the 3D Monte Carlo simulations. The sky view factor (V_d) represents the portion of visible sky limited by surrounding terrain (Zakšek et al., 2011), while the terrain configuration factor (C_T), the counterpart of the sky view factor, represents the portion of surrounding terrain which is visible to the ground target (Dozier and Frew, 1990). For an unobstructed infinite slope with the slope of α and aspect of β and a given solar illumination geometry (i.e., solar zenith angle, SZA, and solar azimuth angle, SAA), the cosine of the local solar incident angle (μ) can be calculated by

$$\begin{array}{}\text{(6)}& \begin{array}{rl}\mathit{\mu }& =\cos \left(\mathrm{SZA}\right)\cdot \cos \left(\mathit{\alpha }\right)+\sin \left(\mathrm{SZA}\right)\cdot \sin \left(\mathit{\alpha }\right)\\ & \cdot \cos (\mathrm{SAA}-\mathit{\beta }).\end{array}\end{array}$$

The SZA and SAA are assumed to be constant within a grid cell, but α and β vary within a grid cell. The grid cell average solar incident angle, $\stackrel{\mathrm{‾}}{\mathit{\mu }}$, can be expressed as follows:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}\stackrel{\mathrm{‾}}{\mathit{\mu }}& =\stackrel{\mathrm{‾}}{\cos \left(\mathrm{SZA}\right)\cdot \cos \left(\mathit{\alpha }\right)}\\ & +\stackrel{\mathrm{‾}}{\sin \left(\mathrm{SZA}\right)\cdot \sin \left(\mathit{\alpha }\right)\cdot \cos (\mathrm{SAA}-\mathit{\beta })}\\ & =\cos \left(\mathrm{SZA}\right)\cdot \stackrel{\mathrm{‾}}{\cos \left(\mathit{\alpha }\right)}+\sin \left(\mathrm{SZA}\right)\cdot \cos \left(\mathrm{SAA}\right)\\ & \cdot \stackrel{\mathrm{‾}}{\sin \left(\mathit{\alpha }\right)\cdot \cos \left(\mathit{\beta }\right)}+\sin \left(\mathrm{SZA}\right)\cdot \sin \left(\mathrm{SAA}\right)\\ & \cdot \stackrel{\mathrm{‾}}{\sin \left(\mathit{\alpha }\right)\cdot \sin \left(\mathit{\beta }\right)},\end{array}\end{array}$$

where overlines represent grid-averaged values. To further improve the regression parameterization, μ, V_d, and C_T are normalized by cos (α). The land surface albedo is adjusted, instead of modifying incoming solar radiation, to maintain the surface energy conservation and the consistency between the surface and the first levels of atmosphere above the surface (Lee et al., 2015). Specifically, to keep the absorbed solar radiation of the ground surface unchanged, Lee et al. (2015) built the relationship between direct (${\mathit{\alpha }}_{\mathrm{dir}}^{\mathrm{TOP}}$) and diffuse (${\mathit{\alpha }}_{\mathrm{dif}}^{\mathrm{TOP}}$) albedo over mountains and those (${\mathit{\alpha }}_{\mathrm{dir}}^{\mathrm{PP}}$ and ${\mathit{\alpha }}_{\mathrm{dif}}^{\mathrm{PP}}$) over flat surfaces as follows:

$$\begin{array}{}\text{(8)}& {\displaystyle }{F}_{\mathrm{dir}}^{\mathrm{PP}}\cdot \left(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{dir}}^{\mathrm{TOP}}\right)=\left(F_{\mathrm{dir}}^{\mathrm{TOP}}+F_{\mathrm{rdir}}^{\mathrm{TOP}}\right)\cdot \left(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{dir}}^{\mathrm{PP}}\right),\text{(9)}& {\displaystyle }{F}_{\mathrm{dif}}^{\mathrm{PP}}\cdot \left(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{dif}}^{\mathrm{TOP}}\right)=\left(F_{\mathrm{dif}}^{\mathrm{TOP}}+F_{\mathrm{rdif}}^{\mathrm{TOP}}\right)\cdot \left(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{dif}}^{\mathrm{PP}}\right).\end{array}$$

Substituting Eqs. (1)–(4) into Eqs. (8)–(9) leads to

$$\begin{array}{}\text{(10)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{dir}}^{\mathrm{TOP}}=\mathrm{1}-\left(\mathrm{1}+f_{\mathrm{dir}}+f_{\mathrm{rdir}}\right)\cdot \left(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{dir}}^{\mathrm{PP}}\right),\text{(11)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{dif}}^{\mathrm{TOP}}=\mathrm{1}-\left(\mathrm{1}+f_{\mathrm{dif}}+f_{\mathrm{rdif}}\right)\cdot \left(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{dif}}^{\mathrm{PP}}\right).\end{array}$$

The parameterizations represented by Eqs. (5), (10), and (11) were implemented in ELM to account for the sub-grid topographic effects on solar radiation fluxes. In this study, the $F_{\mathrm{couple}}^{\mathrm{TOP}}$ is neglected due to its limited impacts and nonlinear relationship with land surface albedo (Lee et al., 2011) and will be further considered. Equation (5) was pre-calculated for different SZAs using high-resolution DEMs (see Sect. 2.5), which are shown in Tables S1 and S2 in the Supplement. These parameterizations were developed at a 10 km × 10 km spatial scale, and Lee et al. (2013) demonstrated that they can be applied to various spatial resolutions larger than 10 km × 10 km.

Figure 1 A schematic diagram of the five components of incoming solar radiation over mountains adapted from Lee et al. (2011). The components of incoming solar radiation included are (1) direct flux, (2) diffuse flux, (3) direct-reflected flux, (4) diffuse-reflected flux, and (5) coupled flux. Local solar zenith angle (θ), sky view factor (V_d), and terrain configuration factor (C_T) are also marked.

2.3 Model setup and experiment design

The Tibetan Plateau (TP), also known as the Third Pole, plays an important role in regulating the Earth climate system (Lu et al., 2020; Yang et al., 2009). The TP has complex topographic features, where the central part is relatively flat and the western and southern regions have remarkable terrain undulations (Fig. 2). Figure S1 shows the heterogeneous spatial variations of the topographic factors used in the solar radiation parameterizations. Therefore, TP is an ideal region to study topography-related land surface processes in ELM.

Figure 2 Geographical distributions of the (a) mean and (b) standard deviations of elevation derived from 90 m DEM at 0.125^∘ × 0.125^∘ spatial resolution over the TP. The black line represents the contour line of 1.5 km.

Offline ELM simulations over the TP were performed for the period of 2000 to 2010 with and without the topographic parameterization, and the simulations are denoted as TOP and PP, respectively. The simulations were performed in the prescribed satellite vegetation phenology mode in which seasonally varying leaf area index is prescribed based on the MODIS data (Lawrence and Chase, 2007). The 3-hourly Global Soil Wetness Project meteorological forcing dataset version 1 (GSWP3v1) (Dirmeyer et al., 2006; Yoshimura and Kanamitsu, 2013) with 0.5^∘ × 0.5^∘ spatial resolution was used to drive all the model simulations. The bilinear interpolation technique was used to downscale the GSWP3v1 data to the required spatial resolution, and the “coszen” (i.e., the cosine of the solar zenith angle)-based, nearest-neighbor, and linear-interpolation methods were used to downscale the solar, precipitation, and other data to the half-hourly temporal resolution, respectively. ELM was configured to run over the TP at five different spatial resolutions: r0125 (0.125^∘), r025 (0.25^∘), r05 (0.5^∘), f09 (about 1^∘), and f19 (about 2^∘). The model outputs were archived at half-hourly frequency. The impact of initial conditions on subsequent analysis was avoided by discarding the results of the first year.

2.4 Model analysis

The ELM-based simulations, TOP and PP, at r0125 resolution were used to analyze the sub-grid topographic effects on surface energy budget (i.e., land surface albedo, net solar radiation, sensible heat flux, and latent heat flux), snow cover fraction, snow depth, and surface temperature. Surface temperature was calculated from the emitted longwave radiation using the Stefan–Boltzmann law, with the assumption that surface emissivity is equal to 1. The seasonally averaged values were computed from the half-hourly ELM outputs for different seasons: winter (DJF), spring (MAM), summer (JJA), and autumn (SON). Both the absolute differences (i.e., TOP − PP) and relative differences (i.e., (TOP − PP) $/$ PP) were used to analyze the sub-grid topographic effects and their spatial patterns.

The relationship between sub-grid topographic effects and elevations was analyzed by dividing the elevations into four intervals: 1.5–2.5, 2.5–3.5, 3.5–4.5, and > 4.5 km, which account for about 11 %, 9 %, 14 %, and 23 % of the rectangular region shown in Fig. 2, respectively. Regions with elevations lower than 1.5 km were not included in the analysis due to their flat topography (Fig. 2). Grid cells with a mean slope of zero were also excluded from this analysis. Additionally, grid cells with zero snow cover fraction were excluded when analyzing results for snow cover and snow depth.

A random forest model was used to quantitatively analyze the relative contribution of different factors to the topography-driven differences. The random forest model is a regression tree-based bootstrapped non-parametric machine learning model, which allows the calculation of the variable importance by estimating the out-of-bag (OOB) errors (Breiman, 2001). The OOB error represents the averaged prediction error for each sample z_i, calculated by only using the prediction trees that do not include z_i in the bootstrap samples. The variable importance can be measured by permutating the driving variables and then calculating the averaged decrement of OOB errors after permutation. Specifically in this study, based on Eqs. (5) and (7), we selected the quantities $\stackrel{\mathrm{‾}}{\sin \left(\mathit{\alpha }\right)\cdot \cos \left(\mathit{\beta }\right)}$, $\stackrel{\mathrm{‾}}{\sin \left(\mathit{\alpha }\right)\cdot \sin \left(\mathit{\beta }\right)}$, σ_h, $\stackrel{\mathrm{‾}}{V_{\mathrm{d}}}$, $\stackrel{\mathrm{‾}}{C_{\mathrm{T}}}$, and the PP-simulated land surface albedo as the driving variables. Combined with the driving variables, all the ELM-derived seasonally averaged data was used to train the random forest model to measure the relative importance of different factors in controlling the sub-grid topographic effects.

ELM outputs for TOP and PP at the remaining four spatial resolutions (r025 to f19) were processed to derive seasonally averaged values for studying the sensitivity of the sub-grid topographic effects to spatial scales.

2.5 Remote sensing data

The Shuttle Radar Topography DEM (SRTM) data at 90 m spatial resolution were used to derive the topographic factors required for the TOP simulations. The spatial mean and standard deviations of elevation, slope, aspect, sky view factor, and terrain configuration factor were computed for each ELM grid cell at all five spatial resolutions.

Table 1 Specifications of the remote sensing data used in the study.

The MODIS data from 2001–2010 were used to compare with the ELM simulations. All MODIS data listed in Table 1 were downloaded from the Google Earth Engine Platform (Gorelick et al., 2017). Specifically, these data included both direct (i.e., black-sky) and diffuse (i.e., white-sky) albedo data from the daily MCD43A3 v6 products with 500 m spatial resolution (Schaaf et al., 2002), snow cover data from daily MOD10A1 v6 products at 500 m spatial resolution (Hall et al., 2002), both daytime and nighttime surface temperature data from the daily MOD11A1 v6 products with 1 km spatial resolution (Wan, 2014), and latent heat flux data from the 8 d MOD16A2 v6 products with 500 m spatial resolution (Mu et al., 2007, 2011). Only the MODIS pixels with good quality indicated by the quality assurance flag were used in the analysis. All MODIS data were upscaled using the area-weighted averaging method to conform with the ELM resolutions.

2.6 Comparison with remote sensing data

MODIS data (introduced in Sect. 2.5) were used to compare with both TOP and PP at r0125 resolution. All MODIS data from 2001–2010 was averaged to the seasonal scales. The MODIS instantaneous surface diffuse and direct albedo datasets were derived for the local solar noon, and the MODIS instantaneous surface temperature data were derived for daytime and nighttime corresponding to the MODIS overpass time: 10:30 and 22:30 (local solar time), respectively. The ELM-simulated surface albedo and surface temperature were extracted at the corresponding MODIS time to compute the seasonally averaged values. The consistencies between ELM simulations and MODIS data were evaluated by computing the difference between MODIS data and PP (δ_PP) and TOP (δ_TOP). Furthermore, the change in the bias with respect to the MODIS data was computed as $\left|{\mathit{\delta }}_{\mathrm{TOP}}\right|$–$\left|{\mathit{\delta }}_{\mathrm{PP}}\right|$.

3 Results

3.1 Sub-grid topographic effects on surface energy budget, surface temperature and snow cover/depth

Compared to PP, TOP overall has higher net solar radiation (Fig. 3) and lower land surface albedo (Figs. 4, S2 and S3) for all seasons. The net solar radiation for PP shows an expected and opposite spatial pattern to the land surface albedo. The absolute differences in net solar radiation between TOP and PP can be as large as around 20 W m⁻² for all seasons, and the relative differences in the winter season can be as large as 25 %. In some small portions of the northern regions, TOP also shows lower net solar radiation than PP in winter and autumn, possibly due to the self-shadow or cast-shadow from the surrounding terrain. For PP, the spatial differences in surface albedo between the northwest and southeast of the study region are caused by the spatial differences in snow cover (Fig. 4a and b). In summer, the land surface albedo in the western regions decreases due to snow melt. The land surface albedo for all seasons in the western and southern regions shows large absolute and relative differences between TOP and PP that can be as large as 0.1 and 50 %, respectively, during winter. The spatial pattern of the difference in land surface albedo between TOP and PP is similar to the heterogeneous spatial pattern of topography (Fig. 2).

Figure 3 PP-simulated net solar radiation for different seasons (top row), and absolute (middle row) and relative (bottom row) differences between TOP and PP for different seasons.

Figure 4 Absolute differences between TOP and PP for four seasons in different variables (from top to bottom): land surface albedo, snow cover fraction, snow depth, surface temperature, sensible heat flux, and latent heat flux.

Larger net solar radiation in TOP compared to PP leads to lower snow cover and depth and higher surface temperature (Figs. 4, S2 and S3). TOP has lower snow cover fractions for most western regions in winter and spring (Fig. 4). Compared to PP, the absolute and relative decreases in snow cover fraction of TOP can be larger than 0.1 and 20 %. The absolute difference in snow depth can reach up to 10 cm. Snow albedo feedback may have contributed to the large differences between TOP and PP, as larger net solar radiation in TOP reduces snow cover, which may further increase the net solar radiation. Surface temperature has a similar spatiotemporal pattern as the net solar radiation (Fig. 4). The absolute difference in surface temperature between TOP and PP is generally within 1 K for all seasons. The western regions have large differences in surface temperature and snow cover during winter.

TOP has higher sensible and latent heat fluxes than PP, due to the higher net solar radiation (Figs. 4, S2 and S3). TOP shows higher sensible heat flux than PP for all seasons, and the absolute and relative differences can be as large as 10 W m⁻² and 20 %, respectively (Figs. 4 and S3). The difference in the latent heat flux is smaller compared to the difference in the sensible heat flux (Fig. 4) and is generally within 5 W m⁻². But the relative difference in latent heat flux may be larger than 20 % in winter (Fig. S3). How the partitioning of surface heat flux between sensible and latent heat fluxes responds to the difference in net solar radiation between TOP and PP may vary by season and region depending on the soil moisture, vegetation, and other factors.

3.2 Contribution of different factors

The random forest model can well predict the sub-grid topographic effects on solar radiation with high coefficients of determination (R²) for all seasons (Fig. 5a–d), which demonstrates that the topographic factors can explain the difference between TOP and PP in land surface albedo well. Further variable importance analysis (Fig. 5e–h) shows that the contributions of different factors to the sub-grid topographic effects are different. The first two terms (i.e., $\stackrel{\mathrm{‾}}{\sin \left(\mathit{\alpha }\right)\cdot \cos \left(\mathit{\beta }\right)}$ and $\stackrel{\mathrm{‾}}{\sin \left(\mathit{\alpha }\right)\cdot \sin \left(\mathit{\beta }\right)}$), related to the sub-grid distribution of slope and aspect, can account for 62.5 % of the differences in surface albedo during winter (Fig. 5e). The slope and aspect affect the direct solar radiation, which dominates the total solar radiation under clear-sky conditions. The sky view factor, terrain configuration factor and land surface albedo for PP, which mainly affect the diffuse and reflected radiation, account for 2.7 %, 2.3 %, and 24.7 % in winter, respectively. The dominant factors for the differences between TOP and PP can be different in different seasons (Fig. 5e–h). In summer, the contributions of the first two terms decrease to 47.1 % (Fig. 5g). This is because the solar position (i.e., solar illumination geometry) is different in different seasons. In winter, the solar zenith angle is large over the TP and thus there are strong shadowing effects, while the sun is moves northward and gets closer to the nadir position from spring to summer, which can lead to the reduced shadowing effects. Similar results were obtained for other variables (e.g., net solar radiation and surface temperature) and thus are not shown in this paper.

Figure 5 The performance of the random forest modeling in predicting the relative difference in land surface albedo between TOP and PP for four seasons (a–d). The relative importance of different factors in predicting the differences in surface albedo between TOP and PP for four seasons (e–h). R² is the coefficient of determination and the different factors are described in the text.

3.3 Sensitivity to elevations

The differences of surface energy budget, surface temperature and snow cover and depth between TOP and PP show elevation-dependent patterns (Figs. 6–7). Generally, as the elevation increases, TOP has a lower land surface albedo and snow cover fraction than PP and therefore a higher net solar radiation, surface temperature, and sensible and latent heat flux than PP. These elevation-dependent patterns are similar for all seasons, although the differences between TOP and PP are larger in winter than in summer (Fig. 6). Taking land surface albedo in winter as an example (Fig. 6a), for elevation between 1.5–2.5 km, TOP has smaller values than PP in 53 % of the regions; for elevation between 2.5 and 3.5 km, the area fraction is 57 %; for elevation between 3.5 and 4.5 km, the area fraction is 68 %; and for elevation above 4.5 km, the area fraction is 74 %. At higher elevations, the larger decrease in land surface albedo of TOP leads to a larger increase in surface fluxes (Fig. 7a, e–f) and surface temperature (Fig. 7d), along with a larger decrease in snow cover and snow depth (Fig. 7b–c). In addition, the quantiles in Fig. 7 also show that as the elevation increases, the relative differences of net solar radiation, snow cover fraction, snow depth, surface temperature, and sensible and latent heat flux between TOP and PP can become larger, and the relative differences of land surface albedo can exceed 10 % for all elevation bands.

Figure 6 Boxplots of the relative differences in land surface albedo between TOP and PP for all seasons at four different elevation bands. Red points represent the mean values.

Figure 7 Boxplots of the relative (or absolute for surface temperature) differences in net solar radiation (a), snow cover fraction (b), snow depth (c), surface temperature (d), sensible heat flux (e), and latent heat flux (f) between TOP and PP in winter at four different elevation bands. Red points represent the mean values.

3.4 Sensitivity to spatial scales

The sub-grid topographic effects on surface energy balance, snow cover, and surface temperature are sensitive to the spatial scales. The sub-grid topographic effects on land surface albedo in winter show similar spatial patterns across spatial scales (Fig. 8a–e). There are similar trends of the sub-grid topographic effects on land surface albedo with elevations at different spatial scales (Fig. S4). The relative contributions of different topographic variables are similar at different spatial scales (Fig. S5). Larger spatial heterogeneity in land surface albedo is present at finer spatial scales, but the pattern is smoothed at coarser spatial scales (Figs. 8 and S4). As the spatial resolution becomes coarser, the terrain becomes flatter, and thus the differences between TOP and PP are smaller. However, the relative difference between TOP and PP can still be as large as 15 % at coarse spatial scales (i.e., f19; Fig. 8e). The statistical distributions of the relative differences in land surface albedo over the TP at different spatial scales are similar, with 0 % to −5 % as the frequent value (Fig. 8f). For snow cover, surface temperature, and other energy balance variables, similar results are noted from Figs. S6–S10, and the sub-grid topographic effects are still large even at a spatial resolution as low as around 2^∘. For instance, for the spatial resolution of f19, the relative differences of net solar radiation, sensible heat, and latent heat flux and snow cover fraction can be larger than 8 %, 20 %, 20 %, and 20 %, respectively. The absolute difference in surface temperature for f19 is within 0.1 K, but for f09 it is still as large as 0.5 K.

Figure 8 The relative differences of land surface albedo between TOP and PP in winter at different spatial scales (a–e) and the statistical histogram of their frequent distributions (f).

3.5 Comparison with MODIS data

Overall, TOP shows better consistencies with the MODIS land surface albedo data than PP (Figs. 9 and S11). In the western regions, PP overestimates direct and diffuse albedo in winter and underestimates them in spring (Fig. S11), possibly due to the bias of snow cover in the model simulations (Fig. S11). In most other regions, PP generally overestimates direct albedo for all seasons and underestimates diffuse albedo outside of summer. The bias in PP, δ_PP, for direct and diffuse albedo can exceed 0.2. Compared to PP, direct albedo of TOP overall has a smaller bias relative to the MODIS data in the western regions, except in spring (Fig. 9). The improvement of TOP in direct albedo can be larger than 0.1. However, for diffuse albedo, the performance of TOP in most regions is similar to or even worse than PP (Fig. 9). The difference in diffuse albedo between δ_TOP and δ_PP is within 0.02 in about 86 % of the whole domain in Fig. 9 in winter.

Figure 9 The change in the bias ($\left|{\mathit{\delta }}_{\mathrm{TOP}}\right|$–$\left|{\mathit{\delta }}_{\mathrm{PP}}\right|$) with respect to the MODIS data for four seasons in different variables (from top to bottom): direct albedo, diffuse albedo, snow cover fraction, daytime and nighttime surface temperature, and latent heat flux.

TOP generally outperforms PP in winter, when compared to MODIS snow cover, surface temperature, latent heat flux data (Figs. 9 and S11). In the western regions, PP has higher snow cover fractions than MODIS data in winter but lower snow cover fractions in other seasons. In other regions, PP has lower snow cover fractions in all seasons (Fig. S11). TOP has smaller biases relative to the MODIS data than PP in winter, and the absolute value of $\left|{\mathit{\delta }}_{\mathrm{TOP}}\right|$–$\left|{\mathit{\delta }}_{\mathrm{PP}}\right|$ can be larger than 0.1. TOP has slightly larger biases in spring, but there is no large difference between TOP and PP in summer and autumn due to the low snow cover. The spatial distribution of δ_PP in snow cover fraction is consistent with the pattern of biases in land surface albedo shown in Fig. S11. For daytime surface temperature, there is a larger difference between PP and MODIS, which can exceed 5 K. TOP can reduce the biases by ∼ 0.5–1 K in the central regions, especially in winter. For nighttime surface temperature, PP has systematically higher values than the MODIS data, and the difference between TOP and PP is small in summer and autumn but large in winter and spring (Fig. 9). For latent heat flux, there are big differences between PP and the MODIS data. In contrast, TOP has a slightly better performance than PP in winter, but for other seasons, TOP possibly has worse performance than PP when compared to the MODIS data (Fig. 9).

Figure 10 Boxplots of the differences in bias for TOP and PP ($\left|{\mathit{\delta }}_{\mathrm{TOP}}\right|$–$\left|{\mathit{\delta }}_{\mathrm{PP}}\right|$) with respect to MODIS data for (a) direct albedo, (b) diffuse albedo, (c) snow cover fraction, (d) daytime surface temperature, (e) nighttime surface temperature, and (f) latent heat flux in winter at four different elevation bands. Red points represent the mean values.

As the elevation increases, TOP shows higher consistencies with the MODIS data in winter (Fig. 10). When the elevation is below 3.5 km, TOP and PP have similar performance, but at higher elevations TOP overall has lower biases in direct albedo (Fig. 10a), snow cover fraction (Fig. 10c), daytime surface temperature (Fig. 10d), and latent heat flux (Fig. 10f). The bias in direct albedo is smaller in TOP as compared to PP for 54 % and 63 % of the study region in elevation bands 3.5–4.5 and > 4.5 km, respectively. The difference in the bias for snow cover fraction between TOP and PP remains unchanged for the elevation bands 3.5–4.5 and > 4.5 km. TOP has smaller bias in daytime surface temperature as compared to PP for 57 % of the study region at elevation > 4.5 km. The bias in latent heat flux is smaller for TOP than PP for elevation band 3.5–4.5 and > 4.5 km for ∼ 60 % of the study region. The differences in bias between PP and TOP are small for diffuse albedo and nighttime surface temperature in most of the regions. For example, when the elevation is above 4.5 km, the difference in biases of diffuse albedo is within 0.01 for 73 % of the regions and is within 0.02 for 91 % of the regions. For nighttime surface temperature, the difference in biases increases with elevation and is within 0.1 K in about 61 % of the regions when the elevation is above 4.5 km.

4 Discussion

Sub-grid topographic effects on solar radiation play an important role in surface energy balance, surface temperature, and snowmelt over complex terrain. Simply neglecting the sub-grid topography can lead to large errors in simulating surface energy balance. Compared to flat surfaces, the land surface albedo over the complex terrain of TP generally decreases and net solar radiation increases (Figs. 3, 4 and S3), which increases the surface temperature (Fig. 4). The snow cover fractions and snow depth decrease due to increased snow melt and possibly snow–albedo feedback (Fig. 4), which may alleviate the snow depth overestimation over the TP in ESMs (Wei and Dong, 2015). The effects of sub-grid topography on solar radiation also show seasonal variations, which are more pronounced in winter because larger solar zenith angles in winter over the TP can cause stronger shadowing effects (Hao et al., 2018b) and large snow cover areas in winter can cause stronger reflected radiation from adjacent topography (Helbig et al., 2010). In addition, the sub-grid topographic effects are elevation dependent (Figs. 6–7) because mountain tops with higher elevations tend to receive more solar radiation due to the topographic effects and thinner atmosphere, while valley areas with lower elevations receive relatively less solar radiation due to the shadowing effects (X. Fan et al., 2019; Lee et al., 2015). Compared to PP, TOP produces results more consistent with the MODIS observations, especially in the high-elevation and snow-covered regions over the TP (Figs. 9–10). Generally, direct albedo of TOP shows higher consistencies with MODIS data than PP, when snow cover fraction is larger or the snow cover fraction of TOP has higher consistencies with MODIS (Fig. 11). These demonstrate that accounting for the sub-grid topographic effects over complex terrain improves the performance of ELM. In a high-resolution coupled model, the highly concentrated differences between TOP and PP along the southern edge of the TP could lead to important differences in simulating clouds, convection, terrain-induced circulation. and transport of aerosols, with potentially important implications for modeling the South Asian monsoon and its hydrologic impacts. Future studies including the sub-grid topographic effects in coupled simulations will address their impacts on coupled land–atmosphere processes.

Figure 11 Relationship between the differences in bias for TOP and PP ($\left|{\mathit{\delta }}_{\mathrm{TOP}}\right|$–$\left|{\mathit{\delta }}_{\mathrm{PP}}\right|$) with respect to MODIS data for direct albedo and PP simulated snow cover fraction (a) or the differences in bias for TOP and PP ($\left|{\mathit{\delta }}_{\mathrm{TOP}}\right|$–$\left|{\mathit{\delta }}_{\mathrm{PP}}\right|$) for snow cover fraction (b) in winter. The red line is the regression line, and R is the correlation coefficient.

Sub-grid topographic effects are strongly dependent on spatial scales. The sub-grid topographic effects are more pronounced at the finer resolution (Figs. 8 and S6–S10) and tend to be spatially smoothed at a coarse resolution (Lee et al., 2011, 2013). Therefore, it is necessary to consider the sub-grid topographic effects on solar radiation in high-resolution land surface modeling. However, the relative differences in net solar radiation between TOP and PP can still reach up to 8 % in some regions even at the coarse spatial resolution of 2^∘ (Fig. S6). This demonstrates that the sub-grid topographic effects on solar radiation cannot be neglected even for simulations at coarse spatial resolutions.

Uncertainties in remote sensing data may affect their reliability as ground truth for evaluating the ELM simulations. The MODIS land surface albedo products have shown good consistencies with ground measurements (Moustafa et al., 2017; Wang et al., 2004), but the semi-empirical kernel-driven-model-based algorithms used to derive the MODIS land surface albedo do not account for topography explicitly (Schaaf et al., 2002; Hao et al., 2020), which may lead to large errors over rugged terrain (Hao et al., 2018a, b). MODIS snow cover data have shown relatively poor performance when compared to ground measurements, especially over the regions of TP with higher elevation and shallower snow depth (Pu et al., 2007; Yang et al., 2015; Zhang et al., 2019). The accuracy of MODIS surface temperature products depends on the accuracy of land cover products and the prescribed surface emissivity values (Duan et al., 2019). The MODIS evapotranspiration product is sensitive to the algorithm used to account for the environmental stresses over the TP, as well as the atmospheric forcing data used to generate the product (X. Li et al., 2019). However, the topography-induced differences between TOP and PP can be comparable to the errors of MODIS data. For example, Wang et al. (2004) reported that compared to ground measurements, MODIS albedo had a maximum error of 0.036 in a semi-desert region on the TP, which is smaller than the maximum difference of 0.1 between TOP and PP (Fig. 4). Wang et al. (2007) showed that the mean and maximum errors of MODIS surface temperature were 0.27 and 2.61 K, respectively, at a semi-desert site on the western TP, which is comparable to the maximum difference of 1 K between TOP and PP (Fig. 4). Salomonson and Appel (2004) showed that using the Landsat 30 m observations as the benchmark, the mean error of MODIS snow cover fraction was smaller than 0.1, which is comparable to the difference of 0.1 between TOP and PP (Fig. 4). Mu et al. (2007) showed that the 8 d MODIS latent heat flux had a mean bias from −5.8 to 39.9 W m⁻², possibly larger than the difference between TOP and PP in our study (Fig. 4). In addition, the SRTM data, used to derive the topographic factors for the parameterization, have shown large errors in some regions (Grohmann, 2018; Mukherjee et al., 2013). More accurate topographic factors can be derived using globally consistent, high-quality DEM data such as the Copernicus 30 m global Digital Elevation Model (GLO-30) (https://spacedata.copernicus.eu, last access: 11 October 2021). The quality of remote sensing data needs to be validated comprehensively before its use in evaluation of LSMs.

The inclusion of sub-grid topographic parameterizations for solar radiation in ELM improves the representations of surface energy balance to some degree, but many shortcomings in ELM's existing radiative transfer modeling scheme limit the potential for further improving the ELM simulations. The 1D two-stream approximation method used in ELM represents the vegetation canopy as a homogeneous “big leaf” (Yuan et al., 2017) and neglects the vertical multi-layer structure (Bonan et al., 2018) and the horizontal leaf clumping (Bailey et al., 2020; Braghiere et al., 2020; C. Li et al., 2019). In the snow-covered regions, the ELM parameterizations for the effects of snow impurities (i.e., black carbon and dust mixing) on light scattering and absorption processes need to be refined to account for internal mixing and non-spherical shapes of snow grains (Dang et al., 2019; He et al., 2018). In addition, ELM also does not account for the influence of adjacent terrain on longwave thermal radiation (Yan et al., 2020). These may partly explain the inconsistencies between ELM simulations and MODIS data, especially for diffuse albedo and nighttime surface temperature (Fig. 10).

In this study, the same atmospheric forcings were used in the simulations at different spatial scales, which could be a source of error at a finer resolution (Fiddes and Gruber, 2014; Tesfa et al., 2020). Furthermore, the sub-grid parameterizations neglect the spatial correlation between sub-grid topography and plant functional types. The spatial pattern of vegetation types generally depends on the topographic distribution, which controls terrestrial water, energy, water, and carbon cycle (Reed et al., 2009). These aforementioned simplifications may affect the accurate representations of the sub-grid topographic effects on solar radiation in ELM at a coarse resolution. Combining the sub-grid topographic parameterizations for solar radiation implemented in ELM in this study with ELM's new sub-grid topography structure (Tesfa et al., 2017) and downscaling of atmospheric forcing (Tesfa et al., 2020) is anticipated to further improve the representations of the land surface processes at different spatial scales (Ke et al., 2013). A future study will investigate the impact of sub-grid topographic parameterizations for solar radiation on the land–atmosphere interactions by performing ELM simulations with an active atmospheric model.

5 Conclusions

The computationally efficient sub-grid topographic parameterizations for solar radiation of Lee et al. (2011) were implemented in ELM in this study. The results show that topography has large effects on surface energy budget, snow cover and depth, and surface temperature that cannot be neglected. The absolute differences with and without accounting for sub-grid topography on net solar radiation, sensible heat flux, and latent heat flux exceed 20, 10, and 5 W m⁻², respectively. Similarly, the differences in land surface albedo, snow cover fraction, snow depth, and surface temperature exceed 0.1, 0.1, 10 cm, and 1 K, respectively. Nearly all the relative differences of these variables, except surface temperature, reach up to 20 %. The magnitude of the sub-grid topographic effects on solar radiation is seasonally-dependent and elevation-dependent, and is also sensitive to the spatial scales. Although the sub-grid topographic effects on solar radiation are larger at finer spatial scales, they cannot be simply neglected even at coarse spatial scales. For example, the relative difference in land surface albedo when accounting for sub-grid topography in winter reaches up to 15 % for the coarse spatial scale of 2^∘. ELM simulations with the sub-grid topographic parameterization for solar radiation have better agreement overall with the MODIS data for simulated surface energy balance, snow cover, and surface temperature over the TP. These results highlight the necessity of accounting for the sub-grid topographic effects in LSMs and show that our improvements in ELM are promising to advance understanding and modeling of the role of the surface topography on terrestrial processes.