# An Evapotranspiration Assimilation Method Based on Ensemble Kalman Filter and À Trous Wavelet

Authors: Shaohui Chen; Jianwei Qi; Xiaomin Sun; Xiangzheng Deng; Jing Tian
Journal: Advances in Meteorology
 (2013)

---

## Abstract

It is challenging to assimilate the evapotranspiration product (EP) retrieved from satellite data into land surface models (LSMs). In this paper, a perturbed ensemble Kalman filter (PEKF) and à trous wavelet transform (AWT) integrated method are proposed to implement the evapotranspiration assimilation. In this method, the AWT is used to decompose the EPs into multiple channels since it is very powerful in fusing high frequency spatial information of multisource data, and then the Kalman filter is performed in the AWT domain. The proposed method combines the advantages of the PEKF that is capable of accommodating model error and observation error, and the AWT can effectively perform multiresolution fusion. Assimilation experiment conducted with the Noah model and the EP retrieved from the MODIS data shows that the proposed method performs better than the traditional ensemble Kalman filter (EnKF) and PEKF methods. The analysis results fit well with the evapotranspiration observation at two field sites with different land surface conditions. These indicate that the proposed method is promising for assimilating regional scale satellite retrieved EP into LSMs.


---

## Body

## 1. Introduction

Evapotranspiration (ET) is an important component of the water and energy exchanges between the atmosphere and land surface. It is crucial to accurately estimate ET for studying global or regional water and energy balances. Hence, good quality of spatial and temporal ET production (EP) can help to improve comprehension of water and energy cycle. However, this kind of EP is generally difficult to obtain in both dimensions of space and time because ET is influenced by many factors, such as air and skin temperatures, soil moisture, vegetation fraction, and horizontal advection. Up to now, there are two approaches to estimate the ET. One is site observations or remote sensing retrievals. Site observations have high spatial resolutions, but can only provide the EP for limited spatial locations [1]. Remote sensing retrievals have high spatial resolutions and can cover large range, but can only retrieve the instantaneous EP. The other is land surface models (LSMs). LSMs are probably the most efficient approach for continuously estimating ET on a large range [1]. Because of the imperfection of the physics of LSMs and the uncertainties of input and driving data, the EP of the LSMs may contain significant errors. Hence, data assimilation (DA) has been applied to integrate observational ET into LSMs [2].DA provides a framework for improving the LSMs by updating the state variables of the LSMs (SVLs) with observations and can combine the high spatial resolution of the observation with the high temporal resolution of the LSM. DA can be realized by two kinds of schemes: continuous assimilation and sequential assimilation. In continuous assimilation, the SVLs are modulated to be close to the observations. Continuous assimilation methods may cause the abrupt change of the SVLs before and after the DA, which will make the subsequent simulation of the LSMs to easily produce obvious errors. In sequential assimilation, the SVLs are updated according to some forecast principle. The update is usually completed by adding the product of the gain that is estimated from model errors and observation errors and the difference between the observations and the outputs of the LSM to the SVLs. Sequential assimilation methods can produce a statistically optimal and dynamically consistent state estimate of the land surface by considering observation errors and model errors. Among sequential assimilation methods, Kalman filter- (KF-) derived methods yield accurate and consecutive ET estimate and have been applied widely in recent years [3].Up to now, many KF-derived DA methods have been proposed. The KF method explicitly computes the error covariances through an additional matrix equation that propagates error information from one update time to the next, subject to possibly uncertain model dynamics [4, 5]. When the LSMs are linear, the KF method not only can get the optimal estimate of the SVLs, but also can estimate the error of the optimal estimate. However, the LSMs are usually highly nonlinear, and in order to apply the KF to the nonlinear LSMs, the KF method is extended to extended KF (EKF) method by using the first derivative of Taylor formula to linearize the LSMs. For the complicated, discontinuous, and nonlinear LSMs, the performance of these methods derived from the EKF is not very steady [6]. In order to overcome the limits of the KF- and EKF-based methods, Evensen [7, 8] proposed ensemble KF (EnKF) according to random dynamic forecast theory proposed by Epstein [9]. Burgers et al. [10] improved the EnKF to perturbed EnKF (PEKF) by randomly perturbing the SVLs, the observations, and the forcings using respective uncertainties. Because the DA methods derived from the EnKF and PEKF avoid the high dimensional nonlinear matrix operation, it is easy to integrate the DA methods into the LSMs. Hence, this branch of the DA methods has become an active research front in recent years.Qin et al. [11] developed a variational data assimilation scheme based on the weak constraint concept. Actually, it assimilated surface skin temperature into a simple land surface model for estimation of ET. It is inconvenient to use the automatic differentiation technique which derived the adjoint codes to evaluate the gradient of the cost function when one had the ET retrievals. Pipunic et al. [12] confirmed that assimilating remotely sensed latent and sensible heat flux could potentially produce better heat flux predictions than assimilating soil moisture, or skin temperature observations from which they are derived. However, this experiment is only tested using synthetic data, and it needs further study about how the assimilation of ET using the actual remotely sensed ET is completed. Jang et al. [13] used MODIS data to calculate ET during clear sky conditions, while the MODIS-MM5 four-dimensional data assimilation system provided input variables for the calculation of ET under cloudy sky conditions, which means that the two types of ET are not merged in fact from the meaning of DA. Xu et al. [14] estimated turbulent fluxes through assimilation of geostationary operational environmental satellites data using the EnKF method. Though the use of geostationary operational environmental satellites data tackles the problem of remotely sensing data sparseness, the potential application of the experiment is constrained by its execution at site scale rather than domain scale. French et al. [15] forecasted spatially distributed cotton ET by assimilating remotely sensed and ground-based observations. However, this method could not be feasible to provide continuous ET estimates that are better than can be achieved with either data alone if temporally continuous point observations are not available. However, more research is required to determine if and how well assimilation of ET can improve heat flux, soil moisture, and soil temperature predictions from the LSMs.The main purpose of the DA is to combine the complementary information from measurements and models of the LSMs into an optimal estimate of the geophysical fields of interest [4]. The idea is similar to data fusion. By taking the fusion of low spatial and high spectral resolution multispectral data (LRMD) and high spatial and low spectral resolution panchromatic data (HRPD) as an example, the purpose of data fusion is to inject the details of the HRPD into the LRMD with the spectral properties of the LRMD reliably preserved. From the data fusion viewpoint, the KF, EnKF, and PEKF methods indirectly inject the high spatial information of the observations into the LSMs, and the indirect injection loses many details of the observations that cannot be extracted by the three methods. In data fusion, many studies have showed that à trous wavelet transform (AWT) is very powerful in injecting the details of the HRPD into the LRMD with the minimum influence to the latter.The AWT was introduced by Holdschneider et al. in 2002 using an à trous (holes) algorithm and can preserve the translation invariance; that is, a translation of the original signal necessarily implies a translation of the corresponding wavelet coefficients [16]. The AWT is a computationally easy, dyadic, redundant, undecimated, nonorthogonal, and symmetric decomposition and provides good localization in both frequency and space domains by decomposing the data into multiple channels with the same size and decreasing resolutions [17]. These advantages make the AWT suitable for fusing multisource data. Therefore, it yields a better integration of the spatial and spectral quality than other methods. Recently, many fusion procedures based on the AWT have been proposed [18, 19].The aim of using the AWT to ET assimilation is to improve ET prediction and hence to indirectly correct the heat flux predictions, meanwhile achieving physically correct soil moisture and temperature estimates through DA. Hence, this study tests the application of the AWT to assimilate remote sensed EP for producing better heat fluxes. Seldom papers could be found about the use of the AWT in DA [20–22] which shows promising results. However, more research is required to unearth the potential of the AWT to improve heat flux predications along with soil moisture and temperature predictions from the LSMs.In Section2, these advantages of the AWT are introduced into the DA. A novel assimilation method is proposed for improving the EP of the LSMs by means of the MODIS data-retrieved EP (MDRE), based on the joint use of the AWT and the PEKF. We use the AWT to decompose the MDRE and the EP of the LSMs into multiple channels and then perform the traditional PEKF in the AWT domain. The proposed AWT-PEKF method combines the advantages of the AWT and the PEKF and can fully inject the details of the MDRE into the EP of the LSMs.Experimental results conducted in Section3 illustrate that there is consistency between the improvement analysis and the quality report of the assimilated EP. The performance of the AWT-PEKF method is quantified in improving the LSMs to predict the time varying EP. Intercomparisons are also made to show the advantage of the AWT-PEKF method over the conventional approaches based on the EnKF and PEKF. Estimation results indicate that the AWT-PEKF method can effectively retrieve the EP with a satisfactory accuracy. Section 4 concludes the paper.


## 2. Combined AWT-PEKF Data Assimilation Method

The PEKF can be understood as a purely statistical Monte Carlo method, where the ensemble of the SVLs evolves in state space with considering the mean of the ensemble as the best estimate and the spreading of the ensemble as the error variance [9]. The PEKF has gained popularity because of its simple conceptual formulation and relative ease of implementation; for example, it requires no derivation of a tangent linear or adjoint models and no integrations backward in time. Further, the computational requirements are affordable and comparable with other popular assimilation methods [9]. Hence, in this paper, the PEKF is adopted to assimilate the EP. The following briefly describes the PEKF.Given theI SVLs of the LSMs at time n-1, denoted by  Xn-1=(x1n-1,…,xin-1,…,xIn-1)T where the superscript T denotes transpose, randomly perturb Xn-1K times using Gaussian noise to produce the ensemble X(1,…,K)n-1=(X1n-1,…,Xkn-1,…,XKn-1). The PEKF consists of a forecast step:
(1)Xbkn=M(Xkn-1),1≤k≤K.In (1), M denotes the nonlinear LSM, the subscript b denotes the background, and n denotes the simulation time.The forecast error can be analyzed as(2)BnHT=〈(Xbn-Xbn¯)(H(Xbn)-H(Xbn)¯)T〉,(3)HBnHT=〈(H(Xbn)-H(Xbn)¯)(H(Xbn)-H(Xbn)¯)T〉,(4)Xbn¯=1K∑k=1KXbkn,(5)Xb(1,…,K)n=(M(X1n-1),…,M(Xkn-1),…,M(XKn-1)),(6)H(Xbn)¯=1K∑k=1KH(Xbkn),(7)H(Xb(1,…,K)n)=(H(Xb1n-1),…,H(Xbkn-1),…,H(XbKn-1)),
where H denotes the observation operator. Then the gain matrix G is obtained as
(8)G=BnHT(HBnHT+On)-1.
In (8), O denotes the observation error matrix.Finally, the assimilation is performed as(9)Xan=Xbn+G(Yon-H(Xbn)).In (9), Yon denotes the observation at time n.The PEKF method uses the ensemble to describe the covariance matrix of the SVLs and avoids the explicit calculation of covariance matrix needed in the EKF. A modest misspecification of the initial ensemble normally does not influence the results very much over time. The PEKF allows for a wide range of noise models and one is not restricted to using Gaussian distributed noise.It can be found from (9) that the assimilation is completed based on single grid and does not consider the interrelation between grids. The interrelation between grids is mainly reflected by the detail and textural information. The reason that the spatial resolution of the observation is higher than the simulation of the LSM is that the observation has more detail and textural information than the simulation of the LSM. The assimilation scheme based on single grid cannot fully inject the detail and textural information of the observation into the simulation of the LSM. The main purpose of the DA is to integrate the high accuracy of the observation and the temporal evolution of the LSMs together.From the perspective of data fusion, the assimilation process can be considered as constructing one coefficient with both the same temporal response as the simulation of the LSM and the same spatial response as the observation at a particular grid location. With the development of the AWT, we expect much room for improvement over the traditional single-grid-based assimilation scheme to merge the observation and the simulation of the LSM in the AWT domain since the wavelet filter can consider the interrelation between the grids during the detail extraction of the observation.The AWT can provide good localization in both frequency and space domains in terms of decomposing the data with finite energy into multiple channels, each one of them with a different degree of resolution [23]. Because of its shift invariance property, the AWT has been successfully used for fusing the data with different resolutions [24]. The AWT representation of the data can be obtained by using à trous algorithm [17].This algorithm consists basically in the application of consecutive convolutions between the data under analysis and a scaling function at distinct degradation levels [25]. One of the most widely used scaling functions for the execution of the à trous algorithm is the b3-spline [23] as follows:
(10)1256(1464141624164624362464162416414641).If the original data is represented byD, the wavelet coefficients, Cj, for the level j are obtained by the difference between two consecutive degraded data, Dj+1and Dj, as shown in the following equation:
(11)Cj=Dj-Dj+1.To carry out the data synthesis from a degradation levelJ, an additive criterion should be applied in which all the coefficients obtained are added to the last degradation level of the original data, as shown in the following equation:
(12)D=DJ+∑j=1JCj.By manipulating independently the approximation components and the wavelet planes of the data, the AWT shows excellent performance in integrating the spectral information of the multispectral data and the spatial information of the panchromatic data. For the DA, the purpose is to integrate the high spatial resolution information of the observation into the simulation of the LSMs. Hence, the purposes of the data fusion and the DA have similarity.It can be seen from (9) that the update is completed by adding the product of the gain and the difference to the SVLs. The high spatial information of the observation used to improve the SVLs is contained in the difference and is modulated into the SVLs by the gain. From the perspective of the AWT-based fusion, the additive DA cannot consider the difference between the wavelet planes. For different wavelet planes, the gain may be different and is varied with the wavelet plane. It will be beneficial to complete the DA in the AWT domain. What is more, the additive DA cannot consider the interrelation between the grids, which leads to the under injection of the detail information of the observation. With the AWT, much room for improvement over the traditional DA can be expected.After the forecast of the LSM using the ensemble, see (1). Apply the AWT to the forecast ensemble and the observation as
(13)Xbkn=∑j=0JXbkjn,(14)Yon=∑j=0JYojn.In (13), Xbk0n denotes the approximation component of   Xbkn, and {Xbkjn}j=1J denotes the J wavelet planes of Xbkn. In (14), Yo0n denotes the approximation component of Yon, and {Yojn}j=1J denotes the J wavelet planes of  Yon.Apply the AWT toH as
(15)H(Xbkn)=∑j=0JH(Xbkn)j.In (15), H(Xbkn)0 denotes the approximation component of H(Xbkn), and {H(Xbkn)j}j=1J denotes the J wavelet planes of H(Xbkn).Then, (2)–(4) and (6) can be, respectively, rewritten as:
(16)BjnHjT=〈(Xbjn-Xbjn¯)(H(Xbn)j-H(Xbn)j¯)T〉,(17)HjBjnHjT=〈(H(Xbn)j-H(Xbn)j¯)(H(Xbn)j-H(Xbn)j¯)T〉,(18)Xbjn¯=1K∑k=1KXbkjn,(19)H(Xbn)j¯=1K∑k=1KH(Xbkn)j.Then the gain matrixGj at decomposition level j is obtained as
(20)Gj=BjnHjT(HjBjnHjT+Ojn)-1.In (20), Ojn denotes the observation error covariance matrix at decomposition level j at time n. Ojn can be obtained using the perturbed observations [10].Finally, the assimilation is performed as(21)Xan=Xbn+∑j=0JGj(Yojn-H(Xbn)j).Compared with (9), equation (21) has two obvious advantages. First, in (9), the execution of the assimilation is completed grid by grid, and the relationship between grids is not considered; in (21), the interrelation between the local grids is considered and is quantified by filtering the grids equivalent to the length of the filter during the decomposition process. Second, in (9), the same gain is used to extract the new information absent in the SVLs from the observation; in (21), the gain is calculated in the AWT domain and is obtained according to the wavelet planes. It is more effective to update the SVLs using wavelet plane varying gains according to the observation.For assimilating the MDRE into the EP of the LSM, the DA procedure based on the AWT-PEKF method comprises the following steps (Figure1):(1)
run the LSM in open time loop;(2)
given the Gaussian probability density function, produce the ensemble of the SVLs according to their probability characteristic. Run the LSM during ensemble loop;(3)
integrate the LSM to get the simulation result over the simulation range, namely, complete range loop;(4)
end ensemble loop. Obtain the mean EP by averaging the simulated EPs as(22)EP0n=∑k=1KEPkn.
In (22), K denotes the size of the ensemble of the SVLs. EP0n denotes the mean EP at the current simulation time n;(5)
perturb the MDREK times and produce the ensemble of the MDRE, denoted by {MDREkn}k=0K [9]. MDRE0n denotes the original MDRE;(6)
apply the AWT on the MDRE and EP. Respectively, obtain the decomposition results as(23)EPkn=∑j=0JEPkjn,0≤k≤K,(24)MDREkn=∑j=0JMDREkjn,0≤k≤K.
In (23) and (24), J is the number of the AWT decomposition levels;(7)
calculate the gain at decomposition levelj, Gj, using the decomposition results as (25):
(25)Gj=BjnHjT(HjBjnHjT+Ojn)-1,(26)BjnHjT=〈(EPjn-EP0jn)(EPjn-EP0jn)T〉,(27)HjBjnHjT=〈(EPjn-EP0jn)(EPjn-EP0jn)T〉.
In (25)–(27), the observation operator, H, is 1 because it does not need transformation between the assimilated variable and the observation;(8)
innovate the EP of the LSM as(28)EPn=EP0n+∑j=0JGj(MDREjn-EPjn);(9)
integrate the LSM until ending open time loop.Figure 1
The schematic flowchart of the AWT-PEKF method.


## 3. Experimental Section

In this section, the AWT-PEKF method is tested through assimilating the EP of the Noah model with the MDRE retrieved using one operational two-layer method (OTLM) [26], and the assimilated result is validated with the site observations from two field stations in Chinese Ecosystem Research Network (Cern), each of which has a different underlying land cover condition. The Noah model is developed by the environmental modeling center of national center for environmental prediction [27]. In addition, we also compare the assimilated results with those of the EnKF and PEKF methods.


### 3.1. Site Observations

Cern was established in 1988 for studying the problems related with Chinese ecology. The network provides continuous observations of ecosystem level exchanges of biogeochemical, water, energy, and momentum at diurnal, monthly, seasonal, and annual time scales. Cern currently comprises forty-six sites, which distribute over the whole of China. The data of the site observations used in this study are obtained from the Yucheng and Lawn sites. At the two sites, ET has been continuously measured by eddy correlation system which is composed of a 3D sonic anemometer and an open path CO2/H2O analyzer since 2001. The ET data that were obtained by averaging the original data of 10 Hz over 30 min are used in the experiment as the validating data.The Yucheng site was installed on Huabei plain of China. The field site is located near the Yucheng county, Shandong province, China, and the geographical coordinate is approximate 36.96°N, 116.63°E, and the altitude is 20 m above sea level. The climate of the area is warm and semi-humid with annual mean temperature 13.1°C, mean precipitation 600 mm per year, and the average monthly relative humidity 66.44%. The soil is mainly classified as sandy-clay loam and sandy loam. The characteristic of this site is typical of agricultural plot. The farm has been continuously alternated between wheat and corn. Lawn site was temporally installed about 54 kilometers from the Yucheng site in 2005. The latitude and longitude are 36.46°N and 116.13°E. The underlying vegetation of this site is grass. The leaf area index is about 2.0. The average canopy height is about 20–40 cm. The area is 10 km2. Figure 2 shows the schematic positions of the two sites, and the research range is mentioned in the following section.Figure 2
The schematic positions of Yucheng site (⚫), Lawn site (■), and research range (□).The two sites represent the main land covers of the studied area. The underlying surfaces of the two sites are relatively homogeneous. Therefore, the observation data of the ET from the two sites can represent the ET statuses of the two homogeneous areas around the two sites. The two areas correspond at least to nine grids of the MDRE and the EP of the Noah model at the spatial resolution of 1 km. Hence, it is sufficient to use these observations as the validating data.


### 3.2. The MDRE

In order to retrieve the MDRE, the OTLM proposed by Zhang et al. [26] is used. In the OTLM, pixel component arranging and comparing algorithm (PCACA) and layered energy-separating algorithm (LESA) are two key procedures. The two procedures are used to partition mixed surface temperature and mixed surface albedo according to vegetation fraction to obtain those of bare soil and vegetation. Then, net radiation of bare soil and vegetation is calculated using the partitioned results. Finally, the MDRE is retrieved with a spatial resolution of 1 km.The OTLM has several advantages. First, the PCACA is based on vegetation fraction and ground temperature trapezoid relation theory and initiates a new method of decomposing mixed pixels. Second, it is very convenient because only single angle remote sensing data are required which can be obtained from most of the satellite sensors. Third, the core of the LESA is Bowen-ratio energy balance method, which reduces the uncertainties in surface energy partition based on the Beer law. Fourth, key nonremote sensing parameters that influence regional ET can also be obtained by using the OTLM. The OTLM has been successfully applied to Huabei area, China, [28]. Hence, the OTLM is used to retrieve the MDRE. Detailed descriptions about the OTLM can be found in [26, 28].


### 3.3. The Temporal Extension of the MDRE

Because the MODIS is easily influenced by the clouds, the data are not good enough to retrieve the MDRE every day. Figure3 shows the available MDREs retrieved from the Aqua data at 14:30 local solar time in May, 2005. The unit for Figure 3, the following mentioned figures, and the site observation is W/m2.The available MDREs retrieved from the Aqua data at 14:30 local solar time in May, 2005. (a) The MDRE on May 2, 2005; (b) the MDRE on May 7, 2005; (c) the MDRE on May 14, 2005; (d) the MDRE on May 19, 2005; (e) the MDRE on May 25, 2005; (f) the MDRE on May 31, 2005.
(a)
(b)
(c)
(d)
(e)
(f)On the contrary, the better the assimilated results are, the more there are the MDREs. Hence, it first requires to temporally extend the MDRE by the available MDREs and the EPs of the Noah model. The series of the extended MDRE is then used as the synthetic observations. In order to extend the MDRE, an intensity modulation method is employed because this method is very simple and can be easily implemented into the Noah model to perform high speed real-time MDRE extension for a long temporal span. This method can interpolate the absent MDRE with the fidelity to the evolvement trend of the Noah EP. This unavailableMDREn at the simulation time n is obtained as
(29)MDREn=∑s=1S(|s-n|∑t=1S|t-n|EPnEPsMDREs).In (29), s, n, and t denote the time. {MDREs}1≤s≤S is the set of the available good MDREs, and S is the size of the set. EPs is the Noah EP corresponding to the MDREs, and EPn is the Noah EP at the simulation time n. s-n describes the departure of the EPs from the EPn and characterizes the temporal difference of the EPs from the EPn, while t-n describes the departure of the MDREt from the MDREn, and quantifies the contribution of the MDREt to the MDREn. By using a ratio between EPn and EPs, spatial details of MDREs are modulated to the predicted MDREn without altering its spatial properties and contrast.


### 3.4. The Noah Model

The Noah model is originated from a physically based land surface-vegetation-atmosphere-transfer scheme. During the past ten years, it underwent substantial upgrades, including modifications to the formulations of canopy conductance, bare soil evaporation, vegetation phenology, ground heat flux, and so forth. These model enhancements significantly improve its performance, and are physically more faithful to nature and thus most likely the route for more improvements in the future [29, 30].The Noah model was chosen for assimilation test for several reasons. It can simulate many states of the land surface, including soil temperature, skin temperature, and the energy and water fluxes of the land surface energy and water balances. In various coupled and uncoupled assessments, it has been proven to have the ability to reproduce the observed land surface energy, and water budgets effectively [29, 30]. Because of its comprehensive nature, the model yields more output variables. The model is updated periodically on the NCEP website (ftp://ftp.emc.ncep.noaa.gov/mmb/gcp/ldas/noahlsm/, accessed on April 8, 2013) and version 2.7.1 is used in this experiment (release date: February 9, 2005).The Noah model contains four vertical soil layers: a thin 10 cm top layer, a second root zone layer of 30 cm, a deep root zone of 60 cm, and a subroot zone of 100 cm. It can be run for 13 vegetation covers (2 of which use the same parameter values) and nine different soil types (two of which also use the same parameters). It has 33 parameters: 10 related to vegetation and 23 that describe soil properties. It also has 16 initial states (when run with four root layers). The model uses a local greenness fraction from the normalized difference vegetation index (NDVI) to establish seasonality in the model for each of the 13 vegetation types.In this assimilation experiment, the Noah model was configured to have dimensions of 250 × 250 at 1 km × 1 km spanning a domain bounded by 35.77°N to 38.26°N, 114.81°E to 117.3°E. Time step is 10800 seconds. Assimilation time span is between 1 January, 2005 and 31 May, 2005 after the Noah model is run from 1 June, 2004 because the MODIS data in this span are less affected by clouds, and it is beneficial for retrieving the MDRE. On this grid, the elevation was derived from the 1 km digital elevation of the GTOPO30 database [31]. The vegetation classification was derived from the global, 1 km, Advanced Very High Resolution Radiometer (AVHRR) based, 13 class vegetation database [32]. The 9-class soil texture data were derived from the top layer of the 1 km, 11-layer soil dataset, state soil geographic database [33]. The monthly vegetation cover fraction data are taken from the satellite-based AVHRR 5-year global monthly climatology of green vegetation fraction [34]. Six EPs of the Noah model at 15:00 solar time corresponding to the nearest time of the six MDREs are shown in Figure 4.The EPs of the Noah model at 15:00 solar time closest to the available MDREs in Figure3. (a) The EP on May 2, 2005; (b) the EP on May 7, 2005; (c) the EP on May 14, 2005; (d) the EP on May 19, 2005; (e) the EP on May 25, 2005; (f) the EP on May 31, 2005.
(a)
(b)
(c)
(d)
(e)
(f)


### 3.5. Results

The data from the two sites are first used to evaluate the ET values of the available MDREs and the Noah model at the corresponding locations before DA. For validation, we utilize root-mean-square error (RMSE) as the estimation index. RMSE is defined as(30)RMSE=1M∑m=1M(o-d)2.In (30), o and d, respectively, denote the observation and validated data. M denotes the number of the data. The unit for the RMSE is W/m2 in this paper. The estimation results are, respectively, shown in Table 1.Table 1
The RMSE results of the MDRE and the Noah EP using the two site observations.
Yucheng site
Lawn site
MDRE
EP
MDRE
EP
RMSE
18.74
51.58
23.47
68.69It can be found from Table1 that the Noah EP has large RMSE values than the MDREs compared with the two site observations. The main reason is that the spatial resolution of the MDRE is higher than that of the Noah EP, and more details can be found from the MDRE by comparing Figures 3 and 4. The configuration of the Noah model to the spatial resolution of the MDRE did not produce new information in the EP because of the coarse resolution of the driving data. It can be seen from Table 1 that the EP can be improved significantly with the MDREs through DA.When performing the AWT-PEKF method, the SVLs including skin temperature, soil temperature and volumetric liquid soil moisture for the first soil layer, canopy water content, and the MDRE are perturbed fifty times according to their error Gaussian distributions, respectively. The four SVLs are perturbed because they are the key variables in the calculation of the ET in the Noah model. For comparisons, the EnKF, and PEKF methods are also performed. The assimilated EPs from the EnKF, PEKF, and proposed methods are, respectively, shown in Figures5, 6, and 7.The assimilated EPs from the EnKF method at 15:00 solar time. (a) The assimilated EP on May 2, 2005; (b) the assimilated EP on May 7, 2005; (c) the assimilated EP on May 14, 2005; (d) the assimilated EP on May 19, 2005; (e) the assimilated EP on May 25, 2005; (f) the assimilated EP on May 31, 2005.
(a)
(b)
(c)
(d)
(e)
(f)The assimilated EPs from the PEKF method at 15:00 solar time. (a) The assimilated EP on May 2, 2005; (b) the assimilated EP on May 7, 2005; (c) the assimilated EP on May 14, 2005; (d) the assimilated EP on May 19, 2005; (e) the assimilated EP on May 25, 2005; (f) the assimilated EP on May 31, 2005.
(a)
(b)
(c)
(d)
(e)
(f)The assimilated EPs from the AWT-PEKF method at 15:00 solar time. (a) The assimilated EP on May 2, 2005; (b) the assimilated EP on May 7, 2005; (c) the assimilated EP on May 14, 2005; (d) the assimilated EP on May 19, 2005; (e) the assimilated EP on May 25, 2005; (f) the assimilated EP on May 31, 2005.
(a)
(b)
(c)
(d)
(e)
(f)In order to evaluate the three methods, the two site observations are employed to validate the DA results by the RMSE. Table2 and Figure 8 show the results.Table 2
The RMSE results of the three methods using the two site observations.
EnKF
PEKF
AWT-PEKF
Yucheng site
28.75
22.54
16.81
Lawn site
25.18
21.97
18.16Figure 8
The RMSE histograms of the three methods using the two site observations.The RMSE reveals the accurate degree of the EPs produced by each method. The lower the RMSE is, the better the assimilation effect is, and vice versa. It can be seen from Table2 and Figure 8 that all methods improve significantly the accuracy of the EPs compared with the EPs without assimilation: the RMSE values are reduced from 51.58 and 68.69 in Table 1 to 28.75, 22.54, and 16.81 and 25.18, 21.97, and 18.16, respectively, for the Yucheng and Lawn sites. The three methods allow an obvious accuracy improvement of the EPs when increasing the spatial details of the EPs. Table 2 and Figure 8 show that the RMSE values of the AWT-PEKF method are the lowest among the three methods. Therefore, we can draw a conclusion that the EPs of the AWT-PEKF method have the least error. The EPs of the PEKF method have slight error. The EPs of the EnKF method have the largest error. This is expected because the ET equation used in the Noah model is a highly nonlinear equation and the perturbation only for the SVLs used in the EnKF can only propagate analytically with finite precision, which introduces large errors in the EP estimation and leads to suboptimal performance.As seen from Tables1 and 2, the AWT-PEKF method produces the best assimilation effect. This is probably because the AWT-based assimilation scheme does good to improve the EP simulation of the Noah model. The AWT-PEKF method gets the advantage of the EnKF and PEKF methods, because the novel assimilation scheme injects the spatial details of the MDRE into the Noah EP by taking into account both the interrelationship between grids, as is the case of the AWT decomposition process, and ensemble-based innovation, as is the case of the PEKF forecast. The two procedures allow the AWT-PEKF method to produce the EPs closer to the ET observation. One can therefore conclude that the Noah model can predict the EP relatively accurate using the AWT-PEKF method.By combining the quantitative estimation results and the intercomparison, it can be seen that the AWT-PEKF method gives the EPs closer to the measured EPs than the EnKF and PEKF methods when the assimilated EPs are compared with the observations from the Yucheng and Lawn sites.


### 3.6. Discussion

Though the AWT-PEKF method outperforms the EnKF and PEKF methods in the experiment, three points are needed to be studied further. The first is that the validating data are sparse. In validation, there are only two available field sites. Because it is difficult to get the ET observation, it needs to confirm whether or not the assimilated EPs obtained by the AWT-PEKF method are also close to the ET observation at other locations. The second is that it needs to test whether or not the AWT-PEKF method is also effective in assimilating other variables. Other variables, such as skin temperature and soil moisture, are forecasted differently from the ET and are also influenced by many factors. Hence, the extension of the AWT-PEKF method to other variables is a big job. The third is that the assimilation of the ET cannot influence the consecutive simulation of the Noah model.As for the first point, the potential solution is to perform the AWT-PEKF method in the area where the field sites are available, meanwhile the MDREs can be retrieved from cloud-free MODIS data. As for the second point, the AWT-PEKF method can be easily extended to assimilate skin temperature and soil moisture products only if the two products retrieved from the MODIS data are prepared in advance. As for the third point, in the following study, we will introduce another novel assimilation method, in which the SVLs relative with the ET are simultaneously updated in order to transfer the assimilation effect into the consecutive simulation of the LSM. The assimilation idea will be presented in another research paper.


## 3.1. Site Observations

Cern was established in 1988 for studying the problems related with Chinese ecology. The network provides continuous observations of ecosystem level exchanges of biogeochemical, water, energy, and momentum at diurnal, monthly, seasonal, and annual time scales. Cern currently comprises forty-six sites, which distribute over the whole of China. The data of the site observations used in this study are obtained from the Yucheng and Lawn sites. At the two sites, ET has been continuously measured by eddy correlation system which is composed of a 3D sonic anemometer and an open path CO2/H2O analyzer since 2001. The ET data that were obtained by averaging the original data of 10 Hz over 30 min are used in the experiment as the validating data.The Yucheng site was installed on Huabei plain of China. The field site is located near the Yucheng county, Shandong province, China, and the geographical coordinate is approximate 36.96°N, 116.63°E, and the altitude is 20 m above sea level. The climate of the area is warm and semi-humid with annual mean temperature 13.1°C, mean precipitation 600 mm per year, and the average monthly relative humidity 66.44%. The soil is mainly classified as sandy-clay loam and sandy loam. The characteristic of this site is typical of agricultural plot. The farm has been continuously alternated between wheat and corn. Lawn site was temporally installed about 54 kilometers from the Yucheng site in 2005. The latitude and longitude are 36.46°N and 116.13°E. The underlying vegetation of this site is grass. The leaf area index is about 2.0. The average canopy height is about 20–40 cm. The area is 10 km2. Figure 2 shows the schematic positions of the two sites, and the research range is mentioned in the following section.Figure 2
The schematic positions of Yucheng site (⚫), Lawn site (■), and research range (□).The two sites represent the main land covers of the studied area. The underlying surfaces of the two sites are relatively homogeneous. Therefore, the observation data of the ET from the two sites can represent the ET statuses of the two homogeneous areas around the two sites. The two areas correspond at least to nine grids of the MDRE and the EP of the Noah model at the spatial resolution of 1 km. Hence, it is sufficient to use these observations as the validating data.


## 3.2. The MDRE

In order to retrieve the MDRE, the OTLM proposed by Zhang et al. [26] is used. In the OTLM, pixel component arranging and comparing algorithm (PCACA) and layered energy-separating algorithm (LESA) are two key procedures. The two procedures are used to partition mixed surface temperature and mixed surface albedo according to vegetation fraction to obtain those of bare soil and vegetation. Then, net radiation of bare soil and vegetation is calculated using the partitioned results. Finally, the MDRE is retrieved with a spatial resolution of 1 km.The OTLM has several advantages. First, the PCACA is based on vegetation fraction and ground temperature trapezoid relation theory and initiates a new method of decomposing mixed pixels. Second, it is very convenient because only single angle remote sensing data are required which can be obtained from most of the satellite sensors. Third, the core of the LESA is Bowen-ratio energy balance method, which reduces the uncertainties in surface energy partition based on the Beer law. Fourth, key nonremote sensing parameters that influence regional ET can also be obtained by using the OTLM. The OTLM has been successfully applied to Huabei area, China, [28]. Hence, the OTLM is used to retrieve the MDRE. Detailed descriptions about the OTLM can be found in [26, 28].


## 3.3. The Temporal Extension of the MDRE

Because the MODIS is easily influenced by the clouds, the data are not good enough to retrieve the MDRE every day. Figure3 shows the available MDREs retrieved from the Aqua data at 14:30 local solar time in May, 2005. The unit for Figure 3, the following mentioned figures, and the site observation is W/m2.The available MDREs retrieved from the Aqua data at 14:30 local solar time in May, 2005. (a) The MDRE on May 2, 2005; (b) the MDRE on May 7, 2005; (c) the MDRE on May 14, 2005; (d) the MDRE on May 19, 2005; (e) the MDRE on May 25, 2005; (f) the MDRE on May 31, 2005.
(a)
(b)
(c)
(d)
(e)
(f)On the contrary, the better the assimilated results are, the more there are the MDREs. Hence, it first requires to temporally extend the MDRE by the available MDREs and the EPs of the Noah model. The series of the extended MDRE is then used as the synthetic observations. In order to extend the MDRE, an intensity modulation method is employed because this method is very simple and can be easily implemented into the Noah model to perform high speed real-time MDRE extension for a long temporal span. This method can interpolate the absent MDRE with the fidelity to the evolvement trend of the Noah EP. This unavailableMDREn at the simulation time n is obtained as
(29)MDREn=∑s=1S(|s-n|∑t=1S|t-n|EPnEPsMDREs).In (29), s, n, and t denote the time. {MDREs}1≤s≤S is the set of the available good MDREs, and S is the size of the set. EPs is the Noah EP corresponding to the MDREs, and EPn is the Noah EP at the simulation time n. s-n describes the departure of the EPs from the EPn and characterizes the temporal difference of the EPs from the EPn, while t-n describes the departure of the MDREt from the MDREn, and quantifies the contribution of the MDREt to the MDREn. By using a ratio between EPn and EPs, spatial details of MDREs are modulated to the predicted MDREn without altering its spatial properties and contrast.


## 3.4. The Noah Model

The Noah model is originated from a physically based land surface-vegetation-atmosphere-transfer scheme. During the past ten years, it underwent substantial upgrades, including modifications to the formulations of canopy conductance, bare soil evaporation, vegetation phenology, ground heat flux, and so forth. These model enhancements significantly improve its performance, and are physically more faithful to nature and thus most likely the route for more improvements in the future [29, 30].The Noah model was chosen for assimilation test for several reasons. It can simulate many states of the land surface, including soil temperature, skin temperature, and the energy and water fluxes of the land surface energy and water balances. In various coupled and uncoupled assessments, it has been proven to have the ability to reproduce the observed land surface energy, and water budgets effectively [29, 30]. Because of its comprehensive nature, the model yields more output variables. The model is updated periodically on the NCEP website (ftp://ftp.emc.ncep.noaa.gov/mmb/gcp/ldas/noahlsm/, accessed on April 8, 2013) and version 2.7.1 is used in this experiment (release date: February 9, 2005).The Noah model contains four vertical soil layers: a thin 10 cm top layer, a second root zone layer of 30 cm, a deep root zone of 60 cm, and a subroot zone of 100 cm. It can be run for 13 vegetation covers (2 of which use the same parameter values) and nine different soil types (two of which also use the same parameters). It has 33 parameters: 10 related to vegetation and 23 that describe soil properties. It also has 16 initial states (when run with four root layers). The model uses a local greenness fraction from the normalized difference vegetation index (NDVI) to establish seasonality in the model for each of the 13 vegetation types.In this assimilation experiment, the Noah model was configured to have dimensions of 250 × 250 at 1 km × 1 km spanning a domain bounded by 35.77°N to 38.26°N, 114.81°E to 117.3°E. Time step is 10800 seconds. Assimilation time span is between 1 January, 2005 and 31 May, 2005 after the Noah model is run from 1 June, 2004 because the MODIS data in this span are less affected by clouds, and it is beneficial for retrieving the MDRE. On this grid, the elevation was derived from the 1 km digital elevation of the GTOPO30 database [31]. The vegetation classification was derived from the global, 1 km, Advanced Very High Resolution Radiometer (AVHRR) based, 13 class vegetation database [32]. The 9-class soil texture data were derived from the top layer of the 1 km, 11-layer soil dataset, state soil geographic database [33]. The monthly vegetation cover fraction data are taken from the satellite-based AVHRR 5-year global monthly climatology of green vegetation fraction [34]. Six EPs of the Noah model at 15:00 solar time corresponding to the nearest time of the six MDREs are shown in Figure 4.The EPs of the Noah model at 15:00 solar time closest to the available MDREs in Figure3. (a) The EP on May 2, 2005; (b) the EP on May 7, 2005; (c) the EP on May 14, 2005; (d) the EP on May 19, 2005; (e) the EP on May 25, 2005; (f) the EP on May 31, 2005.
(a)
(b)
(c)
(d)
(e)
(f)


## 3.5. Results

The data from the two sites are first used to evaluate the ET values of the available MDREs and the Noah model at the corresponding locations before DA. For validation, we utilize root-mean-square error (RMSE) as the estimation index. RMSE is defined as(30)RMSE=1M∑m=1M(o-d)2.In (30), o and d, respectively, denote the observation and validated data. M denotes the number of the data. The unit for the RMSE is W/m2 in this paper. The estimation results are, respectively, shown in Table 1.Table 1
The RMSE results of the MDRE and the Noah EP using the two site observations.
Yucheng site
Lawn site
MDRE
EP
MDRE
EP
RMSE
18.74
51.58
23.47
68.69It can be found from Table1 that the Noah EP has large RMSE values than the MDREs compared with the two site observations. The main reason is that the spatial resolution of the MDRE is higher than that of the Noah EP, and more details can be found from the MDRE by comparing Figures 3 and 4. The configuration of the Noah model to the spatial resolution of the MDRE did not produce new information in the EP because of the coarse resolution of the driving data. It can be seen from Table 1 that the EP can be improved significantly with the MDREs through DA.When performing the AWT-PEKF method, the SVLs including skin temperature, soil temperature and volumetric liquid soil moisture for the first soil layer, canopy water content, and the MDRE are perturbed fifty times according to their error Gaussian distributions, respectively. The four SVLs are perturbed because they are the key variables in the calculation of the ET in the Noah model. For comparisons, the EnKF, and PEKF methods are also performed. The assimilated EPs from the EnKF, PEKF, and proposed methods are, respectively, shown in Figures5, 6, and 7.The assimilated EPs from the EnKF method at 15:00 solar time. (a) The assimilated EP on May 2, 2005; (b) the assimilated EP on May 7, 2005; (c) the assimilated EP on May 14, 2005; (d) the assimilated EP on May 19, 2005; (e) the assimilated EP on May 25, 2005; (f) the assimilated EP on May 31, 2005.
(a)
(b)
(c)
(d)
(e)
(f)The assimilated EPs from the PEKF method at 15:00 solar time. (a) The assimilated EP on May 2, 2005; (b) the assimilated EP on May 7, 2005; (c) the assimilated EP on May 14, 2005; (d) the assimilated EP on May 19, 2005; (e) the assimilated EP on May 25, 2005; (f) the assimilated EP on May 31, 2005.
(a)
(b)
(c)
(d)
(e)
(f)The assimilated EPs from the AWT-PEKF method at 15:00 solar time. (a) The assimilated EP on May 2, 2005; (b) the assimilated EP on May 7, 2005; (c) the assimilated EP on May 14, 2005; (d) the assimilated EP on May 19, 2005; (e) the assimilated EP on May 25, 2005; (f) the assimilated EP on May 31, 2005.
(a)
(b)
(c)
(d)
(e)
(f)In order to evaluate the three methods, the two site observations are employed to validate the DA results by the RMSE. Table2 and Figure 8 show the results.Table 2
The RMSE results of the three methods using the two site observations.
EnKF
PEKF
AWT-PEKF
Yucheng site
28.75
22.54
16.81
Lawn site
25.18
21.97
18.16Figure 8
The RMSE histograms of the three methods using the two site observations.The RMSE reveals the accurate degree of the EPs produced by each method. The lower the RMSE is, the better the assimilation effect is, and vice versa. It can be seen from Table2 and Figure 8 that all methods improve significantly the accuracy of the EPs compared with the EPs without assimilation: the RMSE values are reduced from 51.58 and 68.69 in Table 1 to 28.75, 22.54, and 16.81 and 25.18, 21.97, and 18.16, respectively, for the Yucheng and Lawn sites. The three methods allow an obvious accuracy improvement of the EPs when increasing the spatial details of the EPs. Table 2 and Figure 8 show that the RMSE values of the AWT-PEKF method are the lowest among the three methods. Therefore, we can draw a conclusion that the EPs of the AWT-PEKF method have the least error. The EPs of the PEKF method have slight error. The EPs of the EnKF method have the largest error. This is expected because the ET equation used in the Noah model is a highly nonlinear equation and the perturbation only for the SVLs used in the EnKF can only propagate analytically with finite precision, which introduces large errors in the EP estimation and leads to suboptimal performance.As seen from Tables1 and 2, the AWT-PEKF method produces the best assimilation effect. This is probably because the AWT-based assimilation scheme does good to improve the EP simulation of the Noah model. The AWT-PEKF method gets the advantage of the EnKF and PEKF methods, because the novel assimilation scheme injects the spatial details of the MDRE into the Noah EP by taking into account both the interrelationship between grids, as is the case of the AWT decomposition process, and ensemble-based innovation, as is the case of the PEKF forecast. The two procedures allow the AWT-PEKF method to produce the EPs closer to the ET observation. One can therefore conclude that the Noah model can predict the EP relatively accurate using the AWT-PEKF method.By combining the quantitative estimation results and the intercomparison, it can be seen that the AWT-PEKF method gives the EPs closer to the measured EPs than the EnKF and PEKF methods when the assimilated EPs are compared with the observations from the Yucheng and Lawn sites.


## 3.6. Discussion

Though the AWT-PEKF method outperforms the EnKF and PEKF methods in the experiment, three points are needed to be studied further. The first is that the validating data are sparse. In validation, there are only two available field sites. Because it is difficult to get the ET observation, it needs to confirm whether or not the assimilated EPs obtained by the AWT-PEKF method are also close to the ET observation at other locations. The second is that it needs to test whether or not the AWT-PEKF method is also effective in assimilating other variables. Other variables, such as skin temperature and soil moisture, are forecasted differently from the ET and are also influenced by many factors. Hence, the extension of the AWT-PEKF method to other variables is a big job. The third is that the assimilation of the ET cannot influence the consecutive simulation of the Noah model.As for the first point, the potential solution is to perform the AWT-PEKF method in the area where the field sites are available, meanwhile the MDREs can be retrieved from cloud-free MODIS data. As for the second point, the AWT-PEKF method can be easily extended to assimilate skin temperature and soil moisture products only if the two products retrieved from the MODIS data are prepared in advance. As for the third point, in the following study, we will introduce another novel assimilation method, in which the SVLs relative with the ET are simultaneously updated in order to transfer the assimilation effect into the consecutive simulation of the LSM. The assimilation idea will be presented in another research paper.


## 4. Conclusions

In this paper, we study the hybrid use of the AWT and PEKF methods for assimilating the MDRE into the EP of the Noah model in order to improve the consecutive simulation of the Noah model. The AWT is used to decompose the MDRE for injecting its detail information represented by wavelet planes, while the PEKF is used to complete the assimilation by the model and observation uncertainties. The AWT-PEKF method retains the respective advantages of the AWT and PEKF. Firstly, it is based on multigrid, and the interrelation between grids is considered using the wavelet filter during the filtering procedure. Secondly, according to the gain derived using the PEKF from the model and observation uncertainties, the details and textures of the MDRE are modulated into the EPs of the Noah model using the AWT from image fusion viewpoint.The performance of the proposed method is compared with those of the EnKF and PEKF methods using one assimilation experiment. Intercomparison results of the RMSE confirm the effectiveness of the AWT-PEKF method in improving the spatial accuracy of the EPs. Overall evaluation shows that the AWT-PEKF method is promising and superior to the traditional EnKF and PEKF methods. Several issues are unresolved, such as the validation of the assimilated results, the effects of ensemble size, initial perturbation fields on assimilated results, and the actual performance of this new method in real other variable assimilations. These aspects require further investigation.



---

*Source: 531810-2013-07-22.xml*