The temporal evolution of radiative transfer parameters vegetation opacity and soil surface roughness is estimated by assimilating L-band brightness temperatures. The L-MEB model is used within a particle filter in synthetic applications with increasing complexity. Finally, a real-world experiment is conducted with multiangular SMOS observations covering a SCAN site.

Abstract

ESA’s Soil Moisture and Ocean Salinity (SMOS) mission has been designed to extend our knowledge of the Earth’s water cycle. Soil Moisture and Ocean Salinity records brightness temperatures at the L-band, which over land are sensitive to soil and vegetation parameters. On the basis of these measurements, soil moisture and vegetation opacity data sets have been derived operationally since 2009 for applications comprising hydrology, numerical weather prediction (NWP), and drought monitoring. We present a method to enhance the knowledge about the temporal evolution of radiative transfer parameters. The radiative transfer model L-Band Microwave Emission of the Biosphere (L-MEB) is used within a data assimilation framework to retrieve vegetation opacity and soil surface roughness. To analyze the ability of the data assimilation approach to track the temporal evolution of these parameters, scenario analyses were performed with increasing complexity. First, the HYDRUS-1D code was used to generate soil moisture and soil temperature time series. On the basis of these data, the L-MEB forward model was run to simulate brightness temperature observations. Finally, the coupled model system HYDRUS-1D and L-MEB were integrated into a data assimilation framework using a particle filter, which is able to update L-MEB states as well as L-MEB parameters. Time invariant and time variable radiative transfer parameters were estimated. Moreover, it was possible to estimate a “bias” term when model simulations show a systematic difference as compared to observations. An application to a USDA-NRCS Soil Climate Analysis Network (SCAN) site showed the good performance of the proposed approach under real conditions.

The objectives of ESA’s SMOS mission over land are to provide soil moisture observations for weather forecasting, climate monitoring, and investigating the global freshwater cycle (Kerr et al., 2001, Mecklenburg et al., 2012). Soil Moisture and Ocean Salinity has been providing two-dimensional brightness temperature images of the Earth surface at a frequency of 1.4 GHz (L-band) since it was successfully launched in November 2009 (Mecklenburg et al., 2012). Raw measurement performances vary between 30 km and more than 50 km for ground resolution. The radiometric sensitivity (1.2 s integration time for each polarization) ranges between 0.8 to 3 K, whereas temporal sampling varies between 1 and 3 d, depending on latitude, nature of the target, and location within the instrument field of view (Kerr et al., 2001). Finally, the images are translated into soil moisture maps using the Level 2 soil moisture retrieval algorithm (Kerr et al., 2012) whereby the target accuracy in soil water content is 0.04 m3 m−3.

The SMOS Level 2 processor generates soil moisture from SMOS brightness temperatures and is based on the microwave forward model L-MEB proposed by Wigneron et al. (2007). In general, L-MEB is used to compute top-of-the-atmosphere brightness temperatures based on a set of geophysical state variables (e.g., soil moisture, soil temperature, vegetation water content, vegetation temperature), static state parameters (e.g., soil texture described through sand and clay fraction), and parameters describing the radiative transfer (e.g., single scattering albedo, polarization coupling factor, effective surface roughness). The retrieval as outlined in Kerr et al. (2012) is based on a Bayesian approach where modeled brightness temperatures are compared against the observations through a cost function.

The cost function is then minimized by changing soil moisture and vegetation opacity. The additional state variables and parameters are either obtained through auxiliary data products, for example, forecast and analyses fields from ECMWF, or from the peer-reviewed literature (Kerr et al., 2010). Prior to launch, the optimization of L-MEB for different surfaces has been addressed by numerous studies, most of which were based on the analysis of ground-based and airborne L-band data, for example, de Rosnay et al. (2006), Grant et al. (2007), Guglielmetti et al. (2008), Juglea et al. (2010), Panciera et al. (2008), Saleh et al. (2009), Wigneron et al. (2000), and Merlin et al. (2009). During the life time of the mission, the processor has been continuously improved based on extensive validation work. The retrieval is state-of-the-art in operational ground segments and is already very complex in that (i) land surface heterogeneity within the SMOS footprints is taken into account and (ii) the cost function uses brightness temperature observations at multiple polarizations and incidence angles to make full use of the observations available.

The accuracy of the retrieval can be limited through uncertainties (i) in the observations, (ii) in the microwave emission model, (iii) in the prescribed state variables, and (iv) in the parameter fields. These errors can be systematic and static (e.g., a wrong soil classification at a given location), systematic and temporarily varying (e.g., a seasonal bias in the soil temperature data), or random (e.g., instrument noise). These uncertainties in the observations, the forward model, and the auxiliary data will manifest themselves in systematic and random errors in the retrieved state variables, that is, soil moisture and vegetation water content. Random uncertainties are partly considered by the SMOS Level 2 processor in that the accuracies of the observations and the retrieved state variables can be specified. However, systematic uncertainties can hardly be addressed in current retrievals.

We employed a sampling importance resampling–particle filter (SIR-PF) to demonstrate the potential of data assimilation techniques for retrieval applications. The approach was tested by using a synthetic case study with three different scenarios. The advantage of performing a synthetic study is that the true system is exactly known and is therefore not prone to model inaccuracies or random observation errors. Although we were working with synthetic data mimicking the SMOS mission, the methodology is transferable to other satellite observations.

The first objective was to demonstrate the simultaneous retrieval of a state variable and a parameter. The use of data assimilation for state estimation is a common technique in NWP and terrestrial hydrology where the accuracy of the initial conditions partly determine the forecast skill (Palmer and Hagedorn, 2007). On the other hand, data assimilation has been successfully used for parameter estimation in subsurface hydrology (Evensen, 2009, Han et al., 2012, 2013, Kurtz et al., 2012, Montzka et al., 2012, Moradkhani et al., 2005b). We focused on the retrieval of vegetation opacity and the effective surface roughness parameter. Next to soil moisture, both parameters were found to have the biggest influence on the brightness temperature at low microwave frequencies (Jones et al., 2004, Vereecken et al., 2012).

Soil moisture—among the other input variables—is a prescribed state variable with a given uncertainty that is taken into account as well as the uncertainties of the retrieval variables and the observations. Since global soil moisture estimates and their error estimates are also available from other satellite missions and NWP forecasting systems (Balsamo et al., 2009), this setup could be of potential interest for the vegetation and carbon modeling communities (Kaminski et al., 2012).

Systematic differences between the modeled first guess brightness temperatures and the observations are often referred to as “biases,” although they can vary in space and time. In data assimilation applications in NWP, “bias corrections” have received a lot of attention in the past (Eyre, 1992, Harris and Kelly, 2001, Watts and McNally, 2004). Recently, adaptive “bias” correction schemes have been introduced, which update the systematic differences inside the assimilation scheme (Auligne et al., 2007; de Lannoy et al., 2007). We present a similar approach in the SIR-PF framework.

In summary, the objectives of this study were to (i) introduce a retrieval algorithm based on the SIR-PF, (ii) test its applicability by retrieving simultaneously temporarily constant parameters (vegetation opacity and soil surface roughness), (iii) analyze the ability of the system to track the temporal changes of geophysical state variables, (iv) account for systematic differences between the modeled first guess and observations, and (v) demonstrate the applicability to real-world data. For this last objective the version 5.04 reprocessed SMOS L1C data for 2010 were used.

Methodology

We developed a one-dimensional analysis system to assess the performance of the data assimilation system and the parameter retrieval accuracy. The integrated analysis system comprises the hydrological model HYDRUS-1D, the L-MEB forward operator, and the SIR-PF. In this section, we provide a short overview of these key components for the analysis system. Details on the generation of the synthetic data scenario are given in the subsequent part of the paper.

HYDRUS-1D

For the hydrological simulator, the HYDRUS-1D model (Šimůnek et al., 1999, 2008) was selected to generate reference soil moisture and soil temperature replicates. This model solves the Richards equation for unsaturated water flow: 
graphic
where θ is the volumetric soil moisture (m3 m−3), t is time (day), h is the matric potential in head units (cm), z is the depth (cm), and K(h) is the hydraulic conductivity as a function of matric potential (cm d−1). Soil water retention θ(h) and hydraulic conductivity K(h) functions are described by the Mualem–van Genuchten approach (van Genuchten, 1980). Heat transport is calculated according to Sophocleous (1979) by: 
graphic
where T is the soil temperature (°C), λ is the thermal conductivity of the soil (kg cm−3 d−1 °C−1), Cw and C are the volumetric heat capacities (kg d−2 cm−1 °C−1) of the liquid phase and porous medium, respectively, and Jw is the water flux density (cm d−1). The hydrological model was driven with perturbed precipitation and air temperature data to generate ensembles of soil moisture and temperatures as inputs for the land surface emission model L-MEB.

The L-Band Microwave Emission of the Biosphere Model

The L-MEB is a so called τ model, where the radiative transfer in the vegetation layer is described through its optical depth τ and single scattering albedo ω. A full derivation of the τ-ω model from the general vector radiative transfer theory can be found in Drusch and Crewell (2005). L-Band Microwave Emission of the Biosphere can be used to compute top-of-atmosphere brightness temperatures Tb over various land covers at L-band (Wigneron et al., 2007): 
graphic
where p represents the measured polarization (horizontal = H or vertical = V), and inc represents the incidence angle. TV is the vegetation physical temperature and TG is the soil effective physical temperature. The reflectivity of a rough soil, Γ, which is also sensitive to the incidence angle and measured polarization, is described by the model soil parameters hs and Np following Wang and Choudhury (1981) and Wigneron et al. (2011): 
graphic
where Γ* is the reflectivity of a smooth soil, which strongly depends on soil moisture and temperature. For the computations presented in this paper, Np = 0 is assumed [Wigneron et al., 2001, 2011]. The vegetation transmissivity γ is determined by the vegetation opacity at nadir, τ, and the parameter ttp that corrects the optical depth for non-nadir views at each polarization p: 
graphic

Detailed descriptions as well as applications of L-MEB can be found in several publications (Cano et al., 2010, Panciera et al., 2009, Saleh et al., 2007, 2009, Wigneron et al., 2007, 2008). We estimated the vegetation opacity τ as well as the roughness parameter hs by using the Tb ensemble simulations in the data assimilation framework.

Sampling Importance Resampling–Particle Filter

The SIR-PF, used here for Tb data assimilation, is a commonly used particle filter that was introduced by Gordon et al. (1993) and clarified by Moradkhani et al. (2005a). The procedure is based on the sequential importance sampling algorithm as explained by Arulampalam et al. (2002). It has been successfully used in other studies (Montzka et al., 2012) with different applications including soil moisture (Montzka et al., 2011, 2013), snow water equivalent (Leisenring and Moradkhani, 2011), flood inundation mapping (Matgen et al., 2010), radiance data assimilation for operational snow and rainfall runoff modeling (DeChant and Moradkhani, 2011b), seasonal streamflow forecasting with initial condition characterized by data assimilation (DeChant and Moradkhani, 2011a), streamflow forecasting using combined SIR-PF and Bayesian model averaging (Parrish et al., 2012), and suspended sediment load prediction (Leisenring and Moradkhani, 2012).

The posterior (analysis) probability is obtained using Bayes’ law; a priori information as expressed by prior probability is combined with the knowledge from new information through a likelihood function. This Bayesian presentation uses the batch of data to estimate the posterior distribution given all the information available through a time series of data (Moradkhani et al., 2005a). However, this form makes no attempt to include information gained from new observations when becoming available. The flexibility required to use the new information is provided by a sequential Bayesian scheme. Moradkhani et al. (2005a) demonstrated that the methods based on sequential Bayesian estimation seem better able to benefit from the temporal organization and structure of information, achieving better conformity of the model output with observations. Following Moradkhani and Sorooshian (2008), if xt represents the state variable to be estimated within the Bayesian framework, owing to its stochastic nature, the pertinent information about it at time t can be obtained from the observations Yt= [y1, y2,…yt] through the recursive application of Bayes law where the posterior density of state variables can be obtained as follows: 
graphic
where the forecast density of can be computed using Chapman Kolmogorov equation (Moradkhani et al., 2005a). One issue that may arise while using the recursive Bayes law in Eq. [6] is that the multidimensional integration of forecast density is analytically impossible, making the closed form solution of Eq. [6] intractable. Therefore, employing a Monte Carlo approach is a practical solution to this problem. The most rigorous Monte Carlo procedure in obtaining the state variables conditional distribution is to use the particle filtering.
Two sequential estimation operations are recognized in filtering applications: first is the forecasting step, which is the transition of state variables from one observation time to the next as shown in Eq. [7]: 
graphic
where M represents the radiative transfer model, χt the model states, par the model parameters, and εt the model error. The second step is the analysis (updating) of the forecasted (a priori) states with the new observations. Ensemble procedures offer a practical alternative to an exact Bayesian solution by relying on discrete estimation of forecast and update densities through a set of random variables and corresponding weights: 
graphic
 
graphic
 
graphic
where the forecast and update densities are presented by summation of N Dirac delta functions, where N is the number of particles, and xi and wi denote the state and weight assigned to the ith particle, respectively. The weights before and after updating are indicated by negative and positive signs in the superscript.
The forecasting step in particle filtering is similar to other filtering methods where the model is evolved for each particle with equal weight: 
graphic
where ut are uncertain model inputs. Therefore, the forecast density will be written as: 
graphic
Until resampling, the updated (particles) are kept the same as only the weights, not the forecast values, are updated. Therefore, 
graphic
and from Eq. [12], the filtering posterior, forumla is calculated as follows: 
graphic
where C is the normalizing constant. The likelihood function forumla, which measures the likelihood of a given model M through its misfit e(.), can be computed using: 
graphic

The SIR-PF radiative transfer parameter estimation is performed as follows: HYDRUS-1D is used to calculate 1000 soil moisture and soil temperature realizations. Moreover, the samples of air temperature are used as vegetation physical temperature TV, which is valid for the calculation of the emissivity in the morning hours (SMOS overpass ∼6:00 h). The L-MEB parameters are sampled uniformly in their feasible range (see Table 2). Here the Latin hypercube sampling (LHS) (Iman et al., 1981) without any further constraints is used to generate the plausible initial parameter sets. Then model states replicates and the parameters ensemble are incorporated in the L-MEB model and forwarded in time. After calculating the weight vector for the model states, the resampling is applied not only for the states but also for the parameters. We adopted the procedure by Moradkhani et al. (2005a), who recommended a minor parameter particle perturbation after each assimilation step to avoid sample impoverishment. Here we used a minor perturbation value of 2% of the prior parameter range.

Synthetic Data Scenario Analysis

We developed a one-dimensional analysis system to assess the performance of the data assimilation system and the parameter retrieval accuracy. As described above, the integrated analysis system comprises the hydrological model HYDRUS-1D, the L-MEB forward operator, and the particle filter. The first three scenarios use synthetic data providing idealized conditions, for example, well-defined errors; synthetic data assimilation experiments are often done before applications with real observations (Seuffert et al., 2003). All synthetic data experiments share a common baseline that is described in the following paragraphs.

The hydraulic properties for the simulation were taken from field measurements of the Terrestrial Environmental Observatories (TERENO) test site Selhausen close to the city of Jülich, North Rhine–Westphalia, Germany and are listed in Table 1. According to USDA classification, the soil is a luvic cambrisol with three horizons ranging from 0 to 0.33, 0.33 to 0.57, and > 0.57 m (Bauer et al., 2012). The soil texture is silt. A detailed description of the test site is given by Weihermüller et al. (2007). As proposed by Bauer et al. (2008), the thermal conductivity of the soil was taken from Chung and Horton (1987).

The overall simulation domain was 4 m deep with a permanent groundwater table at the lower end. The upper boundary was set to atmospheric conditions, whereby the climatic data were taken from the station in the city of Aachen (50.78°N, 6.09°E, 202 m asl) for the year 2010. The potential evapotranspiration was calculated according to Penman–Monteith equation as formulated in Allen (1998).

The precipitation and air temperature forcing data were perturbed to generate ensembles with 1000 members. Given the multiplicative nature of errors, precipitation data were perturbed with a lognormal distribution with a coefficient of variation of 0.1, whereas the air temperature was perturbed with a normal distribution with a coefficient of variation of 0.001. The true precipitation and evapotranspiration and their perturbed values are shown in Fig. 1.

The volumetric soil moisture simulated for the top soil (−2 cm) as well as simulated soil temperature at top (−2 cm) and −50 cm with their realizations are plotted in Fig. 2.

The resulting uncertainty in the modeled soil moisture data is comparable to the envisaged accuracy of 0.04 m3 m−3 for the SMOS and Soil Moisture Active Passive (SMAP) missions.

In general, the soil penetration depth δp (m) of microwave sensors such as SMOS is a function of dielectric permeability and electrical conductivity that will be mainly influenced by soil moisture and is therefore not constant in time and space. Nevertheless, Kerr (2007) stated that the penetration depth can be assumed to be approximated as the average over the first 5 cm. To simulate this condition, the computed soil moisture from the HYDRUS outputs was averaged over the first 5 cm. The parameters used for the L-MEB radiative transfer calculations are listed in Table 2. To mimic the typical temporal sampling of a spaceborne instrument operating in a sun-synchronous polar orbit we selected data every 3 d, in this case simulations for 6:00 AM local time. In addition, a random error of 2 K has been added to the brightness temperature computations to mimic instrument noise.

In all three synthetic data scenarios, brightness temperature measurements were assimilated for state update in L-MEB to retrieve vegetation opacity and soil surface roughness:

  • Scenario 1: retrieval of time invariant vegetation opacity and soil surface roughness,

  • Scenario 2: retrieval of temporally changing vegetation opacity and constant soil surface roughness,

  • Scenario 3: retrieval of temporally changing vegetation opacity and constant soil surface roughness under the presence of a systematic difference in the simulated brightness temperatures

The parameters used for this simulation are given in Table 2.

Results for Scenario 1

The hs and τ were fixed to 0.1 (i.e., the global standard value in the SMOS Level 2 processor) and 0.3, respectively. The prior uncertainty range, sampled with LHS, was given with 0 to 0.4 for hs and 0 to 0.9 for τ, respectively (see Table 3).

In Fig. 3, the simulated true Tb time series for horizontal and vertical polarization are shown as well as the corresponding filtered simulations mimicking observations from a spaceborne radiometer. Because of the changing physical temperature of the target, Tb also exhibited a strong seasonal variability. Similarly, by using simulated morning data only, their curve occurred at the lower boundary of the diurnal cycle of the simulated truth Tb.

In winter, that is, from Day of Year (DoY) 0 to DoY 80 as well as DoY 330 to DoY 365, several jumps in the Tb curves were detectable, which could be explained by freezing and thawing processes especially within the top soil layer. These freeze–thaw cycles altered the dielectric properties of the soil substantially because of water phase changes and therefore caused rapid changes in daily and seasonal Tb of 20 to 30 K. These findings were already reported by Rautiainen et al. (2012) and Entekhabi et al. (2010). As a consequence, all data points with values of top soil temperature below 273.15 K (0°C) were rejected from the assimilation to avoid misinterpretation.

In general, the particle filtering result plotted in Fig. 3 represent well the simulated observations for both H and V polarization because of the simultaneous state and parameter estimation. With 1000 model realizations at every assimilation step, the state variations between these 1000 realizations were completely driven by the parameter variations. The parameter changes over time estimated by the data assimilation are shown in Fig. 4. At t = 0, both hs and τ showed larger uncertainty in correspondence with their prior uncertainty range. In spite of the uncertainties of the soil moisture, soil temperature, and the instrument noise in the observations, the simulated true parameters were well represented after a few assimilation steps. This was especially true for τ, where the uncertainty range became relatively small. On the other hand, the uncertainty of hs was larger, but the mean of all particles was close to the simulated real value.

The results of Scenario 1 showed in general good performance; therefore, the approach seems promising for the estimation of radiative transfer parameters. This is important to gain larger knowledge about the characteristics of these parameters, for example, their temporal evolution. Moreover, with the number of realizations used (1000 particles), it is possible to represent the probability density function of states and parameters. On the basis of the results, it can be concluded that it is possible to estimate vegetation opacity and soil surface roughness using data assimilation within the model framework introduced.

Results for Scenario 2

The vegetation opacity parameter of most land cover types has a strong annual variability due to changes in the phenological stages and the water content stored within the canopy. On the basis of Scenario 1, a second scenario was developed with temporal variable vegetation opacity. By using a simple Gompertz growth function (Eq. [16]) a dynamic τ could be generated, starting with bare soil conditions (τ = 0): 
graphic
where τ increases for a hypothetical crop during the growing season (beginning 1 March = DoY 59) to maturity with a maximum of a = 0.3, and a rapid reduction of τ from 0.3 to 0 after harvest on 30 June (DoY 180). A second cycle of growth, maturity, and harvest started 29 October (DoY 301); b is the displacement in time of −100, and c is the growth rate of −0.005 during these cycles. The “true” τ calculated this way was the target value of the experiment.

In Fig. 5, the Tb for Scenario 2 is presented. Especially in vertical polarization, the two growth cycles were visible by higher Tb, which are rapidly reduced during harvest. In general, the SIR-PF was able to correctly follow the vegetation growth, but after harvest and the following period a systematic discrepancy occurred in terms of a time delay of assimilated Tb compared to the simulated “truth.” This behavior could be explained by the uncertainty evolution of the vegetation opacity and roughness parameter (Fig. 6). While τ was able to follow the dynamic growth of the vegetation, harvest and corresponding large change of opacity were too abrupt to be estimated in time. One reason for this temporal mismatch is that the assimilation was performed only every 72 h. The second reason is that after several resamplings during the stable growth season the uncertainty of τ was strongly reduced with all particles having a value of τ ≅ 0.3 and without particles with τ ≅ 0. The minor parameter perturbation, which was implemented to account for sample impoverishment, helped to change τ to a correct level. It took approximately 10 d until the correct τ was retrieved. This type of delay was expected; an immediate response would only be possible if the model predictions were very uncertain while the measurements were characterized by very small uncertainty. To a certain extent, hs compensated for inadequately estimated τ during harvest.

The results of Scenario 2 clearly showed that the data assimilation approach is able to track the temporal variability of the sensitive radiative transfer parameters that is vegetation opacity under the assumption of a constant soil surface roughness. On the other hand, extreme events such as harvest cannot be reproduced on time but will be captured after new data become available. For this specific scenario, the error in Tb for the next observation after the rapid change in τ is about 25 K in horizontal and 10 K in vertical polarization, respectively. Nevertheless, the occurrence of such events on typical radiometer pixel scales of 302 to 502 km2 is not very probable on a single day. Only for some weather events and only in few occasions the radiative transfer parameters might change within hours, and then additional data sources may help increase the parameter estimation accuracy. In general, the uncertainty reported can be accepted.

Figure 6 shows that the uncertainty of the hs retrieval was increasing with increasing τ. This is in line with Eq. [3–5], where the contribution of the soil to the microwave emission was decreasing with increasing τ. When the soil contribution to brightness temperatures observed above the canopy was low, soil parameter estimation became more difficult even at L-band. This indicated that variable environmental conditions helped to separate different parameters during the estimation procedure, which in this case were a vegetation parameter and a soil parameter.

Scenario 3: Parameter Retrieval under the Presence of a Systematic Difference

As previously discussed, the relative difference between the model simulation and the measured brightness temperatures is considered as “bias” (Auligne et al., 2007). To investigate the ability of the proposed data assimilation approach to account for this error, in Scenario 3 an artificial systematic difference in L-MEB Scenario 3 was introduced. Hereby, Scenario 3 was based on Scenario 2, but a persistent overestimation of simulated Tb by a factor of 5% was assumed. The prior “bias” uncertainty range is ± 20%, that is, Tb biased = Tb true ± 0.2Tb true.

Figure 7 shows the “true” Tb of Scenario 3. Moreover, it presents the 5% “bias” of simulated L-band Tb. In contrast to the simulated Tb of Scenario 2, these were placed at the top border of the diurnal Tb cycle, which was not plausible for SMOS ascending node data for the simulated region. Nevertheless, the SIR-PF result accounted for the systematic difference in simulated Tb observations and represented colder day times that occurred in the early mornings. This behavior was still overlaid by the deviating Tb during harvest of Scenario 2.

Similarly, the parameter development of hs and τ showed the same behavior as already plotted for Scenario 2 (Fig. 8). To a certain extent, this resulted from the very accurate estimation of the systematic difference. This is remarkable because additional tests (results not presented here) indicated that the systematic differences in Tb could also be compensated by changing the hs or τ parameters when sudden changes occurred. This hypothesis was supported by the behavior during the harvest periods, where the “bias” irregularly compensated for the τ estimation momentum, that is, one simulated measurement alone did not contain the information for an adequate estimation of radiative transfer parameters and “bias” at such an event. Only a combination of several different observations provided valuable input for parameter estimation. For example, if sudden large changes in observed brightness temperatures occur, the approach may not be able to estimate the real parameter uncertainty. In light of the reason above, we excluded the frozen soil conditions from the analysis. To a certain extent, this problem could be solved by assimilation of multi-incidence angle Tb observations in real SMOS data.

Scenario 4: Real Data Application

The Scenarios 1 to 3 showed that the performance of the approach was quite good for a perfectly parameterized data assimilation system. With Scenario 4, the approach would be analyzed with real-world data. A new challenge of this scenario was to deal with the multiple incidence angles of SMOS for one overpass. This is a large advantage of SMOS as compared to other passive microwave sensors. To this end, SMOS Level 1C data (reprocessed using version 5.04) for 2010 were selected for a site belonging to the USDA-NRCS SCAN (Schaefer et al., 2007). The SCAN site Nunn (40°51′21″ N, 104°44′30″ W, No. 2017, ID: 91CO123003) is located in the United States (Colorado), where the effects of L-band radio frequency interference (RFI) are known to be low. Here top soil moisture (−2 and −4 cm), top soil temperature (−2 cm), soil temperature in depth (−40 cm) in hourly resolution as forcing data for L-MEB (Fig. 9), as well as soil parameters are available. The soil fractions for the A horizon are given in the pedon information with sand = 70.6%, silt = 12.9%, and clay = 16.5%. The soil bulk density is given with 1.41 g cm−3. The site and its surrounding areas were covered by low vegetation, mainly a mix of grassland and agriculture.

For the site, the nearest most representative SMOS pixel (discrete global grid [DGG] 191282 for Nunn) was selected. Once the choice of DGG was made, the SMOS data for this DGG were selected. Only observations from within the narrow swath were selected to maximize the number of observations available for retrieval. Furthermore, if a “bias” occurred, it might change with the incidence angle because of limitations of the radiative transfer model. Therefore, only observations at incidence angles between 40° and 45° were taken into account, which means a surface area of approximately 65 km diameter was covered (representing the spatial resolution at −3dB antenna gain). In addition, focusing on observations taken under a limited range of viewing angles reduces uncertainties introduced through different spatial resolutions of the SMOS observations. The approach followed in this study did not simulate the instrument and the antenna, therefore the resulting values of hs and τ were unweighted values for the whole 65 km diameter area.

The L1C brightness temperatures were transformed from the antenna reference frame into the Earth reference frame. This included a geometric transformation and a correction for Faraday rotation. Observations for which the resulting radiometric accuracy was > 5 K were deleted. No correction was done for atmospheric effects, that is, differences between the signal at the top of the atmosphere and at the land surface; however, atmospheric effects are very small at L-band (Drusch et al., 2001) and can generally be neglected. Final filtering was then done using the RFI and border flags available in the L1C data products. Observations affected by RFI and observations for which the pixel was close to the border of the extended alias free zone were deleted. Similar to the Scenarios 1 to 3, in this real-world application only ascending node SMOS data were used for assimilation because of small discrepancies between ascending and descending orbit (Dente et al., 2012, Oliva et al., 2013) and the possibility of large thermal gradients between top soil and deeper soil layers causing inaccuracies in the soil moisture retrieval.

To exploit the full potential of the multimeasurement data set, the data assimilation scheme was modified. For each particle, initializing L-MEB with a specific “bias,” hs, and τ combination, all inc combinations were driven. To calculate the weight for each particle, the difference between multiangular SMOS observations and multiangular L-MEB simulations was averaged. The likelihood function given in Eq. [13] was modified accordingly: 
graphic
with ninc number of incidence angles available per observation for the specific site. Here it is ranging from 1 to 3. For the parameter estimation, the period with frozen soil conditions at the beginning of 2010 up to DoY 60 was neglected. This was necessary because frozen conditions were not well captured by the model. This caused large deviations between simulated and observed Tb, which made a parameter estimation impossible for this period. A first analysis showed that soil moisture in 2 cm depth performs better than 4 cm depth to capture the temporal variability of Tb, which is in line with the findings of Escorihuela et al. (2010).

Figure 10 shows the Tb evolution for the Nunn site for the year 2010 for ascending node. Simulated Tb were calculated by forcing L-MEB with measured soil moisture in 2 cm depth and measured soil temperature in 2 and 40 cm depth. Moreover, they were generated based on the temporally constant parameters hs = 0.1 and τ = 0.2 for comparison. Here the general increase in Tb in spring and the decrease in Tb in fall in correspondence to the temperature changes were altered by several Tb reducing rain events. In general, the SMOS observations, presented here for incidence angles between 42° and 43° showed a similar behavior as the simulations, but the variability was relatively larger. In general, observed Tb were about 25 to 30 K warmer than simulated Tb. This was changed rapidly during rain events when observed Tb were close to simulated Tb. Additionally, SMOS measurements showed a noise of ∼ 10 K during the dry period from DoY 200 to DoY 280. The result of the data assimilation procedure for a 42.5° incidence angle is presented in Fig. 10 as well. It was able to follow the SMOS observations time series, even with a relaxation of the aforementioned Tb drops, because of the fact that not only the presented (Fig. 10) SMOS observations between 42° and 43° incidence angles are assimilated but also a larger spectrum. The combination of several incidence angle observations resulted in a robust simulation. Moreover, the observed noise during the dry period from DoY 200 to DoY 280 was decreased.

The estimated parameters for Nunn for the year 2010 are presented in Fig. 11. The average estimated “bias” was ranging between 5 and 15% (average of 9.48% for the whole year 2010), confirming the assumptions of Scenario 3 of a positive “bias.” The annual mean systematic difference in absolute values was 26 K for horizontal and 30 K for vertical polarization. The soil surface roughness remained relatively constant throughout the year at ∼ 0.3 (average of 0.32 for the whole year 2010). However, similar to Scenarios 1 to 3, the uncertainty was relatively high. After the cold period, the vegetation opacity evolution began with very small values of < 0.1. Then from DoY 150 it slightly increased. At DoY 180, it dropped downto ∼ 0.05, which corresponded to a rain event. This gave a subtle hint to the observation of interception water on the plant canopy. Soil Moisture and Ocean Salinity Mission observed the reduced emission, and it was interpreted as wetter conditions than the soil moisture measurement indicated. This was not considered in the model and implicated problems during the parameter estimation procedure. This could be solved by further constraints. After that, vegetation opacity continuously increased to 0.3 at DoY 260 and later decreased to < 0.1 from DoY 295 on. Compared to the findings of Saleh et al. (2006, 2007), values of τ ranging from 0.1 to 0.3 were very reasonable for low vegetation areas like agriculture. This vegetation opacity evolution gave an indication for the natural vegetation conditions at the Nunn site.

To value the impact of the proposed data assimilation approach for radiative transfer parameter retrieval, a comparison to the SMOS Level 2 soil moisture product (version 5.01) was performed. The Level 2 2010 data over the Nunn site (DGG = 191282) were filtered for the correct landcover and retrieval model and snow, interception, and “external events” (e.g., rainfall). The data flags did give a warning that the fit found for the Level 2 retrievals was of rather poor quality. On the other hand, the RFI probability was indicated to be very low for this site.

The comparison is shown in Fig. 12. On the left, the SMOS Level 2 soil moisture is compared to the top soil moisture of the SCAN site Nunn. A systematic difference is clearly visible, and the coefficient of determination (R2) is 0.39. This is in line with current studies, for example, Al Bitar et al. (2012). As the SCAN data have been used to estimate the radiative transfer parameters, it is obvious that the retrieved soil moisture based on the data assimilation algorithm (Fig. 12, right) fits very well to the reference without any systematic differences. Just a few outliers are visible, referring to the beginning of the assimilation period, where the parameters are not well-estimated. This resulted in a coefficient of determination (R2) of 0.95. Note that the SMOS L1C and Level 2 data sets have been processed with different versions and are not directly linked. Even by using the yearly averaged estimated systematic difference of 9.48% and the averaged estimated soil surface roughness of 0.32 in correspondence with the simultaneously estimated variable vegetation opacity, the coefficient of determination (R2) was not reduced significantly. This indicates that once “bias” and soil surface roughness are estimated, improved estimation of vegetation opacity or further parameters may be possible.

Discussion

The results show that the developed data assimilation system based on L-MEB and a sequential importance resampling particle filter was able to simultaneously estimate the parameters hs and τ, even if τ changes over time. It is illustrated that the uncertainty for hs was relatively larger than for τ. This is because of the fact that the sensitivity of the model to hs was smaller than that of τ. Moreover, in addition to hs and τ, the data assimilation framework was able to estimate a systematic difference in the synthetic Tb observations. In Scenario 3, it was demonstrated that a systematic model observation difference could be well estimated with a very small uncertainty, even during simultaneous estimation of hs and τ. As it was expected, the uncertainties in the parameter estimation in Scenario 4, the real-world application, were larger. Nevertheless, it was still possible to estimate plausible parameter values for hs and τ, as well as for a “bias” between point-scale brightness temperature simulations and regional scale brightness temperature observations. Note that the systematic difference originated from different sources, for example, from fixed and not estimated parameters in the retrieval model and from discrepancies between the assumptions made for the point-scale simulations and the given conditions at the regional scale.

In general, the proposed approach performed well for the selected experimental conditions. However, the simultaneous parameter estimation bore risks for compensatory effects between different parameters. This included inadequate parameter estimation as well as an inadequate description of its uncertainty statistics. The rapid changes in vegetation opacity and the corresponding reactions of soil surface roughness simulated in Scenario 2 revealed this appraisement. It could be a problem that during a transition period of 10 d, L-MEB parameters were affected by such errors. Nonetheless, in real-world applications, a sudden change in these parameters over an area representing a passive microwave radiometer footprint may not be likely.

The real-world application gives valuable information for the selected site. For example, for Nunn, as compared to the simulation, Tb simulations and SMOS Tb observations showed a “bias” of 30 K for horizontal and 26 K for vertical polarization. This was in correspondence to the SMOS Level 2 soil moisture product, where the model parameters were different to those estimated in this study and where the “bias” was not considered. Soil surface roughness can be treated as a constant of 0.32, which is in line with findings of, for example, Peischl et al. (2012). Vegetation opacity increases from < 0.1 to 0.3 during late summer and decreases to 0.1 during fall.

A possibility to enhance the quality of the results is to enhance the filtering of SMOS observations to reduce their noise, for example, using data with a higher radiometric accuracy only. Another chance is to constrain the evolution of the parameters. Especially, a constrained variability dynamic of the systematic difference may result in a decrease of its compensatory effects as reported in Scenario 3 and enhanced vegetation opacity retrieval. If the systematic difference is constant over time, once estimated for a site, an estimation might not be necessary in future studies.

Conclusions and Outlook

We investigated parameter retrieval for the L-MEB model by means of four scenarios. These scenarios were defined with increasing complexity, whereby Scenarios 1 to 3 were synthetic simulations and Scenario 4 was an application to real measured data for the SCAN site Nunn in Colorado.

With Scenario 1 we analyzed the ability of the data assimilation framework to estimate temporally constant vegetation opacity and soil surface roughness. Scenario 2 tested whether temporally variable parameters could be estimated. The results showed that the estimation of constant as well as time-varying parameters by means of a particle filter-based data assimilation framework containing L-MEB is possible. A third scenario assumed a hypothetical systematic difference of +5% between Tb simulations and synthetic Tb observations. The analysis showed a good performance to simultaneously estimate vegetation opacity, soil surface roughness, as well as this “bias.” Further enhancements may be obtained by constraining the “bias” estimation. Moreover, further developed assimilation techniques like a particle smoother (Chin et al., 2007), taking into account a time window instead of a single time set, may improve the estimation of radiative transfer parameters under fast changing environmental conditions.

In Scenario 4, real SMOS Tb observations in ascending node were assimilated to L-MEB and vegetation opacity, soil surface roughness, and the systematic difference are estimated. For the SCAN site Nunn, the vegetation opacity showed a seasonal variability between < 0.1 and 0.3, the soil surface roughness remained constant at ∼ 0.3, and the systematic difference between point-scale Tb simulations and regional-scale Tb observations by SMOS was estimated with 30 K for horizontal and 26 K for vertical polarization. The differences are usually a combination of several sources of uncertainty: (i) erroneous calibration of the instrument or adverse environmental effects and aging, (ii) the radiative transfer model, (iii) erroneous input parameters and state variables, and (iv) differences in spatial support of footprint-scale observations and point-scale models. Moreover, the estimation of a soil moisture-related roughness parameter as published in Escorihuela et al. (2007) and Saleh et al. (2007) may improve the results. An issue of future research will be the analysis of different data assimilation and resampling approaches, which will give further information about the quality of the resulting parameter estimation. Moreover, it is recommended to develop and apply enhanced methods for the assimilation of multiple-incidence angle SMOS Tb observations.

As SMOS is able to retrieve vegetation opacity directly, a comparison of this product to vegetation opacity retrieved by the presented approach might give further information about their uncertainties. A first analysis indicates a rather good agreement of the mean values for the study site. Moreover, correlations of the estimated vegetation opacity to vegetation indices of multispectral remote sensing data (such as MODIS) will be identified.

The extension of these experiments to a larger spatial scale will be another issue of future research. To reach this goal, a network of “checkpoints,” where all necessary in situ data are available, may help to estimate global parameter characteristics. Moreover, at these sites the estimated parameters may validate the generic products, which are used for the SMOS Level 2 processing procedure.

These “checkpoints” can be obtained by The International Soil Moisture Network (Dorigo et al., 2011), the U.S. SCAN Network (Schaefer et al., 2007), the German TERENO (Zacharias et al., 2011), or further networks (Zhao et al., 2013). Moreover, modeling results can be used by the SMOS L2 processor as a first guess for soil moisture, which is available, for example, from ECMWF. The methodology itself does not require in situ data although it was demonstrated for its applicability on the point scale. By this approach, local parameter as well as systematic difference characteristics can be taken into account for the processing of an enhanced SMOS soil moisture product.

This study was supported by the German Ministry of Economics and Technology through the German Aerospace Center (50EE1040), by the European Space Agency (Support to Science Element (STSE) Program: SMASPARES), and by the German Research Foundation DFG (Transregional Collaborative Research Centre 32—Patterns in Soil–Vegetation–Atmosphere Systems: Monitoring, modeling and data assimilation). The authors would like to thank the Natural Resources Conservation Service (NRCS) for making available the SCAN data.

All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher.
Open access