*2.4. Albedo Upscaling*

According to Liang [10], unless the surface area is large and sufficiently homogeneous or a sufficient number of measurements is obtained during the period of satellite overpass, "point" measurements are not feasible to validate coarse-scale albedo products through a direct comparison. The in situ-based measurements are considered to be spatially representatively at homogenous land surfaces, and; therefore, can be directly used [11]. However, the in situ measurements and coarse-scale albedo are generally at an unmatched scale over heterogeneous land surfaces [24]. In order to validate satellite-derived coarse resolution albedo products from in situ measurements, upscaling from ground "point" data to the coarse resolutions is essential. This upscaling process consists of two steps. First, a "calibration factor" needs to be calculated based on the "point" measurements and the high-resolution EO pixels within the projected field-of-view (FoV) of the tower albedometer. Second, the high-resolution albedo products need to be upscaled to coarse-resolution using the pre-calculated "calibration factor".

High-resolution (<100 m) albedo products are not currently available from EO satellites. However, they can be derived using the methods described in [25], and in [26] for the Chinese HuanJing (HJ) satellite sensor. Their method used early machine learning (neural networks) based on numerous radiative transfer simulations. In this study, atmospherically corrected bi-directional reflectance factors (BRFs), from high-resolution EO, alongside coarse-resolution albedos, predicted from a MODIS BRDF climatology over a larger area, were employed as inputs. Firstly, pure endmembers for the corresponding area of the reference albedometer measurements were extracted based on the high-resolution spectral reflectance data. This was followed by aggregating the derived high-resolution spectral BRFs to coarse resolutions. In the next step, a linear regression was established between the endmember abundance and MODIS BRDF-derived albedo-to-nadir-reflectance ratios. This linear regression was established based on the data at coarse resolutions, but could be used to assign each high-resolution pixel around the tower with a corresponding albedo value. This was based on the assumption that over a larger area, coarse-resolution BRDFs can be represented by high-resolution mosaics of pixel BRDFs, weighted by their coverage proportions. Then, the high-resolution spectral reflectance values were converted to shortwave reflectance values through the use of narrowband to broadband conversion coefficients. In the end, a "calibration factor" was obtained by calculating the ratio of shortwave albedo derived from the albedometer and that from the high-resolution EO. This factor could then be applied to upscale the shortwave albedo from high resolutions to coarser resolution.

#### 2.4.1. Retrieval of High-Resolution Shortwave Albedos

We modified the method proposed by Shuai et al. [27] and used it to retrieve high-resolution shortwave albedos from Landsat-8 high-resolution BRFs and coarse-resolution MODIS BRDF climatology. The idea of Shuai's method is to first classify the spectral features of the land surface through unsupervised classifications, using 6 non-thermal bands of 30-m Landsat data. The Landsat data are then re-projected and aggregated to MODIS resolution, such that the albedo-to-nadir-reflectance ratios can be calculated for the "pure" pixels defined from the classification. Finally, the high-resolution albedo is produced for each of the Landsat pixel using the derived ratios. In Shuai's method, the detection of "pure pixels" from coarse-resolution pixels is an essential requirement for the subsequent calculation of albedo-to-nadir-reflectance ratios. However, this approach is unlikely to find "pure pixels" for each of the classes at a coarse resolution image, especially for heterogeneous land surfaces. In the modified retrieval method, the albedo-to-nadir-reflectance ratios still play a key role in deriving the high-resolution albedos, but the existence of "pure pixels" in coarse resolutions is not necessarily required.

The 1-km MODIS BRDF climatology parameters were produced from the 500-m MCD43A1 products [28]. The surfaces reflectances derived from this kernel-driven BRDF model are described as:

$$R(\lambda, \theta\_{\rm in}, \theta\_{\rm out}, \phi) = f\_{\rm iso}(\lambda) + f\_{\rm vol}(\lambda)k\_{\rm vol}(\theta\_{\rm in}, \theta\_{\rm out}, \phi) + f\_{\rm gen}(\lambda)k\_{\rm gas}(\theta\_{\rm in}, \theta\_{\rm out}, \phi) \tag{4}$$

where *λ* is the bandpass of a given spectral channel; *θin*, *θout* and *φ* are the solar zenith, view zenith and relative azimuth angles, respectively. *k* is the BRDF RossThick–LiSparse–Reciprocal (RTLSR) kernel and *f* is the spectrally-dependent kernel weighting, with subscripts *iso*, *vol* and *geo* representing the isotropic, volumetric and geometric-optical components, respectively. Integration of the BRFs over all view angles results in a DHR, and a further integration over all illumination angles results in a BHR:

$$\mathrm{DHR}\_{\mathrm{M}}(\boldsymbol{\theta}\_{\mathrm{in}}(\mathrm{L}\_{\mathrm{\mathsf{S}}})) = \frac{1}{\pi} \int\_{0}^{2\pi} \mathrm{d}\varphi \int\_{0}^{1} \mathrm{R}\_{\mathrm{M}}(\lambda, \,\theta\_{\mathrm{in}}(\mathrm{L}\_{\mathrm{\mathsf{S}}}), \,\theta\_{\mathrm{out}}, \,\varphi) \mathrm{u}\_{\mathrm{V}} \mathrm{d}\mathrm{u}\_{\mathrm{V}} \tag{5}$$

$$\text{BHR}\_{\text{M}} = \frac{1}{\pi} \int\_0^{2\pi} \text{d}\,\text{q} \int\_0^1 \text{DHR}\_{\text{M}}(\Theta\_{\text{in}}(\text{L}\_8)) \text{u}\_{\text{s}} \text{d}\mathbf{u}\_{\text{s}} \tag{6}$$

where *uv*(=sin *θout*) and *us*(=sin *θin*) are the variables of integration. The shortwave BRFs at MODIS 1-km resolution for the Landsat-8 solar zenith and view zenith angle are given by

$$\mathcal{R}\_{\mathcal{M}}(\Omega(\mathcal{L}\_8)) = \mathcal{R}\_{\mathcal{M}}(\lambda, \: \theta\_{\text{in}} = \theta\_{\text{in}}(\mathcal{L}\_8), \: \theta\_{\text{out}} = \theta\_{\text{out}}(\mathcal{L}\_8), \: \\*\*p = \varphi(\mathcal{L}\_8)) \tag{7}$$

where <sup>Ω</sup>(*<sup>L</sup>*8) is the Landsat-8 sun and sensor geometry. Then, the ratios between the shortwave albedo and shortwave reflectance values can be computed for all the pixels at 1-km MODIS resolution:

$$a\_D = \frac{DHR\_M(\theta\_{in}(L\_8))}{R\_M(\Omega\_L)};\ a\_B = \frac{BHR\_M}{R\_M(\Omega\_L)}\tag{8}$$

An endmember is defined as a land "type" that is assumed to have a unique spectral signature. Here, the N-FINDR endmember extraction algorithm [29] was adopted to extract the pure endmembers from the 30-m Landsat-8 spectral reflectance data. This was followed by re-projecting the Landsat-8 data from Universal Transverse Mercator (UTM) to MODIS (sinusoidal) projection, and aggregating the pixels from 30-m to 1-km resolution. Then the proportion of each endmember was calculated for each of the aggregated 1-km pixel using a fully constrained least squares (FCLS) linear un-mixing method [30]. Then, the following equation was established:

$$\mathbf{A}\,\,\mathbf{W} = \mathbf{R} \tag{9}$$

where **A** is a (m, n) matrix with m being the number aggregated pixels, and n being the number of endmembers. Each row of the matrix **A** contains the proportions of derived endmembers. **W** is a (n, 1) matrix containing the weighting parameters. **R** is a (m, 1) matrix with the elements representing the MODIS BRDF climatology albedo-to-nadir-reflectance ratios. The weighting function **W** is solved as,

$$\mathbf{W} = \left(\mathbf{A}^{\mathrm{T}} \,\mathbf{A}\right)^{-1} \,\mathbf{A}^{\mathrm{T}} \,\mathbf{R} \tag{10}$$

where the superscript T refers to matrix transpose. Given the Landsat-8 spectral reflectance values, the shortwave reflectance values can be calculated through the use of a set of narrowband to broadband conversion coefficients [18]:

$$
\alpha^{L\_8} = 0.356 \mathfrak{a}\_2 + 0.13 \mathfrak{a}\_4 + 0.37 \mathfrak{a}\_5 + 0.085 \mathfrak{a}\_6 + 0.072 \mathfrak{a}\_7 - 0.0018 \tag{11}
$$

where *α*2, *α*4, *α*5, *α*6 and *α*7 represent the Landsat-8 blue, red, NIR, SWIR-1 and SWIR-2 narrow bands, respectively. The shortwave reflectance is transformed to shortwave albedo through the following formula:

$$\mathbf{B} = (\mathbf{C} \ \mathbf{W}) \circ \mathbf{L} \tag{12}$$

where **C** is (k, n) matrix with k being the number of 30-m Landsat-8 pixels. Each row of **C** contains the derived endmember proportions of Landsat-8 pixels. **L** is a (k, 1) matrix that contains the shortwave reflectance derived from Equation (11) for each pixel, and ◦ is the Hadamard product. The processing chain for generating high-resolution albedo using Landsat-8, as an example, is illustrated in Figure 5.

#### 2.4.2. Upscaling of Albedo from Tower to Coarse Resolutions

The pyranometers that measure downwelling and upwelling radiation are mounted on towers with a fixed height of 10 m at the SURFRAD sites, while the BSRN and FLUXNET sites utilize towers at different heights, usually dependent on the height of the canopy-top. The reference albedo was located by assuming the albedometer measures a circular area from the top of the tower [31]. The diameter of this circular area, which represents the effective projected FoV of the tower albedometer, was estimated as:

$$\text{ID} = 2 \tan(FoV^\circ / 2) \cdot (h\_{\text{tower}} - h\_{\text{ToC}}) \tag{13}$$

where *htower* and *hToC* are the height of tower and averaged height of vegetation, respectively. *FoV*/2 is half of the effective field of view in degrees, which is 81◦ [32]. For a pyranometer mounted on a 10-m tower, the projected FoV on the surface is 126 m. The values for all the sites are given in Table 1. The reference albedo was approximated by averaging the corresponding albedo values of pixels within a high-resolution EO shortwave albedo product. Then, a "calibration factor" could be derived from the ratio between the in situ albedo and the reference albedo. To produce the coarse-resolution albedos, the high-resolution albedo product needed to be aggregated to a coarse resolution first, and then modified with this "calibration factor". The above-mentioned process for producing coarse-resolution albedo product is illustrated in Figure 6.

**Figure 5.** Processing chain for calculating high-resolution albedo from Landsat BRFs and MODIS BRDFs.

**Figure 6.** Steps of upscaling high-resolution albedo to coarse-resolution pixel size. The process of calculating high-resolution EO shortwave albedo is illustrated in Figure 5.
