## Abstract

This study evaluates the performance of a new stationarity-based method for parameter estimation of a simple coupled water and energy balance model using in situ and remotely sensed surface soil moisture [from Advanced Microwave Scanning Radiometer for Earth Observing System (EOS; AMSR-E)] and surface temperature [from a combined Moderate Resolution Imaging Spectroradiometer (MODIS) and AMSR-E product]. Parameter estimation is carried out using both direct calibration to measured surface fluxes (latent, sensible, and ground heat) and a recently published method based on enforcing stationarity of model-predicted surface state tendency terms. The latter stationarity-based method was developed for parameter estimation without knowledge of observed fluxes—that is, it requires only forcing terms (e.g., radiation, wind speed, air temperature) and surface states (moisture and temperature). In addition, the stationarity-based method can easily handle gaps in atmospheric forcing and surface state data, as it does not integrate over time to simulate fluxes. The evaluation is conducted at three AmeriFlux sites. Changing the data sources of surface states (in situ measured and remotely sensed data) leads to little degradation in estimating turbulent fluxes despite the relatively poor quality of satellite data at some sites. In all cases, direct calibration yields marginally better predictions than the stationarity-based method, with site-averaged root-mean-square errors for daily total energy fluxes approximately 5–6 W m^{−2} lower. However, direct calibration requires observed fluxes in the objective function, which imposes a strong limit on its application.

## 1. Introduction and background

### a. Introduction

Many methods for estimating turbulent fluxes using remotely sensed data have been developed over the past few decades. Remotely sensed data may be directly introduced in semiempirical models or through energy budget methods that combine empirical and physically based relationships in order to compute the energy fluxes. For example, Bastiaanssen et al. (1998a,b) developed a Surface Energy Balance Algorithm for Land (SEBAL) to estimate the spatially distributed surface energy balance using remote sensing information (e.g., surface temperature, surface reflectivity, and vegetation index), limited ground data, and a set of empirical relationships. Su (2002) proposed a single-source Surface Energy Balance System (SEBS) to estimate the atmospheric turbulent fluxes from satellite earth observation data in combination with meteorological information. SEBS consists of several modules to estimate net radiation and ground heat flux, and the partitioning of available energy into latent heat and sensible heat fluxes.

Remotely sensed data can also be assimilated into physically based land surface scheme (LSS) models. Compared with the semiempirical approaches, LSS models provide a detailed description of soil and vegetation canopy processes, fluxes, and states, instead of a limited number of final fluxes. Moreover, it is possible to implement procedures to assimilate data acquired by a large range of satellite-borne instruments, differing in wavelength domains, acquisition time, and geometry (Olioso et al. 1999). One category of physically based methods involves directly inputting into LSS models quantities derived from remote sensing such as surface temperature and soil moisture content. For example, Boulet et al. (2000) built a simple LSS model designed for scaling applications and remote sensing utilization, and also indicated that the analytical simplicity of the model is suitable for the investigation of spatial variability of the water and energy balance. Another category consists of updating state variables and/or adjusting parameters in LSS models through minimizing the difference between the model predictions and information derived from remote sensing. For example, Ottlé and Vidal-Madjar (1994) assimilated soil moisture inferred from infrared remote sensing into a hydrological model, which improved simulations of the energy budget. Caparrini et al. (2004) developed an approach to estimate the surface turbulent fluxes through assimilation of radiometric surface temperature sequences.

The empirical methods for estimating fluxes usually need local calibration. The methods using LSS models usually require a large set of input parameters and initial conditions describing surface properties that must be correctly specified in order to provide accurate assessment of the energy and water fluxes. Most of these properties vary in time and space, and are often assessed through in situ experiments (Demarty et al. 2005). As stated by Courault et al. (2005), one of the main problems when using LSS models is that the spatial resolution of the remotely sensed data leads to difficulties in determining the values of the “effective” parameters corresponding to the composite surface. Cleugh et al. (2007) pointed out that one of the biggest impediments to global, multitemporal satellite-based evapotranspiration monitoring is the conflicting requirement for algorithms that are biophysically realistic yet simple enough for global parameterization and implementation.

### b. Background

A new method for parameter estimation based on stationarity constraints of soil moisture and temperature has been developed and tested with field measurements at point scales (Sun et al. 2011). Salvucci (2001) initially demonstrated that the water storage change during an interval of time, conditionally averaged on soil moisture storage during that interval, tends to zero because of equilibrating tendencies in the water balance. This property can be used to estimate the moisture dependence of net drainage and evaporation from precipitation measurements (after disappearance of the soil moisture storage change term in the water balance). Applying this property, the parameters of a water balance model can be estimated by matching the expectation of modeled water fluxes conditioned on soil moisture to the expectation of precipitation conditioned on soil moisture. However, while the method robustly estimated the combined sensitivity of evapotranspiration and net drainage to the soil moisture, it had some difficulty partitioning these two terms, in particular for cases of relative dry surface soils.

Sun et al. (2011) further extended the idea to energy balance, and formulated an objective function for estimating parameters of the coupled water and energy balance system. By coupling water and energy balance, the authors expected more accurate parameter estimation and flux partitioning. After applying this stationarity-based method to a conceptually simple coupled water and energy balance model (with 16 parameters) and testing it at two AmeriFlux sites (Vaira Ranch, California, and Kendall Grassland, Arizona), their results indicated that the stationarity-based method worked quite well in estimating and partitioning turbulent fluxes, with the root-mean-square errors (RMSEs) for daily sensible heat and latent heat fluxes of approximately 15–20 W m^{−2} at both sites. Sun et al. (2012) further applied the stationarity-based method with gridded atmospheric forcing and remotely sensed surface states over the Southern Great Plains, assessing the evapotranspiration estimates and the parameter behavior in the larger domain. Salvucci and Entekhabi (2011) extended the methodology to apply it to more complex, multiple soil-layer land surface models, and tested the proposed method using the Noah land surface model (Ek et al. 2003).

The novelty of the stationarity-based method is that it does not require the measured fluxes; that is, the measured fluxes do not appear in the objective function. Compared with direct calibration parameter estimation methods (*direct calibration* hereafter denotes calibration to measured fluxes in this study), this unique characteristic of the stationarity-based method releases the burden of knowing field-measured fluxes. In addition, continuous forcing and surface state data are not required and it is not sensitive to initialization, as the method is based on statistical conditional sampling instead of integrating the model in time. This is another important feature when applying the method with remote sensing, as remotely sensed data are usually sparse because of coverage and atmosphere influences such as clouds.

In this paper, we test the stationarity-based parameter estimation method with remotely sensed soil moisture and temperature, discuss the model performance switching in situ surface states to remote sensing data, and compare the performance of the stationarity-based method and direct calibration (to measured fluxes) at three AmeriFlux sites. The paper is organized as follows: the stationarity-based method is briefly introduced in section 2; the model formulations and study sites are described in section 3; the data used are described in section 4; the model test results are presented in section 5; the results are discussed in section 6; and a summary is presented in section 7.

## 2. Stationarity-based method for parameter estimation

The method, developed based on the stationarity characteristic of soil moisture and temperature, is briefly described below. Readers are referred to Sun et al. (2011) for the complete derivation.

For the water balance (WB), the conditionally averaged rate of change of estimated moisture storage () can be written as

In (1), is the model-estimated rate of change of storage (using an estimate of the parameter vector denoted ), is the true storage change rate (using the true parameter vector ), represents error in the estimated rate of storage change due to misspecification of the parameter vector, and represents error associated with the model formulation.

Salvucci (2001) demonstrated (through mathematical derivation, Monte Carlo studies and observations) that . Given that vanishes, and under the assumption that —that is, assuming the models are flexible enough with respect to moisture that model structural error will vanish—Eq. (1) becomes

In (2), *P* is precipitation, ET is evapotranspiration, and *Q* is the combined water losses resulting from surface runoff and drainage. ET and *Q* are modeled through the atmospheric forcing, surface states, and an estimate of the parameter vector (). Model structure and the vector of free parameters are listed in the appendix.

For the energy balance (EB), a similar derivation applying stationarity of surface temperature () yields

In (3), is the net radiation, *H* is the sensible heat flux, *LE* is the latent heat flux, and *G* is the ground heat flux; is the effective heat capacity of the medium. Note that is now augmented to include the parameters specific to the energy balance flux components (see the appendix).

Combining the error terms from both the water [Eq. (2)] and energy balance [Eq. (3)], we set up an objective function in terms of the two stationarity terms and as follows:

In Eq. (4), the two tendency terms that approximate the moisture- and temperature-dependent errors have been expressed in common units of watts per square meter by multiplying the moisture term by the product , where is latent heat of vaporization (taken as J kg^{−1}) and is density of water (taken as 1000 kg m^{−3}).

Using Eqs. (2) and (3), the parameters of the land surface energy and water balance model can be estimated by minimizing the weighted objective function [Eq. (4)]. The novel characteristic of this approach is that the parameters governing the individual water and energy fluxes can be estimated without measuring each flux. Instead, only “forcing” terms (net radiation, precipitation, wind speed, humidity, and air temperature) and state variables (moisture and temperature) are required.

We briefly introduce the optimization process here. More details are described by Sun et al. (2011). Taking the water balance [Eq. (2)] as an example, we first divide the range of the conditioning variable (the observed soil moisture, S, in this example) into eight bins with the same increment. We then calculate the arithmetic mean of all the values [i.e., , , ] for some set of parameters, for which *S* falls in a given bin. In this way we reduce the entire time series of , found by evaluating Eq. (2) for the specific model, to eight conditionally averaged values, which are further weighted by the number of observations in each bin. To arrive at Eq. (4), we simply follow the same steps for the energy balance [Eq. (3)] and then combine the squared conditional means of the tendency terms [as shown in Eq. (4)]. To minimize the objective function, we select the generalized pattern search (GPS) algorithm, found in the Genetic Algorithm and Direct Search Toolbox in MATLAB, as our main optimization strategy. In fact, Eq. (4) can be further decomposed into seasonal and perturbation terms, allowing for more robust estimation of parameters by highlighting disparate temporal scales of variability in the error terms. See Sun et al. (2011) for a complete description of how to evaluate the terms in Eq. (4).

## 3. Models and sites description

### a. Model

The model used in this paper is that used by Sun et al. (2011). A parallel resistance system with a Jarvis-type stomatal conductance is used to model land surface fluxes in terms of moisture and temperature gradients. The Richardson number is used for atmospheric stability correction. The Reynolds number and a free parameter (*C*) are used to account for the difference between the aerodynamic and skin temperature. See the appendix and refer to Sun et al. (2011) for the details of this model formulation. There are in total 16 parameters (listed in Table 1) to be estimated in the coupled water and energy balance models. Note that we do not simulate or propagate the surface state variables. Instead, the observations (in situ or remotely sensed) of the surface states are used to drive the model.

### b. Study Sites

#### 1) Vaira Ranch (VR)

The Vaira Ranch site (38.4067°N, 120.9507°W, elevation 129 m), located in California, is a grazed C3 grassland opening in a region of oak/grass savanna. The soil type is rocky silt loam, and the climate is Mediterranean. Forcing and surface state data from 2004 to 2007 are used in the analysis. The annual mean precipitation and air temperature were 565.8 mm and 16.0°C, respectively, during the study period. The eddy covariance technique is used to measure the fluxes at Vaira Ranch and the other two sites.

#### 2) Kendall Grassland (KG)

The Kendall Grassland site (31.7365°N, 109.9419°W, elevation 1531 m) is in Arizona. The tower is located in a small intensively studied area within the Walnut Gulch Experimental Watershed. The soil is coarse loamy, mixed, and contains some limestone fragments. The climate is temperate and semiarid, and the landscape is a warm-season C4 grassland with a few shrubs interspersed. Forcing and surface state data from 2006 to 2007 are available and used for testing, with the annual mean precipitation and air temperature at 293.6 mm and 17.3°C, respectively, during the period.

#### 3) ARM SGP main (ARM)

The Atmospheric Radiation Measurement Program (ARM) Southern Great Plains (SGP) main site (36.6058°N, 97.4888°W, elevation of 314 m) is in Oklahoma, with the surrounding landscape used solely for agricultural purposes, which routinely undergoes four stages of development: till, plant, fertilize, and harvest. The soil type is a silty clay loam. The climate is temperate continental. Forcing and surface state data from 2005 to 2006 are chosen according to data continuity in the time series. The annual mean precipitation and air temperature were 630.2 mm and 15.5°C, respectively, during the study period.

## 4. Data

We use field-measured forcing (precipitation, net radiation, air temperature, wind speed, and relative humidity) and in situ or remotely sensed surface states (soil moisture and temperature) as inputs to the model. Forcing data are half hourly originally and averaged to daily values, creating time series for each site that range from 2004–07 at Vaira Ranch, 2006–07 at Kendall Grassland, and 2005–06 at ARM SGP main. At Vaira Ranch, the field-measured leaf area index (LAI) estimates are available and linearly interpolated to daily values. At the Kendall Grassland and ARM SGP main sites, Moderate Resolution Imaging Spectroradiometer (MODIS)-derived LAI estimates (MOD15A2, http://daac.ornl.gov/MODIS/) are used to represent the seasonality of the vegetation cover, with the original 8-day composite products disaggregated to daily values to be consistent with the daily forcing and surface state data. At all three sites, the first layer volumetric moisture content () (representing the top 3 cm at Vaira Ranch, the top 5 cm at Kendall Grassland and ARM SGP main) is transformed and used to represent the dimensionless soil moisture index *S* (i.e., ) in tests comparing the impact of remotely sensed and in situ moisture on flux prediction.

For remotely sensed surface states, we use land surface moisture data retrieved from the Advanced Microwave Scanning Radiometer for Earth Observing System (EOS; AMSR-E), which provides global soil moisture with a high temporal (day, night) and low (0.25°) spatial resolution. Many methods have been proposed to retrieve soil moisture from AMSR-E. We chose the dataset produced through the Vrije University Amsterdam in collaboration with National Aeronautics and Space Administration (VUA–NASA) (Owe et al. 2001, 2008). We chose to use land surface temperature (LST) products derived MODIS and AMSR-E (VUA–NASA) on the *Aqua* platform in this study, to be consistent with the observation time of the soil moisture data. The diurnal cycle of LST does not need to be resolved, since the stationarity-based method is driven with daily averaged forcing and surface states in this study. MODIS temperature observations have a higher spatial resolution, and potentially more accuracy, than AMSR-E (Holmes et al. 2009). On the other hand, AMSR-E temperature observations have better performance in cloudy conditions than MODIS, which leads to fewer gaps in the data. A promising method for combining the temperature observations from MODIS and AMSR-E on board *Aqua* based on the linear relationship between these two temperature products was suggested by Parinussa et al. (2008). Following their work, we developed a strategy for combining MODIS (MYD11A1, V005, http://lpdaac.usgs.gov) and AMSR-E (VUA–NASA) LST as follows: upscale the MODIS 1-km product to 0.25° to match the AMSR-E product; use upscaled MODIS as a base dataset; if a MODIS observation is missing at a time when an AMSR-E observation exist, then replace the missing MODIS data with a linear regression estimate based on the relationship between MODIS and AMSR-E derived from all days when both exist. For AMSR-E and MODIS products, the descending (~0130 LT ) and ascending (~1330 LT ) data are used for calculating a single daily (average) value. By combining the estimates from MODIS and AMSR-E, the data coverage was significantly increased at all three study sites.

Figure 1 shows the comparison of the scaled AMSR-E soil moisture and combined MODIS and AMSR-E LST with the field-measured surface states. At all three sites, the combined MODIS and AMSR-E LST appear to match the field measurement relatively well, with the coefficient of determination values above 0.9 and RMSE values of about 3°C. Regarding AMSR-E soil moisture, the best correlation and accuracy is at the Varia Ranch site, with an of 0.7 and an RMSE of 0.14. The Kendall Grassland site shows the worst accuracy (RMSE of 0.26) and the ARM SGP site shows the worst correlation ( of 0.17).

When testing the models with remotely sensed data, the effective surface temperature used in the model calculations (see the appendix) is calculated with the combined satellite LST, and the satellite broadband emissivity () is assumed to be 1 for simplicity. However, effective emissivity is still considered a free (i.e., searchable) parameter in the model. The impact of assuming that the satellite broadband emissivity is 1, but treating the emissivity in the model as free, can be seen by equating the two resulting estimates of longwave outgoing radiation and solving for an effective temperature ():

Now, assume = 1,

is the bulk surface temperature used through the models ( in section 2 and the appendix), and is a searchable free parameter. In essence, is a factor for adjusting the combined LST from the two satellite estimates to account for estimation bias, especially due to representing daily average LST from the average of two instantaneous values.

## 5. Results

By minimizing the stationarity-based objective function [Eq. (4) in section 2], the parameters are estimated for each model at each study site. The estimated parameters (reported in Table 2) are then used to calculate the water and energy fluxes. For comparison, the results of running the model with field-measured surface states and forcing are presented here as well. In this case, surface temperature is calculated from the measured longwave outgoing radiation.

In addition, in order to compare with the standard model calibration methods, we estimate model parameters (and calculate the resulting fluxes) by minimizing a mean squared error between the three measured and predicted energy fluxes (denoted below as a best-fit objective function) using

In both cases {parameter estimation via stationarity [Eq. (4)] or direct calibration to measured fluxes [Eq. (7)]}, after the parameters are estimated, we re-estimate the model fluxes using an implicit land surface temperature found by imposing the energy balance on each day, similar to the standard Penman and Penman–Monteith approach. This implicit land surface temperature closely tracks the measured values but is less noisy, thereby resulting in better flux estimates (Sun et al. 2012). The reestimated model fluxes are still evaluated using measured (in situ or remotely sensed) soil moisture data.

### a. Vaira Ranch Site

The predicted fluxes are first compared with the observed fluxes as conditionally averaged values (Fig. 2), in which the average dependence of water and energy fluxes on soil moisture and surface temperature can be visually observed. The top row (Figs. 2a,b) displays the results of estimating the model parameters with field-measured surface states, and the bottom row (Figs. 2c,d) displays the results using remotely sensed surface states. Solid lines are conditional averages of the model predictions, and dashed lines are conditional averages of the measured fluxes.

In both cases, the energy (denoted EB) and water (denoted WB) fluxes are partitioned quite well and the dependence of each flux on the state variables is reasonable. For example, the latent heat flux *LE* (red line in the EB panel, Fig. 2) increases along with surface temperature in cold days, reaches a peak around 18°C, and then decreases with surface temperature until approaching zero because of water stress. There is some mismatch for the sensible heat flux *H* (green line in the EB panel, Fig. 2) in the cold season. Comparing Figs. 2a and 2c, *H* (in the colder bins) appears better estimated when driving the models with remotely sensed data. This could indicate that the remotely sensed estimate of surface temperature is more spatially representative of the flux footprint than the in situ estimate. Note that in order to balance the net radiation for stationarity of temperature, *LE* is underestimated in the cold bins to offset the overestimate of *H*.

While it is informative to display the behavior of the fluxes in the conditional average plots, a good match on conditional means in each bin does not necessarily correspond to good predictive analysis criterion (e.g., as measured by the and/or the RMSE). Thus, the turbulent fluxes are also compared in scatterplots with and RMSE calculated (see Fig. 3) and time series plots (see Fig. 4).

In Fig. 3, the top two rows show results from stationarity-based parameter estimation [using Eq. (4)] and direct calibration [using Eq. (7)], with field-measured soil moisture and temperature. The bottom two rows show the same cases, but using remotely sensed soil moisture and temperature. The plots are arranged in the same way in Fig. 4. Testing with field-measured surface states, the model performance (in estimating turbulent fluxes) are similar in both cases, with differences of less than 0.1 and a difference of RMSE of about 3 W m^{−2}. See Sun et al. (2011) for discussions on the possible reasons of the differences. Testing with remotely sensed surface states, the for turbulent fluxes is almost the same in both cases, and the difference of RMSE for *LE* and *H* is nearly null.

For the stationarity case, when switching from field-measured to remotely sensed surface states, the model performance is not degraded at all. Instead, while remains about the same, the RMSE for *LE* and *H* is approximately 2 W m^{−2} smaller with remotely sensed surface states than those with field-measured ones. This is consistent with what we see in the conditional average plots of Fig. 2. Table 3 lists RMSEs for each energy flux in the stationarity case and the calibration case, as well as difference between the two cases.

### b. Kendall Grassland Site

Figure 5 displays the same analysis as Fig. 2 for the Kendall Grassland site. Again, the fluxes are partitioned and estimated quite well with both field-measured and remotely sensed surface states. There are some mismatches in the extreme bins (very hot, very cold, very wet, very dry), but usually there are far fewer data points in these bins. At this site, 3% of the data fall in the last two wet bins in the water balance. In the energy balance, 6% of the data fall in the first two cool bins. The dependence of the water and energy fluxes on soil moisture and temperature is again reasonable. For example, *LE* remains small and insensitive to temperature on cold days, probably because of limited soil moisture and available energy. It increases and peaks (in the warm season) with increasing temperature coinciding with increasing soil moisture (coming from monsoon precipitation). At very high temperatures, it decreases again because of limited soil moisture. To a certain degree, the temperature axis in Figs. 2 and 5 can be considered to represent seasonal changes (as we describe here), but care should be taken in such interpretations as the data that get averaged to form the conditionally averaged fluxes are determined solely based on temperature and not day of year. Thus, a hotter-than-average day in the cold season could be averaged with an abnormally cool day in the warm season.

Values of and RMSE for the turbulent fluxes are calculated and shown in the scatterplots for each case in Fig. 6 and time series plots in Fig. 7. Testing with field-measured soil moisture and temperature, the model performs quite well, with near 0.9 and RMSE values under 14 W m^{−2} for *LE* and *H* for both stationarity- and calibration-based estimates. Switching to remotely sensed surface states, the decreases to approximately 0.8 and the RMSE increases to 14–17 W m^{−2} for *LE* and *H*. There are minimal differences in and RMSE between the stationarity case and the corresponding calibration case, which indicate that the stationarity-based method works about as well as calibration, despite the fact that it does not require measurements of *LE*, *H*, or *G*.

### c. ARM SGP site

The estimated fluxes at the ARM SGP main site have been compared with the observed ones in conditional average plots (Fig. 8), scatterplots (Fig. 9), and time series plots (Fig. 10). From Fig. 8, the behaviors of *H* and *LE* are similar to that found at Kendall Grassland, but the magnitudes of *H* and *LE* are much closer at this site (i.e., the Bowen ratio is closer to one on average). The model captures the dependence of the water and energy fluxes on moisture and temperature, though there is a consistent bias of small magnitude in each bin. Similar to Kendall, only about 5% of the data fall in the hottest bin and 24% of the data fall in the two hottest bins. Switching from field-measured surface states to remotely sensed ones, the situation gets worse; that is, the magnitude of bias, in overestimating *H* (underestimating *LE*) in the cold end and underestimating *H* (overestimating *LE*) in the hot end, becomes larger.

From Fig. 9, we see that the and RMSE between the measured and modeled turbulent fluxes are generally not very good for this site. Testing with field-measured surface states, is near 0.5–0.7 and the RMSE is around 17–23 W m^{−2}. Switching to remotely sensed surface states, the decreases to 0.4–0.6 and the RMSE increases to 19–25 W m^{−2}. Again, the differences in and RMSE between the stationarity case and the corresponding calibration case are generally small, with the difference of near 0.1 for *LE* and 0.2 for *H*, and the difference of RMSE around 3–4 W m^{−2} for both *LE* and *H*. These results indicate that the relatively poor performance is more likely due to model structure and/or data issues, rather than robustness issues of the stationarity-based method.

## 6. Discussion

We use the error criterion, RMSE, as a general guideline for evaluating the model performance at each site. To better view the change of model performance by replacing field-measured surface states with remotely sensed data, we compare the RMSEs for the energy fluxes (for parameters estimated using the stationarity-based method) and display them in Fig. 11. Evaluating *G* (bottom box in Fig. 8), the model provides generally consistent performance at three sites. The difference between the RMSEs, when switching field-measured to remotely sensed surface states, is near or less than 1 W m^{−2} at any site. There are relatively larger differences in the RMSE values associated with *H* (middle panel) and *LE* (top panel). The difference of RMSE for *LE* between cases with field-measured and remotely sensed surface states is small (less than 1 W m^{−2}) at the ARM SGP site, but larger (about 2–4 W m^{−2}) at the Varia Ranch and Kendall Grassland sites. Regarding the RMSE for *H*, the differences are around 2 W m^{−2} through the three sites.

The degradation in model performance can be explored by considering the surface measurement errors in the context of the model structure and complexity. At the Varia Ranch site, replacing the field-measured surface states with remotely sensed data leads to a small decrease of the RMSE values for *LE* and *H*. This could possibly be due to the better representation of the remotely sensed estimate of surface temperature in the flux footprint than the in situ estimate, especially in the lower temperature end (see Fig. 2, EB panels; Fig. 1, VR site). At the Kendall Grassland site, the model performance degrades by switching surface states from field measurement to remotely sensed data, with the RMSE for *LE* increased by 4 W m^{−2} and that for *H* increased by 2 W m^{−2}. This relatively larger degradation of model performance for *LE* than *H* is likely due to the poor accuracy of remotely sensed soil moisture data (largest RMSE through all three sites, see Fig. 1) as well as the poor satellite coverage at this site. AMSR-E observations exist for only 57% of the time through the 2-yr time series at the Kendall Grassland site, while about 80% of AMSR-E soil moisture data are available in the 4-yr time series at the Vaira Ranch site, and 75% data in the 2-yr time series at the ARM SGP site. At the ARM SGP site, the model performance degrades by using remotely sensed surface states instead of field-measured ones as well. But the RMSEs for turbulent fluxes increase in much smaller magnitude, by 0.5 W m^{−2} for *LE* and 1.5 W m^{−2} for *H*. The performance degradation for *H* is larger than that for *LE*. This is possibly because the error introduced by remotely sensed data can be compensated for by the parameters in the *LE* model (in total 11 parameters involved), while more errors remain in *H* because of less flexibility (a simpler model formation with three searchable free parameters). Overall, however, the degradation of model performance using remotely sensed surface states is not substantial.

Comparing the stationarity case and calibration case, as shown in Table 3, the sum of the RMSE difference for energy fluxes is around 6 W m^{−2} using field-measured surface states and 5 W m^{−2} using remotely sensed data averaged through all three sites. Considering the substantial advantage of the stationarity-based method versus calibration method—that is, the stationarity-based method does not require measured turbulent fluxes for calibration [see Eqs. (4) and (7)]—we regard these increases in the RMSE values as acceptable. Furthermore, the stationarity-based method works quite well using remotely sensed surface states, with the highest RMSE of about 22 W m^{−2} for *LE* and 25 W m^{−2} for *H* at the ARM SGP site. For perspective, an assessment of some 30 published validations of remote sensing–based estimated flux against ground-based measurements of evapotranspiration shows an average RMSE value of about 50 W m^{−2} for instantaneous fluxes (Kalma et al. 2008). In the context of daily fluxes, the order of 20 W m^{−2} for RMSEs is considered to be acceptable.

Even though the presumed physical meaning of the model parameters are unchanging, their influence on model behavior is not necessary similar across the sites for the same model or across the models for the same site (Hogue et al. 2006). In our study, we apply the same model to different study sites with different datasets. From Table 2, we see that the optimal parameters at each site are different when using in situ versus remotely sensed surface states. For example, the soil hydraulic conductivity (*k*) increases about 30% at the ARM SGP main site and nearly doubles at the Vaira Ranch and Kendall Grassland sites, when switching from field-measured surface states to remotely sensed data. This indicates that pre-assigning the parameter values from lookup tables (which are based on soil and biome types) or even through field experiments might not be always suitable, especially when applying the model with remote sensed data in larger scales. Possible reasons include heterogeneity of the land surface, data issues related to remote sensing, and, in general (but not in our study), different model formulations.

In some cases in Table 2, estimated values of individual parameters differ from typical values. For example, the estimated value of the vegetation height () is the upper boundary of the value range (100 cm) in most cases. However, this parameter-hitting boundary issue is not due to the stationarity-based methodology, but rather the highly simplified model itself and/or measurement errors of the forcing data, because the issue also happens in the case of direct calibration by simply matching the modeled fluxes to observations.

To further explore the issue, we reran the model [matching modeled fluxes to the observations as in the direct calibration case through minimizing Eq. (7)] with field-measured surface states, but this time we extended the upper boundary of (from 100 to 300 cm) while maintaining other parameter boundaries. We did this for the Kendall Grassland and ARM SGP main sites. The estimated parameter values and the RMSE of the fluxes are shown in Table 4. There is nearly no change in the RMSE of the fluxes at both sites after extending the upper boundary of . At Kendall Grassland site, the chosen value of is close to the upper boundary and nearly triples the original chosen value. This large difference in the chosen value, however, does not affect the RMSE much (Table 4). In fact, the effect is compensated through adjusting other parameters, such as *A*, the optimum air temperature (), and parameter *C* (associated with ). The situation is similar at the ARM SGP main site. The chosen value of hits the upper boundary again at 300 cm, but the effect of increasing is compensated by the change of chosen (influencing soil evaporation). The fact that a huge change of value does not impact the fluxes also illustrates that is not a sensitive parameter in the model. The upper boundary of 100 cm still seems high for . But considering that the role of in the model is mostly in setting the aerodynamic roughness, larger values could be accounting for roughness because of rolling hills and interspersed shrubs (Sun et al. 2011). Thus, we think 100 cm should be a reasonable value for the upper boundary of .

To examine and control the well-known equifinality issue in parameter estimation problems, we minimize the objective function though different optimization methods, such as grid search, gradient-based method (e.g., the fmincon function in MATLAB), the generalized pattern search (GPS) algorithm, and then check the parameters and objective function values in each case. In our testing, results from a single application of the GPS algorithm were essentially identical to those from applying fmincon with hundreds of random initial points; see Sun et al. (2011) for more discussion. Sun et al. (2012) presented the optimized parameter values (applying the same optimization procedure grid by grid) in a larger domain. The smooth pattern of the parameter fields also indicates that we achieved the global minimum, and no major equifinality issue appeared in our results.

## 7. Summary

In this paper, we tested the stationarity-based method of parameter estimation with remotely sensed surface states (AMSR-E soil moisture, combined MODIS and AMSR-E LST), and analyzed the degradation of the model performance by using the field-measured surface states and remotely sensed states at three AmeriFlux sites. In addition, we have contrasted what we classify as “direct-calibration-based parameter estimation” with what we term “stationarity-based parameter estimation.” The former we define here as searching a parameter space in order to minimize a mean square prediction error [e.g., as in Eq. (7)]. The latter minimizes an indirect estimate of mean square prediction error [ in Eq. (4)]. The stationarity-based parameter estimation could, therefore, also be classified as calibration. However, we feel that the distinction between Eq. (7), which requires measured fluxes, and Eq. (4), which does not, is of such practical importance that it is useful to classify the methods differently.

Summarizing the results from the analysis using field-measured and remotely sensed surface states under stationarity-based parameter estimation, there is limited degradation of model performance in estimating turbulent fluxes because of changing data sources of surface states. Furthermore, the performance degradation is minimal when using the stationarity-based parameter estimation method instead of calibration to measure fluxes. Compared with the direct calibration method, the significant advantage of the stationarity-based method is that it releases the burden of knowing the measured fluxes (measured fluxes do not appear in the objective function). Because we are not simulating over time, our proposed method has other distinct advantages—for example, it is not sensitive to initialization, it can easily handle gaps in forcing and surface states data (a common occurrence), and it can be applied using daily averaged data. The stationarity-based method is derived only from stationarity characteristic of surface states and conservation statements (energy and water) at the land surface, and thus it can be applied at the scale appropriate to the input data and the model.

## Acknowledgments

This study was funded by the NASA Terrestrial Hydrology Program under Grants NNG06GE48G, NAG5-11695, and NNX12AP78G.

### APPENDIX

#### Model Formulation

The model is derived, in part, from components of the Noah land surface model (Chen et al. 1997), the Mosaic model (Koster and Suarez 1996), and the Simple Biosphere Model (Sellers et al. 1996).

For energy balance, is a forcing input taken from field measurement.

The *H* and *LE* heat fluxes are modeled using bulk difference equations as follows (e.g., Garratt 1992):

In Eqs. (A1) and (A2), is the density of air (approximated here as 1.24 kg m^{−3}); is the heat capacity of the air (1005 J kg^{−1} K^{−1}); *P*_{a} is the surface atmosphere pressure; is the wind speed measured at the same height *z* (typically 2 m); is the bulk surface temperature, where we use the skin (i.e., radiometric) temperature in the model, which is estimated from outgoing longwave radiation—that is, , where is the observed outgoing longwave radiation, is the Stefan–Boltzmann constant ( W m^{−2} K^{−4}), and is the emissivity parameter to be estimated. In addition, is the air temperature at height *z*; and are the vapor pressure at the same levels as and , respectively; and and are the corresponding specific humidity values, respectively. The surface vapor pressure () is assumed saturated at the surface temperature ().

The and are the bulk transfer coefficients for heat and water vapor, respectively. Here we use the modified Louis model suggested by Chen et al. (1997): is parameterized following the modified Louis approach given by Mahrt (1987) for the stable boundary layer case and Holtslag and Beljaars (1989) for the unstable case—that is,

In Eq. (A3), the function is given by

for stable conditions and by

for unstable conditions.

For the calculation of Richardson number (Ri), we use (Dingman 2002)

In Eqs. (A4)–(A6), *a* is a constant (set equal to 1 in the model); is the von Kármán constant, set to 0.4; and is the height at which the wind speed and air temperature are measured. Following common practice (Dingman 2002), the zero-plane displacement is approximated as and the roughness height for momentum is approximated as . In our model, we estimated vegetation height as a free constant parameter.

Following Chen et al. (1997) we specify the relation between and as

In Eq. (A8), is the kinematic molecular viscosity ( m^{2} s^{−1}), Re is the roughness Reynolds number, is the surface friction velocity, approximated (without stability correction) as , and *C* is a free parameter.

In the equation, is further modified based on the formulation of . For *LE*, we need to consider canopy and soil resistance to water vapor and how these resistances depend on soil moisture, vapor pressure deficit (VPD), and temperature conditions. To keep the model as simple as possible, we define another surface conductance , which is composed of canopy conductance and soil conductance for the bulk surface area:

In Eq. (A9), represents canopy conductance, and we apply a Jarvis-type stomatal conductance (Jarvis 1976) that is then scaled by LAI and is associated with three environmental stress terms, defined as

Here *A* is a free parameter accounting for maximum leaf conductance. The first environmental stress function characterizes the influence of soil moisture, defined as

which is a transpiration efficiency meant to represent stomatal conductance limitations under dry soil conditions. Two parameters ( and ) are used to represent as a piecewise linear function. The soil moisture stress function is evaluated relative to the dimensionless soil moisture index , where represents the measured soil moisture data and is the maximum of soil moisture data in the studied time series. This two-parameter model for soil moisture control on stomatal conductance is flexible, allowing for control of the slope (i.e., sensitivity of conductance to moisture and for control of the range of moisture over which limits conductance).

The second stress term representing the influence of VPD is defined as the same formation in Noah land surface model (Ek et al. 2003; Alfieri et al. 2008),

where *g* is a species-dependent empirical parameter. In Noah LSM, *g* is taken as a constant. Here it is a free parameter.

The third environmental stress term is related to temperature, following Jarvis (1976), Matsumoto et al. (2005), and Samanta et al. (2008), defined as

where , and representing the optimal, lowest, and highest temperatures for transpiration, respectively; and are fixed as 0° and 50°C, respectively, while is a free parameter.

In Eq. (A9), represents soil conductance, defined as

In Eq. (A13), *B* is a free parameter accounting for maximum soil conductance; is scaled by effective soil area index (SAI), which is fixed to be equal to 1; and is an evaporation efficiency meant to represent bare soil evaporation limitations under dry soil conditions. Similar to the soil moisture stress function in canopy conductance, is dependent on the dimensionless soil moisture index (*S*), but with two distinct parameters and to be estimated—that is,

For the combined soil and vegetation area, *LE* can now be written as

The following simple model is adopted for the ground heat flux:

where is deep soil temperature, which is estimated as the time average of ; and is the apparent conductivity of the soil and is treated as a searchable parameter.

In the water balance, *P* is a forcing input. ET is converted from *LE* through . The net drainage and runoff term (*Q*) is represented as follows:

where *k*, *w*, *c*, and *n* are parameters that account for both percolation when the soil is wet [first term of Eq. (A17) dominates second term], and capillaries rise from deeper layers when the soil is dry (second term dominates first).

## REFERENCES

*Physical Hydrology*. 2nd ed. Prentice Hall, 646 pp.

*The Atmospheric Boundary Layer*. Cambridge University Press, 316 pp.

*ECMWF Workshop on Parameterization of Fluxes and Land Surface,*Reading, United Kingdom, ECMWF, 121–147.