Estimating daily surface downward shortwave radiation over rugged terrain without bright surface at 30 m on clear-sky days using CERES data

ABSTRACT

In this study, the authors propose a model, called the Daily Downward Shortwave Radiation Random Forest Model over Rugged Terrain (DSRMT), to accurately calculate the downward shortwave radiation over a terrain without bright surface on clear days at a daily scale (DSR_{daily−rugged}). It was built by using the random forest method based on the comprehensive samples from CERES4_SYN1deg_Ed4A within 17 typical mountainous regions. DSRMT could directly estimate DSR_daily-rugged from the instantaneous direct and diffuse solar radiation on a flat surface during 10:30–14:30hrs on each day by comparing with the terrain factors from a digital elevation model, broadband albedo from the Global Land Surface Satellite, and ancillary information. The in-situ validation results showed that it generally delivered superior performance in estimating DSR_daily-rugged at any time during 10:30–14:30hrs, especially at noon, yielding a validated root mean-squared error (RMSE) of 24.90–29.22 Wm⁻² and mean absolute error (MAE) of 19.16–22.94 Wm⁻², and the average weighted DSR_daily-rugged were usually more accurate with the RMSE and MAE of 21.63 and 17.14 Wm⁻². Overall, DSRMT was found to deliver satisfactory performance because of its high accuracy, robustness, ease of implementation, and efficiency, so it has the strong potential to be widely used in practice.

KEYWORDS:

• Daily downward shortwave radiation
• rugged terrain
• estimation
• modeling
• DSRMT
• CERES4
• remote sensing

Nomenclature

Parameter	=	Meaning
DSR	=	Downward shortwave radiation
DSRrugged	=	DSR over rugged terrain
Ddir-rugged	=	Direct solar radiation over rugged surface
Ddif-rugged	=	Diffuse radiation over rugged surface
Dref-rugged	=	The reflected radiation of the surrounding terrain
Ddir-flat	=	Direct solar radiation over flat surface
Ddif-flat	=	Diffuse solar radiation over flat surface
Dins-dir-flat	=	The instantaneous Ddir-flat
Dins-dif-flat	=	The instantaneous Ddif-flat
DSRins-flat	=	The instantaneous DSR values over flat surface
Dins-dir-rugged	=	The instantaneous Ddir-rugged
Dins-dif-rugged	=	The instantaneous Ddif-rugged
Dins-ref-rugged	=	The instantaneous Dref-rugged
DSRins-rugged	=	The instantaneous DSR values over rugged surface
DSRdaily-rugged	=	The daily DSR over rugged surface
DSRdaily−flat	=	The daily DSR over flat surface
$\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$	=	The final daily mean value
${\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}}_{\mathrm{DCF}}$	=	The final daily mean value using direct correction formula (DCF)
φ	=	Latitude
α	=	The surface broadband albedo
θ0	=	Solar zenith angle
φ0	=	Solar azimuth angle
θi	=	The angles between the incident ray and the slope surface normal (solar illumination angle)
u0	=	The cosine of θ0
us	=	The cosine of θi
Φ	=	A binary function that is Φ = 0 if obstructed and Φ = 1 if not
S	=	Slope
A	=	Aspect of slope
Vd	=	Sky view factors
Vc	=	Terrain view factors
CI	=	Clear index
Zd	=	The elevation of adjacent pixels
Z0	=	The elevation of the target pixel
d	=	The resolution of DEM
Gsc	=	The solar constant (0.0820 MJm−2·min−1)
dr	=	The inverse relative distance from the Earth to the Sun
ws	=	The sunset hour angle
δ	=	Sun declination

1. Introduction

The downward shortwave radiation (DSR; spectral range of 0.3–3.0 μm) is the part of solar radiation that persists at the top of the atmosphere (Letu et al. 2020). As an essential radiative component, the DSR plays an important role in driving the climate system and the hydrological cycle, which are important for agriculture and forestry, meteorology, weather forecasting, and climate monitoring (Letu et al. 2020; Yang et al. 2010; Wild 2016; Chen et al. 2012). With technological developments, the demands on the high spatial resolution of the DSR are growing in various applications, particularly at the regional scale. Examples include investigations of the heterogeneity of urban landscapes and local agricultural management (Ma et al. 2016; Gupta et al. 2004; Yan et al. 2018; Zhang et al. 2020; Gowda et al. 2007; Zhao et al. 2020). However, accurately calculating the DSR at a high spatial resolution, especially that those below 5 km, requires considering the significant influence of rugged terrains, which cover about 24% of the Earth' s surface, on it (Liou, Lee, and Hall 2007; Wang, Zhou, and Liiu 2004; Wen et al. 2018; Yan et al. 2020). Otherwise, the uncertainty in the instantaneously estimated DSR may be as large as 600 Wm⁻² (Hao et al. 2021; Wang et al. 2018; Ma et al. 2023). It is thus important to obtain an accurate DSR over rugged terrains (DSR_rugged).

However, DSR_rugged cannot be easily measured on site due to the scarcity of measuring stations in mountainous regions, and almost no product can directly provide its value. Several methods have been proposed to estimate DSR_rugged. According to previous studies (Dubayah and Loechel 1997; Chen et al. 2013; Wang et al. 2014; Yan et al. 2018), DSR_rugged consists of three radiative components: direct solar radiation (D_dir-rugged), diffuse radiation (D_dif-rugged), and reflected radiation from the surrounding terrain (D_ref-rugged). Hence, DSR_rugged is usually calculated by summing these three components, which are calculated separately. These three radiative components (D_dir-rugged, D_dif-rugged, and D_ref-rugged) can either be simulated by using radiative transfer models (RTMs) and considering the influence of the terrain on the interactions between the atmosphere and the land surface (e.g. absorption, diffusion, and emission), or they can be corrected through parameterization by using the direct correction formula (DCF). The discrete anisotropic radiative transfer (DART) model (Gastellu-Etchegorry, Martin, and Gascon 2010; Wang et .al 2020a) and LargE-Scale remote sensing data and image Simulation Framework (LESS) (Qi et al. 2019; Chu et al. 2021; Ma et al. 2022a; Yan et al. 2020) are the most popular incarnations of the RTM for simulations. The two models differ in terms of their radiative components and efficiency, but are not recommended for simulating large-scale scenes (especially the DART model) (Malbéteau et al. 2017; Yan et al. 2020). The DCF is used more often to obtain DSR_rugged due to its simplicity and ease of implementation. D_dir-rugged is corrected based on the corresponding direct solar radiation over a flat surface (D_dir-flat) by using the solar illumination angle (θ_i) of the inclined and unshielded surface considered. D_dif-rugged is corrected based on the corresponding diffuse solar radiation over a flat surface (D_dif-flat) by using sky view factors (V_d), and D_ref-rugged is directly calculated by using the surface broadband albedo (α) and terrain view factors (V_c) (Giles 2001; Dozier, Bruno, and Downey 1981; Li et al. 2012; Olyphant 1986; Wu et al. 2018). DCF has been successfully used to obtain the instantaneous DSR_rugged (DSR_ins-rugged) from satellite data (e.g. the Landsat Thematic Mapper < TM > and the Moderate-Resolution Imaging Spectroradiometer < MODIS>) in areas with a rugged terrain in China, such as the Dayekou watershed of the Heihe River Basin and the Chengde area (Wang et al. 2018) as well as the Tibetan Plateau (TP) (Zhang et al. 2015a; Wang et al. 2018). Yan et al. (2020) found that DSR_ins-rugged calculated based on the DCF was in a good agreement with ground measurements, and was comparable to and even better than that the results obtained from the LESS. However, the performance of the DCF is highly reliant on the accuracy of inputs to it, especially the atmospheric parameters and factors of the terrain (Zhang et al. 2015a). Moreover, it can be used only at the instantaneous scale, and cannot be used to calculate the daily DSR over a rugged terrain (DSR_daily-rugged), where this is needed for applications. Some researchers have used methods of temporal expansion that are commonly used on flat surfaces, including the direct averaging method and the sinusoidal model (Roupioz et al. 2016), which is built by assuming that variations in the DSR conform to the sine curve, to expand the calculated DSR_ins-rugged to obtain DSR_daily-rugged. However, Yan et al. (2018) noted that the uncertainty in the value of DSR_daily-rugged based on these methods can be as large as 60 Wm⁻². Therefore, a reasonable and effective method is needed to accurately expand the calculated DSR_ins-rugged to DSR_daily-rugged, especially when data on the former are limited.

To address the above issues, the relationship between DSR_ins-rugged and DSR_daily-rugged needs to first be thoroughly explored, and this requires a large number of comprehensive samples. However, ground measurements over rugged terrains are scarce, and it is difficult to comprehensively simulated such regions based on data from LESS or DART (Yan et al. 2020; Malbéteau et al. 2017; Wang et al. 2018). As proposed in one past study (Wang et al. 2020c), DSR_rugged can be obtained from DSR products by using the DCF method. In this study, we use global hourly data on solar radiation, including D_dir-flat and D_dif-flat at 1° from CERES4_SYN1deg_Ed4A (CERES4, for short) (Rutan et al. 2015) to generate the required samples (see Section 3.1). Note that hourly data were considered to be at the instantaneous scale in this study. We then analyze these samples to develop an empirical model, called the Daily Downward Shortwave Radiation Random Forest Model over Rugged Terrain (DSRMT), to estimate the DSR_daily-rugged on clear days by using the instantaneous D_dir-flat and D_dif-flat (D_ins-dir-flat and D_ins-dif-flat) as well as ancillary information as inputs. Note that rugged terrains with bright surface which was defined as the surface with albedo > 0.7 were not considered in this study. The remainder of this paper is organized as follows: Section 2 introduces the data and methods used here, Section 3 details the development of the proposed DSRMT, including sample preparation, the exploration of the relationship between DSR_ins-rugged and DSR_daily-rugged, and model building. The results of validation of the DSRMT are provided and analyzed in Section 4, and conclusions of this study are summarized in Section 5.

2. Materials

2.1. Data and pre-processing

2.1.1. In situ measurements

We collected DSR measurements for validation from six sites (Figure 1) located in mountainous areas, including one site (called GMN-PG3) in the Sierra Nevada in Spain and five sites in the Saihanba Forest Park of Chengde city in the Hebei province (Chengde Experimental Area, for short) of China (Zhou et al. 2018; Yan et al. 2018).

Figure 1. Spatial distribution of the six mountainous sites: (a) All sites over the globe. (b) One site, GMN-PG3, in the Sierra Nevada in Spain, and (c) The other five sites (sites 2–4 and site 6–7) in the Chengde Experimental Area. The legend in (a) represents the classes of mountains as defined by the United Nations Environment Program World Conservation Monitoring Center (UNEP-WCMC), and the legends in (b) and (c) show the elevations obtained from the Shuttle Radar Topography Mission (SRTM).

Table 1 gives detailed information on these sites. Data on all six sites were obtained from two observation networks/programs. The Guadalfeo Monitoring Network (GMN; https://doi.org/10.1594/PANGAEA.895236) was built in Sierra Nevada to study the dynamics of snow and hydrological processes in this mountainous area (Polo et al. 2019). It contained four sites, but only a few of them could provide radiative measurement (Aguilar, Pimentel, and Polo 2021). By considering data from matching years, only hourly DSR measurements from 2009, 2011, and 2018 at GMN-PG3 (Figure 1b) were used in this study. The other five sites were located in the Chengde Experimental Area (Figure 1c), which was built to investigate variations in radiation over mountainous areas. The DSR was measured at these sites by using the EKO radiation table or the CNR4 radiation table at one-minute intervals (Yan et al. 2018; Yan et al. 2020; Zhou et al. 2018). Following quality control, the DSR measurements at sites 2, 3, 4, 6, and 7 from 2018 and 2019 were used in this study.

Table 1. Detailed information of the six sites used in this study.

Site Name	Latitude (°)	Longitude(°)	Elevation (m)	S * (°)	A * (°)
site2	42.397	117.398	1848.427	19	269
site3	42.393	117.397	1852.700	26	196
site4	42.393	117.395	1810.066	22	285
site6	42.396	117.390	1756.811	29	189
site7	42.387	117.400	1838.169	30	138
GMN-PG3	36.872	−3.236	1332	20.81	15.26

*S and A represent slope and aspect.

We first chose the quality-controlled in-situ measurements under clear-sky conditions. Only measurements with a corresponding clearness index (CI, see Section 2.1.3) higher than 0.7 were used from GMN-PG3. Only measurements that were flagged as having been made under a clear sky were used from the Chengde Experimental Area. All clear-sky DSR measurements were then transferred to the local time format. The measurements from the Chengde Experimental Area were aggregated into their hourly means as long as 50 samples were available within the given hour, and all hourly measurements, including those from GMN-PG3, were aggregated into their corresponding daily means without any missing hourly data.

According to Yan et al. (2018), only radiative measurements made parallel to an inclined surface, as shown in Figure 2b, can be used to characterize the actual DSR_rugged. Thus, we needed to correct measurements from GMN-PN3 for which the instrument had been deployed parallel to the horizontal surface, as shown in Figure 2a. We used the simple method of correction developed by Wang et al. (2020c) based on the solar zenith angle (θ₀) and θ_i (Figure 2c): (1) $DS{R}_{ins-rugged}={\displaystyle \frac{DS{R}_{ins-flat}\times u_s}{u_0}}$(1) (1a) $u_s=cos{\theta }_i=u_{0}cosS+sin{\theta }_{0}sinScos({\phi }_0-A)$(1a) (1b) $u_0=cos{\theta }_0$(1b) where DSR_ins-flat (Wm⁻²) represents instantaneous DSR values over a flat surface, θ_i (the angle with respect to the normal of the terrain) was calculated based on the solar azimuth angle φ₀ (from the north), θ₀ (angle with respect to the vertical direction), S (angle with respect to the horizontal) and A (angle with respect to the north); φ₀ and θ₀ were calculated by using the Solar Position Calculator (Mikofski 2022).

Figure 2. The different ways to deploy the measuring instrument on a sloped surface: (a) parallel to the horizontal surface, (b) parallel to the inclined surface, and (c) diagram of the corresponding θ_i and solar zenith angle θ₀.

Finally, we obtained 175 daily samples from the Chengde Experimental Area in 2018 and 2019, and 119 daily samples from GMN-PG3 in 2009, 2011, and 2018.

2.1.2. Remotely sensed data

Several kinds of remotely sensed data were used in this study, including D_ins-dir-flat, D_ins-dif-flat, and DSR_ins-flat from CERES4 and MCD18, the albedo from the Global Land Surface Satellite (GLASS), and the DEM from the SRTM.

2.1.2.1. D_ins-dir-flat, D_ins-dif-flat, and DSR_ins-flat from CERES4

According to previous studies (Wang et al. 2020b; Zhang et al. 2015b), the DSR from the SYN1deg series of products of CERES, published by the National Aeronautics and Space Administration (NASA; https://ceres.larc.nasa.gov), is one of the best global radiative products. Its values of D_ins-dir-flat, D_ins-dif-flat, and DSR_ins-flat were calculated from an RTM by considering the parameters of clouds and aerosols from MODIS, the geostationary satellite (GEO), and the Global Modeling and Assimilation Office (GMAO) (Rutan et al. 2015) at a spatial size of 1° and at hourly intervals from 2000 onward as inputs. The values of D_ins-dir-flat, D_ins-dif-flat, and DSR_ins-flat were taken from CERES4, from 8:30–17:30hrs local time in 2015, to construct the samples over mountainous areas in light of its outstanding performance, fine temporal resolution, and global coverage. The matching data from in-situ measurements in 2009, 2011, 2018, and 2019 were used to validate the model. The details of samples preparation are provided in Section 3.1.

2.1.2.2. D_ins-dir-flat, D_ins-dif-flat, and DSR_ins-flat from MCD18

MCD18 (https://lpdaac.usgs.gov) was released by NASA, was generated by using the look-up-table (LUT) method based on the top-of-atmosphere reflectivity acquired from MODIS and auxiliary data, including geographical locations, surface albedo from the MCD43A3, vapor pressure from the Modern-Era Retrospective Analysis for Research and Applications Version 2 (MERRA2), DEM from the Global 30 Arc-Second Elevation Dataset (GTOPO30), and data on the surface reflectance climatology. In contrast to CERES4, MCD18 provided the data on D_ins-dir-flat, D_ins-dif-flat, and DSR_ins-flat at a spatial size of 5 km at the satellite overpass time. The creators of MCD18 used the direct averaging method to obtain the DSR daily values, and the results indicated that the daily DSR from MCD18 was superior to that obtained from GLASS but was slightly inferior to that obtained from CERES4. We matched the data on D_ins-dir-flat, D_ins-dif-flat, and DSR_ins-flat obtained from MCD18 with those acquired through in-situ measurements in 2009, 2011, 2018, and 2019 to assess the feasibility of the model.

2.1.2.3. GLASS broadband albedo

The GLASS broadband albedo (α, http://glass-product.bnu.edu.cn) product (Liu et al. 2013a; Liu et al. 2013b; Qu et al. 2014) was directly estimated by using satellite reflectance data from MODIS and the Advanced Very High Resolution Radiometer (AVHRR) through the LUT method. The latter was established based on data measured at the sites and simulations of the model under different conditions in terms of region, type of land cover, and season at a global scale. The GLASS broadband albedo had a high accuracy than MCD43 (Liu et al. 2013b; He et al. 2013; He, Liang, and Song 2014), especially in certain areas, such as TP (An et al. 2021b). We used the α at a spatial resolution of 250 m and a temporal resolution of four days from 2015 for modeling as same as the time of CERES4. The albedo at the six mountainous sites over five years (2009, 2011, 2015, 2018 and 2019) were used for model validation in this study. The days on which the albedo was unavailable were assigned the albedo on the nearest day to it (within four days) for which this value was available by assuming that the variations in the albedo were minor in this period. For the sake of convenience, the albedo values of the samples were assumed to be equal to the average value of the adjacent terrain as Ma et al. (2023) for a 3 × 3 window.

2.1.2.4. SRTM DEM

DEM data from the SRTM (https://lpdaac.usgs.gov/products/srtmgl1v003/) were generated by NASA, the National Geospatial Intelligence Agency (NGA), and the German and the Italian space agencies by using radar interferometry. They cover about 80% of the global land surface between 60° N and 56° S in the system of coordinates of the World Geodetic System 1984 (WGS84) at a spatial resolution of 1 rad/s (∼30 m) (Rodríguez, Morris, and Belz 2006). Tang et al. (2021) have claimed that it is the best DEM product available. We used DEM from the SRTM to calculate factors of the terrain, including S, A, V_d, and V_c, and to determine whether the target pixel was blocked by neighboring pixels (see Equation (5(2a) $S=arctan\sqrt{f_x^2+f_y^2}$(2a) )). Figure 3 shows an example. Information on the terrain of adjacent pixels in a 3 × 3 window is usually applied to obtain the topographic factors of the target pixel Z₅ (Wu et al. 2018; Zhou and Liu 2004).

Figure 3. Diagram of calculation of factors of the terrain of the target pixel Z₅. Z represents the elevation of the pixel.

S, A, V_d, and V_c were calculated as follows (Wu et al. 2018): (2) $f_x={\displaystyle \frac{z_8-z_2}{2\times d},f_y={\displaystyle \frac{z_{6-}{z}_4}{2\times d}}}$(2) (2a) $S=arctan\sqrt{f_x^2+f_y^2}$(2a) (2b) $A={270}^{\circ }+\arctan \left({\displaystyle \frac{f_y}{f_x}}\right)-{90}^{\circ }{\displaystyle \frac{f_y}{|f_x|}}$(2b) (2c) $V_d={\displaystyle \frac{(1+cosS)}{2}}$(2c) (2d) $V_c=1-V_d$(2d) where d represents the resolution of the DEM data (30 m in this study).

2.1.3. Other parameters

Clouds are a key factor influencing the process of atmospheric radiative transfer, especially as they determine the ratios of D_dir and D_dif in the DSR. Clouds are more likely to be generated in mountainous areas than in flat areas (Yan et al. 2018). Hence, their influence was considered to be non-negligible in estimating DSR_rugged. The CI is defined as the ratio of the DSR to extraterrestrial radiation (R_se, Wm⁻²): (3) $CI={\displaystyle \frac{DS{R}_{ins-flat}}{R_{se}}}$(3) (3a) $R_{se}={\displaystyle \frac{1440G_{sc}{d}_r}{\pi }(w_s\sin (\phi )\sin (\delta )+\cos (\delta )\sin (w_s))}$(3a) (3b) $d_r=1+0.033\cos \left({\displaystyle \frac{2\pi doy}{365}}\right)$(3b)

(3c) $\delta =0.409\sin \left({\displaystyle \frac{2\pi doy}{365}-1.39}\right)$(3c) (3d) $w_s=\arccos (-\tan (\phi )\tan \delta )$(3d) (3e) $\mathrm{rad}={\displaystyle \frac{\pi }{180}\ast (decimaldeg)}$(3e) where G_sc is the solar constant (0.0820 MJm⁻²min⁻¹), d_r is the inverse relative distance from the Earth to the Sun (unitless), w_s is the sunset hour angle (rad), φ is the latitude (rad), δ is the solar declination (rad), and doy is the day of the year (Irmak et al. 2003). The CI is used to represent the transmittance of the sky, and many studies have used it to classify conditions of the sky. We define the sky as clear in this study when CI > 0.7, as suggested by Ruiz-Arias et al. (2009).

Moreover, the map of mountains released by the UNEP-WCMC (http://www.unep-wcmc.org) was used to determine the mountain regions over the globe as a reference. This map divides the world’s mountains into six categories based on S, the elevation, and local range of elevation according to Table 2 (Kapos 2000; Kapos et al. 2000). This standard of classification can be used to identify large mountain ranges that can be mapped and avoid the inclusion of mid-elevation plateaus. This map has been used by the UNEP-WCMC to determine the world's mountain forests, and has been widely applied in related research (Ehrlich, Melchiorri, and Capitani 2021; An et al. 2021; Umuhoza et al. 2021). We chose 17 typical mountain regions across the world for sample collection based on this map.

Table 2. Classes of Mountains defined in UNEP-WCMC.

Class (elevation in meters)	Additional Criterion	% of the Earth's Surface
>4,500		1.2
3,500 ∼ 4,499		1.8
2,500 ∼ 3,499		4.7
1,500 ∼ 2,499	S > 2°	3.6
1,000 ∼ 1,499	S > 5° or LER > 300m	4.2
300 ∼ 999	LER > 300m	8.8
Total		24.3

LER (Local Elevation Range) represents the elevation range in an area assessing over a seven-mile radius.

2.2. Method

The flowchart of this study is given in Figure 4. First, comprehensive samples for DSR_ins-rugged and the corresponding DSR_daily-rugged were prepared by using data on D_ins-dir-flat, D_ins-dif-flat, and DSR_ins-flat from CERES4 that were randomly chosen from 17 global mountainous regions selected from the UNEP-WCMC map. The details are provided in Section 3.1. Second, we explored the relationship between DSR_ins-rugged and DSR_daily-rugged based on the samples under various conditions of the factors of the terrain (S, A, V_d, and V_c), location, and time. Third, these results were used to build the DSRMT. Finally, the performance of the DSRMT was evaluated on the validation samples and in-situ measurements. Its applicability was assessed by using the DSR-related variables from MCD18 as input.

Figure 4. The flowchart of this study.

2.2.1. Introduction to DCF

As mentioned above, the DSR_ins-rugged samples over mountainous areas were calculated from CERES4 by using the DCF method as below: (4) $\begin{array}{rl}& DS{R}_{ins-rugged}=D_{ins-dir-rugged}+D_{ins-dif-rugged}+D_{ins-ref-rugged}\\ & =D_{ins-dir-flat}\times \left({\displaystyle \frac{u_s}{u_0}}\right)\times \mathrm{\Phi }+D_{ins-dif-flat}\times V_d+V_c\times \alpha \times DSR\end{array}$(4) where D_{ins-dir-rugged} (Wm⁻²) and D_{ins-dif-rugged} (Wm⁻²) represent the instantaneous direct and diffuse solar radiations over a rugged terrain for a given sample, respectively, and D_{ins-ref-rugged} (Wm⁻²) is the reflected solar radiation from the adjacent region (3 × 3 window). The average DSR and α (unitless) were assumed to be the same as those of this sample, u_s and u₀ represent the cosine values of θ_i and θ₀, respectively, (see Equation (1(1) $DS{R}_{ins-rugged}={\displaystyle \frac{DS{R}_{ins-flat}\times u_s}{u_0}}$(1) )), Φ is a binary shadow function (unitless) indicating whether the given pixel was sheltered by a shadow (in which case its values was zero, and was one otherwise), V_d (unitless) and V_c (unitless) were calculated from the DEM of the SRTM according to Equation (2(1a) $u_s=cos{\theta }_i=u_{0}cosS+sin{\theta }_{0}sinScos({\phi }_0-A)$(1a) ), and α was obtained from the GLASS broadband albedo product.

With respect to Φ, a shadow in a pixel of an image of a rugged terrain is commonly defined as a self-generated shadow (u_s< 0), or that formed by the sheltering of the pixel by surrounding terrain, as shown in Figure 5. We thus needed to determine whether a sheltered pixel was blocked by adjacent pixels.

Figure 5. Schematic diagram of determining whether a pixel was blocked by neighboring pixels.

According to Equation (5(2a) $S=arctan\sqrt{f_x^2+f_y^2}$(2a) ) (Li et al. 2012), if the elevation of pixels adjacent to the given one (Z_d) is higher than the theoretical value of the target pixel (h) when it is sheltered, then the target point is determined to be sheltered. Thus, Φ = 0, and vice versa: (5) $Z_d\ge h=Z_0+d\times tan({90}^{\circ }-{\theta }_0)$(5) where Z₀ is the elevation of the target pixel.

2.2.2. Random forest (RF)

We built the DSRMT by using the random forest (RF) method. Developed by Breiman (1996), the RF is one of the most popular non-linear machine learning methods because of its high accuracy, ease of implementation, low computational cost, and high speed, especially in predictive analysis (Babar et al. 2020). RF regression (shown in Figure 6) is performed as follows (Hou et al. 2020): (1) Extract data either randomly or at a fixed ratio to construct a regression tree. (2) Set the number of nodes, and randomly select a subset of the input variables on each node. Split each decision tree according to the minimum Gini coefficient until it meets the threshold set by the user. (3) Convert the generated regression trees into a random forest, and calculate the final prediction by using the mean of the values predicted by multiple trees. Four hyper-parameters—n_estimaators, max_depth, min_samples_split, and min_samples_leaf—need to be tuned to obtain the optimal model. To determine the optimal hyper-parameters in the RF model, we followed Li et al. (2022) in applying a circular approach by using the minimum root mean-squared error (RMSE) between the training and the testing procedures. We implemented the RF by using the Scikit-learn toolbox (Pedregosa et al. 2012) on the Python platform in a Microsoft Windows 10 system with 32 GB of memory.

Figure 6. The structure of random forest regression.

2.2.3. Validation of accuracy

Four statistical measures—the RMSE, the mean absolute error (MAE), the square of the correlation coefficient (R²), and the relative RMSE (rRSME) – were used to assess the accuracy of the model. (6) $\mathrm{MAE}={\displaystyle \frac{{\sum }_{i=1}^N|(Y_i-X_i)|}{N}}$(6) (7) $\mathrm{RMSE}=\sqrt{{\displaystyle \frac{1}{N}\sum _{i=1}^N{(X_i-Y_i)}^2}}$(7) (8) ${\mathrm{R}}^2=1-{\displaystyle \frac{{\sum }_{i=1}^N{(X_i-Y_i)}^2}{{\sum }_{i=1}^N{(X_i-\overline{X})}^2}}$(8) (9) $\mathrm{rRMSE}={\displaystyle \frac{RMSE}{\overline{X}}}$(9) where Y_i and X_i are the estimation and the measurement of the ith group of samples, and N represents the number of samples. R² was also used to represent the correlation between DSR_ins-rugged and DSR_daily-rugged.

3. Development of DSRMT model

3.1. Sample preparation

As described above, we randomly chose 17 study areas in typical mountainous regions across the world within 60° S−60° N from which to collect the samples, as shown in Figure 7a. These 17 regions were mainly covered by forests, deserts, grasslands, barren land, shrubs, and snow, and their elevations ranged from 1000 to 10000 m, except the 17th region (Figure 7a17). This was located in a desert with a relatively flat surface. Seven of the regions were located in typical mountainous areas of Asia: the Himalaya Mountains (Figure 7a7), Tianshan Mountains (Figure 7a8), Altai Mountains (Figure 7a9), Qilian Mountains (Figure 7a10), Jablonov Mountains (Figure 7a11), Qinling Mountains (Figure 7a12), and the Arabian Plateau (Figure 7a16). Five regions (Figure 7a2−a6) were located in the Andes Mountains of South America, two (Figure 7a13−a14) were in the Great Dividing Range of Australia, and three study regions were located in the Alps in Europe (Figure 7a15), the Rocky Mountains in North America (Figure 7a1), and the Sahara Desert in Africa (Figure 7a17). Note that the mountainous areas at high latitudes (>60° S/N) were not considered because of the limited coverage of the SRTM DEM, and only five study regions were in the Southern Hemisphere (Figure 7a4−a6 and a13−a14) because of its smaller land surface.

Figure 7. Spatial distributions of (a) the 17 study regions within typical mountainous areas (see Table 2) across the world, and the training and validation samples acquired from them (a1)–(a17). The legend in (a) represents the classes of mountains and that in (b) represents the land cover type as defined by the International Geosphere–Biosphere Program (IGBP) based on MCD12 C1 in 2015.

The pixels for each study area obtained from CERES4 were randomly sampled as shown in Figure 7a1−a17. To ensure the comprehensiveness of the samples, the pixels were sampled by considering the latitude, doy, and four factors of the terrain (S, A, V_d, and elevation) were calculated by using SRTM DEM. We also made sure that they did not contain shadows (see Equation (5(2a) $S=arctan\sqrt{f_x^2+f_y^2}$(2a) )). The slopes of these samples were solo slopes because of the fine resolution of DEM (Xian et al. 2023; Wen et al. 2018). Following this, values of DSR_ins-rugged were calculated from those of D_ins-dir-flat, D_ins-dif-flat, and DSR_ins-flat from 8:30−17:30hrs local time in 2015 from CERES4. The corresponding DSR_daily-rugged at these pixels on each day were obtained by aggregating all values of DSR_ins-rugged on the day. To satisfy the requirements of this study, only values of DSR_daily-rugged under a clear sky (CI > 0.7) without bright surface (the corresponding GLASS α>0.7) were used as samples. Finally, all samples in each study region were divided into training (black dots) and validation (red dots) sets of nearly the same size, except for the data from four study regions (Figure 7a7, a12, a13, and a16) that were located in the desert and the TP. Samples from these regions were used only for validation.

The statistics on all DSR_daily-rugged samples based on various factors are presented in Figure 8. The samples were uniformly distributed in terms of A, S, and doy (Figure 8a, b, and e). Most of the samples were from the Northern Hemisphere, while the smallest numbers of samples were from within 40°–90° S (Figure 8c) and high-elevation zones (Figure 8f). The number of samples gradually increased with V_d (Figure 8d), possibly because of the more frequent occurrence of clouds over higher slopes with smaller values of V_d (Yan et al. 2018).

Figure 8. The statistics on all training (orange bar) and validation samples (yellow bar) of DSR_daily-rugged based on various factors: (a) A (0−360°) with an interval of 45°, (b) S (0−90°) with an interval of 10°, (c) φ (latitude) within 40° S−90° S, 20° S−40° S, 0°−20° S, 0°−20° N, 20°−40° N, and 40°−90° N, (d) V_d (0.5∼1) with an interval of 0.1, (e) doy (1∼365) with an interval of 60, and (f) elevation as determined by the six mountain classes defined in Table 2.

In total, we prepared 698,309 samples on DSR_daily-rugged, including 382,220 for training and 316,089 for validation, across various elevations, types of land cover, and factors of the terrain to develop our model. The corresponding 1,802,619 samples of DSR_ins-rugged were obtained from 8:30−17:30hrs local time, and were applied to all training samples of DSR_daily-rugged.

3.2. Relationship between DSR_ins-rugged and DSR_daily-rugged

The relationship between DSR_ins-rugged and DSR_daily-rugged at different times from 8:30–17:30hrs was preliminary explored based on all the training samples. The results are shown in Figure 9. The scatter plots indicate that DSR_ins-rugged was more linearly related to DSR_daily-rugged closer to noon and at 10:30–14:30hrs (Figure 9c–g), while the relationship between them was poor at all other times (e.g. 8:30, 9:30, 15:30, 16:30 and 17:30hrs, Figure 9 a, b, h–j). These results coincide with those of Yan et al. (2018) who claimed that the influence of the terrain on the DSR decreases close to noon, when the solar altitude increases. Moreover, smaller values (<100 Wm⁻²) of DSR_ins-rugged imply that diffuse and reflected solar radiation from the adjacent terrain accounted for a relatively large portion of the overall radiation. Hence, the high density (yellow color) of the low values, even at noon, indicates that the influence of the adjacent terrain on the DSR should be considered during the day.

Figure 9. Scatter plots between DSR_ins-rugged and DSR_daily-rugged at (a) 8:30, (b) 9:30, (c) 10:30, (d) 11:30, (e) 12:30, (f) 13:30, (g) 14:30, (h) 15:30, (i) 16:30, and (j)17:30hrs, local time.

By considering the results in Figure 9 and the times at which most polar-orbiting satellites passed over the study regions, we chose the data obtained at 10:30, 11:30, 12:30, 13:30, and 14:30hrs, when the relation between DSR_ins-rugged and DSR_daily-rugged was relatively sound, for further investigation. The variations in the relationship between DSR_ins-rugged and DSR_daily-rugged with the values of S and A at the five times are shown in Figure 10a−e, respectively, Figure 10f gives the corresponding sample size in each case.

Figure 10. The correlation between DSR_ins-rugged and DSR_daily-rugged as represented by R² under different combinations of S and A at five times: (a) 10:30, (b) 11:30, (c) 12:30, (d) 13:30, and (e) 14:30hrs. (f) The corresponding sample size at the five times, and the different colors indicate the numbers of different samples.

The above shows that DSR_ins-rugged and DSR_daily-rugged were relative closely correlated, with values of R² ranging from 0.6–1 when S was smaller than 50° at the five times. When S was larger than 50°, the values of R² decreased with the increase in S and A, especially when A > 180° before noon and A < 180° after noon. The lowest values of R² (<0.3) were obtained when A∈(160°−320°) in the morning (10:30 and 11:30hrs, Figure 10a and b) and A∈(40°−160°) in the afternoon (13:30 and 14:30hrs, Figure 10d and e). The influence of S and A on R² declined close to noon (11:30 and 12:30hrs, Figure 10b and c) compared with that at the other times considered. This coincides with the results in Figure 9, but the values of R² were still lower than 0.4 when S was larger than 70° and A∈(80°−280°) (Figure 10c). Combining these results with values of φ₀ shows that low values of R² were obtained in shady areas and areas with a steep slope. The terrain had a minor influence on DSR_ins-rugged or DSR_daily-rugged when S was smaller than 10° at all times, for all values of A. To confirm this, we examined the relation between DSR_{daily−rugged} and the corresponding daily DSR over a flat surface (DSR_daily−flat) obtained from CERES4 with different slopes, and the results given in Figure 11 demonstrated this point.

Figure 11. The relationship between DSR_{daily−rugged} and DSR_daily−flat as represented by the normalized RMSE, normalized MAE obtained by normalizing the original RMSE and the MAE through min-max normalization method for clearer presentation, and R² with different slopes. The closer R² to zero, and close the normalized RMSE and the MAE were to one, the worse was the correlation, and vice versa.

We also examined the influence of α, which was related to the value of D_ref-rugged of the adjacent terrain, on the relationship between DSR_ins-rugged and DSR_daily-rugged. Table 3 lists the results at the five times (10:30, 11:30, 12:30, 13:30, and 14:30hrs) as α was varied from 0 to 0.7. It shows that the RMSE generally increased and varied greatly as α increased at all five times. It increased to 10 Wm⁻² before 12:30hrs and decreased to 5 Wm⁻² afterward. Hence, the influence of α needed to be considered in the relation between DSR_ins-rugged to DSR_daily-rugged, especially when D_ref-rugged was very large.

Table 3. The relationship between DSR_ins-rugged and DSR_daily-rugged as represented by the RMSE and R² for different albedo values at an interval of 0.2 from 0−0.7 at five times (10:30, 11:30, 12:30, 13:30, and 14:30hrs).

albedo	10:30	11:30	12:30	13:30	14:30
RMSE (Wm^−2)
0 − 0.2	45.16	39.13	40.02	46.89	56.67
0.2 − 0.4	43.73	37.33	38.64	44.97	53.83
0.4 − 0.6	43.40	41.53	43.80	47.97	56.80
0.6 − 0.7	51.08	47.31	47.42	50.16	58.68
R^2
0 − 0.2	0.70	0.78	0.78	0.70	0.56
0.2 − 0.4	0.70	0.79	0.78	0.70	0.53
0.4 − 0.6	0.75	0.80	0.76	0.72	0.60
0.6 − 0.7	0.64	0.70	0.70	0.67	0.57

In conclusion, DSR_ins-rugged was adequately but varyingly related to DSR_daily-rugged at the five times considered under a clear sky when S > 10°, and their relation was significantly influenced by the values of S and A as well as the incident solar angle and α.

3.3. Building DSRMT model

The results in Section 3.2 shows that the DSR_daily-rugged could be linked to DSR_ins-rugged at different times from 10:30–14:30hrs by combining several factors of the terrain with ancillary information. After multiple experiments, DSR_daily-rugged was estimated as: (10) $DS{R}_{daily-rugged}=f(D_{ins-dir-rugged},D_{ins-dif-flat},DS{R}_{ins-flat},Te{r}_{i=1,2,3,4},\alpha ,CI,d_r,\phi )$(10) where D_{ins-dir-rugged}, D_ins-dif-flat, and DSR_ins-flat were taken from CERES4, and D_{ins-dir-rugged} was calculated from Equation (4(2) $f_x={\displaystyle \frac{z_8-z_2}{2\times d},f_y={\displaystyle \frac{z_{6-}{z}_4}{2\times d}}}$(2) ). Ter_{i = 1,2,3,4} represents the factors of the terrain (S, A, V_d, and V_c) obtained from the SRTM DEM by Equation (2(1a) $u_s=cos{\theta }_i=u_{0}cosS+sin{\theta }_{0}sinScos({\phi }_0-A)$(1a) ), α is the surface daily broadband albedo obtained from GLASS, CI and d_r were calculated by using Equations (3) and (3b), φ is the latitude of the sample, and u_s and u₀ were calculated from Equations (1a) and (1b).

Five models of the DSRMT were trained on data from each of the five times (10:30, 11:30, 12:30, 13:30, and 14:30hrs) by using the RF method for considering the overpass time of most polar-orbiting satellites and the temporal resolutios of their products. They were called mod1030, mod1130, mod1230, mod1330, and mod1430, respectively. After several experiments, the five optimal models were identified and their settings are given in Table 4. The accuracies of training of the five DSRMT models were very similar to one another, with values of R² around 0.96, an RMSE ranging from 11.22–12.06 Wm⁻², and values of MAE close to zero.

Table 4. Hyper-parameter settings used to identify five optimal DSRMT models. The three values in brackets for each hyper-parameter of every model represent the start, interval, and end values, respectively, and the values in parentheses represent the value of the confirming hyper-parameter.

	Hyper-Parameters
	n-estimators	max depth	min samples split	min samples leaf
mod1030	[30,10,70] (40)	[10,5,40] (20)	[2, 2,10] (2)	[2,2, 10] (4)
mod1130	[30,10,70] (60)	[10,5,40] (15)	[2, 2,10] (2)	[2,2, 10] (8)
mod1230	[30,10,70] (60)	[10,5,40] (15)	[2, 2,10] (2)	[2,2, 10] (8)
mod1330	[30,10,70] (40)	[10,5,40] (15)	[2, 2,10] (2)	[2,2, 10] (8)
mod1430	[30,10,70] (40)	[10,5,40] (15)	[2, 2,10] (2)	[2,2, 10] (8)

If more than one estimate of DSR_daily-rugged was obtained on a day by DSRMT, their weighted average was calculated as the final daily mean $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ as shown in Equation (11(3b) $d_r=1+0.033\cos \left({\displaystyle \frac{2\pi doy}{365}}\right)$(3b) ). The weights were determined by the number of DSR_daily-rugged estimates and the corresponding value of u_s of the DSR_ins-rugged obtained by φ₀, θ₀, S, and A (Equation (1(1) $DS{R}_{ins-rugged}={\displaystyle \frac{DS{R}_{ins-flat}\times u_s}{u_0}}$(1) )): (11) $\overline{DS{R}_{daily-rugged}}={\displaystyle \frac{{\sum }_{i=1}^n(u_{s,i}\times DS{R}_{daily-rugged,i})}{{\sum }_{i=1}^n{u}_{s,i}}}$(11) where n represents the number of DSR_daily-rugged estimates of the DSRMT, and ranges from two to five.

The overall accuracy of the value of $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ calculated based on the outputs of the five DSRMT models (n = 5) against all training samples is presented in Figure 12. The accuracy of $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ increased significantly by about 2 Wm⁻² in terms of the RMSE (RMSE = 9.72 Wm⁻²) compared with the accuracy of training of each of the five DSRMT models.

Figure 12. The overall accuracy of the $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$, obtained by the weighted average of the outputs of the five DSRMT models based on all training samples.

Therefore, the DSRMT could thus estimate DSR_daily-rugged at any time from 10:30hrs to 14:30hrs on the same day based on factors of the terrain and other easily obtainable information, as long as the corresponding D_{ins-dir-rugged}, D_ins-dif-flat, and DSR_ins-flat were available. Note that DSR_daily-rugged could be estimated from one DSRMT model, and the DSRMT model could be trained at any moment during daytime (10:30hrs∼14:30hrs were suggested) according to the available training samples besides the five moments discussed in this study, so the application of DSRMT is very flexible. The weighted average of u_s needed to be applied when more than two estimates of DSR_daily-rugged were obtained by the DSRMT. However, as analyzed above, the DSRMT just works on a rugged terrain with a large S (S > 10°) and a value of α not larger than 0.7 under a clear sky.

4. Results and analysis

The performance of the DSRMT was assessed based on validation samples calculated from CERES4 and in-situ measurements (see Section 2.1.1), and then the influence of several factors on its outcomes were analyzed.

4.1. Validating the accuracy of DSRMT

4.1.1. Samples for validation from CERES4

The accuracy of each DSRMT model at the five times (mod1030, mod1130, mod1230, mod1330, and mod1430) was validated by using all validation samples of DSR_daily-rugged from CERES4, and the results are shown in Table 5. Note that all inputs were available at the five times, such that the five models were assessed on the same data. Overall, the RMSE values of the five models were in the range 17−23 Wm⁻², and their R² values were lower than 0.92. mod1230 had the best accuracy. They recorded the smallest values of the RMSE and MAE of 17.41 and 13.10 Wm⁻², respectively, and a highest R² of 0.92. mod1130 was the second best, with slightly larger values of the RMSE (18.17 Wm⁻²) and MAE (13.70 Wm⁻²). The accuracy of the DSRMT models worsened close to 10:30hrs and 14:30hrs, and mod1430 had the worst accuracy, with an RMSE of 22.67 Wm⁻², an MAE of 17.31 Wm⁻², and a value of R² of 0.86. However, all five models performed reasonably well because of their comparable training and validation accuracy.

Table 5. The accuracy of validation of the five DSRMT models using validation samples from CERES4.

Model	Validation accuracy			No. of samples
	RMSE (Wm^−2)	MAE (Wm^−2)	R^2
mod1030	19.51	14.82	0.90	86,733
mod1130	18.17	13.70	0.91	86,733
mod1230	17.41	13.10	0.92	86,733
mod1330	19.82	15.10	0.89	86,733
mod1430	22.67	17.31	0.86	86,733

We then calculated $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$, the weighted average of the outputs of the five DSRMT models (n = 5 in Equation (11(3b) $d_r=1+0.033\cos \left({\displaystyle \frac{2\pi doy}{365}}\right)$(3b) )), and validated the results on all validation samples from CERES4. The results are shown in Figure 13a.

Figure 13. The accuracy of validation of $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ based on the weighted average of five DSR_daily-rugged estimates of the five DSRMT models by using (a) all validation samples from CERES4, (b) validation samples from CERES4 for regions that provided the training samples, and (c) validation samples from CERES4 for regions from which training samples were not chosen.

Figure 13a shows that the overall accuracy of the estimated $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ was satisfactory, with an RMSE of 15.45 Wm⁻², an MAE of 11.62 Wm⁻², and an R² of 0.94. This was much better than the results of any of the five DSRMT models, with a range of reduction of 1.96‒7.22 Wm⁻² in the RMSE and 1.48‒5.69 Wm⁻² in the MAE. Although this is slightly different from the accuracy of training shown in Figure 12 (with an increase of 5.73 Wm⁻² in the RMSE), it is acceptable due to differences in distribution between the samples used for training and validation. To better illustrate this, we split all validation samples into those from the regions that had also provided the training samples and those from regions from which training samples had not been chosen. Their accuracies of estimation of $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ are shown in Figure 13b and c, respectively. The accuracy in Figure 13b is slightly higher than that in Figure 13c, with RMSEs of 14.62 and 17.04 Wm⁻², MAEs of 10.87 and 13.18 Wm⁻², and R² values of 0.95 and 0.93, respectively. A significant overestimation can be observed at low values (<100 Wm⁻²) in Figure 13c. The performance of statistical models relies heavily on the samples used for training. Thus, the models possibly performed worse in situations where the training samples were not representative. Four of the regions from which the validation samples were randomly selected (Figure 7a7, a12, a13, and a16) featured relatively large differences in their values of DSR_rugged owing to their special natural and atmospheric conditions.

We examined and compared the accuracies of calculation of $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ based on the weighted average of u_s, but with different combinations of the five DSR_daily-rugged estimates. The results are displayed in Table 6. On the whole, the accuracy of estimates of $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ improved and became more robust when the number of DSR_daily-rugged estimates was increased from two to five, with RMSE values in the ranges of 16.04−19.78 Wm⁻², 15.61−17.6 Wm⁻², 15.39−16.25, and 15.41 Wm⁻². Similar accuracies of the estimated $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ were obtained even though the combinations were different. The accuracy was higher if the output of a more accurate model, such as mod1130 and mod1230, was added.

Table 6. The accuracy of validation of $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ calculated by using different combinations of the five DSR_{daily−rugged} estimates from the DSRMT models from 10:30–14:30hrs based on all validation samples from CERES4. The letters A, B, C, D, and E represents the DSR_{daily−rugged} estimates of mod1030, mod1130, mod1230, mod1330, and mod1430, respectively.

Combinations	RMSE (Wm^−2)	MAE (Wm^−2)	R^2	combinations	RMSE (Wm^−2)	MAE (Wm^−2)	R^2
A – E	15.45	11.62	0.94	A, B, E	15.89	11.96	0.94
A, B	17.28	13.16	0.92	A, C, D	15.61	11.71	0.94
A, C	16.04	12.06	0.93	A, C, E	15.70	11.79	0.94
A, D	16.32	12.23	0.93	A, D, E	16.51	12.50	0.93
A, E	16.88	12.68	0.93	B, C, D	15.99	12.04	0.93
B, C	16.27	12.24	0.93	B, C, E	16.10	12.16	0.93
B, D	16.49	12.41	0.93	B, D, E	16.85	12.83	0.93
B, E	16.96	12.84	0.93	C, D, E	17.60	13.46	0.92
C, D	17.17	13.04	0.92	A, B, C, D	15.40	11.55	0.94
C, E	17.68	13.49	0.92	A, B, C, E	15.39	11.55	0.94
D, E	19.78	15.20	0.90	A, B, D, E	15.74	11.85	0.94
A, B, C	15.87	11.98	0.93	A, C, D, E	15.86	11.98	0.94
A, B, D	15.75	11.83	0.94	B, C, D, E	16.25	12.32	0.93

The performance of all five DSRMT models was acceptable, but their results varied and improved close to noon. We recommend using the weighted average of u_s to obtain the final DSR_daily-rugged estimates if more than one output is available. Moreover, the accuracy can be improved if the results of the DSRMT at noon are used.

4.1.2. Using in-situ measurements

For a more objective validation, we also verified the performance of the DSRMT against the in-situ measurements in mountainous areas. Table 7 presents the accuracy of calculation of DSR_daily-rugged by the five DSRMT models at six mountainous sites against their in-situ measurements. The accuracies of the five DSRMT models were different but acceptable. However, their performance was generally worse than that obtained when using validation samples from CERES4, with RMSE values ranging from 24.9–29.22 Wm⁻² and the MAE ranging from 19.16–22.94 Wm⁻². We think that this result was obtained due to uncertainty in the DSR-related products of CERES4 used to train the DSRMT. Of the five models, mod1230 delivered the best performance with an RMSE of 24.90 Wm⁻² and an MAE of 19.16 Wm⁻². This coincided with the results obtained by using validation samples from CERES4. Surprisingly, the DSRMT performed much better on data from the afternoon than in the morning in this case, especially on data from 14:30hrs, with an RMSE of 24.34 Wm⁻² and an MAE of 20.18 Wm⁻².

Table 7. Results of validation of the five DSRMT models in terms of estimating DSR_daily-rugged in comparison with in-situ measurements from six mountainous areas.

Model	Validation accuracy
	RMSE (Wm^−2)	MAE (Wm^−2)	R^2
mod1030	26.39	19.28	0.92
mod1130	29.22	22.94	0.92
mod1230	24.90	19.16	0.92
mod1330	25.67	20.82	0.92
mod1430	24.34	20.18	0.91

We then obtained the accuracy of estimation of $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ by using the weighted average of the estimates of the five DSRMT models according to u_s, and compared them with ground measurements. The accuracies of $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ obtained from different combinations of models were different, with values of RMSE in the range of 21.63‒25.84 Wm⁻² and MAE in the range of 17.14‒19.72 Wm⁻². These results were better than or comparable to those of a single model (see Table 7). We also calculated the corresponding $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ at these sites by using the DCF method (${\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}}_{\mathrm{DCF}}$) and averaging all DSR_ins-rugged estimates from CERES4 during the day. Figure 14a and b present the results of validation of the optimal values of $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ and ${\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}}_{\mathrm{DCF}}$ against the ground measurements.

Figure 14. The overall accuracy of validation against in-situ measurements of (a) the optimal $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ estimates according to u_s based on the weighted average of the outputs of the DSRMT models, and (b) ${\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}}_{\mathrm{DCF}}$ estimates obtained by averaging all estimates DSR_ins-rugged by using the DCF method during the day.

Figure 14 shows that the accuracy of estimation of the optimal $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ based on the DSRMT (Figure 14a) was much higher than that of ${\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}}_{\mathrm{DCF}}$ when using the DCF method (Figure 14b), with an RMSE that was greater by ∼11 Wm⁻² and an MAE larger by ∼6 Wm⁻². However, the accuracy of the estimated DSR_daily-rugged obtained by any one of the five DSRMT models (see Table 7) was even better than that of ${\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}}_{\mathrm{DCF}}$. $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ was closer to the 1:1 line in the plots, while significant overestimations at large values (>300 Wm⁻²) and underestimations at small values (<150 Wm⁻²) were observed for ${\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}}_{\mathrm{DCF}}$, especially at the GMN-PG3 site, site 2, and site 4. Taking these sites as examples, we present the temporal variations in $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ and ${\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}}_{\mathrm{DCF}}$ in Figure 15 as well as their scatter plots against the in-situ measurements. The values of $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ (red line) and ${\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}}_{\mathrm{DCF}}$ (green line) generally captured the variations in the in-situ measurements (black dots) at these three sites on most days (Figure 15a1, b1, and c1), but the values of $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ were closer to the in-situ measurements. Specifically, the two estimates were generally smaller than the corresponding in-situ measurements at sites 2 and 4, especially in case of large values, and the values of ${\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}}_{\mathrm{DCF}}$ were smaller. The estimated $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ was slightly smaller than the ground measurements at site GMN-PG3 but ${\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}}_{\mathrm{DCF}}$ was larger than the in-situ measurements at this site. The scatter plots in Figure 15 show that the accuracy of estimation of $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ (Figure 15a2, b2, and c2) improved significantly at the three sites compared with that of ${\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}}_{\mathrm{DCF}}$ (Figure 15a3, b3, and c3), with reductions in the RMSE of 22.23, 13.86, and 10.01 Wm⁻², and those in MAE of 14.9, 14.8, and 9.4 Wm⁻², respectively.

Figure 15. Temporal variations in the values of $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ based on the DSRMT (red line), ${\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}}_{\mathrm{DCF}}$ based on the DCF method (green line), and in-situ measurements (black dot) at (a1) GMN-PG3, (b1) site 2, and (c1) site 4. (a2–a3), (b2–b3), and (c2–c3) show their corresponding scatter plots against the in-situ measurements. Note that the time given on the abscissa is not continuous in (a1)–(c1).

Furthermore, the accuracy of estimation of DSR_daily-rugged, which was expanded from the data on DSR_ins-rugged obtained from CERES4 at each of the five times (10:30–14:30hrs) by using sinusoidal model (Bisht et al. 2005) designed for flat surface but widely used, was validated against the same in-situ measurements for comparison. As shown in Table 8, the estimation accuracy yielded values of RMSE of 24.04 – 45.31Wm⁻², MAE of 17.30 – 39.24 Wm⁻², and R² of 0.84 – 0.94. It delivered the best performance at 12:30hrs, following which its accuracy declined. Combined with the data in Table 7, these results indicate that sinusoidal model can perform comparably to the proposed DSRMT model only at 12:30hrs, when topographic effects are thought to be the weakest, but its performance is not robust, and relies heavily on the location of the site and the availability of DSR_ins-rugged. Thus, the proposed DSRMT is superior to the sinusoidal model.

Table 8. The validation results of the sinusoidal method by using data on DSR_ins-rugged from CERES4 at each of the five times as the input against the same in-situ measurements.

Time	Validation accuracy
	RMSE (Wm^−2)	MAE (Wm^−2)	R^2
10:30	45.31	39.24	0.86
11:30	32.10	26.04	0.92
12:30	24.04	17.30	0.94
13:30	28.25	20.35	0.90
14:30	32.96	25.29	0.84

The above results demonstrate the superiority of the proposed DSRMT model in terms of estimating DSR_daily-rugged over the average means over the traditional DCF method and the sinusoidal method, particularly in terms of robustness. The accuracy of estimation of DSR_daily-rugged can be improved by using the weighted average of μ_s.

4.2. Factors influencing the performance of DSRMT

To assess the factors influencing the performance of the DSRMT and further better understand the mechanisms of the topographical effect on DSR_rugged, we performed a sensitivity analysis of all variables of the five DSRMT models on the Python platform, and results are shown in Figure 16.

Figure 16. The contributions of all variables of the five DSRMT models.

As shown in Figure 16a–k, the contribution of each variable to the DSRMT from 10:30–14:30hrs (represented using bars of different colors) was ranked from large to small. Different variables made similar contributions to the five DSRMT models. The five most influential variables were D_{ins-dir-rugged}, DSR_ins-flat, A, D_ins-dif-flat, and φ, but significant differences were observed in contributions of the first four variables (D_{ins-dir-rugged}, DSR_ins-flat, A, and D_ins-dif-flat) at the five times considered. Specifically, D_{ins-dir-rugged} (Figure 16a), which represents the terrain-calibrated values of D_ins-dir-flat, made the largest contribution in all five models (41.54% to mod1030, 51.07% to mod1130, 52.41% to mod1230, and 37.19% to mod1330). Its influence increased from 10:30hrs, reached its maximum value (over 50%) near noon (11:30 and 12:30hrs), and then declined until 14:30hrs. DSR_ins-flat (Figure 16b) came in the second place in all five models, contributing ∼20% to each except mod1030 (29%). Variations in its contributions were relatively stable, where this coincides with the result mentioned above: that the radiation reflected from the region adjacent to a given one could not be ignored, especially when θ₀ was large. The contribution of D_ins-dif-flat (Figure 16d) varied, as did that of D_{ins-dir-rugged}, but it accounted for only 5.97%−9.80% in the five models. This result was obtained mainly because we considered only clear-sky conditions in this study. Compared with the other factors of the terrain, A played an important role in the DSRMT (Figure 16c), especially when θ₀ was large, and contributed 18.15% to mod1430. S (Figure 16i) and V_d (Figure 16g) had similar and minor effects on the models.

In summary, D_{ins-dir-rugged} made the largest contribution to the DSRMT at all the times of day considered, while the impact of factors of the terrain declined close to noon (11:30 and 12:30hrs). Of the latter, the results were the most sensitive to A. We now further analyze the influence of the most sensitive factors as well as shadows on the performance of the DSRMT.

4.2.1. DSR-related variables

The sensitivity analysis above shows that the accuracy of the three DSR-related variables is crucial to the performance of the DSRMT. We thus examined its performance by using D_{ins-dir-rugged}, D_ins-dif-flat, and DSR_ins-flat from MCD18 as inputs. As described in Section 2.1.2, MCD18 provides only DSR-related parameters at the satellite overpass time. We hence used the DSRMT model for the time closest to that at which the MCD18 passed over the study region to calculate DSR_daily-rugged. If more than one value of DSR_daily-rugged was estimated in one day, $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ was calculated according to Equation (11(3b) $d_r=1+0.033\cos \left({\displaystyle \frac{2\pi doy}{365}}\right)$(3b) ). We used this procedure to validate values of $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ estimated by the proposed model based on data from MCD18 at the six mountainous sites, and compared them with the corresponding $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ obtained from CERES4 data at these sites (see Section 4.1.2). The comparison is shown in Figure 17.

Figure 17. Accuracy of $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ estimates obtained by the DSRMT against in-situ measurements by using DSR-related variables from (a) MCD18 and (b) CERES4 as inputs.

Only 165 samples from the in-situ measurements were used for comparison after matching the data. The resulting accuracy of estimation of $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ based on data from MCD18 was comparable to, but slightly better than, that obtained from data from CERES4, with RMSEs of 23.60 and 24.96 Wm⁻², MAEs of 18.09 and 19.25 Wm⁻², and R² values of 0.91 and 0.93, respectively. This result might have been obtained due to discrepancies in accuracy and spatial resolution between the DSR-related variables provided by CERES4 and MCD18. Therefore, the proposed DSRMT works well as long as accurate DSR-related parameters, D_ins-dir-flat, D_ins-dif-flat, and DSR_ins-flat, are available.

4.2.2. Terrain-related factors

We explored the influence of factors of the terrain, including A, S and the elevation, on the performance of the DSRMT. Table 9 presents the accuracy of validation of $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ estimates generated by the five DSRMT models under seven ranges of values of S (<20°, 20−30°, 30−40°, 40−50°, 50−60°, 60−70°, and >70°), nine ranges of values of A (<40°, 40−80°, 80−120°, 120−160°, 160−200°, 200−240°, 240−280°, 280−320°, and >320°), and seven ranges of elevation (<300, 300 m−1000, 1000 m−1500, 1500 m−2500, 2500 m−3500, 3500 m−4500 m, and >4500 m).

Table 9. The accuracy of validation of estimates of $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ generated by the five DSRMT models under various values of S, A, and elevation.

	RMSE (Wm^−2)	rRMSE (%)	MAE (Wm^−2)	R^2	No. of validation samples
	S
<20°	12.32	4.62	9.40	0.97	17,048
20° − 30°	13.34	5.14	10.21	0.97	17,279
30° − 40°	14.93	5.99	11.47	0.95	16,640
40° − 50°	16.52	6.75	12.51	0.91	14,913
50° − 60°	17.52	7.31	13.29	0.85	12,510
60° − 70°	18.67	8.28	14.39	0.79	5,543
>70°	23.05	10.93	16.98	0.65	2,800
	A
<40°	16.73	7.02	12.79	0.95	9,741
40° − 80°	16.49	6.65	12.96	0.92	7,964
80° − 120°	16.42	6.46	12.01	0.92	8,104
120° − 160°	16.11	6.19	11.79	0.91	10,031
160° − 200°	13.74	5.32	10.7	0.95	10,766
200° − 240°	15.34	6.04	11.27	0.93	12,732
240° − 280°	14.67	5.89	10.78	0.94	8,230
280° − 320°	15.72	6.42	11.86	0.94	9,609
>320°	13.90	5.83	10.74	0.97	9,556
	elevation
<300m	14.25	5.86	11.04	0.96	6,688
300 − 1000m	13.64	5.31	10.40	0.95	42,248
1000 − 1500m	15.87	6.56	12.27	0.95	12,378
1500 − 2500m	17.71	7.13	13.16	0.91	7,446
2500 − 3500m	19.64	8.06	14.51	0.89	5,402
3500 − 4500m	16.58	6.93	12.21	0.93	9,722
>4500m	21.06	8.63	16.55	0.87	2,849

The results showed that uncertainty in the estimated $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ was the lowest when S < 20° and the elevation <1000 m, and increased with an increase in S and the elevation until it reached its maximum value at S > 70° and elevation > 4500 m. The RMSEs of the models increased to 23.05 Wm⁻² at S > 70° and 21.06 Wm⁻² at elevation higher than 4500 m (rRMSE = 10.93%/8.63%). Regarding of S, the results were coincident with that in Figure 10 that DSR_daily-rugged related to DSR_ins-rugged more and more weakly when S increased. However, in terms of A, the uncertainty in the estimated $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ at western (>160°) were generally better than that at eastern (<160°), and was the smallest on south (160°−200°) with the smallest RMSE/rRSME values (13.74 Wm⁻²/5.32%).

Previous studies have noted that significant discrepancies occur in the accuracy of factors of the terrain when a variety of DEM data are used in various methods (Wu et al. 2018; Ma et al. 2022b). We thus investigated the influence of the DEM data and the method of calculation used to obtain the factors of the terrain on the performance of the DSRMT. The SRTM is one of the best sources of DEM data according to Tang et al. (2021). Ma et al. (2022b) found that the accuracy of the simulated DSR_ins-rugged and DSR_daily-rugged based on the DART and the solar spectrum vector code (6S) in the Chengde Experimental Area can be significantly improved if the SRTM DEM is used. Therefore, the SRTM DEM is appropriate for calculating factors of the terrain in this study. In the context of the methods of calculation, Wu et al. (2018) found that the discrepancy in V_d obtained when using different methods is more significant than those in S and A. The value of V_d calculated in this study was larger by about 0.04−0.18 than that obtained by using the method reported by Wu et al. (2018). The maximum D_ins-dif-flat under a clear sky was usually less than ∼60 Wm⁻², because of which D_{ins-dif-rugged} as calculated from D_ins-dif-flat through V_d was overestimated by 2.4–10.8 Wm⁻² at most, which is negligibly small if the data are expanded to the daily scale. Therefore, the influence of uncertainty in the factors of the terrain on the performance of the DSRMT, either due to the DEM product used or the method used for calculation, can be neglected.

4.2.3. Shadow

Shadows, which are usually caused by the surrounding terrain or the target object itself, are an inevitable factor influencing the estimation of radiation over a rugged terrain. To explore the applicability of the proposed DSRMT in cases involving shadows, we compared in-situ measurements recorded in the presence of shadows (no. of samples = 131) over the Chengde Experimental Area with its estimations of $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ and ${\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}}_{\mathrm{DCF}}$ based on CERES4 data, as shown in Figure 18.

Figure 18. The results of validation based on CERES4 in comparison with in-situ measurements in the presence of shadows. (a) $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ estimates of the proposed model and (b) ${\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}}_{\mathrm{DCF}}$ estimates of the DCF method.

The results indicate that shadows worsened the performance of both the DSRMT and the DCF, but the accuracy of estimation of $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ by the former was still higher than that of ${\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}}_{\mathrm{DCF}}$ by the latter, with an MAE that was smaller by ∼4.5 Wm⁻² and an RMSE smaller by ∼4.9 Wm⁻². The accuracy of estimation of $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ by the proposed method was still acceptable.

In summary, the DSRMT delivered relatively robust performance, but is easily influenced by DSR-related inputs.

5. Conclusions

To obtain accurate estimations of DSR_daily-rugged over rugged terrain, we proposed a conventional, easy-to-use, and highly efficient model called the DSRMT in this study. It uses comprehensive samples in DSR-relevant data from CERES4 based on the RF method. The DSRMT links DSR_daily-rugged non-linearly to three instantaneous DSR-relevant parameters on flat surfaces – D_{ins-dir-rugged}, calculated from D_ins-dir-flat, D_ins-dif-flat, and DSR_ins-flat – four factors of the terrain (S, A, V_d, and V_c) calculated from SRTM DEM data, and ancillary information (e.g. the GLASS α and geometric information). We built five DSRMT models for five times of the day (10:30, 11:30, 12:30, 13:30, and 14:30hrs) by considering the relationships between DSR_ins-rugged and DSR_daily-rugged. DSR_daily-rugged was easily obtained by the DSRMT as long as one group of inputs was provided at one of the five times, and the final $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$ was calculated by the weighted average of u_s if more than one DSR_daily-rugged estimate was obtained in one day. However, the proposed DSRMT works only over rugged terrain without bright surface under clear sky.

We used validation samples from CERES4 and ground measurements at six mountainous sites to validate the performance of the DSRMT, and compared it with the DCF method and the sinusoidal model. It delivered satisfactory performance at all of the five times considered, with the best results at noon. Moreover, $\overline{\mathrm{DS}{\mathrm{R}}_{\mathrm{daily}\text{-}\mathrm{rugged}}}$, obtained by the proposed model based on the weight average, yielded even better values of the RMSE and MAE, ranging from 21.63–25.84 Wm⁻² and from 17.14–19.72 Wm⁻², respectively. This was superior to the results of the averaged values obtained by the DCF method (RMSE = 32.13 Wm⁻², MAE = 23.75 Wm⁻², R²= 0.90) and most results of sinusoidal model (RMSE: 24.04–45.31 Wm⁻², MAE: 17.30–39.24 Wm⁻², R²: 0.84–0.94), and it delivered more robust performance. Several conclusions can be drawn based on the results of this study: (1) Of all the variables considered, the three DSR-related parameters were the most influential for the performance of the DSRMT, and its accuracy can be improved if optimal values of these parameters are used (e.g. from MCD18). (2) The solar radiation reflected from regions adjacent to the one of interest cannot be ignored when estimating DSR_daily-rugged, especially when θ₀ is large. (3) Aspect is the most crucial terrain-related factor influencing the DSR over a rugged terrain, especially in the west and the south. Its topographic effect on the DSR was minor when S was smaller than 10°. (4) Shadows, from the surrounding terrain within a 3 × 3 window of SRTM DEM, worsened the accuracy of estimation of the DSRMT, although its results were still acceptably accurate.

Overall, the proposed DSRMT model can be used obtain accurate values of DSR_daily-rugged as long as data on the instantaneous, direct, and diffuse solar radiation on clear days are available during 10:30‒14:30hrs. It thus offers significant promise for widespread use. Efforts are underway to expand its scope of application, such as over rugged terrain with α>0.7 and under cloudy skies. In addition, we plan to further validate the DSRMT based on more ground measurements.

Authorship contribution statement

Hui Liang: Investigation, Methodology, Software, Data curation, Writing-original draft. Bo Jiang: Conceptualization, Methodology, Resources, Investigation, Data curation, Writing-review & editing, Funding acquisition, Supervision. Jianghai Peng, Shaopeng Li: Software support, Discussion. Jiakun Han, Xiuwan Yin: Discussion.

Acknowledgments

The authors thank NASA and GLASS for providing remotely sensed data. We also thank Prof. María J. Polo for help with using the GMN data. The authors are also grateful to Prof. Guangjian Yan, Prof. Xihan Mu, Prof. Donghui Xie, and Mr. Yingji Zhou for providing valuable in-situ measurements from Chengde. This dataset was provided by the Chengde Remote Sensing Test Site of the State Key Laboratory of Remote Sensing Sciences.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was funded by the National Key Research and Development Program of China [2020YFA0608704] and the National Natural Science Foundation of China [grant no 41971291].