Streamflow data assimilation for soil moisture analysis

Abstract. Streamflow depends on the soil moisture of a river catchment and can be measured with relatively high accuracy. The soil moisture in the root zone influences the latent heat flux and, hence, the quantity and spatial distribution of atmospheric water vapour and precipitation. As numerical weather forecast and climate models require a proper soil moisture initialization for their land surface models, we enhanced an Ensemble Kalman Filter to assimilate streamflow time series into the multi-layer land surface model TERRA-ML of the regional weather forecast model COSMO. The impact of streamflow assimilation was studied by an observing system simulation experiment in the Enz River catchment (located at the downwind side of the northern Black Forest in Germany). The results demonstrate a clear improvement of the soil moisture field in the catchment. We illustrate the potential of streamflow data assimilation for weather forecasting and discuss its spatial and temporal requirements for a corresponding, automated river gauging network.


Introduction
Quantitative precipitation forecasting (QPF) one of the most complex challenges in numerical weather prediction (NWP) (e.g.Rotach et al., 2009;Wulfmeyer et al., 2008).QPF failures can be due to errors in numerics, limited spatial resolution of the model, erroneous model physics, incorrect initial conditions, and limited predictability.The skill of QPF, particularly on the mesoscale, is still strongly limited by uncertainties in initial conditions.Particularly, dynamics in complex terrain and the inhomogeneous distribution of water vapour are considered the most important unknowns in the initial fields.The water vapour field of the continental lower troposphere and therefore cloud formation and precipitation is influenced by the interaction of the atmosphere with the land surface through the energy and water fluxes.Corresponding studies show the soil moisture influence on quantity and spatial distribution of precipitation (e.g.Sch är  et al., 1999;Hohenegger et al., 2008;Trier et al., 2004).Particularly, in summertime, continental QPF depends on the initialization of root zone soil moisture and other land surface states (Reichle et al., 2002;Hohenegger et al., 2009).Soil moisture not only depends on the weather but also on the local land surface characteristics (soil texture, vegetation, orography).But for this highly heterogeneous quantity only scarce representative measurements are available at point locations (e.g.B árdossy and Lehmann, 1998; Grayson and Western, 1998).Multiple efforts to apply remote sensing to regions of scarce or shallow vegetation to obtain the skin layer soil moisture are currently under way (e.g.Crow and Wood, 2003;Dunne and Entekhabi, 2006;Drusch and Viterbo, 2007;Gao et al., 2007).But these techniques so far do not provide data for soil moisture estimates under dense vegetation and within the total soil profile.
Hence the knowledge of the soil moisture distribution is a key issue in NWP.
As the lower boundary of weather forecast and climate models, land surface models (LSM) calculate the coupled water and energy balance at each grid cell of the atmospheric model.On these scales (≥1 km 2 ), soil texture, topography and vegetation and, therefore, water and energy fluxes, soil moisture, runoff and soil temperature are highly heterogeneous (e.g.Kabat et al., 1997).This heterogeneity can neither be measured nor modelled explicitly at acceptable cost.For each grid cell the precipitation is balanced by the sum of evapotranspiration, runoff and soil moisture change.Evapotranspiration and soil moisture cannot be measured at this scale over the large areas an atmospheric model is applied to (e.g.Beven, 2001;Pitman et al., 2004).Also climate simulations rely on a proper root zone soil moisture initialization (e.g.Conil et al., 2009;Seneviratne et al., 2006).
Remotely sensed land surface data and air-temperature are currently assimilated to overcome errors in soil moisture and temperature simulation in NWP models (e.g.Hess, 2001;Seuffert et al., 2004;Crow and Wood, 2003;Gao et al., 2007).
Still unresolved problems are the soil moisture analysis in densely vegetated areas and in the root zone.Recently, various approaches of data assimilation were set up and analyzed to retrieve the root zone soil moisture at the regional scale in hydrological Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion models.They mainly use Kalman Filter techniques and their modifications, which are outlined in detail e.g. by Evensen (2006).Further, Evensen (2003) gives a detailed description and literature review of the Ensemble Kalman Filter (EnKF, Evensen, 1994).Walker et al. (2002) apply a modified Kalman Filter technique with a distributed hydrological model to retrieve the three-dimensional soil moisture from surface soil moisture measurements.This is a valuable approach in hydrology but due to the intense computational cost of a distributed hydrological model not a tool currently suitable for NWP.
At the German Weather Service (DWD) Hess ( 2001) implemented a method based on the EKF (Extended Kalman Filter) technique into COSMO that adjusts the soil state to meet the observed atmospheric state.However, in his approach the soil moisture and soil temperature do not necessarily reflect the reality, i.e., its usage is not consistent with the hydrologic interaction of the land surface and lower atmosphere.This is proven by Drusch and Viterbo (2007), who assimilated screen-level variables in ECMWF's Integrated Forecast System.It becomes a problem when improving parameterizations that rely on quantities that changed by data assimilation, e.g.soil moisture and temperature or atmospheric fields like the water vapour.Moradkhani et al. (2005) and Dunne and Enthekhabi (2006) e.g.use the Ensemble Kalman Smoother for root zone soil moisture analysis assimilating L-band radiobrightness temperatures in an area of the Southern Great Plains (USA) whose vegetation is mainly wheat and grasses (Drusch et al., 2001).
A data source that has only received attention in the past couple of years is streamflow from operational river gauging networks.Streamflow is a quantity that can be measured at relatively high accuracy (about >90%, LfU, 2002).If the runoff is transported to and within the river network, it can be compared to measured streamflow at gauging stations.plied the EnKF for soil moisture update in real time flood forecasting in a 622 km 2 catchment in Austria.However, they use a simple soil moisture model which is not complex enough as land surface model for atmospheric models.Clark et al. (2008) demonstrate that the standard implementation of EnKF is inappropriate and show the improved performance when streamflow is transformed into log space before applying EnKF with the distributed hydrological model TopNet.This is due to the large ranges in streamflow between peak flow and low flow, which can be 2 orders of magnitude or more.Streamflow analyses also allow for an evaluation of the model performance (e.g.Lohmann et al., 2004;Warrach-Sagi et al., 2008).In this study, we go a step further and study the potential of streamflow data assimilation for soil moisture analysis in a catchment, namely for initialisation of numerical weather prediction and climate models.We followed the most recent development in EnKF and applied it to the streamflow data assimilation for soil moisture initialization in a land surface model of the numerical weather predication model COSMO.
In southern Germany a network of automated river and precipitation gauges has been installed in the past couple of years by the federal services for flood monitoring.
The federal state Baden-W ürttemberg has implemented a flood forecast centre, which is able to provide half-hourly updates of streamflow measurements at approximately 140 gauges at the rivers Rhein, Neckar, Donau, Main and their main contributeries.Similar warning systems are available in the federal state Bayern and Rheinland-Pfalz.
Such automated networks provide a valuable source for operational streamflow data assimilation.
The square root algorithm for the EnKF (Evensen, 2004) is set up to assimilate streamflow data in TERRA-ML to analyse the soil water content of the soil profile down to 2.43 m soil depth simulated by TERRA-ML.By mean of an observing system simulation experiment (OSSE) in the Enz River catchment (Germany) we studied the potential of streamflow data assimilation and its spatial and temporal requirements for an automated river gauging network.The Enz catchment is on the downwind side of the Black Forest, i.e.QPF by the weather forecast model is often underestimating, making it a Introduction

Conclusions References
Tables Figures

Back Close
Full The study is carried out exemplarily with the land surface model TERRA-ML coupled to a river routing model (Warrach-Sagi et al., 2008).The multi-layer soil and vegetation model TERRA-ML serves as the lower boundary of the operational non-hydrostatic mesoscale weather forecast model COSMO (Doms et al., 2005).(COSMO is the acronym for the Consortium for Small-scale Modelling (http://www.cosmo-model.org/.)However, the data assimilation system can be set up for any land surface model that includes a river routing model to simulate streamflow.

Description of TERRA-ML and the river routing scheme
This study applies the stand alone version of TERRA-ML in the framework set up by Ament and Simmer (2006).The model configuration and parameters of TERRA-ML are taken from the German Weather Service's COSMO.In the framework, TERRA-ML is set up as if it is called by the COSMO, with the exception that the meteorology is read from a file instead of forecasted at the time step by the COSMO.This framework has the advantage that it allows the simulation of a gridded area (e.g.watershed) per time step mimicking a simulation with a weather forecast model.An important modification of TERRA-ML in this study is the parameterization of the hydraulic conductivity and diffusivity following Campbell (1974) instead of Rijtema (1969) due to the results of Graßelt et al. (2008) and Warrach-Sagi et al. (2008).TERRA-ML and the river routing model are set up as described in detail by Warrach- els.For each grid cell the one-dimensional vertical land surface model TERRA-ML is applied.The locally generated runoff of the LSM needs to be transported into and along the river system to compare it to streamflow measurements at gauging stations and to calculate the streamflow at various locations of the river.Based on the routing scheme described in detail by Lohmann et al. (1996), Lohmann et al. (2004) present a lumped optimized linear routing model, which Warrach-Sagi et al. ( 2008) coupled to TERRA-ML.The routing scheme describes the time runoff takes to reach the outlet of a grid cell and the water transport in the river network.It is assumed that water flows uni-directionally from grid cell to grid cell with eight possible directions through each side and corner of the grid cell.

The streamflow data assimilation system
The Ensemble Kalman Filter (EnKF) has been reviewed by many authors recently (e.g.Evensen, 2003Evensen, , 2006;;Pauwels and DeLannoy, 2006;Clark et al., 2008) and therefore here only a short description of its implementation for the streamflow data assimilation is given.
Both, model results and observation, deviate from the true state.The goal of data assimilation is to find the best estimate of the state (e.g.soil moisture) from model simulations and measurements (e.g.streamflow).One method is to estimate the mean state and the "maximum likelihood" including its covariance as uncertainty measure as it is provided e.g. by the EnKF.
Various algorithms solve the EnKF equations (see e.g.Evensen, 2006).For this study we chose the square root algorithm for EnKF (http://enkf.nersc.no)from Evensen (2004) due to the following aspects: it is stable, needs relatively little computing time, requires relatively little memory and it is straight-forward to implement.The following base line equations describe the EnKF as it is implemented: x a e,n = x b e,n + K e,n (y e,n − H(x b e,n )) ( 5) Bold letters represent matrices, x and y are the vectors for the model state and observation.b is the background (i.e. initial state), a is the analysis, e is the ensemble member, n is the time step, T is the transpose and K is the Kalman gain matrix.A, B and R are error covariance matrices of the analysis, background and observation, I is the identity matrix, H is the observation operator (in this case the river routing model), which transforms the variable from model space to observation space, H is the tangent linear observation operator matrix of H and M is the model operator (in this case TERRA-ML).
Depending on the location within the catchment the water needs more or less time to travel as streamflow through the river network.Water far away from the gauge arrives later than the runoff from grid cells close by.The travel time depends on the river itself and the form and orography of the catchment.This means that the streamflow measured at a gauge depends on the soil moisture distribution in the catchment for a time window from time t=0 to t=m * d t. m is the time step, d t is the time interval of one time step.This period of streamflow data needs to be assimilated.By this the EnKF Introduction

Conclusions References
Tables Figures

Back Close
Full on the catchments' time window to obtain the soil moisture (see e.g.Sect.4.2).Following Clark et al. (2008), the streamflow is transformed into log space before computing the error covariances since they demonstrated that this improves the filter performance.

Study area: the Enz
The Black Forest is a mountain range that reaches from 47.5 • S to 49 The 260 km 2 Enz catchment upstream of Pforzheim (upstream of the Nagold confluence) is on the downwind side of the Black Forest, i.e. precipitation is often underestimated by weather forecast models.Therefore the Enz catchment (Fig. 1) was chosen for the streamflow data assimilation study.No water reservoirs interrupt the river system.Elevation of the catchment ranges between 350 and 930 m above sea level.The catchment is characterized by forested (mixed deciduous and evergreen coniferous trees) upland areas and agriculturally used lowlands.Sandy and loamy soils dominate the upper Enz area (Fig. 2).Between 1997 and 2002 annual precipitation in the catchment ranged from 1088 to 1451 mm.

Conclusions References
Tables Figures

Back Close
Full  2008) applied the coupled TERRA-ML-routing model to the Enz catchment upstream of Pforzheim and compared it to simulations of the flood forecast centre Baden-W ürttemberg and to observations.They showed that the model results and observations agree reasonably well.However, as is always the case, model results and observations both include errors and both differ from the true state.To assess the potential and requirements for streamflow data assimilation, an OSSE is set up.The results of the TERRA-ML-routing model for 1997 in the Enz river catchment (Warrach-Sagi et al., 2008) are assumed to be the "true" state, named "CONTROL" hereafter.
The data assimilation experiment starts on the 5 May 1997 with an ensemble of initial soil moisture fields in the catchment and an ensemble of streamflow at various locations in the river network.The CONTROL streamflow serves as "observation" which is assimilated for the soil moisture analysis.The analysis is then compared to the "true" state, i.e. the CONTROL soil moisture.
A flow duration check is carried out to obtain the assimilation time window for the whole basin at Pforzheim (upstream of Nagold confluence) and the sub basins Große Enz (90 km 2 ), Kleine Enz (71 km 2 ), Eyach (43 km 2 ) and upstream of H öfen (222 km 2 ), downstream of the confluence of the Eyach into the Enz (Fig. 1).For the flow duration check at the initial time step 0.002 kg/m 2 runoff are assumed for each grid cell.No more runoff is assumed afterwards.The routing model calculates the streamflow for each catchment (Fig. 3).Depending on the size and structure of the catchment, the time window until all water has left the catchment varies between 25 and 62 h.Experiments showed that in most cases an assimilation window of 90% of the time window lead to the best results in soil moisture distribution and catchments' mean soil moisture.This is due to the fact that towards the end of the time window too long streamflow time series are assimilated into the grid cells close to the observation location.This will be discussed in Sect.4.4 in more detail.Introduction

Conclusions References
Tables Figures

Back Close
Full

Ensemble preparation
For the OSSE a period is chosen, which does not include extreme events (such as flooding or drought or strong precipitation).A period in spring was chosen, when not only soil texture but also vegetation and weather control the soil moisture.Furthermore, in spring and summer soil moisture impacts the development of convection in the atmosphere.This study starts on the 5 May 1997 (day 125).The initial soil moisture of the CONTROL simulation is perturbed applying the 2-D-pseudorandom sampling method and algorithm (http://enkf.nersc.no) of Evensen (2004) to obtain 100 ensemble members of initial soil moisture fields, which include no step-functions within the 2-D-area.
(See Evensen (2004) for more details on this approach.)The soil moisture of each grid cell is chosen to vary between +10% and −40% of the CONTROL soil moisture.This is to account for the typical underestimation of precipitation in NWP simulations in this area and to account for the fact, that the precipitation might have been simulated in the wrong location within the catchment.The 2-D-pseudorandom fields vary between up to d =±1 and examples are shown for 2 ensemble members in Fig. 4. According to the random number of each grid cell (i , j, k) of each ensemble member e, the soil moisture η in the grid cell is perturbed to i , j and k are the indices of the grid cell in eastward, northward and downward direction, c is the control state.Like in nature, soil moisture in the ensemble for each grid cell is always limited between saturation and air dryness point.
The CONTROL streamflow is perturbed by adding Gaussian noise.The 1-D-pseudorandom sampling method and algorithm (http://enkf.nersc.no) of Evensen (2004) to obtain 100 ensemble members is applied and streamflow perturbed by up to ±15%, assuming that the error might be occasionally larger than the <10% assumed by LfU (2002).Introduction

Conclusions References
Tables Figures

Back Close
Full

Soil moisture analysis
TERRA-ML's soil column is 2.43 m deep (see Sect. 2).The soil water content (SWC) of each grid cell depends on its soil depth and its soil moisture η with the density of water ρ w , z(k) is the depth of the lower boundary of soil layer k, l is the lowest soil layer.Streamflow data is assumed to be available at an half hourly time step like the observations made by the automated gauges in Baden-W ürttemberg.
Figure 5 shows the ensemble spread of the catchments' mean SWC at the initial time t=0 for the background and the analysis for the catchment upstream of Pforzheim assimilating streamflow data from Pforzheim.The analysis ensemble has a lower spread and is closer to the CONTROL SWC.The ensemble mean SWC at t=0 is 525 kg/m 2 for the background, 544 kg/m 2 for the analysis and 557 kg/m 2 for the CONTROL.Fig- ure 6 shows the spatial distribution of the SWC at time t=0.Note that single cells show larger SWC mainly due to different soil texture (peat and loam, see Fig. 2).While the ensemble mean of the background SWC is everywhere 5-6% lower than the CON-TROL SWC, the ensemble mean analysis SWC shows an improvement (Fig. 7).The analysis differs in more than half of the catchment by ±4% from the CONTROL SWC.
Only in a few upstream grid cells it is worse (8%) than the background SWC.In about 1/3 of the catchment the analysis only differs ±2% from the CONTROL.The promising results from the 260 km 2 catchment led to a study about the potential impact of a denser network of gauges for the soil moisture analysis.Gauges were assumed to be at the outlet of the Große Enz (90 km 2 ), the outlet of the Kleine Enz (71 km 2 ) and the outlet of the Eyach (43 km 2 ).Little impact was reached for the Eyach, Introduction

Conclusions References
Tables Figures

Back Close
Full mean SWC at t=0 is 535 kg/m 2 for the background, 571 kg/m 2 for the analysis and 568 kg/m 2 for the CONTROL.Figures 9 and 10 show that only close to the outlet the analysis leads to worse SWCs than the background.Figure 11 shows the impact of assimilating streamflow from the CONTROL model simulation from the Große Enz outlet, Kleine Enz outlet and Eyach with a 0.5 hourly time step.Most areas show a positive or neutral impact of the data assimilation.
All in all the simulations show a gradient in the impact of the data assimilation.Close to the gauge location of the assimilated streamflow and at the furthest upstream grid cells the data assimilation shows worse results than in the middle areas.This is due to flow duration in the river network and the assimilation window.The grid cells close to the gauge would need shorter assimilation windows.However, the OSSE shows that the streamflow data assimilation has the potential to improve the soil moisture throughout the catchment and that a more dense gauging network would help to improve this even further.

Conclusions
Numerical weather forecasting and climate modeling require an accurate soil moisture initialization for their land surface models.So far the areal distribution of root zone soil moisture cannot be measured.Streamflow depends on the soil moisture of a river catchment and is measured at gauging stations of the rivers at relatively high accuracy.
A retrospective EnKF was set up to assimilate streamflow into the multi-layer land surface model TERRA-ML of the regional weather forecast model COSMO.An OSSE was performed in the Enz River catchment located at the downwind side of the northern Black Forest (Germany).The results confirm the potential of streamflow data assimilation for improving soil moisture analyses.Further, we discussed the spatial and temporal requirements for an automated river gauging network.Half-hourly streamflow data is available from the automated gauges of the flood forecast centre of Baden-W ürttemberg (Germany) for approximately 140 gauges.Half-hourly resolution Introduction

Conclusions References
Tables Figures

Back Close
Full of streamflow data is sufficient for its assimilation for soil moisture analysis.In the upper Enz an automated gauge is operational at H öfen.The OSSE shows that streamflow from this location can already improve SWC in the Enz catchment upstream of H öfen, but that a denser network would improve the SWC even more.Namely at the outlets of smaller sub-catchments, like the Große Enz this would be valuable, since the sub-catchments show a differently structured river network (Fig. 1) and flow duration (Fig. 3).Since the necessary assimilation window depends on the catchment size (e.g.48 h for Pforzheim and 21 h for the Kleine Enz), a denser gauging network would shorten the assimilation time making it even more valuable for initialisation in numerical weather forecast models.
All in all the retrospective EnKF is a powerful method to assimilate streamflow data into a land surface model for root zone soil moisture analysis.The implementation of the square root algorithm for EnKF from Evensen (2004)         579 Pauwels and De Lannoy (2006) published the application of a retrospective EnKF to assimilate streamflow data for soil moisture retrieval.Their synthetic tests show promising results for a 1000 km 2 catchment in Belgium and indicate improvements especially in case of precipitation underestimation.They apply it to the high resolution hydrological model TOPLATS.Komma et al. (2008) successfully ap- for streamflow data assimilation(Pauwels and De Lannoy, 2006).
Sagi et al. (2008), therefore here only a summary is given.In COSMO the model TERRA-ML has got 6 hydrological layers (layer depths from the surface: 0.01 m, 0.03 m, 0.09 m, 0.27 m, 0.81 m, 2.43 m) and 8 thermal layers (layer depths from the surface: 0.01 m, 0.03 m, 0.09 m, 0.27 m, 0.81 m, 2.43 m, 7.29 m, 21.87 m).The lower boundary condition is given by free drainage at 2.43 m depth and a constant climatological temperature below 7.29 m depth.For model simulations, watersheds are divided into grid cells as in atmospheric mod- Figure 1

Fig. 1 .Figure 3 Fig. 3 .
Fig. 1.The orography (based on the 90 m-orographic data from the SRTM) and the river network for the Enz catchment upstream of Pforzheim on the rotated spherical coordinate system of the COSMO on a grid resolution of 0.01 • (approx. 1 km).

Fig. 4 .Fig. 5 .Fig. 6 .
Fig. 4. The initial soil moisture of the CONTROL simulation is perturbed applying the 2-Dpseudorandom sampling method and algorithm (http://enkf.nersc.no) of Evensen (2004) to obtain 100 ensemble members of initial soil moisture fields.The 2-D-pseudorandom fields vary between up to d =±1 and examples are shown for 2 ensemble members.

Fig. 7 .
Fig. 7.The difference in soil water content (SWC) of each grid cell relative to the CONTROL SWC of each grid cell for the initial time t=0 (5 May 1997) for the background and analysis assimilating streamflow from the CONTROL model simulation from Pforzheim from t=0 to t=48 h with a 0.5 hourly timestep.

Fig. 8 .Fig. 9 .
Fig.8.TERRA-ML's soil column is 2.43 m deep.The soil water content (SWC) of each grid cell depends on its soil depth and its soil moisture (Eq.7).Distribution of initial mean SWC (t=0) in the Große Enz catchment (90 km 2 ) between the ensemble members.CONTROL mean SWC is 568 kg/m 2 .

Fig. 10 .
Fig. 10.The difference in soil water content (SWC) of each grid cell relative to the CONTROL SWC of each grid cell for the initial time t=0 (5 May 1997) for the background and analysis assimilating streamflow from the CONTROL model simulation from Große Enz outlet from t=0 to t=30 h with a 0.5 hourly timestep.

Fig. 11 .
Fig. 11.The difference in soil water content (SWC) of each grid cell relative to the CONTROL SWC of each grid cell for the initial time t=0 (5 May 1997) for the background and analysis assimilating streamflow from the CONTROL model simulation from the Große Enz outlet, Kleine Enz outlet and Eyach outlet with a 0.5 hourly timestep.
• N at a width of approximately 50 km in Baden-W ürttemberg (Germany).Reaching from North to South, the Black Forest modifies significantly most frontal systems arriving from the Atlantic.In spite of its relatively low height (largest mountain Feldberg 1493 m a.s.l.), orographic lifting of unstable and moist air masses in this region results in the largest amount of precipitation in Germany except the northern front range of the Alps.In summer, the Black Forest is characterized by strong convection, thunderstorms, and the development of extreme precipitation events.The eastern Black Forest hosts about half of the contributories of the Neckar, a major contributory of the Rhine.The western Black Forest drains directly to the Rhine.Most rivers in Baden-W ürttemberg contain automated gauging stations from the flood forecast center and streamflow data are available every 30 min.