**3. Methods**

Our data assimilation methodology included an EnKF which recursively updated the albedo by coupling the direct estimation approach with a dynamic model. A flow chart of this process is shown in Figure 2. The existing high quality multi-year MODIS albedo product data is averaged to compute albedo climatology, then a simple dynamic model is established based on the climatology to evolve albedo over time and to forecast short-range albedos. A direct estimation approach is used to generate high-resolution TM albedo data. The EnKF technique is used to estimate real-time albedos by combining the predictions from the dynamic model and the high-resolution TM data.

The proposed method was essentially a three-step process:

(1) Obtain MODIS albedo data for the historical period of interest, correct them geometrically, convert black-sky albedo (BSA) and white-sky albedo (WSA) to blue-sky albedo using the scatter ratio, resample to a 30 m resolution, and average the albedo of the historical period to obtain the time-series shortwave albedo climatology;

(2) Transform the cloudless TM image into a 30 m spatial resolution TM shortwave albedo via direct estimation algorithm; and

(3) Use the EnKF method with the MODIS time series albedo as the background field and input the TM high spatial resolution albedo data to estimate the high spatial resolution albedo.

This method effectively integrated the time change information from MODIS data and spatial Landsat data to resolve issues with the low spatial resolution, sparsity, and long return cycle of MODIS data, as well as the sparsity of Landsat data to ultimately obtain high spatial and temporal resolution albedo data.

**Figure 2.** Real-time albedo inversion based on ensemble Kalman Filter (EnKF).

#### *3.1. Validation Site Land Surface Heterogeneity*

Previous researchers have verified the consistency of global land surface albedo products against ground measurement data at MODIS 500 m spatial scales at SURFRAD, Atmospheric Radiation Measurement Southern Great Plains (ARM/SGP) [34], and Cloud and Radiation Testbed-Southern Great Plains (CART/SGP) sites [35]. The MODIS albedo product has high accuracy with RMSE of 0.013–0.018. However, these verifications are based on homogeneous sites using point observations. In heterogeneous surface areas, site observations do not represent the MODIS pixels area, and thus the remote sensing products cannot be fully verified. As such, we must analyze site representativeness before validation with remote sensing data [36]. The heterogeneity of surface space describes the extent of surface heterogeneity within a certain range and reflects the heterogeneity in the area around FluxNet sites [37,38]. In this study, we used variogram model parameters and the relative coefficient of variation (RCV) to analyze the effect of heterogeneity in the area around FluxNet sites according to Roman [39].

The isotropic spherical variogram model was used to evaluate landscape heterogeneity with 1.0 km2, 1.5 km2, and 2.0 km<sup>2</sup> TM albedo subsets. The TM albedo was obtained using the direct estimation approach proposed by He [18]. The spherical model formula is:

$$\gamma\_{\rm sph}(h) = \begin{cases} \mathbf{c}\_0 + \mathbf{c} \cdot \left( 1.5 \cdot \frac{\mathbf{h}}{\mathbf{a}} - 0.5 \left( \frac{\mathbf{h}}{\mathbf{a}} \right)^3 \right) & \text{for } 0 \le h \le \mathbf{a} \\\mathbf{c}\_0 + \mathbf{c} & \text{for } h > \mathbf{a} \end{cases} \tag{3}$$

The range a defines the distance from a point beyond which there is no further correlation of a biophysical property associated with that point. It has also been described as the average patch size of the landscape when the curve reaches the high-level distance [40]. The data are correlated within the range of a. Otherwise, the data are not correlated with each other; that is, the observations outside the range do not affect the estimation results. For a 1.0 km<sup>2</sup> subset amax = 690 m, for a 1.5 km<sup>2</sup> subset amax = 1050 m, and for a 2.0 km<sup>2</sup> subset amax = 1410 m. The sill c is the variation of ordinates when the abscissa is greater than the variable range. The nugge<sup>t</sup> c0 is the variation of the abscissa at 0, which describes the variability of the data at the microcosmic level.

The coefficient of variation (CV) is the ratio of the standard deviation to the mean. It is independent of spatial scale and can provide an estimate of the overall variability of data. The relative CV is:

$$\mathrm{R\_{CV}} = \frac{\mathrm{CV\_{1.5\chi} - \mathrm{CV\_X}}}{\mathrm{CV\_X}} ; \mathrm{x} = 1.0 \text{ km}^2 \tag{4}$$

where Rcv is the relative CV, CV1.5x is the CV calculated from a 1.5 km<sup>2</sup> TM albedo subset at the center of a given observation station, and CVx is the CV calculated from a 1.0 km<sup>2</sup> TM albedo subset at the center of the given observation station. If a site is more representative, it has a weaker spatial heterogeneity, a similar landscape around the site, and Rcv closer to 0.

#### *3.2. Albedo Dynamic Model Construction*

Albedo climatology is determined by the average of albedo over historical period. In this study, we calculated albedo climatology from high-quality multi-year MODIS albedo data as follows:

$$ALB\_t^{\prime} = \frac{1}{n} \sum\_{k=1}^{n} ALB\_t(k) \tag{5}$$

where *ALBt* is the albedo value corresponding to time t on the background field curve, *ALBt*(k) is the MODIS albedo value at time t, and n is the year.

The albedo climatology describes the general change tendency of albedo in a one-year period. There may be some deviation from the ground measurements due to precipitation, anthropogenic activities, and other reasons. Albedo climatology simulates the general trend of land surface albedo and supplies the albedo estimation background. When new observations are available, the EnKF updates the estimation to produce a final value. We constructed our dynamic model based on the climatology used to forecast the short-range albedo:

$$ALB\_t = F\_t \times ALB\_{t-1} \tag{6}$$

where *ALBt* represents the current estimated albedo, *ALBt*−<sup>1</sup> represents the estimated albedo at the preceding time step, and *Ft* is:

$$F\_t = 1 + \frac{1}{ALB\_t + \varepsilon} \times \frac{dALB\_t}{dt} \tag{7}$$

where ε = 10−<sup>3</sup> prevents negative denominators and *dALBt dt* is the growth rate of albedo at time t. The dynamic model was adopted from Samain et al. [41] and Xiao et al. [42].

## *3.3. High-Resolution Albedo Estimation*

Surface albedo can be obtained directly from Top of Atmosphere (TOA) reflectance without requiring atmospheric correction [43]. The apparent surface albedo was separated from the inherent surface albedo, and the empirical relationship between the simulated surface albedo and TOA reflectance was built based on extensive radiative transfer simulations under a variety of atmospheric conditions. To mitigate errors from nonlinearities, a new statistical relationship was extended to generate albedo measurements from the spectral bands and three broadbands [18], including the visible (300–700 nm), near-infrared (NIR; 700–3000 nm), and total shortwave ranges (300–3000 nm), based on the following linear equation:

$$\mathbf{x}\_{\lambda} = \sum \mathbf{p}\_{i}^{TOA} \cdot \mathbf{c}\_{i} + \mathbf{c}\_{0} \tag{8}$$

where *αλ* is the surface albedo for the spectral range of *λ*, ρTOA i is the TOA reflectance for spectral band i, and ci and c0 are the regression coefficients.

The direct estimation approach simulates TOA reflectance from radiative transfer information. Here, the sensor spectral response functions obtained from the USGS and the MODIS BRDF database were used as inputs for the 6S radiative transfer code, aerosol types, water vapor, ozone, and CO2. The linear regression coefficients for each of the geometrical combinations were precalculated and stored in LUTs for operational use. TOA reflectance was simulated using the regression models in the LUT. Cloud-free observations were selected from the quality control document, then the TM surface albedo was obtained from TOA reflectance data. The direct estimation approach has been successfully applied to Landsat Multispectral Scanner (MSS), Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+), and Operational Land Imager (OLI) sensors [18]; extensive validations have also been carried out at SURFRAD, AmeriFlux, BSRN, and GC-Net sites. The direct estimation approach can generate reliable surface albedo estimates with accuracy of 0.022 to 0.034 in terms of RMSEs over snow-free surfaces [44,45]. The direct estimation approach has also been used in OLI and Gaofen-1 satellite (GF-1) surface albedo inversion [46,47].

### *3.4. Ensemble Kalman Filter*

In the 1960s, Kalman [48,49] and Bucy [48] developed the Kalman filter algorithm for optimal estimation of system states according to a linear system state equation and combined input and output system observation data. The traditional Kalman filter algorithm is often applied to linear problems [50,51]. It can estimate the state of a dynamic system from a series of data with measurement noise when the measurement variance is known, and is a common component in communication, navigation, guidance, and control applications. The dynamic model and observer in a Kalman filter are typically expressed as follows:

$$
\omega\_{k+1} = M\_k \mathbf{x}\_k + \omega\_{k+1} \tag{9}
$$

$$z\_k = H\_k \mathbf{x}\_k + \upsilon\_k \tag{10}$$

where *xk* is the dynamic model state vector; *Mk* ∈ *Rn*×*n* is the system matrix, which changes with time; *zk* represents observed values, and *Hk* is a time-varying observation system which transforms the state variable into the same value as the observed value. In this study, the observed values and state variable were land surface albedos. *ωk* is the model error and *υk* is observed noise. The albedo has different levels of inversion accuracy over time.

The predicted value of the state variable of the dynamic model at time k is defined as *x fk* . The data assimilation step serves to calculate the optimal estimate *xak* of the state variables at moment k by combining the predicted value *x fk* and the observed value *zk* of the model state variables at time k, *<sup>x</sup>ak*, the optimal estimate of *xk*, is a linear function of *x fk*and *zk*.

$$\mathbf{x}\_k^a = \mathbf{x}\_k^f + \mathbf{K}\_k[\mathbf{z}\_k - H\_k\mathbf{x}\_k^f] \tag{11}$$

where *zk* − *Hkx fk* is observed incrementally (innovation) and *Kk* is the Kalman gain matrix. Data assimilation is conducted to minimize the variance of *xak* by determining *Kk*. The update equations for state variables and covariance are:

$$\mathbf{x}\_k^a = \mathbf{x}\_k^f + \mathbf{K}\_k[\mathbf{z}\_k - H\_k\mathbf{x}\_k^f] = \mathbf{x}\_k^f + \mathbf{P}\_k^f H\_k^T (H\_k \mathbf{P}\_k^f H\_k^T + \mathbf{R}\_k)^{-1} (\mathbf{z}\_k - H\_k \mathbf{x}\_k^f) \tag{12}$$

$$\mathbf{P}\_k^d = \mathbf{P}\_k^f - \mathbf{K}\_k \mathbf{H}\_k \mathbf{P}\_k^f = \mathbf{P}\_k^f - \mathbf{P}\_k^f \mathbf{H}\_k^T \left(\mathbf{H}\_k \mathbf{P}\_k^f \mathbf{H}\_k^T + \mathbf{R}\_k\right)^{-1} \mathbf{H}\_k \mathbf{P}\_k^f \tag{13}$$

where *xak* is the updated state variable value, *x fk* is the predicted state variable, *Pfk* is the predicted covariance, and *Pak*is the posterior covariance of state variables.

The ensemble Kalman filter algorithm is a complex sequential data assimilation method which produces non-linear models within a Kalman gain scheme. In the ensemble Kalman filter algorithm, x is defined as an n-dimensional model state vector and *A* = (*<sup>x</sup>*1, *x*2, ... *xN*) ∈ *n*×*<sup>N</sup>* is composed of a collection of N model state vectors. The ensemble mean is stored in each column of *A* ∈ *n*×*N*. The ensemble perturbation matrix is defined as:

$$A' = A - \overline{A} \tag{14}$$

The ensemble covariance matrix is:

$$P\_{\mathfrak{c}} = \frac{A' \left(A'\right)^{T}}{N - 1} \tag{15}$$

Given a vector of measurements *d* ∈ *<sup>m</sup>*, where m is the number of measurements, the observation matrix is expressed as follows:

$$D = (d\_1, d\_2, \dots, d\_N) \in \mathfrak{R}^{m \times N} \tag{16}$$

Then, the standard analysis equation of the Kalman filter is:

$$A^a = A + A'A'^T H^T \left( HA'A'^T H + \mathbb{R} \right)^{-1} (D - HA) \tag{17}$$

where *H* ∈ *n*×*<sup>N</sup>* (m is the number of measurements) is the measurement operator relating the model state to the observations; *R* ∈ *m*×*m* is the observational error covariance matrix, and *D* ∈ *m*×*<sup>N</sup>* is the disturbance observation matrix set. H is a linear operator which is not suitable in cases when the operator is nonlinear. By augmenting the model state vector, *x<sup>T</sup>* = [*xT*, *<sup>h</sup><sup>T</sup>*(*x*)]. *Y* ∈ *n*×*<sup>N</sup>* is defined as the set matrix of augmented state vectors and *Y* ∈ *n*×*<sup>N</sup>* is the set disturbance matrix of state variables. As such, the analytical equation is:

$$A^a = A + A' \overline{Y'}^T \overline{H'}^T \left(\overline{H'} \overline{Y'} \overline{Y'}^T \overline{H'}^T + R\right)^{-1} \left(D - \overline{H} \overline{Y}\right) \tag{18}$$

where *H* is a new observation operator. The analytical equation can solve the data assimilation problem for which the observation operator is nonlinear.

In this study, we calculated the background error covariance from albedo background and field observations. We calculated observation error covariance values from TM albedo and field observations over the time series described above. At each time point, we generated N random noises with 0 as the mean and 0.01 as the variance. N is the ensemble number, which affects estimation efficiency and accuracy: A larger N has a lower calculation efficiency and higher accuracy. We set the ensemble number to 100 and added random noise values to the background field values of the current date to obtain N new background field values [52]. Finally, we calculated the standard deviation between the N background field values and the ground field observations as the background field errors of the current date. Observation errors were simulated in the same way. Figure 3 shows the error setting for the background field, TM albedo, and field observation for the AU-DaP site in 2008. The "background error" and "TM error" labels represent simulated background errors and TM albedo errors, respectively. For regional use, the error of the area was set based on the error of the center pixel; a normal distribution noise with a mean of 0 and a variance of 0.01 was added for each pixel.

**Figure 3.** Error setting for albedo background and TM data.

Surface albedo is the only model state here, so the size of the model state vector n is equal to 1. We used the albedo value derived by the direct estimation approach as the albedo observation data, so the sum of the size of the model state vector n and the number of measurement equivalents added to the original model state vector was equal to 2. The simulated albedo was included in the state variables of the dynamic model at the given observation time. We emphasize the observation error covariance R here as per its significant influence on the Aa value. In our experiments, all the errors were simulated at site AU-DaP.
