Simulation and Analysis of the Topographic Effects on Snow-Free Albedo over Rugged Terrain

Abstract

Topography complicates the modeling and retrieval of land surface albedo due to shadow effects and the redistribution of incident radiation. Neglecting topographic effects may lead to a significant bias when estimating land surface albedo over a single slope. However, for rugged terrain, a comprehensive and systematic investigation of topographic effects on land surface albedo is currently ongoing. Accurately estimating topographic effects on land surface albedo over a rugged terrain presents a challenge in remote sensing modeling and applications. In this paper, we focused on the development of a simplified estimation method for snow-free albedo over a rugged terrain at a 1-km scale based on a 30-m fine-scale digital elevation model (DEM). The proposed method was compared with the radiosity approach based on simulated and real DEMs. The results of the comparison showed that the proposed method provided adequate computational efficiency and satisfactory accuracy simultaneously. Then, the topographic effects on snow-free albedo were quantitatively investigated and interpreted by considering the mean slope, subpixel aspect distribution, solar zenith angle, and solar azimuth angle. The results showed that the more rugged the terrain and the larger the solar illumination angle, the more intense the topographic effects were on black-sky albedo (BSA). The maximum absolute deviation (MAD) and the maximum relative deviation (MRD) of the BSA over a rugged terrain reached 0.28 and 85%, respectively, when the SZA was 60° for different terrains. Topographic effects varied with the mean slope, subpixel aspect distribution, SZA and SAA, which should not be neglected when modeling albedo.

1. Introduction

Land surface albedo, defined as the fraction of incident solar radiation (0.3–3 μm) reflected by land surfaces [1,2], is one of the most significant geophysical variables affecting the Earth’s climate and controlling the surface radiation budget. It plays a crucial role in a variety of models, including general circulation models, land surface climate models, energy balance models, hydrology models, and biosphere models. Land surface albedo under ambient light conditions, also known as blue-sky albedo, is a combination of directional hemispheric reflectance, known as black-sky albedo (BSA), and bihemispherical reflectance, known as white-sky albedo (WSA); blue sky albedo takes into account the proportion of diffused skylight illumination and the solar zenith angle [3,4,5,6]. Over recent decades, land surface albedo remote sensing estimation algorithms have been developed, demonstrating that a bihemispherical integration method using the bidirectional reflectance distribution function (BRDF) [7,8,9] has a robust performance and is widely used in albedo estimation. This method for estimating albedo generally assumes that the land surface terrain is flat and homogeneous [10,11], and albedo products have been created with this method using different satellite datasets, such as the Moderate Resolution Imaging Spectroradiometer (MODIS) [12], the polarization and directionality of Earth reflectances (POLDER) [13], the multi-angle imaging spectroradiometer (MISR) [14], and the Clouds and the Earth’s Radiant Energy System (CERES) [15].

Spatial heterogeneities, which are comprised of complex topography and heterogeneous land cover, complicate the estimation of land surface albedo [16]. Directly applying albedo estimation methods that are suitable for flat terrain to a rugged terrain leads to large errors [17,18]. Actually, topography plays different roles in albedo estimation at different spatial scales [19]. For a single slope, physical-based BRDF models, such as the improved four-scale geometric-optical model for sloping terrain (GOST) [20], the improved soil-leaf-canopy radiative transfer model for sloping terrain (SLCT) [21], the improved Li-Strahler geometric-optical canopy model for sloping terrain (GOMST) [8], the vertical vegetation model (VVM) [7], the path length correction (PLC) model [22], have been developed based on the radiative transfer principle and the geometric optical theory, which depend on the orientation of the plant stand and the particular configuration of the sun direction and the terrain slope [7,8]. By integrating the BRDF over the exitance hemisphere for a single irradiance direction, the BSA for a single slope can be obtained; by integrating the BRDF over all viewing and irradiance directions, the WSA for a single slope can be calculated. Investigations have shown that the albedo for a single slope is related to the slope and aspect of the single slope [7,8,23,24]. With an increase in slope, the albedo becomes sensitive to the aspect of the slope, and the slopes facing away from the sun may display larger albedos than those of the sunward facing slopes due to increased mutual shadowing [8]. In addition, the terrain shadowing and diffused radiation from the adjacent slopes significantly influences the single slope albedo [8]. For a rugged terrain, the topographic effects on albedo generally focus on the integrated effects caused by subpixel slopes within one remote sensing pixel [18,19]. Neglecting the subpixel topography variability in albedo estimation over a rugged terrain leads to significant deviations [17,19,25], which can reach a relative error of 33% for a mean slope of 40° [19]. It has been shown that MODIS albedo retrievals are also highly sensitive to subpixel topography, and the MODIS albedo over a rugged terrain can change up to 100% for spruce vegetation in winter [26].

Albedo depends on both land surface characteristics and the atmosphere. Topography alters land surface characteristics and solar illumination geometry. Thus, the topographic effects on albedo over a rugged terrain are related to the spatial distribution characteristics of the subpixel topography, the solar zenith angle (SZA), and the solar azimuth angle (SAA) [17,19,25]. Compared to a single slope, it is difficult to investigate the variation in albedo over a rugged terrain by integrating the effects of subpixel slopes. The lack of a rigorous and effective physical BRDF/albedo model for rugged terrains has contributed to this challenge. Wen et al. [19] developed a land surface albedo estimation and scale correction method over a rugged terrain by the hemispheric integration of surface reflectance over a rugged terrain, where the rugged terrain reflectance was estimated based on the subpixel reflectance. Considering that the analysis of topographic effects on land surface albedo required huge amounts of typical albedo data over the rugged terrain to ensure the reliability of the analysis results, this method was inconvenient for the large-scale simulation of surface albedo because the rugged terrain reflectances under the entire hemispheric view space were required during each albedo calculation under different SZAs and SAAs. Therefore, we dedicated ourselves to developing a simplified method to estimate the albedo over a rugged terrain directly by the subpixel albedo in this paper based on the same idea in Wen et al. [19].

In this paper, we focused on quantitatively investigating topographic effects on snow-free albedo over a rugged terrain based on the developed method. The BSA was related to both the topography and the solar illumination geometry. The shadows induced by adjacent topographies immensely affected the BSA. However, the WSA was independent of solar illumination geometry and shadow effects. Therefore, to emphasize the topographic effects on albedo, BSA was selected instead of WSA in this paper. This paper was organized as follows. First, a BSA estimation method for rugged terrain was proposed in Section 2. Section 3 described the DEM datasets and the generation method of the reference BSA dataset. The proposed BSA estimation method was validated in Section 4.1; Section 4.2 quantitatively investigated topographic effects on snow-free BSA given variations in mean slope, subpixel aspect distribution, SZA, and SAA. Finally, a conclusion was provided in Section 5.

2. Methods

2.1. BSA over a Rugged Terrain

The rugged terrain, which is comprised of a number of subpixel slopes with different slopes and aspects, is shown in Figure 1a, where it is assumed that the entire terrain is horizontal. The SZA and SAA are denoted by ${\theta }_s$ and ${\varphi }_s$, respectively, with respect to the horizontal plane.

According to the definition of albedo [1], the BSA ${\rho }_{BSA\_low}({\theta }_s,{\varphi }_s)$ over a rugged terrain under a clear sky is calculated by:

$${\rho }_{BSA\_low}({\theta }_s,{\varphi }_s)=\frac{{\mathsf{\Phi }}_e^b}{{\mathsf{\Phi }}_i^b}$$

(1)

where ${\mathsf{\Phi }}_e^b$ and ${\mathsf{\Phi }}_i^b$ represent the reflected and incident radiation flux over the whole pixel under a clear sky, respectively. The incident radiation flux over a rugged terrain is equal to the incident radiation flux received by the corresponding projected horizontal pixel. Therefore, ${\mathsf{\Phi }}_i^b$ is equal to

$${\mathsf{\Phi }}_i^b=E_s\cos {\theta }_s{\displaystyle \int \cos a_{j}d{A}_{tj}}$$

(2)

where $E_s$ represents the direct solar irradiance independent of topography, ${\displaystyle \int \cos a_{j}d{A}_{tj}}$ represents the area of the projected horizontal pixel, $a_j$ denotes the slope of the subpixel slope, $d{A}_{tj}$ denotes the incremental surface area of the subpixel slope and the subscript $tj$ is the jth subpixel slope.

2.2. BSA Estimation Method Derivation

The subpixel slope BSA upscaling approach, based on the radiosity theory [27], is a feasible scheme to estimate the BSA over a rugged terrain. However, with the upscaling approach, many parameters are required, including reflectance characteristics, vegetation structural parameters, illumination conditions, and the radiation of each subpixel slope, which results in more uncertainties and computational complexities. The parameterized idea adopted by Wen [19] is an alternative method to estimate BSA due to its simplicity and efficiency.

This idea assumes that a virtual slope with slope $\alpha$ and aspect $\beta$ exists, where incoming and outgoing radiation are the same as the sum of those over a rugged terrain, as shown in Figure 1. Thus, ${\mathsf{\Phi }}_e^b$ can be expressed as:

$${\mathsf{\Phi }}_e^b=E_s\cos i_s^e{\rho }_{BSA\_\mathrm{eq}}(i_s^e,{\phi }_s^e)A_e({\theta }_s,{\varphi }_s)$$

(3)

where ${\rho }_{BSA\_\mathrm{eq}}(i_s^e,{\phi }_s^e)$ represents the BSA of the virtual slope; $i_s^e$ and ${\phi }_s^e$ represent the relative SZA and SAA corresponding to the virtual slope, respectively; and $A_e$ denotes the area of the virtual slope surface, which depends on both the subpixel topography distribution and the solar illumination geometry.

According to the principle from the mountain radiative transfer theory [28], under a specific solar illumination geometry and by neglecting the multi-scattering effects from an adjacent terrain, the incremental reflected radiation flux of the jth subpixel slope, $d{\mathsf{\Phi }}_{ej}^b$, is given by

$$d{\mathsf{\Phi }}_{ej}^b={\mathsf{\Theta }}_{sj}{E}_s\cos i_{sj}{\rho }_{BSA\_high}(i_{sj},{\phi }_{sj})d{A}_{tj}$$

(4)

where $\mathsf{\Theta }$ represents the shadow factor, which is set to 1 for an illuminated slope and 0 otherwise [29,30,31] and can be calculated using the ray-tracing method. Thus, ${\mathsf{\Theta }}_{sj}$ indicates whether or not the subpixel slope is sunlit. Variables $i_{sj}$ and ${\phi }_{sj}$ represent the relative SZA and SAA, respectively, corresponding to the subpixel slope, and ${\rho }_{BSA\_high}(i_{sj},{\phi }_{sj})$ represents the BSA of the subpixel slope.

${\rho }_{BSA\_high}(i_{sj},{\phi }_{sj})$ is affected by topographic obstructions because several parts of the reflected radiation from the outgoing hemisphere are obstructed. Variable $d{\mathsf{\Phi }}_{ej}^b$ can be approximately calculated as follows:

$$d{\mathsf{\Phi }}_{ej}^b=d{\tilde{\mathsf{\Phi }}}_{ej}^b{V}_j$$

(5)

where $d{\tilde{\mathsf{\Phi }}}_{ej}^b$ denotes the reflected radiation flux neglecting obstructions from the adjacent terrain, and $V_j$ represents the sky view factor, which can be calculated using the DEM [32]. Sky view factor $V_j$ represents the unobstructed portion of the sky at a given location and ranges between 0 and 1. A value of $V$ close to 1 indicates that almost the entire hemisphere is unobstructed and visible, which is the case for exposed features, such as planes and peaks; values close to 0 are present in deep sinks and lower regions of deep valleys, where almost no sky is visible [33].

Therefore, ${\rho }_{BSA\_high}(i_{sj},{\phi }_{sj})$ is equal to:

$${\rho }_{BSA\_high}(i_{sj},{\phi }_{sj})=\frac{d{\tilde{\mathsf{\Phi }}}_{ej}^b{V}_j}{d{\mathsf{\Phi }}_{ij}^b}={\tilde{\rho }}_{BSA\_high}(i_{sj},{\phi }_{sj})V_j$$

(6)

where ${\tilde{\rho }}_{BSA\_high}(i_{sj},{\phi }_{sj})$ denotes the BSA neglecting obstructions from an adjacent terrain. Thus, the sum of the subpixel slope radiation fluxes over a rugged terrain, ${\mathsf{\Phi }}_e^b$ is as follows:

$${\mathsf{\Phi }}_e^b=E_s{\displaystyle \underset{A(s)}{\int }\cos i_{sj}{\tilde{\rho }}_{BSA\_high}(i_{sj},{\phi }_{sj})V_{j}d{A}_{tj}}$$

(7)

where $A_{(s)}$ denotes the subpixel slopes that are illuminated by the sun. Combining Equations (3) and (7), we obtain

$${\displaystyle \underset{A(s)}{\int }\cos i_{sj}{\tilde{\rho }}_{BSA\_high}(i_{sj},{\phi }_{sj})V_{j}d{A}_{tj}}=\cos i_s^e{\rho }_{BSA\_\mathrm{eq}}(i_s^e,{\phi }_s^e)A_e({\theta }_s,{\varphi }_s)$$

(8)

To focus on the topographic effects on BSA, land cover type within the rugged terrain is assumed to be homogeneous, and the differences among ${\tilde{\rho }}_{BSA\_high}(i_{sj},{\phi }_{sj})$ for different subpixel slopes are only caused by different solar illumination geometries and DEM characteristics. ${\tilde{\rho }}_{BSA\_high}(i_{sj},{\phi }_{sj})$ has an identical function form with that of ${\rho }_{BSA\_\mathrm{eq}}(i_s^e,{\phi }_s^e)$, except that the latter input angle parameters are different.

To unify the symbols and simplify the deduction, Equation (8) is substituted with

$${\displaystyle \underset{A(s)}{\int }Y(u_j,v_j)V_{j}d{A}_{tj}}=Y(u_e,v_e)A_e(u,v)$$

(9)

where $Y(u_j,v_j)=\cos i_{sj}{\tilde{\rho }}_{BSA\_high}(i_{sj},{\phi }_{sj})$, $Y(u_e,v_e)=\cos i_s^e{\rho }_{BSA\_\mathrm{eq}}(i_s^e,{\phi }_s^e)$, $u_j=\cos i_{sj}$, $v_j=\cos {\phi }_{sj}$, $u_e=\cos i_s^e$, $v_e=\cos {\phi }_s^e$, $u=\cos {\theta }_s$, and $v=\cos {\varphi }_s$. To construct the virtual slope, we obtain the specific formulas for $u_e$, $v_e$ and $A_e$, which determine if the virtual slope exists. An alternative method based on the Taylor expansion, which is similar to the derivation strategy of the Hapke shadow function [34,35], was used to solve Equation (9). Specifically, $Y$ is assumed to be mathematically well behaved, and $Y$ is expanded on both sides of Equation (9) in a Taylor series about $u$ and $\mathrm{v}$:

$$\begin{array}{l}{A}_e(u,v)[Y(u,v)+\frac{\partial Y}{\partial u}(u,v)(u_e-u)+\frac{\partial Y}{\partial v}(u,v)(v_e-v)+\dots ]\\ =Y(u,v){\displaystyle \underset{A(s)}{\int }{V}_{j}d{A}_{tj}}+\frac{\partial Y}{\partial u}(u,v){\displaystyle \underset{A(s)}{\int }(u_j-u)V_{j}d{A}_{tj}}+\frac{\partial Y}{\partial v}(u,v){\displaystyle \underset{A(s)}{\int }(v_j-v)V_{j}d{A}_{tj}}+\dots \end{array}$$

(10)

Since $u$ and $v$ are independent variables, and $Y$ is an arbitrary function of $u$ and $v$, Equation (10) is satisfied only if the coefficients of $Y$ and its partial derivatives are separately equal on both sides of the equation. Neglecting the higher order terms of the Taylor expansion for $Y$, we obtain Equation (11):

$$\begin{array}{l}{A}_e(u,v)={\displaystyle \underset{A(s)}{\int }{V}_{j}d{A}_{tj}}\\ A_e(u,v)={\displaystyle \underset{A(s)}{\int }(u_j-u)V_{j}d{A}_{tj}}/(u_e-u)\\ A_e(u,v)={\displaystyle \underset{A(s)}{\int }(v_j-v)V_{j}d{A}_{tj}}/(v_e-v)\end{array}$$

(11)

Solving Equation (11), $u_e$, $v_e$, and $A_e$ can be specifically and respectively be formulated as:

$$\begin{array}{l}{A}_e(u,v)={\displaystyle \underset{A(s)}{\int }{V}_{j}d{A}_{tj}}\\ u_e=\frac{{\displaystyle \underset{A(s)}{\int }{u}_j{V}_{j}d{A}_{tj}}}{{\displaystyle \underset{A(s)}{\int }{V}_{j}d{A}_{tj}}}\\ v_e=\frac{{\displaystyle \underset{A(s)}{\int }{v}_j{V}_{j}d{A}_{tj}}}{{\displaystyle \underset{A(s)}{\int }{V}_{j}d{A}_{tj}}}\end{array}$$

(12)

By combining Equations (1), (2), and (12), we obtain:

$${\rho }_{BSA\_low}({\theta }_s,{\varphi }_s)=\frac{\cos i_s^e{\rho }_{BSA\_\mathrm{eq}}(i_s^e,{\phi }_s^e)A_e({\theta }_s,{\varphi }_s)}{\cos {\theta }_s{\displaystyle \int \cos a_{j}d{A}_{tj}}}$$

(13)

By introducing Equation (12) into Equation (13), we obtain

$${\rho }_{BSA\_low}({\theta }_s,{\varphi }_s)={\rho }_{BSA\_\mathrm{eq}}(i_s^e,{\phi }_s^e)\frac{{\displaystyle \underset{A(s)}{\int }\cos i_{sj}{V}_{j}d{A}_{tj}}}{\cos {\theta }_s{\displaystyle \int \cos a_{j}d{A}_{tj}}}$$

(14)

Equation (14) shows that the BSA over a rugged terrain can be obtained by the product between the BSA of the virtual slope and a specified factor. Its discrete formula is written as

$${\rho }_{BSA\_low}({\theta }_s,{\varphi }_s)={\rho }_{BSA\_\mathrm{eq}}(i_s^e,{\phi }_s^e)\frac{{\displaystyle \sum _{k=1}^N{\mathsf{\Theta }}_{sk}\cos i_{sk}{V}_k/\cos a_k}}{N\cos {\theta }_s}$$

(15)

where $N$ represents the number of the subpixel slopes within the pixel, and $k$ denotes the kth subpixel slope.

Equations (13) and (14) indicate that this method combines well with any single slope BSA estimation method to estimate the BSA over a rugged terrain. In this paper, the combined PROSPECT leaf optical property model and the SAIL canopy bidirectional reflectance model, also referred to as PROSAIL, is selected as a the single slope surface reflectance and the BSA estimation method. The PROSAIL model has been widely used to study plant canopy spectral reflectance [36] and is relatively mature and efficient. For an infinitely inclined homogeneous vegetation canopy, the topographic influences on the anisotropic reflectance simulated by the PROSAIL model can be categorized into two aspects: photon path length alteration inside the vegetation layer and the adjustment of the extinction coefficient. The first effect is handled by a simple geometric correction of the solar-terrain sensor and the consideration of vertical tree growth [18,24]. Second, a spherical leaf inclination distribution function (LIDF) is assumed, which allows the effective extinction coefficient for a unit path length to be fixed in all directions, where the topographic effects on the extinction coefficient can be neglected. Given the specific vegetation parameter and the terrain configuration, we can use Equations (13) and (14) and the PROSOIL model under topographic considerations to estimate the BSA over a rugged terrain.

2.3. Topographic Effect Analysis Methods

The topographic effects on BSA are influenced by the spatial distributions of the subpixel slopes, SZA and SAA. A local sensitivity analysis method is used in this paper to analyze the effects of these different factors. This method estimates the effect of a single factor on the outputs while maintaining the other factors at their nominal values [37]. The maximum absolute deviation (MAD) and the maximum relative deviation (MRD) are used to quantitatively analyze these sensitivities:

$$MAD=\underset{i=1}{\overset{N}{Max}}\left\{a_{model}^i\right\}-\underset{i=1}{\overset{N}{Min}}\left\{a_{model}^i\right\}$$

(16)

$$MRD=\left[\underset{i=1}{\overset{N}{Max}}\left\{a_{model}^i\right\}-\underset{i=1}{\overset{N}{Min}}\left\{a_{model}^i\right\}\right]/\underset{i=1}{\overset{N}{Max}}\left\{a_{model}^i\right\}$$

(17)

where $N$ represents the amount of estimated BSAs, and $i$ represents the BSA index. Variable $a_{model}^i$ denotes the ith estimated BSA calculated by the proposed method in this paper.

3. Datasets

3.1. Simulated DEM Dataset

To simulate and evaluate the BSA over a rugged terrain, nine simulated DEMs with 30 m spatial resolution were generated to provide various magnitudes of roughness [19]. The statistical mean values of the elevation, slope and sky view factor are listed in

Table 1. Basic parameters for nine simulated DEMs.

Filter	Exaggeration	Mean Elevation (m)	Mean Slope (°)	Mean Sky View Factor
5 × 1	1	29.49	1.33	1.00
3 × 1	1	29.49	1.88	1.00
1 × 1	1	29.49	2.51	1.00
5 × 1	10	294.89	12.84	0.95
3 × 1	10	294.89	17.59	0.90
1 × 1	10	294.89	22.62	0.85
5 × 1	20	589.78	23.70	0.83
3 × 1	20	589.78	30.90	0.74
1 × 1	20	589.78	37.72	0.64

. For this research, only the central 1 × 1 km area (i.e., 33 × 33 grid cells) of the aforementioned DEMs was used to avoid calculation errors at the edge of the DEM.

3.2. Global Digital Elevation Model (GDEM)

With an average elevation exceeding 4500 m, the Tibetan Plateau (Figure 2) is the largest and highest plateau in the world. This plateau is characterized by high spatial heterogeneity due to the presence of mountainous areas [38]. Therefore, it is an ideal region to compare and validate the developed method in this paper. A 100 km × 100 km test area, with elevations ranging from 2492 m to approximately 6769 m, is identified from the plateau. The DEM shown in the upper-right part of Figure 2, with a 30-m spatial resolution of the study area, was collected by the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) of the global digital elevation model, version 2 (ASTER GDEM2). The terrain over the southern parts of the study area was relatively complex and rugged, whereas the northern terrain was relatively gentle. This study area, with a high spatial heterogeneity, was suitable for the study of topographic effects on BSA.

To effectively investigate the topographic effects on BSA, three typical categories for 1-km real DEMs, with typical surface fluctuations, were selected from the study area: gentle slope terrain (10°), moderate slope terrain (20°), and steep slope terrain (30°). Each DEM category was comprised of six representative DEMs with different subpixel aspect distributions. Their slopes and aspect distributions are shown in Figure 3, where the code names of the DEMs (e.g., dem-10-1) are marked in the legend. The subpixel aspect distributions of the six DEMs from 1-km DEM category with a mean slope of 10° had the following characteristics, subsequently: dominant southwest-oriented distribution, relatively uniform distribution, dominant southwest-oriented distribution, relatively uniform distribution, dominant eastward-oriented distribution, and dominant northward-oriented distribution. The DEMs with a mean slope of 20° had the following characteristics, subsequently: dominant northward-oriented distribution, dominant eastward-oriented distribution, relatively uniform distribution, dominant southwest-oriented distribution, dominant northward-oriented distribution and dominant southward-oriented distribution. The DEMs with a mean slope of 30° had the following characteristics, subsequently: dominant northwest-oriented distribution, dominant southwest-oriented distribution, dominant southwest-oriented distribution, dominant southeast-oriented distribution, relatively uniform distribution and dominant eastward-oriented distribution.

3.3. Reference BSA Dataset Simulation Based on the Radiative Approach

Considering that it is difficult to obtain the albedo reference over a rugged terrain, and current land albedo products have poor accuracies over a rugged terrain, the radiosity approach [39,40,41], which is a widely used computer simulation model, was used to generate the reference BSA to evaluate the performance of the proposed method. In this paper, the terrain was described by the DEM, and the land cover was assumed to be homogeneous in the simulation scenario.

Specifically, the procedure included three steps: the anisotropic reference reflectance simulation based on the radiosity approach, the spectral BSA calculation by integrating the reflectance over the hemispheric exitance given an illumination direction, and the broad-band shortwave BSA calculation by integrating the spectral BSA weighted by the incident radiation. In the first step, the reflectance characteristics of each subpixel slope were acquired based on the PROSAIL model under topographic consideration and upscaled to reflectance values over a rugged terrain based on the radiosity approach. In this paper, considering that the BSA of soil has a similar variation with topography as that of vegetation, only the vegetation BSA variation was analyzed; the input parameter specifications in the PROSAIL model are shown in

Table 2. Specification of input parameters in the PROSAIL model.

Model Parameters		Unit	Range
Leaf parameters	Leaf structure index	unitless	1.5
	Leaf chlorophyll content	[μg/cm2]	40
	Leaf dry matter content	[g/cm2]	0.009
	Leaf water content	[cm]	0.01
	Leaf brown pigment	[g/cm2]	0.0
Soil parameters	Reflectance	------	Shown in Figure 4a
Canopy structure parameters	LAI	[m2/m2]	3
	Leaf inclination distribution function	------	Spherical
	Hot spot size parameter	[m/m]	0.01
Atmospheric condition	Incoming radiation	------	Shown in Figure 4b
Illumination view geometry	Solar zenith angle	[°]	0–90 at a 5° interval
	Solar azimuth angle	[°]	0–360 at a 5° interval
	View zenith angle	[°]	0–90 at a 5° interval
	View azimuth angle	[°]	0–360 at a 5° interval
Terrain	DEM	------	Nine simulated DEMs and real DEMs

and Figure 4.

4. Results and Discussion

4.1. Modeled BSA Accuracy Assessment

To analyze the accuracy of the modeled BSA with the proposed method, a simulation scenario was constructed based on both the nine simulated DEMs and three real DEM categories. For the simulation scenarios, the realistic tree shape parameters, component signatures, and other parameters are listed in Table 2. Then, the modeled BSAs are estimated by the proposed method in Section 2.2, and the reference BSAs are generated by the method in Section 3.3. The determination coefficient (R²), root mean square errors (RMSE), mean absolute percentage error (MAPE) and mean bias (Bias) are adopted to evaluate the accuracy of the proposed method.

Figure 5a shows the comparison results of the modeled and reference BSAs over nine simulated DEMs. It is found that these two types of BSAs are close in magnitude, indicating that a majority of data points are distributed around the 1:1 line with a small RMSE (0.0060), MAPE (0.0038) and Bias (0.0038). This demonstrates that the modeled BSA is consistent with the reference BSA.

Table 3. Accuracy statistics of the modeled BSA over different simulated DEMs.

Mean Slope (°)	R2	RMSE	MAPE	Bias
1.33	1.0000	0.0002	0.0001	0.0001
1.88	1.0000	0.0003	0.0002	0.0002
2.51	1.0000	0.0006	0.0003	0.0003
12.84	0.9985	0.0059	0.0037	0.0037
17.59	0.9961	0.0071	0.0050	0.0050
22.62	0.9873	0.0078	0.0059	0.0059
23.70	0.9893	0.0078	0.0059	0.0059
30.90	0.9782	0.0078	0.0065	0.0065
37.72	0.9734	0.0074	0.0066	0.0066

presents the error statistics of the BSA for different simulated terrains. The RMSE varies from 0.0002 to 0.0074, and the MAPE increases from 0.0001 to 0.0066 as the mean slope increases from 1.33° to 37.72°, which indicates that with an increase in mean slope, the accuracy of the proposed method gradually decreases but is still acceptable. This can be explained by the increase in mean slope, which causes the multi-scattering effects and the terrain obstruction effects for the reflected radiation to become increasingly obvious; however, the proposed method neglects multi-scattering effects and approximately considers the terrain obstruction effects. These results confirm the capability of the proposed method for estimating BSA with simulated DEMs.

The real DEMs provide a great number of terrains with different subpixel slope distributions, which are closer to the natural characteristics of the terrain than those of the simulated DEMs with normal distributions. Figure 5b indicates that the modeled BSA shows an adequate performance, with a high R² close to 1, a low RMSE of 0.0038, a MAPE of 0.0023 and a Bias of 0.0023.

Table 4. Accuracy statistics of the modeled BSA over different real DEMs.

Mean Slope (°)	R2	RMSE	MAPE	Bias
10	0.9998	0.0020	0.0011	0.0011
20	0.9994	0.0030	0.0023	0.0023
30	0.9963	0.0056	0.0036	0.0036

presents the error statistics under different real terrain conditions. Overall, the results obtained from the comparison with the real DEMs are similar to those of the simulated DEMs. Therefore, the proposed method provides a high-quality BSA dataset and can be applied to investigate the topographic effects on BSA.

4.2. Topographic Effects on BSA

4.2.1. Factors Influencing the BSA over a Rugged Terrain

Real DEMs were used to investigate the topographic effects over a rugged terrain because they provide a sufficient range of topographies that cover a wide variety of natural terrain characteristics. Figure 6 shows the hemispheric BSAs modeled with different solar illuminations over a flat terrain and the three typical real 1-km DEMs (i.e., dem-10-1, dem-20-1, and dem-30-1), which indicate that BSA distributions with solar illumination geometries are intensively affected by topography. Figure 6a shows that the BSA over a flat terrain monotonously increases with SZA regardless of the SAA value. When the terrain is relatively flat in dem-10-1, several minor changes occur in the shape of the BSA distribution, as shown in Figure 6b. The BSA increases with SZA when the SZA is substantially less than 90° and has a weak relationship with the SAA, but the BSA decreases when the SZA is close to 90° due to the terrain block of incident radiation. With an increase in mean slope for dem-20-1, the shape of the BSA distribution evidently changes and becomes asymmetric, as shown in Figure 6c. The north-facing slope terrains are dominant in dem-20-1, which indicates that the BSA is dependent on the SAA. The BSA increases monotonously with the SZA and is relatively large when the SAA is oriented north. In comparison, the BSA first increases then decreases with the SZA and becomes small when the SAA is oriented south due to the existence of shadows. As the terrain becomes more rugged in dem-30-1, the BSA generally decreases, and the shape of the BSA distribution clearly changes. This phenomenon can be explained by the fact that shadow effects are more obvious in dem-30-1 and substantially influence the BSA. Thus, we conclude that the subpixel slope distribution and the solar illumination geometry are two important factors influencing BSA, and BSA variations with varying illumination angles present different trends under different terrain conditions. Specifically, the mean slope, subpixel aspect distribution, SZA, and SAA are four main controlling factors for quantitatively analyzing the topographic effects on BSA.

4.2.2. BSA Variation with Mean Slope

The 1-km DEMs over the study area were used to investigate the subpixel effects of the mean slope on the BSA (shown in Figure 7). When the SZA is 0°, a shadow does not exist, but the spatial variations in the subpixel slope and aspect within the rugged terrain affect the BSA. The MAD and MRD of the BSA reach 0.08 and 36%, respectively. However, when the SZA is 30°, shadows occur, and the topographic effects on the BSA are further enhanced. The MAD and MRD of the BSA increase to 0.12 and 50%, respectively. When the SZA is 60°, the MAD and MRD of the BSA reach 0.28 and approximately 85%, respectively, due to more obvious shadow effects and terrain obstruction effects.

Generally, BSA presents a decreasing trend with an increase in mean slope. When the SZA is 0°, the BSA clearly decreases with an increase in mean slope because several parts of the outgoing reflected radiation in the hemisphere are obstructed even though there is not a shadow present. However, when the SZA is 30°, and the SZA is 60°, the BSA variation with an increase in mean slope is maintained, but it is not as obvious. This is because BSA is affected by various factors, such as SZA, SAA, shadows, terrain obstruction, and the slope and aspect distribution of the subpixel slopes. When the terrain is gently rugged (20°), the shadow and occlusion effects of the terrain may be weak, but the alteration in regional illumination angle caused by the topography may result in an increase in BSA, as shown in Figure 7b,c. The regional SZA of the subpixel slope facing the sun is small, whereas that of the subpixel slope facing away from the sun can be relatively large, which indicates that the BSA is low when sunward subpixel slopes account for the majority of the data. When the terrain is steep (30°), the effects of shadows and terrain obstruction play a dominant role, which means that many subpixel slopes cannot be illuminated by the sun, and many parts of the outgoing reflected radiation in the hemisphere are obstructed. Therefore, the BSA decreases significantly regardless of the SZA, as shown in Figure 7.

4.2.3. BSA Variations with Sub-Pixel Aspect Distributions

Three typical 1-km DEM categories were used to investigate the effects of subpixel aspect distribution on BSA. Without loss of generality, the SAA is set to 150° in this analysis. For each category, the subpixel aspect distributions of the six DEMs differ significantly. Figure 8 shows the BSA variation with the subpixel aspect distribution. When the SZA is close to 0°, the influence of the subpixel aspect distribution on BSA under different terrains is minimal because shadowing does not occur in this scenario. The influence of the subpixel aspect distribution increases with SZA regardless of the mean slope of the terrain. When the mean slope is 30°, the MAD of the BSA increases from 0.01 to approximately 0.15 with an increase in SZA from 0° to 60°. This is because shadows gradually occur as the SZA increases, and the subpixel aspect distribution is related to the shadow ratio and distribution. When the terrain is relatively flat, the influence of the subpixel aspect distribution on the BSA is minimal, as shown in Figure 8a. With an increase in mean slope, the influence of the subpixel aspect distribution on the BSA gradually increases. The MADs of the BSAs over the terrain with mean slopes of 10°, 20°, and 30° are 0.06, 0.08, and 0.15, respectively. Overall, with an increase in SZA and mean slope, the influences of the subpixel aspect distribution gradually increase. In addition, the influence of the subpixel aspect distribution has a strong relationship with the SAA, which will be discussed in Section 4.2.5.

Based on specific solar illumination geometries and mean slopes, the topographic effects on the BSA are dependent on the spatial distribution characteristics of the subpixel aspect. When the proportion of the subpixel slope facing the sun is high, the BSA is relatively large, whereas the BSA becomes small when the proportion of the subpixel slope facing away from the sun is high. For instance, when the SZA is 60°, and the SAA is 150° (i.e., the sun is to the southeast), the BSA of the dominant eastward-oriented dem-10-5 is at a maximum in Figure 8a, the BSA of the dominant southward-oriented dem-20-6 is at maximum in Figure 8b, and the BSA of the dominant eastward-oriented dem-30-6 is at a maximum (and that of the dominant northwest-oriented dem-30-1) is at a minimum in Figure 8c.

4.2.4. BSA Variation with SZA

Both in situ measurements and remote sensing data from satellite and aircraft platforms have shown that BSA is strongly dependent on SZA [42,43,44,45]. Figure 9b,d show the BSA simulations with different SZAs, which exhibits obvious spatial variation characteristics corresponding to the mean slope distribution (shown in Figure 9a). As the SZA increases, the dynamic range of BSA becomes gradually large due to the modulation in the regional illumination angle and shadow effects. When the SZA is 0°, BSA varies from 0.14 to 0.22. In this case, shadows do not exist because the sun is at nadir, and the variation in BSA is mainly caused by the alteration in the regional illumination angle due to the inclined terrain. When the SZA is 30°, BSA varies from 0.12 to 0.24. The range increases slightly due to the existence of shadows and the increase in SZA. When the SZA is 60°, the BSA range, which becomes even larger, varies from 0.05 to 0.33. The low value of BSA (close to 0) can be explained by the fact that shadow effects are serious in some circumstances when a majority of the subpixel slopes face away from the sun. The high BSA value can be attributed to the high proportion of subpixel slopes facing toward the sun, where few shadows exist, even though the SZA is large.

Furthermore, three typical 1-km DEM categories are used to analyze the variations in topographic effects based on SZA. The variations in BSA based on SZA for different terrains are shown in Figure 10. When the terrain is flat (Figure 10a), BSA increases monotonously as the SZA increases. When the terrain is rugged (Figure 10b,d), BSA first increases then decreases with an increase in SZA. When the mean slope is 30°, the MAD and MRD of the BSA over rugged terrains caused by different SZAs exceed 0.22 and approximately 88%, respectively. This is because when the SZA is close to 90°, the shadow effect is obvious. The SZA at the inflection points gradually decreases from 70° to approximately 50° as the mean slope increases from 10° to 30°, respectively. These inflection points can be attributed to the increase in mean slope, which causes shadows to more easily appear when the SZA is relatively small. It is concluded that the increasing trend in BSA with SZA gradually slows as the mean slope increases.

4.2.5. BSA Variation with SAA

Figure 11 shows the variation in BSA based on DEMs with different mean slopes with a specified SAA for different SZAs. These DEMs are dem-10-3, dem-20-3, and dem-30-3, whose subpixel aspect distributions exhibit a dominant southwest-orientation, a relatively uniform distribution, and a dominant southwest-orientation, respectively. The influence of the SAA on BSA increases with SZA. For a mean slope of 30°, the MAD of the BSA with a 30° SZA is approximately 0.03, while the MAD of the BSA with a 60° SZA exceeds 0.07. The SAA has little impact on BSA when the terrain is relatively flat, whereas the SAA considerably affects the BSA when the terrain is rugged. The influence of the SAA on BSA gradually increases with an increase in mean slope. Figure 11d shows that as the mean slope increases from 0° to 30°, the MAD of BSA varies from 0 to approximately 0.07, which can be explained because the shadow effects are obvious as the SZA and the mean slope increase simultaneously.

Furthermore, the variations in BSA with SAA under different subpixel aspect distributions are also analyzed. Figure 12 shows the variation in BSA with SAA over the different terrains with different aspect distributions when the SZA is 30°. It is concluded that there is a parabolic distribution between BSA and the SAA. When the SAA is close to the predominant aspect of the subpixel slopes, the BSA is relatively large, as shown in the variation curves with SAA in Figure 11 and Figure 12. This is because the shadow area is relatively small in this situation. When the SAA is opposite that of the predominant aspect at the end of the parabolic distribution, the BSA reaches a minimum due to the large shadow ratio. For terrains in dem-10-2, dem-10-4, dem-20-3 and dem-30-5, which have relatively uniform aspect distributions, the SAA has little effect on BSA. When the aspect is unevenly distributed, the parabolic distribution characteristics between the SAA and BSA are distinct. For dem-30-1, the northwest-facing slopes comprise the majority; therefore, the maximum and minimum BSAs appear in the northwestern and southeastern directions, respectively, which is where the sun is located. The BSA variations in dem-30-2 and dem-30-3 with SAA coincide due to the analogous aspect distributions between these two DEMs. These results demonstrate that BSA over a rugged terrain is sensitive to SAA, and the subpixel aspect distribution has considerable influence on the relationship between the SAA and the BSA.

4.3. Method Limitations

The accuracy assessment shows that with an increase in mean slope, BSA estimation errors gradually increase due to model simplification and approximations in problems such as multi-scattering and terrain obstruction. However, in general, the multi-scattering effect from adjacent terrains has a weaker contribution to BSAs, unless the slope lies in deep valleys or the adjacent terrains have high reflectance [28,46]. Therefore, considering that multi-scattering effect is negligible for snow-free land surface and the developed method has high precision in the validation experiments, the proposed method can have an adequate performance in vegetation and soil albedo estimations; these problems may have little effect on the analysis of topographic effects on BSA. Unfortunately, the developed method is not applicable to snow covered land surfaces with high reflectance. In addition, due to the model assumption of homogeneous land cover, the developed method faces a challenge in the mixed-pixel region. Subsequent research will focus on solving these problems.

This paper mainly uses modeled BSA data with the proposed method to analyze the topographic effects on BSA. The developed BSA estimation method was validated against the radiosity approach. Therefore, developing an efficient collection method for in situ albedo data over rugged terrains and validating the developed method against in situ data are urgent. In addition, in practical applications, DEM quality and the mismatching between the DEM and the remote sensing imagery also affect the accuracy of albedo estimation with the developed method. These basic data processing issues need to be further discussed and addressed, although they are beyond the scope of this paper.

5. Conclusions

Neglecting topographic effects may lead to significant bias when estimating land surface albedo. In this paper, we presented an efficient snow-free BSA estimation method over a rugged terrain. The proposed method was validated using simulated DEMs and real DEMs. The validation results showed that the modeled BSAs were consistent with the reference BSA, with an RMSE smaller than 0.01, which confirmed the ability of the proposed method to estimate BSA (with acceptable errors). By comparing the modeled BSA using the proposed method over the real DEM scenario, the topographic effects on BSA were investigated in detail.

The BSA over a rugged terrain is influenced by the subpixel slope distribution (mean slope and subpixel aspect distribution) and the solar illumination angle (SZA and SAA). These factors are related to the modulation in the regional illumination angle and shadow effects, which play key roles in the topographic effects on BSA. The more rugged the terrain and the larger the solar illumination angle, the more obvious topographic effects are on BSA. Specifically, for the subpixel slope distribution, the mean slope has a higher influence on BSA than that on the subpixel aspect distribution. For the mean slope, BSA generally presents a decreasing trend with an increase in mean slope. The larger the SZA, the more obvious the decreasing trend in BSA is with an increase in mean slope. When analyzing the subpixel aspect distribution, for an increase in mean slope, the influence of the subpixel aspect distribution on BSA gradually increases. Under specific solar illumination geometries and mean slopes, the influences of the subpixel aspect distribution on BSA are dependent on the proportion of subpixel slopes facing the sun. For the solar illumination angle, the SZA has a greater impact on BSA than that from the SAA. When the terrain is relatively flat, BSA increases monotonously with an increasing SZA. As the terrain becomes rugged, BSA first increases then decreases with an increase in SZA due to the incident radiation terrain block. For the SAA, BSA over a rugged terrain is sensitive to the SAA, and the influence of SAA on BSA increases with an increasing SZA and mean slope. A parabolic distribution is found between BSA and the SAA. When the SAA is close to the predominant aspect of the subpixel slopes, BSA is relatively large.

The motivation and findings in this study can benefit land surface albedo modeling and retrieval in the field of remote sensing. Subsequent research will focus on practical remote sensing applications and method improvements by considering multi-scattering effects.