Estimation of Shortwave Solar Radiation on Clear-Sky Days for a Valley Glacier with Sentinel-2 Time Series

: Downward surface shortwave radiation (DSSR) is the main energy source for most glacial melting, and Moderate Resolution Imaging Spectroradiometer (MODIS) and Landsat Thematic Mapper (TM) data have been used extensively in the inversion of input parameters for estimating DSSR. However, for valley glaciers under complex climatic conditions, the values of MODIS atmospheric products, especially aerosol products, are often invalid, and TM images are always saturated with snow. Furthermore, an estimation model based on optical satellite images must simultaneously consider terrain and atmospheric e ﬀ ects and the transient nature of ice / snow albedo. Based on a high-resolution (12 m) digital elevation model (DEM), the newly launched Sentinel-2 satellites, rather than MODIS and TM, were used to provide input data for our published mountain radiation scheme in a valley glacier. Considering Laohugou Glacier No. 12 as the study area, 62 typical Sentinel-2 scenes were selected and spatiotemporal DSSR variations on the glacier surface were obtained with a 10 m spatial resolution during a mass-balance year from September 2017 to August 2018. Ground-based measurements on 52 clear-sky days were used for validation and the mean bias error (MBE = − 16.0 W / m 2 ) and root-mean-square di ﬀ erence (RMSD = 73.6 W / m 2 ) were relatively low. The results conﬁrm that DSSR is a ﬀ ected mainly by the solar zenith angle and atmospheric attenuation in ﬂat areas of valley glaciers, while in areas with complex terrain, the DSSR received by the glacier surface is a ﬀ ected primarily by the terrain and ice / snow albedo, which exhibits very high spatial heterogeneity.


Introduction
Glaciers provide valuable aid in the monitoring of climate change and downward surface shortwave radiation (DSSR) is a major driving factor in the melting of glaciers [1][2][3][4][5][6]. Since the 1880s, some researchers have noted the contributions of solar radiation to glacial changes and ablation [7][8][9][10]. In rugged terrain, however, numerous factors and processes seriously affect the spatial distribution of the DSSR received by a given surface point. In addition to the well-known sun-surface geometrical relationship, the following aspects are important components in the complexity of estimating the solar radiation received by valley glaciers: (1) local topographic factors, such as the slope and aspect, the obstruction coefficient, the sky view factor, and the topographic configuration [11][12][13][14][15]; (2) the fact that radiation in the atmosphere is attenuated more easily by aerosols and precipitable water (PW) than by Rayleigh scattering and absorption achieved through the ozone column and uniformly mixed gases [16]; and (3) the transient nature of ice/snow cover and the high-reflection characteristics of the snow surface [17,18]. Therefore, accurately estimating the DSSR received by mountain glacier surfaces has become a popular research topic.
The global radiation received by a rugged terrain surface consists of three components: direct, diffuse, and surrounding-reflected radiation. In general, the energy of direct solar radiation is the largest, followed by that of diffuse radiation and finally that of the radiation reflected from the surrounding terrain. However, when the surrounding environment consists of objects covered by ice and snow, the radiation reflected from the surroundings will be large and non-negligible; consequently, the albedo of the surrounding slope ice under the complicated terrain conditions of a valley glacier will be substantial [19]. Furthermore, local topography drastically alters the spatial distribution of mountain glacial DSSR values. Many researchers [20][21][22][23] have developed a variety of complex topographic factors to finely describe the solar radiation energy received by the slope surface: the cosine of the solar illumination angle (cos i s ) on a sloped grid affected by the terrain slope and aspect, solar zenith, and azimuth angles; the obstruction coefficient (V s ), which indicates whether the surface is obscured by the surrounding terrain; the sky view factor (V iso ), which describes the area ratio of the sky dome seen by the slope surface; and the topographic configuration factor (F ij ), which describes the energy ratio of the surrounding visible pixels to the object pixel.
Some studies have quantitatively estimated the contributions of topographic factors to glacial DSSR values using glacier meteorological observations [11][12][13]. Kruss and Hastenrath published three consecutive articles investigating the latitude effect and the geometrical effects of shortwave solar radiation on the climate responses of mountain glaciers in the cases of both horizontal and sloping surfaces [2,24,25]. Klok and Oerlemans [26] found that neglecting the effects of topographic factors resulted in a 37% increase in the estimated annual incoming shortwave radiation received by a glacier. Jiskoot and Mueller [27] confirmed that the variability in solar irradiance caused by terrain exaggerates the heterogeneity of melting across glacier surfaces on clear-sky days. Mögland et al. [8] qualitatively demonstrated that the retreat of vertical ice walls/cliffs at the summit was the secondary energy source (after solar radiation) affecting the recession of glaciers on Kilimanjaro. Based on the Yang broadband atmospheric attenuation model, Zhu et al. [10] found that shortwave radiation provides up to 60% of the melting energy for both the Zhadang glacier and the Parlung No. 4 glacier, and the ice/snow albedo was assumed to vary as a function of the dew-point temperature. However, the above mentioned studies have shortcomings in one or more of the following aspects: only simple topographic factors such as slope were considered; the estimation of atmospheric transmittance was strongly dependent on ground observations; and glacier albedo was assumed to be homogeneous over large areas or was calculated using a simple empirical scheme according to the characteristics of ice/snow measured in the field.
Satellite remote sensing has become an effective means for achieving glacial monitoring without the need to visit the site [28,29]. It is challenging to obtain field measurements of glacial DSSR because glaciers are poorly accessible due to their adverse natural conditions, and the observation instruments used to monitor glacial melting and movement are unstable. Two methods based on remote sensing data are used to estimate DSSR: empirical models, which use regression relationships between the radiance measured by satellites and ground measurements, and theoretical models, which use radiative transfer models or parameterization schemes. At global or regional scales, many scientific teams have released solar radiation products, such as the Global Energy and Water Cycle Experiment (GEWEX). However, for a valley glacier with complex terrain, mountainous radiative transfer models are more conducive than empirical models to accurately calculating the DSSR. Hence, an estimation model for such a valley glacier should consider the effects of the following three factors simultaneously: terrain effects, atmospheric attenuation, and the spatiotemporal variability of ice/snow albedo.
Moderate Resolution Imaging Spectroradiometer (MODIS) and Landsat Thematic Mapper (TM) remote sensing datasets have been used extensively for the estimation of mountain solar radiation. Zhang et al. [30] accurately estimated solar radiation components for complex mountains by integrating the Li Mountain radiation algorithm [23] and the Yang broadband atmospheric attenuation model [16] with MODIS and TM imagery. However, MODIS atmospheric products cannot satisfy the requirements of practical applications for relatively small mountain glaciers. First, the values of MODIS atmospheric products (especially aerosol products) for mountain glaciers are often invalid; second, the spatial resolution is too coarse to describe the atmospheric conditions of valley glaciers in detail. In addition, Landsat TM images are often saturated with snow pixels, especially on sun-facing slopes [23,31,32]. As alternatives, the Sentinel-2A and Sentinel-2B (S2 A/B) satellites were launched by the European Space Agency (ESA) in June 2015 and March 2017, respectively.
Compared with MODIS and Landsat TM, the use of Sentinel-2 (S2) satellite data to estimate valley glacier DSSR has four advantages: (1) S2 can provide both atmospheric and surface information, thus reducing the errors caused by different satellite data; (2) S2 data can improve the DSSR spatiotemporal resolution and provide reliable information for understanding the melting of glaciers; (3) S2 data can almost completely avoid the saturation defect of the snow signal from optical satellites because of its high radiometric resolution; (4) the data processing software Sen2Cor for S2 provided by the official website can use S2 itself for atmospheric and surface signals to retrieve glacier albedo with the aid of digital elevation model (DEM). Therefore, S2 has become a significant source of information for glacier research [29,33,34].
The purpose of this research is to introduce images of the newly launched advanced S2 satellite into the previously published mountain radiation scheme, replacing MODIS atmospheric products and Landsat TM imagery, to obtain the DSSR received by a valley glacier surface with high temporal and spatial resolution. Due to the complex effects of clouds on the solar radiation received by mountain glaciers [6] and the difficulty associated with monitoring clouds in optical remote sensing data, only the characteristics of DSSR on clear-sky days are considered in this study. The remainder of this paper is organized as follows. In Section 2, the study area and data are introduced. In Section 3, the principles of the DSSR estimation algorithm using S2 data are described; the results and verification are shown in Section 4. The discussion is presented in Section 5. Finally, the conclusions are given in Section 6.

Study Area
Laohugou Glacier No. 12 (96.5 • N, 39.5 • E) is located on the northwest edge of the Tibetan Plateau and faces a northern slope with 9.85 km long, with an area of approximately 20.4 km 2 , which is an important water source for the Shulehe River basin. It is the largest valley glacier located within the Qilian Mountains and is the most typical continental glacier in the region [6,35]. Its mean annual air temperature is more than 0 • C and major precipitation occurs in summer. In recent years, Laohugou Glacier has shrunk due to increases in the warmth and moisture of the regional climate [36]. The terrain of Laohugou Glacier No. 12 is complicated, especially in the southern accumulation zone, with glacier altitude differences greater than 1200 m (from 4202 m to 5427 m a.s.l.) according to the latest TanDEM-X DEM dataset (spatial resolution of 12 m) provided by the German Aerospace Center (DLR) (Figure 1).
The map in Figure 1 shows the locations of two automatic weather stations (AWSs) installed in relatively flat areas (slopes with inclinations of less than 3 • ) [6]. AWS1 (5040 m a.s.l.) lies in the accumulation zone, which is covered by snow beginning in late September or early October until late July, and AWS2 is installed near an altitude of 4550 m in the ablation zone. The map in Figure 1 shows the locations of two automatic weather stations (AWSs) installed in relatively flat areas (slopes with inclinations of less than 3°) [6]. AWS1 (5040 m a.s.l.) lies in the accumulation zone, which is covered by snow beginning in late September or early October until late July, and AWS2 is installed near an altitude of 4550 m in the ablation zone.

Data Sources
As shown in Table 1, this research relies on three data sources. (1) Pyranometer datasets for validating the model were obtained at the two AWSs. AWS1 and AWS2 were equipped with Kipp and Zonen CNR1 four-component net radiometers and were connected to a data collector (CR1000, Campbell Scientific Inc., USA) with a low temperature resistance (-55°C), which acquired observations every 30 min by taking the average every 10 seconds. Notably, due to snow cover and rime formation on the pyranometer sensor and the low solar altitude angle, the values of DSSR measurements can be distorted. The detailed data processing and quality control methods of field observations are provided in the literature [37]. Within these data, the clear-sky condition was visually judged from S2 images. (2) DEM datasets with different spatial resolutions were acquired, including a 90 m Shuttle Radar Topography Mission (SRTM) DEM, a 30 m global DEM (GDEM), and the 12 m DEM acquired from the DLR. The SRTM and GDEM datasets were both obtained from the USGS Earth Resources Observation System (EROS) Data Center. (3) Sixty-two S2A/B Level-1C (L1C) Multispectral Imager (MSI) images were collected during one mass-balance year (from September 2017 to August 2018); the images were freely obtained from the ESA Copernicus Open Access Hub. The DEMs and images were all set to the Universal Transverse Mercator/World Geodetic System 1984 (UTM/WGS84) projection/coordinate reference frame.

Data Sources
As shown in Table 1, this research relies on three data sources. (1) Pyranometer datasets for validating the model were obtained at the two AWSs. AWS1 and AWS2 were equipped with Kipp and Zonen CNR1 four-component net radiometers and were connected to a data collector (CR1000, Campbell Scientific Inc., USA) with a low temperature resistance (-55 • C), which acquired observations every 30 min by taking the average every 10 s. Notably, due to snow cover and rime formation on the pyranometer sensor and the low solar altitude angle, the values of DSSR measurements can be distorted. The detailed data processing and quality control methods of field observations are provided in the literature [37]. Within these data, the clear-sky condition was visually judged from S2 images.

Methods
This paper adopted the mountain radiative transfer scheme from Zhang et al. [30] that combines the Li Mountain radiation algorithm [23] with the Yang broadband atmospheric attenuation model [16]. Based on the DEM, the Li mountainous spectral radiation scheme fully considers the effects of terrain factors, such as the slope, aspect, sky view factor, and obstruction coefficient, on solar radiation and calculates the surrounding-reflected radiation by combining Landsat TM images. The Yang model proposes a general formula for estimating atmospheric broadband transmittance. Nevertheless, our published DSSR estimation method can accurately calculate the DSSR received by the slope surface by introducing MODIS aerosol optical depth (AOD) and PW products and Landsat TM images as input parameters. However, in this study, the MODIS atmospheric products and Landsat imagery were replaced by S2 products in estimating the DSSR of a valley glacier for the following three reasons. First, the values of MODIS aerosol products are often invalid for small and complex valley glaciers. Second, Landsat TM imagery is generally saturated with snow. Third, S2 can provide finer glacier albedo estimates (10 m) for small-scale glaciers. A flowchart of the method used in this study is shown in Figure 2.
Remote Sens. 2019, 11, x FOR PEER REVIEW This paper adopted the mountain radiative transfer scheme from Zhang et al. [30] that com the Li Mountain radiation algorithm [23] with the Yang broadband atmospheric attenuation [16]. Based on the DEM, the Li mountainous spectral radiation scheme fully considers the eff terrain factors, such as the slope, aspect, sky view factor, and obstruction coefficient, on radiation and calculates the surrounding-reflected radiation by combining Landsat TM image Yang model proposes a general formula for estimating atmospheric broadband transmi Nevertheless, our published DSSR estimation method can accurately calculate the DSSR recei the slope surface by introducing MODIS aerosol optical depth (AOD) and PW products and L TM images as input parameters. However, in this study, the MODIS atmospheric produc Landsat imagery were replaced by S2 products in estimating the DSSR of a valley glacier following three reasons. First, the values of MODIS aerosol products are often invalid for sma complex valley glaciers. Second, Landsat TM imagery is generally saturated with snow. Third, provide finer glacier albedo estimates (10 m) for small-scale glaciers. A flowchart of the metho in this study is shown in Figure 2.

Shortwave Solar Radiation in Mountainous Areas
The global radiation received by the surface of a mountainous glacier can be decompose three components: direct, diffuse, and surrounding-reflected radiation. Considering the in direction of solar radiation, diffuse radiation can be further divided into isotropic and aniso diffuse radiation components. As shown in Equation (1), global shortwave solar radiation expressed as follows [30]: , and denote global, direct, diffuse, isotropic, aniso and reflected radiation, respectively. Different solar radiation components are affected by different terrain factors. The radiation component is calculated by considering the terrain shading factor, the cosine of the illumination angle of the slope pixel, and the atmospheric transmittance. The factors that inf anisotropic diffuse radiation are similar to those that influence direct solar radiation with the ad of the solar anisotropy index. Isotropic diffuse radiation is calculated by considering the sk factor, anisotropy index, and atmospheric diffuse radiative transmittance. Calculatin surrounding-reflected solar radiation requires the surface albedo, the topographical structural and the sum of the direct solar radiation and the diffuse radiation from the surrounding visible The specific calculation formulas for these four components can be found in the study by Zh al. [30].

Shortwave Solar Radiation in Mountainous Areas
The global radiation received by the surface of a mountainous glacier can be decomposed into three components: direct, diffuse, and surrounding-reflected radiation. Considering the incident direction of solar radiation, diffuse radiation can be further divided into isotropic and anisotropic diffuse radiation components. As shown in Equation (1), global shortwave solar radiation can be expressed as follows [30]: where E, E dir , E di f , E iso_di f , E aniso_di f , and E re f denote global, direct, diffuse, isotropic, anisotropic, and reflected radiation, respectively. Different solar radiation components are affected by different terrain factors. The direct radiation component is calculated by considering the terrain shading factor, the cosine of the actual illumination angle of the slope pixel, and the atmospheric transmittance. The factors that influence anisotropic diffuse radiation are similar to those that influence direct solar radiation with the addition of the solar anisotropy index. Isotropic diffuse radiation is calculated by considering the sky view factor, anisotropy index, and atmospheric diffuse radiative transmittance. Calculating the surrounding-reflected solar radiation requires the surface albedo, the topographical structural factors, and the sum of the direct solar radiation and the diffuse radiation from the surrounding visible pixels. The specific calculation formulas for these four components can be found in the study by Zhang et al. [30].
In short, there are four input parameters for the instantaneous DSSR estimation model received by the undulating glacier surface: the sun-surface geometry, local topography, atmospheric factors and glacier albedo. The sun-surface geometry consists of the solar zenith angle and solar azimuth, which are determined by the S2A/B satellite overpass time. Various local terrain factors are calculated Remote Sens. 2020, 12, 927 6 of 19 from high-resolution DEMs; the details of the corresponding processes can be found in the studies of Li et al. [22,23]. Various atmospheric factors, especially the PW and aerosol concentrations, are further inverted from S2 images. The glacier albedo can be estimated using multiband S2 satellite remote sensing data. Therefore, extracting the atmospheric factors and glacier albedo from S2 satellite data is the key to DSSR estimation in a valley glacier. The following subsection will introduce in detail how to retrieve the PW and aerosol concentrations and the ice/snow albedo from S2A/B imagery.

Inversion of the Input Parameters
S2 satellite data (including S2A and S2B) have a higher radiometric resolution (12 bit) and a higher revisit cycle of 5 days than MODIS and TM satellite data under the same viewing conditions. The top-of-atmosphere (TOA) L1C reflectance products consist of 13 spectral bands with three different resolutions (10 m, 20 m, and 60 m). The MSI instrument has a wide-spectrum optical range (from the visible and near-infrared (VNIR) spectrum to the shortwave infrared (SWIR) spectrum). Reliable information for both the surface and the atmosphere is recorded in these bands [29,33,34]. For instance, bands 2 and 12 can reflect aerosol information; bands 8a and 9 can invert the PW in the atmosphere; band 10 is sensitive to atmospheric cirrus; and bands 10 and 11 are the characteristic bands for the detection of ice, snow, and clouds.
The S2 Level-2A (L2A) products released by the ESA not only provide bottom-of-atmosphere (BOA) reflectance ortho-images that can be used to perform atmospheric, topographic and cirrus corrections (the last two of which are optional) on L1C images but also output atmospheric products, such as PW and AOD, 550 nm) maps [38]. These factors are ideal input parameters for the atmospheric and surface properties of a solar radiation model. However, for China, the current ESA website provides L2A products only after December 2018. Fortunately, the official website provides Sen2Cor software, which can process L1C products to obtain L2A data by configuring parameters such as the DEM, terrain correction and bidirectional reflectance distribution function (BRDF) in ground image processing parameter (GIPP) files (L2A_GIPP.xml) [38].

Retrieval of AOD and PW Products
To retrieve AOD and PW products, the S2 L1C data are converted into TOA radiances by a radiometric calibration, and the resulting data are pre-classified coarsely into land, water, cloud, and ice/snow maps. Moreover, a database of radiative transfer calculations (i.e., look-up tables (LUTs)) is generated using the libRadtran4 code, which is in 6-dimensional parameter space. In addition, the aerosol types are divided into four categories (rural, urban, desert and maritime).
Finding dark reference areas is an important step in aerosol inversion. The SWIR band (band 12, 2.19 µm) is first utilized to search for dense dark vegetation (DDV) pixels, while bands 4 (red) and 2 (blue) are employed to estimate the visibility. The advantage of these data over MODIS [39] is that the threshold of the dark pixel reflectance (i.e., a medium brightness) can be dynamically adjusted based on the scene; as a result, a certain percentage (2%-5%) of dark pixels can be obtained, after which the AOD value can finally be calculated. With a 3 km × 3 km spatial filter, the average AOD value of non-reference pixels can be retrieved. The algorithms within the Sen2Cor toolbox can be used to estimate the AOD over areas with small, bright surfaces (i.e., mountain glaciers) as an alternative to using MODIS aerosol products, which are often incapable of performing an inversion in such areas.
After performing the aerosol inversion, the atmospheric pre-corrected differential absorption (APDA) water vapour inversion algorithm [40] can be implemented with two S2 bands, namely, band 8a (865 nm) and band 9 (945 nm). Details regarding the specific algorithm can be found in the literature [41].

Retrieval of Ice/Snow Albedo in Mountainous Terrain
The surface reflectance inversion accuracy using optical satellite remote sensing data in rugged terrain is restricted simultaneously by atmospheric, topographic and surface BRDF effects [42]. Similarly, ice/snow albedo is retrieved from S2 L2A products that have been corrected by performing topographic normalization based on an empirical BRDF model, which can reduce the atmospheric, topographic, and BRDF effects from optical satellite images using the ACTOR3 radiative transfer model proposed by Richter (1998) [43] in Sen2Cor. Haze and cirrus corrections are also performed.
The principle of topographic normalization is briefly described here. The surface reflectance ρ l of each pixel on a Lambertian surface can be calculated iteratively based on a DEM by the following: where E dir , E di f , and E global represent the direct, diffuse and global (direct plus diffuse) solar radiation components, respectively, received by sloped terrain. L sat , L p , d, τ v , and ρ terrain are the radiance measured at the sensor, the path radiance, the Earth-Sun distance factor, the transmittance at the ground-sensor view angle, and the average terrain reflectance, respectively. Topographic normalization also includes the preparation of various terrain factors, such as the obstruction coefficient, the sky view factor, and the topographic configuration factor. A detailed description and discussion of this algorithm are given in the literature [43].
To reduce the high reflectance values in regions of low illumination caused by the slope surface BRDF effect in ACTOR3, an empirical geometric function G is introduced to obtain the true spectral reflectance ρ t in Equation (3). The G function varies between a specified lower boundary g and 1 (Equation (4)). G values less than the lower boundary g are set to g, while values greater than 1 are set to 1: where i s indicates the illumination angle and i T represents the threshold angle. In addition, the topographic correction has a considerable dependence on the DEM spatial resolution, and the DEM grid size must be less than the scale of the original image [44,45]. Therefore, 12 m DEM data are used to correct for terrain effects to improve the correction accuracy. Accordingly, by modifying some optional parameters in the Sen2Cor toolbox, the L2A products produced following this approach can be improved compared to the standard products produced by the default parameter settings. The parameters established in the L2A_GIPP files are shown in Table 2. After topographic normalization, L2A products, i.e., surface reflectance images, can be considered to represent the spectral albedo of the glacier surface. Deriving the broadband or total albedo of the glacier surface requires a narrow-to-broadband conversion after conducting atmospheric and topographic corrections. For instance, Naegeli et al. [18] adopted the glacier albedo formula proposed by Knap et al. [46] that had been developed for Landsat 5/7 TM images. Hence, considering the bandwidth difference, two corresponding L2A bands are selected to estimate the spatial distribution of the albedo with the following formula: where α is the glacier albedo and b 3 and b 8 represent band 3 and band 8, respectively.

Evaluation against Ground-Based Measurements
To validate the reliability of the proposed DSSR model, in situ observations acquired on 52 clear-sky days were selected from two AWS meteorological stations for the period from January 2018 to August 2018. The ground-based observations at AWS1 and AWS2 were recorded every 30 min. The S2A/B overpass time was approximately 12:30 local time, and the time difference between the two AWS stations was less than 15 min. Thus, the 12:30 measurements of AWS1 and AWS2 were selected. Figure 3 shows the scatterplots and statistical results for the solar radiation produced by comparing the estimated and measured global radiation. Detailed definitions of the statistical indicators used herein can be found in the literature [47]. The results confirm that the DSSR estimates are consistent with the ground-based measurements (R 2 = 0.864). The root-mean-square difference (RMSD) between the instantaneous DSSR model estimates and measurements is 73.6 W/m 2 , and the RMSD percentage (RMSD %) is 7.9%. The modelled radiation estimates are below the 1:1 line; thus, the mean bias error (MBE) is -16.0 W/m 2 .
considering the bandwidth difference, two corresponding L2A bands are selected to estimate the spatial distribution of the albedo with the following formula: = 0.726 − 0.322 − 0.015 + 0.581 , (5) where is the glacier albedo and and represent band 3 and band 8, respectively.

Evaluation against Ground-Based Measurements
To validate the reliability of the proposed DSSR model, in situ observations acquired on 52 clearsky days were selected from two AWS meteorological stations for the period from January 2018 to August 2018. The ground-based observations at AWS1 and AWS2 were recorded every 30 min. The S2A/B overpass time was approximately 12:30 local time, and the time difference between the two AWS stations was less than 15 min. Thus, the 12:30 measurements of AWS1 and AWS2 were selected. Figure 3 shows the scatterplots and statistical results for the solar radiation produced by comparing the estimated and measured global radiation. Detailed definitions of the statistical indicators used herein can be found in the literature [47]. The results confirm that the DSSR estimates are consistent with the ground-based measurements (R 2 = 0.864). The root-mean-square difference (RMSD) between the instantaneous DSSR model estimates and measurements is 73.6 W/m 2 , and the RMSD percentage (RMSD %) is 7.9%. The modelled radiation estimates are below the 1:1 line; thus, the mean bias error (MBE) is -16.0 W/m 2 .

Mapping Solar Radiation
The incoming solar radiation on 62 clear-sky days between September 2017 and August 2018 was simulated at the S2A/B overpass time. Four typical seasonal periods of solar radiation were selected for analysis. Generally, maps of the spatial distribution of radiation show that the lowest values are observed close to the northern slope and that the maximum insolation is located on the southern slope (Figure 4 (a)). Additionally, the incoming solar radiation received by the glacier surface varies strongly over the considered mass-balance year. Figure 4 (a) indicates that the DSSR is weakest in December (mean of 398.2 W/m 2 ) and strongest in May (mean of 948.9 W/m 2 ).

Mapping Solar Radiation
The incoming solar radiation on 62 clear-sky days between September 2017 and August 2018 was simulated at the S2A/B overpass time. Four typical seasonal periods of solar radiation were selected for analysis. Generally, maps of the spatial distribution of radiation show that the lowest values are observed close to the northern slope and that the maximum insolation is located on the southern slope ( Figure 4a). Additionally, the incoming solar radiation received by the glacier surface varies strongly over the considered mass-balance year. Figure 4a indicates that the DSSR is weakest in December (mean of 398.2 W/m 2 ) and strongest in May (mean of 948.9 W/m 2 ).
As shown in Table 3, the seasonal variation in the incoming irradiance is dominated by the solar zenith angle. That is, for flat areas in glacial regions (i.e., areas with small slope inclinations), the solar zenith angle determines the average incoming solar radiation: the greater the solar zenith angle, the smaller the solar irradiance, and vice versa. The global minimum solar radiation is observed in December and the global maximum solar radiation is observed in June. As shown in Table 3, the seasonal variation in the incoming irradiance is dominated by the solar zenith angle. That is, for flat areas in glacial regions (i.e., areas with small slope inclinations), the solar zenith angle determines the average incoming solar radiation: the greater the solar zenith angle, the smaller the solar irradiance, and vice versa. The global minimum solar radiation is observed in December and the global maximum solar radiation is observed in June. In addition, the relationship between the solar radiation and solar zenith angle was explored using ground-based measurements. According to the data, the amount of solar irradiance decreases rapidly as the solar zenith angle increases (R 2 = 0.85), and that these values are negatively correlated ( Figure 5).   In addition, the relationship between the solar radiation and solar zenith angle was explored using ground-based measurements. According to the data, the amount of solar irradiance decreases rapidly as the solar zenith angle increases (R 2 = 0.85), and that these values are negatively correlated ( Figure 5).  As shown in Table 3, the seasonal variation in the incoming irradiance is dominated by the solar zenith angle. That is, for flat areas in glacial regions (i.e., areas with small slope inclinations), the solar zenith angle determines the average incoming solar radiation: the greater the solar zenith angle, the smaller the solar irradiance, and vice versa. The global minimum solar radiation is observed in December and the global maximum solar radiation is observed in June. In addition, the relationship between the solar radiation and solar zenith angle was explored using ground-based measurements. According to the data, the amount of solar irradiance decreases rapidly as the solar zenith angle increases (R 2 = 0.85), and that these values are negatively correlated ( Figure 5). Additionally, as reflected in Figure 4a, distinct irregular fluctuations in the measured global radiation are caused by the presence of rugged terrain. In flat areas, the difference in solar radiation is small; in contrast, the DSSR received by the southern portion of the basin is strongly heterogeneous due to the rugged topography. In addition, the solar irradiance over an ice/snow-covered surface is affected by the reflection characteristics of the surrounding visible terrain, as shown in Figure 4b. Hence, a visual comparison of these slope and aspect maps with the corresponding spatial heterogeneity of the estimated solar radiation shows a near-perfect match.

Spatial Profile Characteristics
Two interpolated lines were extracted to explore the spatial dependence of the glacier DSSR on the terrain (Figure 6a). One line (considered a vertical line) was the midline through the east branch of the glacier, while the other line (considered a horizontal line) stretched across the southern area of the glacier. Beginning from the origin, a point was selected every 10 m along both lines, resulting in 982 points on the vertical line and 599 points on the horizontal line. Figure 6b shows that the variation in the DSSR along the vertical interpolated line is small, and the solar irradiance intensity is approximately 900 W/m 2 , except for some anomalously low points in the southern region due to the large slope. However, along the horizontal interpolated line, the variation in the DSSR is very strong, with a difference exceeding 1000 W/m 2 (334 to 1460 W/m 2 ). These results imply that data collected at several sparsely distributed AWSs cannot be representative of the incoming solar irradiance throughout a mountain glacier [11,24]. Although some scholars have found that the distribution of solar radiation between the glacier margin and the center flow line is highly heterogeneous, some studies have assumed that the DSSR measured at a single site is representative of an entire mountain glacier, as has often been the case in many mass-balance models [27,48,49]. Evidently, this assumption undoubtedly introduces considerable uncertainty into mountain glacier research. Additionally, as reflected in Figure 4 (a), distinct irregular fluctuations in the measured global radiation are caused by the presence of rugged terrain. In flat areas, the difference in solar radiation is small; in contrast, the DSSR received by the southern portion of the basin is strongly heterogeneous due to the rugged topography. In addition, the solar irradiance over an ice/snow-covered surface is affected by the reflection characteristics of the surrounding visible terrain, as shown in Figure 4 (b). Hence, a visual comparison of these slope and aspect maps with the corresponding spatial heterogeneity of the estimated solar radiation shows a near-perfect match.

Spatial Profile Characteristics
Two interpolated lines were extracted to explore the spatial dependence of the glacier DSSR on the terrain (Figure 6 (a)). One line (considered a vertical line) was the midline through the east branch of the glacier, while the other line (considered a horizontal line) stretched across the southern area of the glacier. Beginning from the origin, a point was selected every 10 m along both lines, resulting in 982 points on the vertical line and 599 points on the horizontal line. Figure 6 (b) shows that the variation in the DSSR along the vertical interpolated line is small, and the solar irradiance intensity is approximately 900 W/m 2 , except for some anomalously low points in the southern region due to the large slope. However, along the horizontal interpolated line, the variation in the DSSR is very strong, with a difference exceeding 1000 W/m 2 (334 to 1460 W/m 2 ). These results imply that data collected at several sparsely distributed AWSs cannot be representative of the incoming solar irradiance throughout a mountain glacier [11,24]. Although some scholars have found that the distribution of solar radiation between the glacier margin and the center flow line is highly heterogeneous, some studies have assumed that the DSSR measured at a single site is representative of an entire mountain glacier, as has often been the case in many mass-balance models [27,48,49]. Evidently, this assumption undoubtedly introduces considerable uncertainty into mountain glacier research.    Figure 6c shows the albedo along the same two interpolated lines to further demonstrate the importance of the contribution of the glacier albedo to the incoming solar irradiance. On a flat surface, even a high albedo cannot change the steady variation in the solar radiation. However, on complex terrain, especially along the horizontal interpolated line, ice/snow albedo exacerbates the spatiotemporal heterogeneity of the DSSR. Therefore, the glacier albedo must be retrieved by satellite images rather than assumed to be a certain value before estimating the DSSR received by rugged terrain.

Spatiotemporal Variation
To quantify the spatiotemporal variation in the DSSR across the mountain glacier, 6 representative topographic locations, i.e., points P01, P02, P03, P04, P05, and P06, were selected for analysis. Among them, points P01, P02, P03, and P06 are located in relatively flat areas at different glacial elevations, as shown in Figure 7a. In contrast, point P04 lies on a 38 • slope in the southeast and point P05 lies on a 24 • slope in the north.
Remote Sens. 2019, 11, x FOR PEER REVIEW 11 of 19 even a high albedo cannot change the steady variation in the solar radiation. However, on complex terrain, especially along the horizontal interpolated line, ice/snow albedo exacerbates the spatiotemporal heterogeneity of the DSSR. Therefore, the glacier albedo must be retrieved by satellite images rather than assumed to be a certain value before estimating the DSSR received by rugged terrain.

Spatiotemporal Variation
To quantify the spatiotemporal variation in the DSSR across the mountain glacier, 6 representative topographic locations, i.e., points P01, P02, P03, P04, P05, and P06, were selected for analysis. Among them, points P01, P02, P03, and P06 are located in relatively flat areas at different glacial elevations, as shown in Figure 7 (a). In contrast, point P04 lies on a 38° slope in the southeast and point P05 lies on a 24° slope in the north.  Figure 7 (b) illustrates curves reflecting the spatiotemporal variation in the instantaneous solar irradiance at the six points over a mass-balance year. The seasonal variations in the incoming solar radiation are similar (sinusoidal) at the four points on the flat ice surface, but the other two points display completely different trends. The solar irradiance at point P05 is very low throughout the year, while that at point P04 is very high during the same period. Therefore, curves of the three DSSR components were generated to explore the reason for this disparity.
As shown in Figure 8, the amount of direct radiation received by point P05 is zero from November 20, 2017 to January 22, 2018 due to shade from the terrain. Global solar radiation is derived mainly from diffuse radiation, but this solar radiation component is very weak due to the small sky view factor. During the other seasons, although point P05 is illuminated by the sun, the direct radiation is still low because the cosine of the solar illumination angle at point P05 is smaller than those at the five other points. However, the spatiotemporal distribution of solar radiation reaching point P04 is completely different. Because point P04 is located on a steep southeast-facing slope, the cosine of the solar illumination angle is always very high, and thus, the direct solar radiation received at this point is very intense. Similarly, on this steep southeast-facing slope, due to the high topographic configuration factor and the fact that the surrounding surface is covered by either ice or snow, the surrounding-reflected radiation at P04 is extremely strong compared with that at the points on a horizontal surface. In addition, the diffuse radiation component is always the largest under the same atmospheric conditions due to the relatively high sky view factor. Therefore, the amount of global solar radiation received by point P04 is abnormally high throughout the mass-balance year and even   Figure 7b illustrates curves reflecting the spatiotemporal variation in the instantaneous solar irradiance at the six points over a mass-balance year. The seasonal variations in the incoming solar radiation are similar (sinusoidal) at the four points on the flat ice surface, but the other two points display completely different trends. The solar irradiance at point P05 is very low throughout the year, while that at point P04 is very high during the same period. Therefore, curves of the three DSSR components were generated to explore the reason for this disparity.
As shown in Figure 8, the amount of direct radiation received by point P05 is zero from 20 November 2017 to 22 January 2018 due to shade from the terrain. Global solar radiation is derived mainly from diffuse radiation, but this solar radiation component is very weak due to the small sky view factor. During the other seasons, although point P05 is illuminated by the sun, the direct radiation is still low because the cosine of the solar illumination angle at point P05 is smaller than those at the five other points. However, the spatiotemporal distribution of solar radiation reaching point P04 is completely different. Because point P04 is located on a steep southeast-facing slope, the cosine of the solar illumination angle is always very high, and thus, the direct solar radiation received at this point is very intense. Similarly, on this steep southeast-facing slope, due to the high topographic configuration factor and the fact that the surrounding surface is covered by either ice or snow, the surrounding-reflected radiation at P04 is extremely strong compared with that at the points on a horizontal surface. In addition, the diffuse radiation component is always the largest under the same atmospheric conditions due to the relatively high sky view factor. Therefore, the amount of global solar radiation received by point P04 is abnormally high throughout the mass-balance year and even  As shown in Figure 8 (b), two large troughs are observed in the 6 spatiotemporal variation curves of the diffuse radiation. Further research indicates that these troughs originate from high atmospheric transmittance over four days, which greatly reduces the amount of diffuse radiation received by the surface.

Impact of the DEM Spatial Scale on Topographic Normalization
DEMs are an important yet basic source of data for topographic normalization in areas composed of rugged terrain. However, the DEM grid size affects the topographic factors and ultimately affects the accuracy of the surface reflectivity inversion [44]. Figure 9 depicts the topographic normalization results for the S2 L1C image acquired on October 14, 2017, using DEMs with three resolutions of 90 m, 30 m, and 12 m; the corresponding spectral reflectance curves for point A located on a north-facing slope are also plotted. As shown in Figure 8b, two large troughs are observed in the 6 spatiotemporal variation curves of the diffuse radiation. Further research indicates that these troughs originate from high atmospheric transmittance over four days, which greatly reduces the amount of diffuse radiation received by the surface.

Impact of the DEM Spatial Scale on Topographic Normalization
DEMs are an important yet basic source of data for topographic normalization in areas composed of rugged terrain. However, the DEM grid size affects the topographic factors and ultimately affects the accuracy of the surface reflectivity inversion [44]. Figure 9  A visual inspection indicates that the effects of the terrain are most obvious in the original imagery of the rugged area, expressed as dark shadows in the shade and bright regions on the sunny slope. After removing these terrain effects, ice/snow features can be restored in more detail, and the surface reflectance inversion accuracy can be improved. Moreover, the higher the spatial resolution of the DEM, the better the removal of the terrain effect and the more reliable the surface reflectance inversion. In Sen2Cor, however, a 90 m SRTM DEM is generally used for the topographic normalization of S2 imagery by default. Therefore, the L2A products produced in this paper with a 12 m DEM can better remove the terrain effects in S2 imagery, thereby improving the reliability of an ice/snow surface reflectance inversion in a mass-balance year.

Terrain Impacts on the DSSR Received by Ice/Snow Surfaces
The solar irradiance of glacier surfaces with complex terrain is strongly influenced by various topographic factors. The most influential factor is the obstruction coefficient. When the obstruction coefficient is 1, the direct radiation received from the sun is zero; that is, the slope pixel cannot be irradiated.
The second largest influencing factor is the aspect. As shown in Table 4, at the S2A/B overpass time on September 1, 2017, the minimum DSSR is received by the northwest-, west-, and north-facing slopes on the glacier surface, while the maximum DSSR is received by the south-and southeast-facing slopes. Especially on the steep south-and southeast-facing slopes, due to the cosine of the solar illumination angle, the terrain configuration factor and the influence of the surrounding ice/snow coverage, the radiation reflected from the surrounding terrain is very strong. As a result, on the (2) A visual inspection indicates that the effects of the terrain are most obvious in the original imagery of the rugged area, expressed as dark shadows in the shade and bright regions on the sunny slope. After removing these terrain effects, ice/snow features can be restored in more detail, and the surface reflectance inversion accuracy can be improved. Moreover, the higher the spatial resolution of the DEM, the better the removal of the terrain effect and the more reliable the surface reflectance inversion. In Sen2Cor, however, a 90 m SRTM DEM is generally used for the topographic normalization of S2 imagery by default. Therefore, the L2A products produced in this paper with a 12 m DEM can better remove the terrain effects in S2 imagery, thereby improving the reliability of an ice/snow surface reflectance inversion in a mass-balance year.

Terrain Impacts on the DSSR Received by Ice/Snow Surfaces
The solar irradiance of glacier surfaces with complex terrain is strongly influenced by various topographic factors. The most influential factor is the obstruction coefficient. When the obstruction coefficient is 1, the direct radiation received from the sun is zero; that is, the slope pixel cannot be irradiated.
The second largest influencing factor is the aspect. As shown in Table 4, at the S2A/B overpass time on 1 September 2017, the minimum DSSR is received by the northwest-, west-, and north-facing slopes on the glacier surface, while the maximum DSSR is received by the south-and southeast-facing slopes. Especially on the steep south-and southeast-facing slopes, due to the cosine of the solar illumination angle, the terrain configuration factor and the influence of the surrounding ice/snow coverage, the radiation reflected from the surrounding terrain is very strong. As a result, on the surface of a glacier, steep sunny slopes receive very strong global solar radiation exceeding the solar constant received at the TOA.  Table 5 indicates that the average value of global solar radiation decreases as the slope increases. At the S2A/B overpass time on 1 September 2017, the minimum DSSR is received by slopes steeper than 25.1 • on the glacier surface, while the maximum DSSR is received by slopes with inclinations of 25.1-42.6 • . However, the standard deviations and the differences between the maximum and minimum DSSR indicate that the solar radiation is very heterogeneous within the same slope zone. In addition, the solar radiation received by slopes with inclinations in the range of 9.8-42.6 • on the glacier surface is higher than the solar constant due to the influence of the slope aspect.

Reliability of the Inversion of Atmospheric Transmittance Using S2 Data
In our previous research [30], it was found that MODIS atmospheric products overestimate aerosol and PW concentrations. Accordingly, atmospheric aerosol and PW products retrieved from S2 data have replaced MODIS atmospheric products. Therefore, these two types of atmospheric products need to be compared, and the reliability of the S2-estimated atmospheric transmittance needs to be further discussed. Figure 10a depicts the ground DSSR measurements and estimated atmospheric transmittances at the two AWSs at the S2 overpass time. The two curves are highly consistent, which can explain why the atmospheric transmittance obtained by the S2 inversion can accurately describe the atmospheric conditions. However, the MODIS aerosol products are all invalid in this region, and there are no aerosol ground observations. Hence, there is no way to verify the estimation accuracy of the aerosol contents from the inversion of S2 images, and it is impossible to compare the two aerosol products. Therefore, only the two types of atmospheric PW products are compared. Figure 10b indicates that the S2 atmospheric PW products can describe the atmospheric conditions with more reliable accuracy than the MODIS atmospheric products.

Limitations of the DSSR Estimation Method
In this study, S2 data are introduced into our previously published mountain radiation model [30] to estimate the DSSR on the surface of a glacier. The statistical results (MBE% = -1.7%, R 2 = 0.864) are higher than our published estimation accuracy (MBE% = -6.2%, R 2 = 0.505), which shows that the proposed method can accurately calculate the ice/snow surface irradiance for valley glaciers. However, the RMSD% (7.9%) is still not higher than the published estimation accuracy (RMSD% = 7.5%). Hence, the DSSR estimation approach for mountain glaciers proposed in this paper still suffers from some limitations.
First, the DSSR estimation accuracy depends on the reliability and spatial resolution of the DEM data employed. The DEM is the basic source of data for estimating the DSSR received by a mountain glacier surface. However, glacial DEMs change rapidly, so it is difficult to obtain timely and reliable DEM data. The DEM provided by the DLR cannot accurately describe the glacier between September 2017 and August 2018. In addition, a DEM with a 12 m spatial resolution cannot accurately describe the local terrain of a glacier.
Second, the proposed mountain radiation model cannot fully or accurately consider multiple scattering sources within a single pixel, which may introduce non-negligible uncertainty for glacier surfaces with high reflectivity. Moreover, this model ignores the multi-source scattering of radiation between the surface and atmosphere; nevertheless, incorporating this consideration will be the next research goal.
Third, the accuracy of solar irradiance verification is affected by the inconsistent footprint between the DSSR values retrieved from satellite imagery and in situ measurements caused by unstable ground-based observation instruments and the mismatched spatial scale. Additionally, there is considerable uncertainty in the retrieval of the AOD and this uncertainty ultimately affects the atmospheric transmittance estimation accuracy in high-reflection areas composed of ice/snow. When a scene contains no dark targets or medium-brightness reference pixels, the program Sen2Cor calculates the AOD with a user-defined visibility (e.g., 40 km).
In addition, the solar radiation received by a glacier surface is underestimated and overestimated simultaneously. In flat regions, the global radiation is underestimated; in contrast, in the southern region of the glacier composed of rugged terrain, the radiation reflected from mutually observed pixels becomes important, and the global radiation is overestimated due to the overestimation of the glacier albedo. This difficulty may be affiliated with the uncertainty in the narrow-to-broadband conversion of S2 satellite signals into glacier albedo. For more details, Greuell et al. [50] discussed the importance of six steps in this conversion process. In summary, although this article focuses on topographic normalization to eliminate (or at least reduce) terrain-induced BRDF

Limitations of the DSSR Estimation Method
In this study, S2 data are introduced into our previously published mountain radiation model [30] to estimate the DSSR on the surface of a glacier. The statistical results (MBE% = -1.7%, R 2 = 0.864) are higher than our published estimation accuracy (MBE% = -6.2%, R 2 = 0.505), which shows that the proposed method can accurately calculate the ice/snow surface irradiance for valley glaciers. However, the RMSD% (7.9%) is still not higher than the published estimation accuracy (RMSD% = 7.5%). Hence, the DSSR estimation approach for mountain glaciers proposed in this paper still suffers from some limitations.
First, the DSSR estimation accuracy depends on the reliability and spatial resolution of the DEM data employed. The DEM is the basic source of data for estimating the DSSR received by a mountain glacier surface. However, glacial DEMs change rapidly, so it is difficult to obtain timely and reliable DEM data. The DEM provided by the DLR cannot accurately describe the glacier between September 2017 and August 2018. In addition, a DEM with a 12 m spatial resolution cannot accurately describe the local terrain of a glacier.
Second, the proposed mountain radiation model cannot fully or accurately consider multiple scattering sources within a single pixel, which may introduce non-negligible uncertainty for glacier surfaces with high reflectivity. Moreover, this model ignores the multi-source scattering of radiation between the surface and atmosphere; nevertheless, incorporating this consideration will be the next research goal.
Third, the accuracy of solar irradiance verification is affected by the inconsistent footprint between the DSSR values retrieved from satellite imagery and in situ measurements caused by unstable ground-based observation instruments and the mismatched spatial scale. Additionally, there is considerable uncertainty in the retrieval of the AOD and this uncertainty ultimately affects the atmospheric transmittance estimation accuracy in high-reflection areas composed of ice/snow. When a scene contains no dark targets or medium-brightness reference pixels, the program Sen2Cor calculates the AOD with a user-defined visibility (e.g., 40 km).
In addition, the solar radiation received by a glacier surface is underestimated and overestimated simultaneously. In flat regions, the global radiation is underestimated; in contrast, in the southern region of the glacier composed of rugged terrain, the radiation reflected from mutually observed pixels becomes important, and the global radiation is overestimated due to the overestimation of the glacier albedo. This difficulty may be affiliated with the uncertainty in the narrow-to-broadband conversion of S2 satellite signals into glacier albedo. For more details, Greuell et al. [50] discussed the importance of six steps in this conversion process. In summary, although this article focuses on topographic normalization to eliminate (or at least reduce) terrain-induced BRDF effects in regions with extreme illumination geometries on the inversion of glacier albedo, there are still some shortcomings in the present study, such as the construction of the ice/snow BRDF and the surface classification (such as snow, debris, and dust), which will require further study in the future. Nevertheless, the overestimation bias cannot be assessed accurately due to the lack of ground observations on the steep southern slope.

Conclusions
The method used in this paper can accurately estimate the DSSR received by the surfaces of valley glaciers. The results show that the incoming solar irradiance of the surfaces of mountainous glaciers varies dramatically among different types of terrain and seasons. Meteorological ground-based site observations often do not represent the spatiotemporal distribution of solar radiation throughout a glacial region. Therefore, a mountain radiation model with satellite remote sensing data as the input atmospheric and surface parameters is an important tool for estimating the solar radiation received by the surfaces of mountain glaciers using high-precision DEM data. Additionally, albedo is an important input parameter for estimating the solar radiation of a mountain glacier because of its significant feedback to the reflected radiation.
The main contributions of this article lie in two aspects. One is to introduce the newly launched S2A/B satellite data into a previously published mountain radiation scheme, which is more conducive to the DSSR estimation of mountain glacier surface. The other is to modify the default parameter settings of Sen2Cor toolbox provided by ESA to obtain the improved L2A products than the standard products, thus reliable albedo of ice/snow cover can be retrieved. This study explored the seasonal cycle of solar irradiance in a mass-balance year. The verification results confirm that this mountain radiation scheme with S2 data as the model input parameters can provide reliable solar radiation information for mountain glacier surfaces and is not dependent on ground-based observations; thus, this approach is suitable for estimating the incoming solar irradiance of glaciers in remote and difficult-to-reach areas. This information could be helpful for better understanding the dynamic variations in solar shortwave radiation received across the surfaces of glaciers in a mass-balance year and this technique could provide basic data for research on glacial changes.
However, this study also has a shortcoming. Due to limitations in the acquisition of observation data at ice sites and the special requirements of glacial terrain and areas for quantitative estimation of solar radiation from mountain glaciers, only a single glacier was selected for the research. In the future, we will make efforts to cooperate with additional large-area valley glacier observation stations to improve the universality of the model estimation.
Author Contributions: Y.Z. designed the research, formulated the model, and wrote the paper. X.Q. helped with the experiments and provided site observations. X.L. and J.Z. helped construct the model and significantly contributed to the manuscript's revision. Y.L. organized and processed the site observations. All authors contributed to the discussion of the results. All authors have read and agreed to the published version of the manuscript.  surrounding-reflected irradiance (W m −2 ) cos i s cosine of solar illumination angle on a sloped grid V s obstruction coefficient V iso sky view factor F ij topographic configuration factor