The impact of model and rainfall forcing errors on characterizing soil moisture uncertainty in land surface modeling

. The contribution of rainfall forcing errors relative to model (structural and parameter) uncertainty in the prediction of soil moisture is investigated by integrating the NASA Catchment Land Surface Model (CLSM), forced with hydro-meteorological data, in the Oklahoma region. Rainfall-forcing uncertainty is introduced using a stochastic error model that generates ensemble rainfall ﬁelds from satellite rainfall products. The ensemble satellite rain ﬁelds are propagated through CLSM to produce soil moisture ensembles. Errors in CLSM are modeled with two different approaches: either by perturbing model parameters (repre-senting model parameter uncertainty) or by adding randomly generated noise (representing model structure and parameter uncertainty) to the model prognostic variables. Our ﬁndings highlight that the method currently used in the NASA GEOS-5 Land Data Assimilation System to perturb CLSM variables poorly describes the uncertainty in the predicted soil moisture, even when combined with rainfall model perturbations. On the other hand, by adding model parameter perturbations to rainfall forcing perturbations, a better characterization of uncertainty in soil moisture simulations is observed. Speciﬁcally, an analysis of the rank histograms shows that the most consistent ensemble of soil moisture is obtained by combining rainfall and model parameter perturbations. When rainfall forcing and model prognostic perturbations are added, the rank histogram shows a U-shape at the domain average scale, which corresponds to a lack of variability in the forecast ensemble. The more accurate estimation of the soil moisture prediction uncertainty obtained by combining rainfall and parameter perturbations is encouraging for the application of this approach in ensemble data assimilation systems.

Abstract. The contribution of rainfall forcing errors relative to model (structural and parameter) uncertainty in the prediction of soil moisture is investigated by integrating the NASA Catchment Land Surface Model (CLSM), forced with hydro-meteorological data, in the Oklahoma region. Rainfall-forcing uncertainty is introduced using a stochastic error model that generates ensemble rainfall fields from satellite rainfall products. The ensemble satellite rain fields are propagated through CLSM to produce soil moisture ensembles. Errors in CLSM are modeled with two different approaches: either by perturbing model parameters (representing model parameter uncertainty) or by adding randomly generated noise (representing model structure and parameter uncertainty) to the model prognostic variables. Our findings highlight that the method currently used in the NASA GEOS-5 Land Data Assimilation System to perturb CLSM variables poorly describes the uncertainty in the predicted soil moisture, even when combined with rainfall model perturbations. On the other hand, by adding model parameter perturbations to rainfall forcing perturbations, a better characterization of uncertainty in soil moisture simulations is observed. Specifically, an analysis of the rank histograms shows that the most consistent ensemble of soil moisture is obtained by combining rainfall and model parameter perturbations. When rainfall forcing and model prognostic perturbations are added, the rank histogram shows a U-shape at the domain average scale, which corresponds to a lack of variability in the forecast ensemble. The more accurate estimation of the soil moisture prediction uncertainty obtained by combining rainfall and parameter perturbations is encouraging for the application of this approach in ensemble data assimilation systems.

Introduction
Soil moisture is a key variable of the land surface water budget. It has an impact on water, energy and biogeochemical cycles; thus, it plays a major role in many research fields, such as hydrology, agriculture and ecology. As the availability of in-situ measurements is scarce, global soil moisture data rely on satellite retrievals (e.g., Njoku et al., 2003;Kerr et al., 2010;Entekhabi et al., 2010a) and land surface model (LSM) simulations (e.g., Wood et al., 1992;Koster et al., 2000). LSM soil moisture predictions can be enhanced by assimilating near-surface satellite soil moisture observations through a land data assimilation system (LDAS). Often, such systems utilize ensemble-based techniques to update LSM soil moisture predictions in response to satellite observations of near-surface soil moisture (e.g., Reichle and Koster, 2005). Proper characterization of the uncertainty in the LSM soil moisture predictions is crucial for the optimal assimilation of the observed near-surface soil moisture data. Specifically, the quality of soil moisture assimilation estimates highly depends on the accuracy of observational and model error estimates (Crow and Van Loon, 2006;Reichle et al., 2008). However, the way model errors are handled in standard land data assimilation systems could be improved.
A common technique to introduce model uncertainty in land data assimilation systems is by directly adding randomly generated noise to the model prognostic variables, representing errors in model structure and parameters (in addition to perturbing model forcings; Reichle et al., 2007). The present study investigates the ability of this approach to exhaustively describe uncertainty in the soil moisture fields predicted by the LSM, and proposes an alternative method to 3500 V. Maggioni et al.: Model and rainfall forcing errors on characterizing soil moisture uncertainty introduce model uncertainty, i.e., by applying perturbations directly to the model parameters. This alternative approach defines deviations from an optimal parameter set that provides comparable model performance in terms of simulated surface and root zone soil moisture.
Rainfall is the dominant meteorological forcing input to the land surface model for soil moisture simulations. Therefore, uncertainty in the input rainfall products can have an important impact on the predicted soil moisture fields and their associated uncertainty. Hossain and Anagnostou (2006a) have proposed a stochastic error model (named SREM2D) to generate ground truth rainfall ensembles from satellite rainfall products. Hossain and Anagnostou (2006c), and recently Maggioni et al. (2011), have investigated the implication of using SREM2D in representing temporal and spatial uncertainty of soil moisture prediction in a land surface model forced with satellite rainfall data. They showed that soil moisture ensembles from land surface models forced with SREM2D-generated rainfall adequately capture the soil moisture error characteristics at different spatial scales.
In addition to satellite rainfall-only uncertainty, Hossain and Anagnostou (2005) have explored the combined effect of model parameter and rainfall forcing uncertainty on the simulation of soil moisture. They showed that rainfall uncertainty alone can explain only part (between 20 % and 60 %) of the uncertainty in soil moisture prediction. This demonstrates the need for further investigation of the interaction between rainfall and model errors to optimize the use of satellite rainfall products in land data assimilation systems.
This study aims at investigating the impact of model and rainfall forcing uncertainty on soil moisture fields simulated by the NASA Catchment Land Surface Model (hereinafter CLSM, or Catchment model; Koster et al., 2000;Ducharne et al., 2000). Specifically, it builds upon the recent study by Maggioni et al. (2011) that investigated soil moisture prediction uncertainty associated with errors in rainfall forcing alone. We refer to model uncertainty to describe model structural errors, model calibration errors and parametric errors (i.e., uncertainties related to the model itself). On the other hand, we define prediction uncertainty as the combination of model uncertainty and input forcing data uncertainty. The present study introduces several novelties that address some of the limitations of the earlier work by Hossain and Anagnostou (2005). First, two different model uncertainty frameworks are compared: the direct perturbation of model parameters versus the direct perturbation of model prognostic variables, which is typically used in LDAS. Second, results are presented in terms of both surface and root zone soil moisture. Third, our study uses spatially distributed data and a spatial error model, while the study by Hossain and Anagnostou (2005) was limited to a one-dimensional simulation that used single point data, thus neglecting spatial error characteristics. Finally, we use a different land surface model (i.e., CLSM), which is part of the quasi-operational NASA Goddard Earth Observing System Model, Version 5 (GEOS-5) system (Rienecker et al., 2008) developed at the NASA Global Modeling and Assimilation Office (GMAO).
This work is presented from the perspective of data assimilation and investigates whether the method that is currently used in the NASA GEOS-5 LDAS (Reichle et al., 2009) to perturb model variables of the Catchment model is able to adequately describe the uncertainty in the predicted soil moisture. The manuscript is structured as follows. Section 2 provides a description of the study area and datasets. Section 3 describes the approach followed to study how uncertainty in simulated soil moisture is affected by (a) model uncertainty alone, (b) rainfall forcing uncertainty alone, and (c) the combination of the two sources (prediction uncertainty). Section 4 provides a discussion of results, and Sect. 5 summarizes the major findings.

Study area and data
The study area is the Oklahoma region in the United States. Specifically, we use a 25 km Cartesian modeling grid ranging between 100 • W and 94.5 • W in Longitude and 34.5 • N and 37 • N in Latitude (Fig. 1). The study period includes three continuous years from 1 January 2004 to 31 December 2006. The Oklahoma region offers a good coverage by the Weather Surveillance Radar 88 Doppler (WSR-88D) network (Maddox et al., 2002), multi-year satellite rainfall products and a dense network of hydro-meteorological stations from the Oklahoma Mesonet (Brock et al., 1995). The Oklahoma Mesonet provides observations with high temporal frequency from 115 automated observing stations that record several meteorological parameters (rainfall, wind, radiation, etc.) and soil moisture at depths of 5, 25, 60, and 75 cm (available every 30 min). For the 3-yr study period, soil moisture observations of sufficient quantity and quality at all four measurement depths were available only at 21 Mesonet stations as shown in Fig. 1. We refer the reader to Maggioni et al. (2011) for further details about the study domain.
Two rainfall products -the rain gauge-calibrated WSR-88D radar rainfall and the NOAA CMORPH satellite rainfall -are employed in this study and interpolated to the 25km Cartesian grid shown in Fig. 1. Along with supplemental surface meteorological forcing data, these radar and satellite precipitation products force the land surface model at the 25km grid resolution to generate soil moisture fields. The radar rainfall product is from the Stage IV WSR-88D precipitation estimation algorithm, which consists of a national mosaic of precipitation estimates based on observations from all WSR-88D radars across the continental US (Fulton et al., 1998). Stage IV data represent the best quality WSR-88D rainfall product available at hourly/4-km resolution and include corrections for ground clutter and anomalous propagation, vertical reflectivity profile effects and systematic variations of the reflectivity-rainfall relationship on the basis of bias estimates through comparisons with rain gauge observations  ( Fulton et al., 1998;Lin et al., 2005). The satellite product is the NOAA-climate prediction center morphing (CMORPH) product, which is based on a unique combination of passive microwave (PMW) retrievals and infrared (IR) data (Joyce et al., 2004). In particular, CMORPH uses motion vectors from half-hourly interval IR images to propagate precipitation estimates derived from PMW data. The spatio-temporal resolution of CMORPH is 8 km/half-hour. As stated above both radar and satellite precipitation datasets are gridded to the 25-km grid and aggregated to a 3-h time step to ensure common spatial and temporal scales.

The Catchment land surface model
The NASA CLSM Ducharne et al., 2000) constitutes the modeling scheme to simulate surface and root zone soil moisture in this study. The model uses the hydrological catchment as a fundamental land surface element with boundaries defined by topography, abandoning the traditional approach of quasi-rectangular model units with boundaries defined by the overlying atmospheric grid. Consequently, catchments lying below the same atmospheric grid cell can differ in terms of the topographical parameters and vegetation type, extracted from high-resolution vegetation data sets. Within each catchment, a redistribution of soil water content according to topography is combined with more traditional vegetation and evaporation parameterizations from the Mosaic LSM Suarez, 1992, 1996), a three-layer snow model (Lynch-Stieglitz, 1994), and a linear soil-heat diffusion scheme. The Catchment model framework is based on the well-established TOPMODEL, developed by Beven and Kirkby (1979). Forced with surface meteorological data, CLSM predicts the soil moisture distribution from the catchment morphology and from three bulk soil moisture prognostic variables, representing equilibrium conditions associated with water table distribution and non-equilibrium conditions near the surface. The first of these three prognostic variables is the catchment deficit, which measures the average amount of water, per unit area, needed to bring all of the soil in the basin to saturation, assuming that the unsaturated zone is initially at equilibrium. The second prognostic variable is the root zone excess, defined as the amount by which the moisture in the top 100 cm of the soil exceeds (or is less than) the moisture estimated by the local equilibrium profile. The third prognostic variable, the surface excess, measures how far the water in the 2 cm ground surface layer is from the equilibrium with the root zone soil moisture.
The horizontal distribution of soil moisture in response to topographical information allows the separation of the catchment into three different hydrological regimes that respond dynamically to changes in the three soil moisture prognostic variables: the (area) fraction over which the ground surface is saturated, the fraction over which the ground surface is not saturated but transpiration proceeds, and the fraction over which no transpiration proceeds as the soil moisture is below the wilting point. Within each of the three regimes, different physical processes control runoff and evaporation. At each time step the model computes surface water and energy fluxes for each fraction, and then combines them into a single catchment flux.
CLSM is forced with several surface meteorological variables, such as precipitation, humidity, air temperature and radiation and with prescribed, climatological vegetation parameters (leaf area index and greenness from satellite observations). Precipitation data in this study are from the aforementioned WSR-88D (Stage IV) and CMORPH products. The remaining surface meteorological forcing data (air  temperature and humidity, radiation, wind, and surface pressure) are extracted from the Global Land Data Assimilation Systems (GLDAS) project (Rodell et al., 2003; http://ldas. gsfc.nasa.gov), based on output from the global atmospheric data assimilation system at the NASA GMAO (Bloom et al., 2005). CLSM simulations were initialized from a spin-up integration by forcing the model with the WSR-88D (Stage IV) rainfall fields and by looping three times through the 3-yr time series of forcing data (2004)(2005)(2006).
The Catchment model was demonstrated in several past studies to realistically describe soil moisture dynamics (e.g., Bowling et al., 2003;Boone et al., 2004). In a recent study, Maggioni et al. (2011) showed consistency and high correlations between soil moisture anomaly time series from the OK Mesonet station observations and corresponding simulations from the Catchment model forced with WSR-88D (Stage IV) rainfall.
In this study CLSM is first forced with the WSR-88D (Stage IV) precipitation data to generate the reference soil moisture fields (Fig. 2). Then, CLSM simulations are performed to investigate the contributions of model and rainfall-forcing error in soil moisture prediction uncertainty: (i) model uncertainty alone through parameter perturbations (case M1); (ii) model uncertainty alone through prognostic perturbations (case M2); (iii) rainfall forcing uncertainty alone (case F); (iv) combination of rainfall forcing and model (through parameter perturbation) uncertainty (case M1F); and (v) combination of rainfall forcing and model (perturbing prognostics) uncertainty (case M2F). The uncertainty analysis is carried out in terms of both surface and root zone soil moisture values. Surface soil moisture henceforth refers to the (0-2) cm soil moisture output from CLSM and to the OK Mesonet soil moisture observations at 5 cm depth, whereas "root zone" soil moisture is defined as the (0-100) cm soil moisture CLSM output and the corresponding depth-weighted average over the 5 cm, 25 cm, 60 cm and 75 cm OK Mesonet measurements.
Details about the different experiments are described next. Specifically, Sect. 3.2 describes the methodology used to study the model uncertainty alone (left column of Fig. 2), Sect. 3.3 presents the setup to study rainfall forcing uncertainty alone (middle column of Fig. 2), and Sect. 3.4 illustrates the method used to analyze the combined rainfall and model uncertainty (right column of Fig. 2).

Model uncertainty
The first framework to characterize the model error (case M1 in Fig. 2) is based on model parameter perturbations. Firstly, an analysis is conducted to identify a small subset of  parameters for which CLSM exhibits sensitivity in terms of simulated soil moisture effects. This was computed with respect to the reference soil moisture, obtained by forcing the model with radar rainfall and with the originally calibrated set of parameters, currently used as part of the GEOS-5 LDAS. Each parameter was independently scaled by multiplicative coefficients ranging from 25 % to 400 % of their standard values. Only one parameter at a time was perturbed with the same scaling factor for all grid cells, while keeping the other parameters constant; the scaled parameter was held constant for the entire 3-yr integration.
Surface and root zone soil moisture fields generated from CLSM based on each value of the parameter-scaling factor are evaluated in terms of two performance metrics: efficiency score (ES) and mean relative error (RelativeBias), defined as follows: whereθ is soil moisture obtained from the perturbed simulation, ϑ is the reference soil moisture and N is the length of the time series. Because the number of parameters that we choose affects the number of ensemble members that we will need to use in the combined uncertainty simulation (Sect. 3.4), the subset of perturbed model parameters should be kept as small as possible. Therefore, although the analysis was performed on several model parameters, we only represent model uncertainty with the two parameters to which modeled soil moisture showed significant sensitivity: the Clapp-Hornberger parameter (b) and the soil wilting point wetness (wpwet). Figure 3 shows the above performance metrics (ES and Relative Bias) versus the parameter scaling factor (presented in %) for these two parameters. The figure also shows the metrics for the matric potential at saturation (ψ), which is reported as an example of a parameter to which the model did not show considerable sensitivity. Metrics are computed for each grid cell and shown as boxplots, where the central mark is the median, and the edges of the box are the 25th and 75th percentiles. Results are presented here only for the surface soil moisture, but are very similar for the root zone simulations. Efficiency scores range between 0.88 and 1 for b and from negative values to 1 for wpwet, while the relative bias varies between −0.3 and 0.1 for b, and between −0.3 and 0.6 for wpwet. It is worth noting that in this study we are Combinations of these two parameters (b and wpwet) are then used to produce model simulations of soil moisture fields. Each set of perturbed parameter values is assigned a likelihood value; in this study it was chosen to be the efficiency score (Eq. 1), which indicates the correspondence between the model predictions and the reference integration. For example, the parameter combination where b is perturbed with a multiplicative factor of 80 % of its original calibrated value and wpwet of 90 % of its original value is assigned an efficiency score of 0.94 (0.88) for surface (root zone) soil moisture (Fig. 4). The total sample of simulations is then split into behavioral and non-behavioral parameter combinations, according to the equifinality assumption (Beven and Binley, 1992). In other words, we want to identify the model parameter values that give similar error metrics in terms of simulated soil moisture. In fact, we assume that the existing parameter set for the Catchment model (defined in previous studies) is the most appropriate, and then use the equifinality concept to determine the range of parameter variations around that parameter set that give comparable model performance (relative to model simulations using the optimal parameter set) according to a pre-defined criterion.
This criterion is based on a breaking point in the efficiency score metric, sorted from high to low values (Fig. 4). The slope of the tangent to the curve drops from −0.012 to −0.031 (i.e., by about 3 times) when the efficiency score drops below 0.8 for surface soil moisture (units of the tangent slope are change in ES per experiment after ranking the experiments by ES). A very similar behavior is observed for root zone soil moisture, where the tangent slope drops from a value of −0.020 to a value of −0.056 when the efficiency score falls below 0.7. The combinations of perturbed parameter values that show efficiency scores higher than the ES values corresponding with the two breaking points are selected as behavioral parameter sets. This results in thirteen sets of parameter values, listed in Table 1, which are then used for case M1. (Experiment "b200w115" exceeds the ES threshold only for surface soil moisture and is therefore excluded). For example, the tag "b080w105" on the x-axis corresponds to the combination that applies deviations of 80 % and 105 % to b and w, respectively. The dashed line shows the cutoff threshold of acceptable efficiency scores.
For each simulation, the same parameter set is used for all grid cells in the domain and kept constant during the 3-yr integration period.
Alternatively, model uncertainty can be represented by directly perturbing model prognostic variables (case M2 in Fig. 2). The three bulk soil moisture prognostic variables of CLSM are surface excess, root zone excess, and catchment deficit (Sect. 3.1). However, Reichle et al. (2007) observed that perturbations in the root zone excess might lead to biases between the ensemble mean and the unperturbed reference integration. Therefore, perturbations are limited to the surface excess and the catchment deficit. Each ensemble member is subject to randomly generated noise that is added to the model prognostic variables to represent errors in model structure and parameters. Specifically, normally distributed additive perturbations are applied; ensemble means are constrained to zero, and time series correlations are imposed via a first-order autoregressive model. The perturbation parameter values are identical to the ones listed in Table 2 of Liu et al. (2011), and shown here in Table 2 for completeness. For consistency with case M1, thirteen simulations are run by perturbing the two model prognostic variables. It is worth noting that while case M1 only addresses uncertainty in the model parameters, approach M2 technically addresses both model structural and parameter uncertainty. In addition, M2 requires independent observations to properly characterize the parameter values of the statistical perturbation model, as it makes assumptions about the model error characteristics (i.e., standard deviation, space-time correlation, Gaussian error), which cannot be known a priori. Moreover, the character of the ensembles generated by M1 and M2 is very different. In the first case, each ensemble member will be biased against the reference in terms of long-term mean soil moisture at a given grid cell, whereas in case M2 each ensemble member is generally unbiased with respect to the ensemble mean (or reference).
In Sect. 4, soil moisture estimates obtained from model simulations using the two methods will be compared to standard normal deviates of Mesonet soil moisture observations to study how the model uncertainty captures the ground measurements.

Rainfall forcing uncertainty
The SREM2D rainfall error model (Hossain and Anagnostou, 2006a) is employed here to determine the error propagation from rainfall forcing to the soil moisture prediction. SREM2D models the multi-dimensional error structure of satellite retrievals with space-time stochastic formulations. The important aspect of SREM2D is its ability to model not only the spatial variability of rain rate estimation error, but also the spatial structure of the successful delineation of rainy and non-rainy areas. As real sensor data actually exhibit spatial clusters of false rain and no-rain area delineations, this characteristic makes SREM2D capable of capturing the magnitude of the satellite rainfall error and variability across scales (Hossain and Anagnostou, 2006b). The input parameters for the satellite rainfall error model are shown in Table 3. In the study of Maggioni et al. (2011) it was shown that, compared to a simpler rainfall error model, SREM2D ensembles provide better encapsulation of the reference (WSR-88D) precipitation.
In this study, CMORPH satellite rainfall is perturbed by SREM2D to produce an ensemble of twenty-four equiprobable reference-like rainfall realizations. This ensemble is used to force CLSM to obtain an ensemble of soil moisture fields that represent the rainfall forcing uncertainty (case F in Fig. 2). In the case of rainfall uncertainty alone, CLSM parameters are set to their original values.

Combined uncertainty
The combined uncertainty is studied by merging the rainfall forcing uncertainty with the model uncertainty, estimated with the two different approaches described above, which results into two experiments: M1F, where model parameter perturbations are added to rainfall forcing perturbations, and M2F, where model prognostic perturbations are added to rainfall forcing perturbations (right column of Fig. 2). In both experiments, CMORPH precipitation forcing has been perturbed through the SREM2D error model. Details about the two experiments are provided next.
Firstly, four parameter sets out of the thirteen behavioral parameter combinations were chosen to represent the model parameter uncertainty in the combined uncertainty experiment. These are the experiments on the anti-diagonal of Table 1 (0.50, 1.15), (0.80, 1.05), (1.25, 0.90), and (2.00, 0.80) for b and wpwet parameters, respectively. The chosen experiments provide a possible combination of perturbations that assures a balanced representation of the range of values of both parameters in the behavioral parameter sets. To some extent, this also addresses the issue of potentially correlated parameters: by selecting the experiments along the anti-diagonal in Table 1, we give a slight preference to anticorrelated over correlated model parameters. However, future studies should evaluate the correlation among model parameters and eliminate potential parameter dependences. For each of those four sets of parameters, six ensemble members were integrated with perturbed CMORPH precipitation inputs (using the SREM2D rainfall error model), for a total of 24 soil moisture ensemble members (case M1F in Fig. 2).
Secondly, another 24-member ensemble of CLSM was integrated by perturbing both precipitation (using SREM2D) and prognostic variables at the same time. Specifically, normally distributed additive perturbations were added to the same prognostic variables as in case M2, i.e., surface excess and catchment deficit. The output from these runs is a soil moisture ensemble carrying rainfall forcing and model uncertainty (case M2F).
The combined uncertainty (rainfall forcing + model error) obtained with the two different methods (cases M1F and M2F) is then compared to the uncertainty introduced by only perturbing rainfall forcing (case F). To investigate the uncertainty associated with each experiment, rank histograms, exceedance ratios and uncertainty ratios (defined below) are presented for both surface and root zone soil moisture. Results are described in the next section.

Discussion of results
We first investigate how the model uncertainty determined by the two techniques encapsulates in-situ measurements from Mesonet. Next, we compare the forcing rainfall uncertainty alone to the combined model and rainfall forcing uncertainty (or prediction uncertainty) through soil moisture time series plots and performance metrics.

Model uncertainty analysis
To compare the two model uncertainty approaches (M1 and M2), we focus on anomaly time series, specifically standard normal deviate time series and associated anomaly correlation coefficients, which capture the correspondence in phase between model estimates and ground observations, disregarding potential bias or differences in the variance (Entekhabi et al., 2010b). Standard-normal deviates are computed at the daily scale as differences between the actual values and the monthly climatological average values of the 3-yr time series, normalized by the corresponding standard deviation. Anomaly correlation coefficients are chosen as a performance metric due to our focus on data assimilation applications, which requires unbiased errors. error approaches. Specifically, Figs. 5 and 6 show standardnormal deviate time series at two representative Mesonet stations, located, respectively, in the eastern (wetter) and western (drier) region of the study area and at the corresponding 25 km grid cells indicated in Fig. 1. Figures 5 and 6 show that variations in the model predicted soil moisture values are consistent with the Mesonet measurements as well as the associated rainfall forcing variations. Furthermore, the standard-normal deviate time series are consistent between surface and root zone soil moisture. The domain-average anomaly correlation coefficients between simulated (ensemble mean) soil moisture and the Mesonet observations are 0.78 (0.64) for the surface soil moisture (root zone soil moisture) for both M1 and M2 model uncertainty approaches. At the western station (drier conditions), correlation is slightly higher than the eastern pixel case (wetter conditions). Specifically, at the eastern location, the anomaly correlation coefficient is 0.53 (0.56) when model uncertainty is introduced by directly perturbing parameters, and 0.54 (0.56) when model uncertainty is assessed by perturbing model prognostic variables for surface (root zone) soil moisture. On the other hand, at the western location the anomaly correlation coefficient is 0.58 for the surface soil moisture, and it increases to 0.81 for root zone soil moisture, for both model uncertainty approaches. We note that correlation coefficients are comparable between M1 and M2 for surface and root zone soil moisture and for both rainfall climatological conditions. From the anomaly time series and correlation coefficients analysis, we can conclude that both model uncertainty approaches capture equally well the measurement variability from Mesonet stations at both soil moisture depths. These considerations give us confidence about the viability of the two approaches to introduce model uncertainty in this study.

Model and rainfall forcing uncertainty analysis
In this section, the rainfall forcing uncertainty alone is compared to the combined model and rainfall uncertainty through time series plots and performance metrics. Figures 7 and 8 show time series of surface and root zone soil moisture ensembles during the warm season of 2005 (June to September) for the following three experiments: case F where only rainfall forcing uncertainty is introduced; case M1F which combines rainfall forcing uncertainty with model uncertainty by perturbing model parameters; and case M2F which includes rainfall forcing uncertainty and model uncertainty by perturbing prognostics. Figures 7 and 8 also include time series of reference model soil moisture, which is defined as the output from CLSM forced with unperturbed WSR-88D rainfall fields. In particular, these figures show time series for the same representative grid cells as in Figs. 5 and 6, one located in the wetter eastern region of the study area (Fig. 7) and one in the drier western part (Fig. 8). Figures 7  and 8 show that for both soil moisture depths, the ensemble envelope is wider when considering the combined rainfall  forcing and model uncertainty through parameter perturbations (case M1F) compared to the other two cases. Moreover, soil moisture time series show slightly wider ensemble envelopes for the combined rainfall forcing and prognostic perturbations compared to the rainfall forcing alone uncertainty experiment. By comparing the two pixel time series, it can be inferred that more variability is observed in wetter rainfall conditions for both surface and root zone soil moisture, which is expected as soil water content variability correlates to rainfall variability and consequently exhibits stronger dependence on rainfall forcing error.
As discussed earlier, wider ensemble envelopes derived from the combined rainfall and model uncertainty will increase the probability of encapsulating the reference soil moisture simulations. Two metrics are introduced to quantify this effect. Specifically, we use the exceedance ratio (ER) metric to evaluate the ability of the ensemble integrations to encapsulate the reference predictions and the uncertainty ratio (UR) metric to evaluate the accuracy of the ensemble envelope width. These metrics are presented for the three experiments with respect to the reference modeled soil moisture. The exceedance ratio is defined as where N exceedance is the number of times and locations for which the reference soil moisture falls outside the ensemble envelope and N t is the total number of times and locations. If ER is equal to 1, it means there is a 100 % probability that the reference will not fall between the lower an upper bounds of the ensemble envelope, whereas if ER is close to 0, there is a low probability that the reference exceeds those bounds, or, in other words, there is a perfect encapsulation of the reference inside the ensemble envelope. The uncertainty ratio, on the other hand, represents the ratio of the aggregate width in the simulated uncertainty relative to the corresponding average actual uncertainty (with respect to the reference soil moisture): whereθ i upper andθ i lower are, respectively, the upper and lower bounds of the simulated ensemble,θ i ens mean corresponds to the ensemble mean, andθ i opt represents the reference soil moisture obtained by forcing the model with unperturbed WSR-88D rainfall fields. A perfect ensemble spread has UR Hydrol. Earth Syst. Sci., 16, 3499-3515 values equal to 1. If UR is less (greater) than 1 the ensemble is underestimating (overestimating) the model prediction error spread, which can translate into insufficient (excessive) weights given to observations in an ensemble-based data assimilation system. Exceedance and uncertainty ratios have been demonstrated to be viable metrics to evaluate ensemble prediction performance in several studies (Hossain et al., 2004;Hossain and Anagnostou, 2005;Borga et al., 2006;Moradkhani et al., 2006 among others). In particular, the combination of these two statistics allows to assess two contrasting issues: if the uncertainty limits are too narrow (that is, ER is high), the comparison with the reference fields suggests that prediction uncertainty is underestimated; on the other hand, if the limits are too wide (that is, UR is high), the model may overestimate the reference uncertainty, which translates into a poor predictive capability.
ER and UR statistics are computed for each grid cell of the domain. Figure 9 shows frequency histograms of ER and UR metrics for the surface and root zone soil moisture. The corresponding mean values of these metrics are reported in Table 4. Consistent with the time series discussed above, the lowest ER is observed for the experiment that combines rainfall forcing and model parameter uncertainties (M1F). Histograms for M1F are shifted towards significantly lower ER values where the average for surface (root zone) soil moisture is 0.39 (0.35). The complement to 1 of the exceedance ratio can also be interpreted as the  On the other hand, the UR in case M1F exhibits values that are closer to 1 compared to the other experiments. In particular, the average UR is equal to 0.58 (0.45) for surface (root zone) soil moisture when perturbations are based only on the forcing rainfall, which increases to 0.65 (0.53) and 1.02 (1.08) when adding the model prognostic perturbations and model parameter perturbations, respectively. This demonstrates that in cases F and M2F, the ensemble significantly underestimates the average actual soil moisture error, whereas in case M1F, the ensemble spread is less biased and statistically closer to the average actual error. Furthermore, we note that root zone soil moisture exhibits slightly higher UR values in the case of M1F and lower UR values for cases F and M2F. This is related to the fact that variability in the forcing rainfall more directly influences the upper centimeter soil moisture, whereas variability added to the model parameters affects the prediction of both surface and root zone soil moisture variables.
In summary, the ER and UR metrics indicate that by combining model parametric uncertainty with rainfall forcing uncertainty, variability in soil moisture prediction error can be well described. Specifically, we showed that this combination increases the ability of the ensemble envelope of encapsulating the reference simulation (lower ER), without overestimating the ensemble spread (UR close to 1). On the other hand, using only rainfall perturbations, or adding prognostic perturbations to rainfall forcing perturbations, results in an underestimation of the actual soil moisture errors (UR values significantly lower than 1) and a lower likelihood of encapsulating the reference simulations (higher ER values).
Even though past studies have shown that CLSM could be run with a small number of ensemble members (Reichle et al., 2002;Forman et al., 2012), we verify here the viability of using twenty-four ensemble members in this particular case. Specifically, we performed an analysis of the two statistics described above for different ensemble sizes that range from 4 to 24 (Fig. 10). We observe that both ER and UR converge at about 16 ensemble members for all the three experiments and for both surface and root zone soil moisture. This gives us confidence that the 24 ensemble members used in this study are sufficient to obtain statistically meaningful results.
Another useful and powerful tool to evaluate ensemble predictions is represented by the rank histogram. If both model and perturbations are able to reproduce the forecast distribution (i.e., the ensemble is statistically consistent), the ensemble members should be equally likely scenarios for the reference simulation. Specifically, the number of times that the reference falls within any two adjacent ensemble members should be independent of the position of these members in the ordered ensemble (Siegert et al., 2012). Therefore, if the ensemble is consistent, a histogram of the rate at which the reference falls into each interval (i.e., the rank histogram) should be flat, up to random variations due to the finiteness of the number of samples (Talagrand et al., 1997;Hamill and Colucci, 1997;Hamill, 2001) the predicted ensemble would show up in a U-shaped rank histogram, whereas consistent biases would translate into a sloped rank histogram. Figure 11 shows rank histograms for surface soil moisture, evaluated for the domain average and for two representative grid cells, one located in the eastern (and wetter) region of the study area and the other in the western (and drier) region. When perturbations are added only to rainfall forcing (case F), the rank histograms show a U-shape when the domain average is evaluated, whereas at the single grid cells a strong bias is observed. In the case of the eastern (wetter) pixel, the ensemble appears to underestimate the reference simulations, as most of the times all the ensemble members are below the reference (i.e., rank is equal to 24). On the other hand, at the western (drier) pixel, the prediction shows a consistent overestimation, as most of the ensemble members are usually above the reference simulation, which corresponds to low values of the rank.
When perturbations are added to rainfall forcing and prognostics (case M2F), no remarkable difference is observed, even though at the grid cell scale the bias slightly decreases compared to case F. However, when perturbations are applied to forcing and model parameters (case M1F), at the domain average, the rank histograms show a considerably flatter histogram, which demonstrate a better estimation of the forecast uncertainty. At the single pixel scale, the bias is significantly reduced, again showing flatter rank histograms when compared to the other two experiments. Very similar results are obtained for root zone soil moisture (not shown here).
The analysis of the rank histograms confirms what was shown by the exceedance and uncertainty ratios. In cases F and M2F, the uncertainty in the forecast ensemble is either underestimated (at the domain scale) or shows consistent biases (at the single pixel locations). Specifically, a positive bias is observed at the drier location, whereas a negative bias is observed at the wetter location. For case M1F, the combination of rainfall forcing and model parameter perturbations reproduces more accurate forecast uncertainty and generates a more consistent ensemble of surface and root zone soil moisture at the domain and grid-cell scales.

Conclusions
In this study, we investigated ways to characterize soil moisture prediction uncertainty needed in land data assimilation systems. Specifically, we focused on the impact of rainfall forcing error and model uncertainty (separately and in combination) on the prediction of surface and root zone soil moisture using the Catchment land surface model, which is part of the NASA GEOS-5 LDAS. Firstly, we examined how ensemble prediction uncertainty encapsulates the reference ground   scale (a, b, c), at the grid cell shown in Fig. 5, located in the eastern region of the study area (d, e, f), and at the grid cell shown in Fig. 6, located in the western region of the study area (g, h, i). measurements from Mesonet stations, comparing two model uncertainty approaches. The first approach uses perturbed model parameters (M1). The second approach, which is also the method that is currently adopted in the GEOS-5 LDAS, adds randomly generated noise to the model prognostic variables during the land model integration (M2). Both methods were shown to capture well the soil moisture (temporal) variations measured at in-situ Mesonet stations. However, the first technique (M1) contributed more spread in the soil moisture ensemble time series compared to the prognostic random perturbation method.
Next, we compared soil moisture prediction uncertainty derived from rainfall forcing uncertainty alone (case F) against uncertainty derived from combining model with rainfall forcing uncertainty, using both model uncertainty approaches (cases M1F and M2F). Rainfall forcing perturbations alone provided narrower ensemble envelopes of simulated surface and root zone soil moisture time series compared to the ones obtained by combining forcing and model uncertainties. Prognostic perturbations added only modest variability to the rainfall forcing uncertainty. On the other hand, the combined uncertainty employing perturbations of model parameters (case M1F) was shown to produce the widest soil moisture ensemble envelopes relative to the other two cases and for both soil moisture depths. This was exhibited in the exceedance and uncertainty ratio metrics. The lowest exceedance ratio, a metric that assesses the capability of the ensemble integrations to encapsulate the reference, was observed when combining rainfall forcing and model parameter uncertainties, while the experiment that combined rainfall-forcing uncertainty with prognostic perturbations exhibited only slightly reduced ER values relative to the rainfall forcing uncertainty alone. Moreover, the uncertainty ratio, which measures the accuracy of the ensemble envelope width, was shown to be close to 1 when perturbations were added to the model parameters (case M1F). This demonstrates that uncertainty in soil moisture was better described by M1F, with ensemble spread close to the model prediction uncertainty spread, defined as the difference between the ensemble mean and the reference model simulation. When rainfall forcing perturbations alone (case F) or rainfall forcing and prognostics perturbations (case M2F) were considered, the uncertainty ratios were significantly lower than 1, exhibiting an underestimation of the actual errors.
The analysis of the rank histograms corroborated what the exceedance and uncertainty ratios showed, that is, the most consistent ensemble of soil moisture was obtained by combining rainfall and model parameter perturbations (case M1F). When only the rainfall forcing was perturbed (case F), the rank histogram showed a U-shape at the domain average scale, which corresponds to a lack of variability in the forecast ensemble. At the grid-cell scale, a consistent positive (negative) bias was observed at the drier (wetter) location. No considerable difference was observed when perturbations are added to both rainfall forcing and model prognostics (case M2F) other than a slight decrease of the bias at the pixel scale. On the other hand, method M1F showed significantly flatter rank histograms and results in an ensemble that is better able to characterize the actual prediction uncertainty.
The accurate estimation of the soil moisture prediction uncertainty is encouraging for the application of this approach in ensemble data assimilation systems. In particular, this study contributes to the development of the NASA GEOS-5 LDAS, providing valuable insights about the interaction between rainfall forcing and model uncertainties in case of satellite rainfall application in land data assimilation. Specifically, it shows that the current method to represent model and forcing uncertainties poorly describes the actual prediction uncertainty.
Other methods to investigate forecast ensemble uncertainty could be considered in future studies. In particular, model averaging methods such as BMA (Bayesian Model Averaging) are computationally affordable and able to provide sharp predictive uncertainty intervals (Raftery et al., 2005). Further studies should also investigate potential correlations among the model parameters, as the perturbation of one soil parameter may affect the others, and how this would influence the model performance. An extension of this study will focus on investigating the impact of the herein proposed combined rainfall and model uncertainty on the assimilation of satellite-retrieved soil moisture in a land data assimilation system. Moreover, other land surface models should be considered. Finally, other regions of the world characterized by different hydroclimatic regimes need to be analyzed.