Soil Moisture Data Assimilation in a Hydrological Model: A Case Study in Belgium Using Large-Scale Satellite Data

In the present study, we focus on the assimilation of satellite observations for Surface Soil Moisture (SSM) in a hydrological model. The satellite data are produced in the framework of the EUMETSAT project H-SAF and are based on measurements with the Advanced radar Scatterometer (ASCAT), embarked on the Meteorological Operational satellites (MetOp). The product generated with these measurements has a horizontal resolution of 25 km and represents the upper few centimeters of soil. Our approach is based on the Ensemble Kalman Filter technique (EnKF), where observation and model uncertainties are taken into account, implemented in a conceptual hydrological model. The analysis is carried out in the Demer catchment of the Scheldt River Basin in Belgium, for the period from June 2013–May 2016. In this context, two methodological advances are being proposed. First, the generation of stochastic terms, necessary for the EnKF, of bounded variables like SSM is addressed with the aid of specially-designed probability distributions, so that the bounds are never exceeded. Second, bias due to the assimilation procedure itself is removed using a post-processing technique. Subsequently, the impact of SSM assimilation on the simulated streamflow is estimated using a series of statistical measures based on the ensemble average. The differences from the control simulation are then assessed using a two-dimensional bootstrap sampling on the ensemble generated by the assimilation procedure. Our analysis shows that data assimilation combined with bias correction can improve the streamflow estimations or, at a minimum, produce results statistically indistinguishable from the control run of the hydrological model.


Introduction
Sources of errors in model estimations are generally of two kinds, input data errors (including initialization) and model errors.One method to decrease deviations of the estimated model states from observations is to optimally "assimilate" observed data of relevant model variables into the model.Data assimilation has a long history as a modeling tool, and in the present article, we will focus on applications to hydrology, with emphasis on soil moisture and streamflow.
Soil moisture is a variable that plays an important role in a hydrological system by regulating several water exchanges and energy fluxes.There are many monitoring networks for soil moisture around the world [1,2].The data from in situ SSM measurements are invaluable in calibrating land surface and hydrological models, while long-term time series can reveal trends due to climate or land use changes.However, in situ soil moisture measurements still remain very sparse, especially compared to the high spatial variability of this variable.
Satellite observations can ensure today complete coverage of large areas and provide a representation of the soil state at regular time intervals.Specifically, as an alternative to ground-based measurements, great efforts have been made in the past two decades to develop soil moisture products from microwave signals with the launch of several instruments on board satellite missions [3][4][5][6][7][8].Examples of such satellite sensors are the Special Sensor Microwave Imager (SSM/I) on board a series of Defense Meteorological Satellite Program platforms, the Tropical Rainfall Measuring Mission (TRMM) Microwave Imager (TMI) by NASA, the Advanced Microwave Scanning Radiometer-Earth Observing System (AMSR-E) by NASA, its successor AMSR2, the Soil Moisture and Ocean Salinity (SMOS) by ESA, the Soil Moisture Active Passive (SMAP) by NASA and the ASCAT sensor by EUMETSAT.In particular, in [8], a Land Parameter Retrieval Model (LPRM) is used for surface soil moisture retrievals in order to generate long global series of data.A special algorithm for hydrological applications (HydroAlgo) has been developed in [9] and used to generate maps of snow depth and soil moisture content from AMSR-E and AMSR2 data.
Each satellite mission carrying such instruments has been designed to satisfy certain technical requirements, e.g., regarding spatial coverage, resolution and scanning frequency.Therefore, significant differences can be expected between the technical specifications of satellite data and the corresponding requirements of several applications in Earth systems, e.g., modeling of the atmosphere and of the hydrological cycle.In fact, many studies have investigated the potential to improve river discharge estimations by using satellite data in hydrological models.This is a complex task, and the outcome depends on the model used, the spatiotemporal scales, the method of using the data (direct use or assimilation) and the quality and availability of the satellite data.For example, [10] obtain improvements in flood forecasting by complementing discharge data assimilation with soil moisture data assimilation.
Several other studies exist investigating potential benefits of assimilating soil moisture data, e.g., [11][12][13][14][15], with encouraging or mixed results.In particular, in [13], the authors conclude that the underestimation of the vertical soil coupling in the SWAT model impedes the ability of data assimilation to update deep soil moisture.On the other hand, in [16], soil moisture products have been successfully used in a simplified continuous rainfall-runoff model in order to provide initial wetness conditions prior to a flood event, demonstrating that these products contain useful information, which can be exploited in flood forecasting.Such studies highlight the issues that may arise due to model structure, catchment characteristics and observation errors.
Another frequently-encountered problem is the difference, or bias, in the mean value and variability between different sources of SSM data (satellite retrievals and model integrations in our context).This bias can pose significant obstacles to the use of the information contained in satellite measurements for any kind of application, including data assimilation, and has to be removed.Different methods to remove such bias have been explored in the literature, in particular linear regression, mean and variance matching and Cumulative Distribution Function (CDF) matching [14,[17][18][19][20][21].In [21], the authors show that proper pre-processing of soil moisture observations can be critical for the performance of the data assimilation.
Estimating an appropriate error for soil moisture observations, used in data assimilation, is still an open problem.One choice is to compare the satellite data with ground-based observations (validation) [22].Independent estimations of the error are also possible, e.g., based on triple collocation [23] and Fourier transforms [24].The main difficulties are the higher uncertainty and coarser spatial resolution of the satellite observations, compared to ground-based measurements, and the errors associated with the retrieval algorithms.For an extensive discussion of such issues arising in SSM data assimilation into a hydrological model, we refer the reader to [25] and the references therein.
The ensemble Kalman filter is a data assimilation method used in geosciences [26,27] and in hydrology in particular [28].It is a generalization of the common Kalman filter [29], appropriate for non-linear systems without the need to calculate the adjoint model.However, it was eventually acknowledged that the use of the stochastic terms required by this method may lead to biased estimations if the model includes non-linear processes [28,30].This problem has been considered in other studies too, and a solution has been proposed in [11,31].This solution can be used in cases where the assimilated variable is unconstrained, which is not the case for soil moisture.This is why the approach of [31] required repeated application of the method in order to further reduce the residual bias due to the bounds of soil moisture.A different problem, but still related to bias, is tackled in [32], where the EnKF is used in order to find an optimal method to reduce model biases.
The existence of bounds for the assimilated variable gives rise to another problem as well.In the ensemble Kalman filter, the generation and use (as perturbations) of stochastic terms from a normal distribution are generally required.Obviously, repeated application of such perturbations on a bounded variable, like soil moisture, will often drive out of their bounds the values that are already near them.This is another source of biases in the estimations, especially in the most extreme cases (very dry or very wet conditions, speaking of soil moisture).A technique based on variable standard deviation was introduced in [33] to solve this problem.Although this method reduces the bias due to exceedance of the bounds, it does not completely eliminate it.
In the present study, our goal is to assimilate satellite data for surface soil moisture in the SCHEME hydrological model and to assess the effects of this assimilation in terms of streamflow.The satellite data for soil moisture that we use here are derived from active microwave measurements with the ASCAT sensor on board the MetOp [7].The derivation of soil moisture values from such measurements relies on a change detection method [34].ASCAT data have been already validated using in situ measurements in southwestern France [35] yielding encouraging results; see also [36,37].More extensive validation of such data makes up part of the ongoing EUMETSAT project H-SAF.We focus here on the case of a catchment in Belgium while our methodology aims to address the two problems in EnKF pointed out previously.
In our approach, surface soil moisture is properly treated as a bounded variable.In particular, the perturbations needed for the ensemble Kalman filter are generated from truncated normal distributions.This ensures that no soil moisture value exceeds the bounds after the application of perturbations.We use truncated distributions in order to generate perturbations for the precipitation, as well, which is thought of as the forcing of the hydrological model.In this case, the perturbations are applied multiplicatively; therefore, truncated log-normal distributions are used.Of course, precipitation is not bounded in theory from the top, and simple log-normal sampling would suffice.The reason for this unusual, at first sight, choice is that it will ensure that the precipitation values that will be generated by applying such perturbations will remain within the climatological margins.
We tackle in the following way the issue of bias, induced by non-linear processes in conjunction with the use of stochastic terms.Initially, a full data assimilation run is performed.Then, another run is performed in which only perturbations are applied to precipitation, that is the model states are not updated by the Kalman equation.The streamflow ensemble generated by this run is carrying only the bias of the process because there is no assimilated information.The deviation of its mean from the control simulation for the streamflow is then used in order to remove the bias from the full data assimilation ensemble.The introduction and use of truncated probability distributions in a data assimilation scheme and the post-processing technique explained above are the two main advances proposed here.
The article is organized as follows.In Section 2, we present a description of the hydrological model and the technical details of the assimilation methodology, including the handling of bounded variables and of bias under the presence of non-linear processes.Subsequently, in Section 3, we deal with the study area and the data used in the numerical experiments.The results are the subject of Section 4. We conclude by providing more insight in Section 5.

Model Description
The SCHEME hydrological model (SCHEldt-MEuse, from the names of the two major rivers of Belgium), used for our simulations, is the distributed version of the IRMB hydrological model [38].The IRMB model has been successfully applied to various catchments ranging from about 100-1600 km 2 , representing several hydrologic conditions in Belgium (e.g., [39]).It has been used to study the climate change impacts on the water cycle under given scenarios in Belgium and Switzerland [40][41][42].
The SCHEME model is optimized for river basins up to 20,000 km 2 .The hydrological processes are lumped within grid cells of 49 km 2 , allowing one to describe the heterogeneity of hydrologic conditions and of hydrometeorological input data.With this design, the SCHEME model is able to simulate a variety of basins and hydrological conditions in the river Scheldt and Meuse River Basins in Belgium and upstream in France.It is currently being used in medium-range streamflow forecasts [43,44].
The SCHEME model structure comprises 9 different land covers with a snow accumulation and melting module.The land cover types are represented by the appropriate variable fractions on each grid cell.Sub-grid soil properties are also indirectly taken into account through a regionalization procedure.The actual evapotranspiration is calculated on the basis of the water intercepted by the vegetation and the water content of two soil layers, as well as the Potential Evapotranspiration (PET) according to the Penman formula.The time step of the model is equal to one day.Surface water is simulated with a unit hydrograph, and the underground water is represented with two reservoirs.The streamflow produced on each grid cell is routed to the outlet with a 1D submodel based on the width function of the river network [45].A sketch of the SCHEME model mechanisms is presented in Figure 1.
Let us recall briefly from [38] the soil water balance mechanism in the SCHEME model.The general principle is the same for the two soil layers.Let W j−1 and W j be the water contents of a given soil layer at the beginning and the end of day j, respectively, W 0 the water content when the soil is saturated and R j the difference between recharge and potential evapotranspiration in day j for the same soil layer.For the upper layer, the recharge is the absorption prior to all surface flooding, and for the lower layer, the recharge is the quantity of the water that infiltrates during the flooding phase.If the capacity of the lower layer is exceeded, then the excess water is percolated to one underground reservoir of the model.If R j < 0, then the water content of the soil will be decreased according to the following equation: If instead R j > 0, then the water content will be increased as: Overall, in the SCHEME model, the soil water content and runoff calculation is based on simple water balance considerations.In other studies, while analogous principles are used for soil moisture, the runoff is estimated using the Soil Conservation Service Curve Number [12][13][14].The PDM model [11,46] is more like the SCHEME model in the sense that runoff production is controlled by the storage capacity of the soil in the context of a conceptual approach.However, an important difference is that in the PDM model, the storage capacity is treated as a random variate to account for spatial variability.A different approach is adopted in the TOPMODEL, where the spatial patterns of soil moisture are regulated by shallow groundwater gradients [15].It is assumed that these gradients can be estimated from local topography.Runoff generation in the TOPMODEL is based on saturation and infiltration excess mechanisms.The parameters of the SCHEME model have been calibrated based on data from a variety of catchments in the Scheldt and the Meuse River Basins from the period 1981-1988.The calibration technique combines elements from the approach in [38] and an automatic algorithm, the "Shuffled Complex Evolution of the University of Arizona" (SCE-UA; [47,48]), which has been proven efficient in locating globally-optimal parameters of hydrological models.The objective function to optimize in this case is the daily root mean square error.After calibration, the parameters are regionalized over the corresponding basins with the use of Artificial Neural Networks (ANN, NevProp software, Version 3, [49]) and geographic indices.These indices are derived from the CORINE Land-Cover (http://www.eea.europa.eu/publications/COR0-landcover)and the Soil Map of the European Countries at 1:1,000,000 [50].First, empirical relationships are being sought between parameter values (the predictands) and indices (the predictors) over the set of catchments used in the calibration step.For each parameter, several indices are combined, an ANN model is trained and the correlation between the predicted values by the ANN model and the predictands is calculated.The correlation is corrected for overfitting with a bootstrap sampling [51] among the catchments in the following way.A model based on the same predictors is trained on the bootstrap sample.The correlation between the values predicted with this model and the predictands is calculated on the bootstrap sample.Then, the same model is applied on the full sample (including all catchments), and a new correlation is obtained.This process is repeated, and the average of the difference between the two correlation values provides an estimation of the bias to remove from the correlation obtained with the full model on the full sample.Different models (combinations of indices) are tested, and the selection is based on the unbiased correlation.Then, the model is applied over the whole basin, including areas not covered with the catchments used for calibration.The sample of parameter values provided by the bootstrap procedure is also used to estimate the model error term in Equation ( 5), as will be explained in Section 3.
The capacity of the upper soil layer is one of the parameters of the model that are optimized and regionalized.Its value is in the range of 5-28 mm in the Scheldt River Basin and corresponds to a soil depth of several centimeters.The saturation of this conceptual reservoir is put into correspondence with the soil wetness index of the satellite product.

General Context
Our data assimilation approach is based on the ensemble Kalman filter technique [26].The EnKF is derived from the standard Kalman filter [29] by generating ensembles to account for model and forcing uncertainties in the presence of nonlinear processes.The goal is to optimally merge model predictions with observations, by providing estimations of the error covariance matrix based on the ensemble of the model states generated with this technique.This procedure effectively eliminates the need to linearize the model used in the simulations.
The general principles of the EnKF can be described as follows.We restrict our interest in finite dimensional systems described by ordinary differential equations.If S is the (differentiable) phase space of such a system and y : R → S a trajectory of the system in it, subject to the dynamics defined by the vector field F, then we have: In discrete form, Equation (3) reads: where k denotes the time step, y k the system state at time step k and M k the function (set of algebraic operations) that must be applied on y k in order to calculate y k+1 , arising from the discretization of the initial dynamical system.In our context, M k arises from a model, and we can write: where the term w k represents the errors and the uncertainties in the model.Let z represent the measurements of the physical quantity related to the model variable y.If the true model state is denoted by m, then we introduce the measurement operator H by the following equation: where ε represents the measurement error.In practice, it is often assumed that ε is normally distributed.We extend the discretization notation of Equation ( 4) to the other variables as well.Thus, we will write During assimilation, the model states are adjusted using information from the observations of the variable z.In order to avoid confusion, we denote the variable y, as calculated by the model without any update, by y − (a-priori or forecast) and the same variable updated through assimilation by y + (a posteriori or analysis).Let finally y i denote the ensemble members, i ∈ {1, . . . ,n}.For example, we would write (y i k ) − for the i-th a priori ensemble member at time step k, and therefore, (y i k+1 Under these conventions, the updated state (y i ) + at time step k is given by the equation: where K k is the so-called Kalman gain matrix at time step k, given by: In Equation ( 8), P − is the error covariance matrix for the forecast estimate y − obtained from the ensemble [26] and R is the variance of the noise term ε.As customarily, the exponent T denotes the transpose of the underlying matrix.
In our case, there is only one variable to assimilate, SSM from satellite observations.Therefore, the assimilation problem is one-dimensional, and all of the matrices involved are simply real numbers.We will consider as true model state m the satellite observation z corrected for bias using as a reference the SSM state calculated by the SCHEME model.In this setting, the measurement operator is simply the ratio H = z/m, with m = z − b(z), b : [0, 100] → R being the bias function when we express the SSM as a percentage.Details will be provided in subsequent sections of this article.

Bounded Variables
One particular feature of soil moisture is that it is bounded between the values 0 and 100.Therefore, during assimilation in an EnKF context, applying Gaussian noise on SSM will often violate this constraint, especially for values near the bounds.This is a long-standing problem that was, to our knowledge, not adequately addressed so far.A method using the variable standard deviation for the error distribution is proposed in [33].This solves the problem partially because it reduces the frequency of the bound crossing, but it does not completely eliminate it.Furthermore, the choice of very small standard deviations for variable values near the bounds may not be always the most appropriate.Our solution is based on the concept of the truncated probability distribution [52].This is a conditional distribution resulting from the restriction of the domain of a given probability distribution.Formally, let X be a real random variable with probability density function f and cumulative distribution function F.Then, the restriction of X on the interval (a, b] is described by the distribution r( f ) (a,b] defined by: where: We notice here that r( f ) (a,b] , or simply r( f ) when the interval is well understood, is essentially a rescaling of f over the given interval so that its integral on it (the probability) is equal to 1. Based on Equations ( 9) and (10), one can derive all of the properties of a truncated distribution necessary for the calculations.For the needs of the present article, we implement sampling algorithms for two distributions of this kind, truncated normal and truncated log-normal.
A truncated variable will be used in our context in the following way.Let us consider an SSM value, for example y 0 .We want to add to y 0 a noise term , coming from a normal truncated distribution r(N(0, σ 2 )) with mean value 0 and standard deviation σ.Since the sum y 0 + has to be between 0 and 100, we have −y 0 ≤ ≤ 100 − y 0 .This defines the range in which r(N(0, σ 2 )) has to be restricted.
We need such additive truncated variables in order to generate the variance R of the noise term ε (Equations ( 6) and ( 8)) and to model uncertainty regarding SSM calculations by the SCHEME model (term w k in Equation ( 5)).For the case of the variance R, we use a truncated normal distribution with σ = 8.0%.This value is in the range obtained from the validation of SSM satellite data (ERS and ASCAT) in Northern France (Grand Morin), using in situ soil moisture measurements [53].
We work similarly with truncated log-normal distributions, which are used to generate multiplicative factors.These factors constitute an implementation of the uncertainty in interpolated precipitation, which is one of the input fields in the SCHEME model.In existing literature, this uncertainty is represented as a rule by full range log-normal distributions [12].However, in our case, allowing the full range log-normal terms can generate completely unrealistic amounts at certain rainfall events, very far away from the limits in which the SCHEME model has been calibrated and known to perform well.Using truncated log-normal distributions instead, conditioned by the climatology of the Demer catchment, solves this problem.

Non-Linear Processes
The ensemble Kalman filter can be used, by its design, in non-linear systems without the need for developing a linearization scheme, like for example in the case of the extended Kalman filter.This has the potential to considerably simplify the study.On the other hand, it became at some point clear that the stochastic terms generated by the filter generally induce bias in the estimations when non-linear processes take place [30].
It is possible to implement a correction method for this bias in order to make evident the net effect of the assimilation procedure.Such a method has been proposed in [31].The case of the SCHEME model considered here is however more special, and we deal with biases in a different way.
In the SCHEME model, the calculation of moisture at the top layer of the soil is carried out as follows.Let ssm be the physical surface soil moisture and W the model variable that represents this moisture.The value of W at each time step depends on throughfall and residual potential evapotranspiration.If it is less than or equal to the capacity W 0 of the upper soil layer (Figure 1), then it provides directly the physical soil moisture for the given time step: ssm = W.However, the SCHEME model allows W to exceed W 0 , according to precipitation, which determines throughfall and evaporation conditions.When this happens, W can be written as: W = W 0 + S, where the term S is used by the model to calculate contributions to the streamflow and to the lower soil layer.In this case, the physical soil moisture of the upper soil layer coincides with the capacity: ssm = W 0 .
This simple condition on the range of W is at the origin of the non-linear behavior of the SCHEME model in the calculation of SSM.In order to quantify the effect of the stochastic terms used by the ensemble Kalman filter, we generate two different ensembles for the streamflow.Let {(s a ) i } be the streamflow ensemble generated by running as usual the full assimilation procedure and {(s p ) i }, i ∈ {1, . . ., n}, be the streamflow ensemble generated by running the model while keeping active only the precipitation perturbations.This means in particular that the state updates through the Kalman equation are switched off.Therefore, the members of s p do not carry information from satellite data.Let also s 0 be the control simulation of the model.
is the average of the ensemble s p at each time step k, then the daily value of bias is given by b k = (s p ) k − (s 0 ) k and the corrected streamflow ensemble by: The notation convention is the same for all of the ensembles s a , s p and s b,a .Another option here is to express the bias as a ratio in order to avoid possible negative values (e.g., in models with large biases).However, the results presented here are obtained with Equation (11), since in our case, this correction yields well-defined (nonzero positive) streamflow values.Trying the ratio version of the bias produces slightly different statistical scores in Tables 1-4 and exactly the same statistical significance results in Tables 5-8 (Section 4).
The reason to examine the effect of precipitation perturbations only is that among the three transformations that affect directly the soil moisture variable W, namely the precipitation and soil moisture perturbations and the Kalman equation, only the first one has the potential to drive W out of range and therefore to activate the decomposition W = W 0 + S. The other two preserve the bounds: the soil moisture perturbations by construction (truncated normal distribution sampling); and the Kalman equation because it is a convex linear combination of two variables in the same interval.In other words, the precipitation perturbations are at the root of the non-linear effects in this case.
The bias correction method proposed in this section presents some analogies with the method of [31] as applied in [11].However, the two are fundamentally different.In our approach, the forcing perturbations are applied while the data assimilation is disabled.This generates a new ensemble carrying only the effects of the non-linear processes.The new ensemble is then used in order to remove the bias due to non-linearity.In [11], the data assimilation is applied, but without the forcing perturbations.This results in a unique time series (unperturbed model prediction) used as the baseline to remove the bias.According to our knowledge, the bias correction method proposed here is new and has not been used in previous studies.

Study Area and Data
Our study is carried out for the catchment of Demer, located in the central-eastern part of Belgium.This catchment covers an area of 1775 km 2 and makes part of the Scheldt River Basin (see Figure 2).The main land use in the Demer catchment consists of crops (46%), meadows (29%) and forests (18%).The remaining 7% is inhabited.The elevation ranges between 17 and 173 meters, with a mildly hilly profile in the south.The length of the main river is 63 km.For the period 1971-2000, the average runoff is about 13 m 3 /s (0.642 mm/day) with a daily maximum of 59 m 3 /s (2.915 mm/day) and daily minimum of 0.5 m 3 /s (0.025 mm/day).
A rain gauge network from the Royal Meteorological Institute of Belgium (approximately one station per 100 km 2 ) provides daily rainfall data in the Demer basin.For the period 1966-1995, the average yearly precipitation in the catchment is about 768 mm.The rainiest month is June, with 77 mm of rainfall, and the driest one is April with 50 mm.Therefore, the precipitation is quite uniformly distributed throughout the year, typical for an oceanic climate under some continental influence.The spatial distribution of precipitation is also quite homogeneous due to the almost flat topography of the area.
In the same reference period , the mean temperature is 9.8 The study period is 1 June 2013-31 May 2016.The satellite data we use for this period come from the EUMETSAT Satellite Application Facility on Support to Operational Hydrology and Water Management (project H-SAF, http://hsaf.meteoam.it/).However, older data starting in June 2009 have been used in order to calculate bias correction for soil moisture.
For assimilation purposes, we use the large-scale surface soil moisture product SM-OBS-1(or H07).This product is generated from MetOp scatterometer data (ASCAT) at a coarse resolution (25 km) controlled by the IFOV of the instrument and provides data for the surface layer of the soil (0.5-2 cm).The measurement of soil moisture with the MetOp scatterometer is fairly direct because of the high sensitivity of the microwaves to the water content of the upper soil layer.This is carried out using the so-called TU-Wien change detection method [54].However, the different contributions to the observed total backscatter from the soil, vegetation and soil-vegetation interaction effects cannot be separated.The main assumptions in this approach are linearity in the relation between the backscattering coefficient and surface soil moisture content, stability in time of roughness and land cover and the influence of the vegetation phenology on backscatter at the seasonal time scale only.Snow and frozen soil conditions are not treated in the change detection method, and special care should be taken, using auxiliary information, in order to avoid spurious contributions in winter.In the Demer catchment, forests and topography have a rather limited influence (see Figure 2 and the catchment metrics in the beginning of this section).Snow presence is also quite limited.The main concern is freezing conditions.We address such issues in the present study by deactivating data assimilation when the mean air temperature falls below 2.0 • C.
We also use the H-SAF product SM-DAS-2 (or H14) for comparisons with the SCHEME model data at the root zone level.This product is generated in the ECMWF Land Data Assimilation System by assimilating ASCAT SSM data (essentially the H07 product that we discussed previously).In this case, the data assimilation process is based on the extended Kalman filter [55].The ECMWF model generates soil moisture profiles according to the Hydrology Tiled ECMWF Scheme for Surface Exchanges over Land (HTESSEL, [56]) providing estimations for four soil layers for a total thickness of 289 cm below the surface.SM-DAS-2 is available at a spatial resolution of 25 km with a daily time step.
Before assimilating the SSM data from H07, a bias correction is applied on them in order to determine the true state in Equation (6).The H07 data used for this purpose cover the period from June 2009-May 2013, and the reference data are obtained from the output of the SCHEME model.First, the spatial averages over the Demer catchment of all of the necessary data are calculated to produce two four-year-long times series, one from H07 and another one from the model.The bias correction is based on the difference between the distributions represented by the two time series (CDF matching, [19]).In Figure 3, we can see the bias as a function of the soil moisture value.Obviously, a piecewise linear interpolation can be used in order to obtain a bias function defined everywhere.As mentioned in Section 2, the satellite data are treated in EnKF as random variables (Equation ( 6)).This is also useful in the implementation of the data assimilation procedure.In particular, during data assimilation, each pixel of the SCHEME model grid is matched with the closest satellite pixel.This means that about ten SCHEME model pixels are matched with the same satellite pixel.However, the corresponding satellite value is perturbed independently for each model pixel and natural cover inside that pixel, using truncated normal distributions as explained in Section 2.2.2.This approach yields different perturbed satellite values for each time step, pixel and sub-pixel features of the SCHEME model.
The model error term w k in Equation ( 5) is calculated as follows in our case.By varying a set of parameters of the SCHEME model, we generate an ensemble of SSM time series with 100 members, as the output of the model runs.Each parameter set defines a variant of the model, which is used in order to produce a long time series of 30 years (1966-1995) for surface soil moisture.In this way, we generate the ensemble of SSM model output.We work with this dataset by dividing first the full SSM range into five classes: 0-20%, 20-40%, 40-60%, 60-80% and 80-100%.These data are then used in order to determine the distribution of the SSM values (Figure 4) and, for each SSM class, the distribution of the standard deviations (Figure 5).The former is used in the generation of the initial SSM ensemble at the start day of the assimilation, while the latter in the generation of the error w k (Equation ( 5)) by sampling a truncated normal distribution.For the generation of the term w k of a given SSM value, we first draw a value from the uniform probability distribution on the interval [0,1].Using the data represented in Figure 5, we can then determine the corresponding bin.This bin makes part of the subdivision of the standard deviation interval, which is known and corresponds to the class to which the given SSM value belongs.With this information, we can determine the (randomly drawn) standard deviation σ corresponding to the given SSM value.Subsequently we use σ in order to define the truncated normal distribution r(N(0, σ 2 )) and randomly draw from it the term w k , as described in Section 2.2.1.The probabilistic context for the generation of the standard deviation σ and, subsequently, of the term w k described above is new and has not been used in other studies according to our knowledge.
In an analogous way, the perturbations for precipitation (Section 2.2.2) are generated by sampling a truncated log-normal distribution with the mean equal to one and the standard deviation that depends on the precipitation value.The dependence is determined empirically based on the climatology of the Demer catchment and the ensemble verification measures for streamflow, described in [30].In this way, the perturbed precipitation is lying always within the interval [0.0, 60.0] in mm/day.This interval covers all of the cases of daily rainfall in the calibration period of the SCHEME model.
The truncated perturbation approach in this context guarantees that no unrealistic precipitation events will be generated during the execution of the data assimilation procedure.Exaggerated precipitation can become a problem with the use of unbounded probability distributions (see Section 2.2.2), and it has not been previously addressed in the existing literature.The same principle holds for any model sensitive to certain forcing variables.

Simulations and Results
The precipitation input and the results of simulations with the SCHEME model are graphically represented in Figures 6-8.Each figure corresponds to one of the three years of the study period (June 2013-May 2016).The precipitation and the soil moisture in the two layers of the model are reduced to time series by calculating the areal averages of the corresponding fields over the Demer catchment.
The effect of the assimilation is very obvious at the level of the soil moisture.Regarding the upper layer of the soil, the ensemble generally follows the satellite values.The situation is very different in the lower soil layer, where no assimilation is taking place, but the moisture values are calculated by the SCHEME model after assimilating satellite data at the upper layer.The results of the reference simulation are the only basis for quantitative comparisons that we have for this layer of the soil.In some cases, we observe large deviations between the ensemble members and the reference simulation (e.g., autumn and winter in Figure 8).
For the streamflow, we have two cases, one ensemble obtained by applying the assimilation procedure (DA) and another one obtained by correcting the previous ensemble for bias (DA-BC), according to the method described in Section 2.2.3.
In the corresponding figures, the ensemble averages are also depicted together with the control simulation of the model and the streamflow observations.The differences between these time series are quantified using well-established synoptic statistical measures such as bias, Absolute bias (A-bias), root mean square error (RMSE), correlation (R) and Nash-Sutcliffe Efficiency coefficient (NSE).The calculations are carried out at seasonal and annual time scales.The results of such synoptic statistics for all possible combinations of data and time periods are shown in Tables 2-4 Finally, in order to assess the significance of the differences in the scores, we calculate the 95% confidence intervals of the ensemble average scores by using a bootstrap approach.More precisely, each ensemble is randomly sampled in two ways, once in its members and once more in time.We consider the difference between the scores of the control simulation and of the ensemble average as statistically significant if the control value lies outside of the confidence interval.
The qualitative features of the results for the streamflow, seasonal and annual, as well, depend on the period considered.For example, for the year June 2013-May 2014, the average of the DA ensemble yields better annual scores than the reference simulation.The situation is reversed in the year June 2014-May 2015, while in the final period from June 2014-May 2015, the ensemble average becomes again better than the control simulation.However, the differences between the scores are small, and this motivates the use of the statistical significance test mentioned before.
On the other hand, the average of the DA-BC ensemble exhibits more consistent performance in time.For example, its annual scores are all better than the scores of the control simulation in the three years studied here.Like in the case of the DA ensembles, the differences in the scores are small and their significance is checked with the calculation of confidence intervals.
The results of statistical significance analysis using bootstrap replications are presented in Tables 6-8.The convention we adopt here is to mark by "Y" differences that are statistically significant and by "N" differences that are not.The "Y" symbol is complemented by the sign + if the corresponding score of the ensemble average is better than the score of the control simulation and by the sign − if it is worse.We use the symbol "Y0" for the only case where the difference is significant, but does not exceed the order of magnitude 10 −3 in terms of absolute values (Table 8).
We can see for instance that the decrease in the scores for the DA ensemble in June 2014-May 2015 reported previously is indeed significant.Over the same period, the DA-BC ensemble has bias that is consistently lower than the bias of the control simulation, also in a statistically-significant way.

Discussion
In this implementation of the ensemble Kalman filter, we can see in a clear way the manifestation of bias due to the use of stochastic perturbations in a nonlinear process.This bias can have unpredictable effects on model performance and can mask the net result of the adjustment applied to the model states through the Kalman equation.
In the three-year period that we study in this article, there are cases of large differences between the control simulation and the ensemble averages for the streamflow during shorter periods of time.For example, from 20 February 2014-6 May 2014, the average of the DA ensemble is a much better approximation of the observations than the control simulation.This is clearly seen in Table 9 and in Figure 9.Of particular interest is the very large difference between the NSE coefficients.On the other hand, the average of the DA-BC ensemble is closer to the control simulation in this period with smaller errors, but still a negative NSE coefficient.10 and Figure 10).Over this period, the DA ensemble average yields poor results while the scores of the control simulation appear quite normal, compared to typical SCHEME model performance.Finally, the average of the DA-BC ensemble is slightly better than the control simulation in this case.The previous situations demonstrate the need to reduce the biases arising in EnKF where non-linear processes take place.The general idea here is to identify the sources of bias and, subsequently, to filter out their effect by running again the model while keeping active the perturbations, which are responsible for bias generation, but without assimilating any data.
It is clear also, at least for this case study, that soil moisture assimilation has a limited impact on the streamflow calculated by the model.This is not surprising because streamflow is mainly determined by precipitation.The assimilation of soil moisture would have more potential in cases where the model does not estimate correctly the water content of the soil, especially in cases with significant amounts of rainfall and saturated soils.Such indications exist also for the SCHEME model in the period that we study here.For example, we see in Figure 6 that in March 2014, there are large differences between the soil moisture values calculated by the model and detected by the satellite.The net effect of assimilating such satellite data is a slightly better streamflow output from the model (in terms of average of the DA-BC ensemble).This improvement is statistically significant, as is demonstrated in Table 6.
The moisture variations of the lower soil layer is another point that warrants further investigation.There are systematic differences between the control simulation and the output of the EnKF, which are larger in autumn and the beginning of winter.The origin of such deviations is in part the bias generated by the precipitation perturbations through the nonlinearities of the model and in part the assimilation itself, which transmits the satellite signal from the surface to the lower soil layer through the processes encompassed in the model, with the largest contribution coming from the latter.The lack of any soil moisture measurements for the lower layer prevents at this point a thorough analysis.However, H-SAF provides a root zone soil moisture product, SM-DAS-2 (H14).In order to compare the estimations of H14 and those of the SCHEME model, we can only rescale both in the interval [0,100], because the units are different.Therefore, no proper quantification of the differences could be made in soil moisture units.Nevertheless, visual inspection of the graphs that we provide in Figure 11 reveals striking similarities between H14 and the ensembles for the lower soil layer, generated by the SCHEME model after assimilating ASCAT SSM data (Figures 6-8).This is an example where two completely different modeling approaches (the SCHEME model with EnKF for ASCAT data, on the one hand, and the ECMWF data assimilation system based on extended Kalman filter, on the other) lead to similar conclusions.).The SCHEME model output is given originally as the degree of saturation (%), while H14 is in volumetric units (m 3 /m 3 ).A rescaling has been applied in both time series for this comparison.

Conclusions
The purpose of this study is to explore the possibilities to improve streamflow estimation from a hydrological model by assimilating soil moisture data.Our analysis leads to improvements of certain aspects of the ensemble Kalman filter for bounded variables used in nonlinear processes.
The ensemble version of the Kalman filter relies on the repetitive generation of stochastic terms, which are used in order to perturb the ensemble members and calculate state corrections.In the existing literature, this leads systematically to biases, because the perturbations do not in principle preserve the bounds of a variable like the soil moisture.Here, we propose a natural solution to this problem based on truncated probability distributions, which can be applied to the case of any bounded variable.
Another problem is the generation of bias in the model output when such perturbations are applied to a variable that is used in nonlinear processes (soil moisture in our case).The solution proposed here is to identify the source of bias and to isolate its effect by running the model in ensemble mode without assimilating any data, but keeping the perturbations that cause the bias generation.The ensemble obtained in this way can be used to remove the nonlinear processes bias from the full assimilation ensemble.
The combination of the previous steps leads to very consistent data assimilation results.In summary, the streamflow estimations after data assimilation in the SCHEME hydrological model and bias correction are in many cases better than the control simulation of the model in a statistically-significant way.In the remaining situations, these estimations are statistically identical to the control.
Our results point out that the EnKF could be also used as a diagnostic tool to detect weaknesses in a model and to improve its performance.The case of March 2014 reported previously is indicative of such situations.It suggests that further investigation is required regarding the optimization of the SCHEME model parameters that control the seasonal streamflow contribution, especially in spring.
Extending the domain of hydrological simulations into other catchments will provide new insight regarding the spatial features of our data assimilation approach assisted by bias correction.Ultimately, a more extensive use of satellite data could be considered in this context.

Figure 2 .
Figure 2. Left panel: The Meuse and the Scheldt River Basins, including the Demer catchment.Right panel: Topographic details and elevation data in meters in the Demer catchment.The small red square is the catchment outlet.

Figure 3 .
Figure 3. Bias for the surface soil moisture product H07, based on data from June 2009-May 2013.The original raw data are also shown.

Figure 4 .
Figure 4. Probability distribution of surface soil moisture values, based on simulation data with the SCHEME model, from January 1966-December 1995.The bins subdivide the interval formed by the whole dataset into ten equal parts.

Figure 5 .
Figure 5. Probability distribution of surface soil moisture standard deviation values, based on simulation data with the SCHEME model, from January 1966-December 1995.The bins subdivide the standard deviation intervals corresponding to each soil moisture class into ten equal parts.

Figure 11 .
Figure 11.Areal average for the lower soil layer moisture of the SCHEME model and the H-SAF product H14, over three successive years (1 June 2013-31 May 2016).The SCHEME model output is given originally as the degree of saturation (%), while H14 is in volumetric units (m 3 /m 3 ).A rescaling has been applied in both time series for this comparison.
• C in the Demer catchment.The hottest months are July and August, with 17.7 • C and 17.3 • C, and the coldest ones are January and February with 2.8 • C and 2.9 • C respectively.

Table 1 .
. Synoptic statistics for streamflow with and without bias correction in the ensemble, 1 June 2013-31 May 2016.The unit for bias, Absolute bias (A-bias) and RMSE is mm/day.
Figure 6.Hydrological data and model output for the period June 2013-May 2014.From top to bottom: precipitation, soil moisture in the upper layer, soil moisture in the lower layer, streamflow with data assimilation only and streamflow with data assimilation and bias correction.Figure 7. Hydrological data and model output for the period June 2014-May 2015.From top to bottom: precipitation, soil moisture in the upper layer, soil moisture in the lower layer, streamflow with data assimilation only and streamflow with data assimilation and bias correction.Figure 8. Hydrological data and model output for the period June 2015-May 2016.From top to bottom: precipitation, soil moisture in the upper layer, soil moisture in the lower layer, streamflow with data assimilation only and streamflow with data assimilation and bias correction.

Table 2 .
Synoptic statistics for streamflow with and without bias correction in the ensemble, 1 June 2013-31 May 2014.The unit for bias, A-bias and RMSE is mm/day.

Table 3 .
Synoptic statistics for streamflow with and without bias correction in the ensemble, 1 June 2014-31 May 2015.The unit for bias, A-bias and RMSE is mm/day.

Table 4 .
Synoptic statistics for streamflow with and without bias correction in the ensemble, 1 June 2015-31 May 2016.The unit for bias, A-bias and RMSE is mm/day.

Table 5 .
Statistical significance of the differences between the scores of the ensemble averages and the corresponding scores of the reference simulation for streamflow, 1 June 2013-31 May 2016.

Table 6 .
Statistical significance of the differences between the scores of the ensemble averages and the corresponding scores of the reference simulation for streamflow, 1 June 2013-31 May 2014.

Table 7 .
Statistical significance of the differences between the scores of the ensemble averages and the corresponding scores of the reference simulation for streamflow, 1 June 2014-31 May 2015.

Table 8 .
Statistical significance of the differences between the scores of the ensemble averages and the corresponding scores of the reference simulation for streamflow, 1 June 2015-31 May 2016.

Table 9 .
Synoptic statistics for streamflow with and without SSM assimilation, 20 February 2014-6 May 2014.The unit for bias, A-bias and RMSE is mm/day.

Table 10 .
Synoptic statistics for streamflow with and without SSM assimilation, 20 September 2014-10 November 2014.The unit for bias, A-bias and RMSE is mm/day.