Hydrologic Data Assimilation: Merging Measurements and Models Steve Margulis Assistant Professor

1. Hydrologic Data Assimilation: Merging Measurements and Models Steve Margulis Assistant Professor Dept. of Civil and Environmental Engineering UCLA CENS Technical Seminar Series April 29, 2005

2. Outline • Introduction and Motivation • Data Assimilation and the Ensemble Kalman Filter • Application: Soil Moisture Estimation via Assimilation of Microwave Remote Sensing Observations into a Hydrologic Model • Application to embedded sensor networks (?): Palmdale Wastewater Irrigation Site (future collaboration with Tom Harmon)

3. Introduction • Hydrologic Observations: • Benefits: • Provide important diagnostic information about real conditions • Yield important validation and model forcing databases • Limitations: • Measurements generally sparse in time and/or space • interpolation/extrapolation • downscaling/upscaling • Contain measurement error • Often measuring states not fluxes • Hydrologic Models: • Benefits: • Representation of our knowledge of physical processes (dynamics) • Physical relationships between observables and states/fluxes of interest • Numerical tool for prediction • Limitations: • Simplified abstractions of reality • Subject to uncertainties in time-varying inputs/time-invariant parameters How can we combine the benefits of both in an optimal framework?

4. What is data assimilation? • Goal: Data assimilation seeks to characterize the true state of an environmental system by combining information from measurementsandmodels. • Typical measurementsfor hydrologic applications: • Ground-based hydrologic and geological measurements (stream flow, soil moisture, soil properties, canopy properties, etc.) • Ground-based meteorological measurements (precipitation, air temperature, humidity, wind speed, etc.) • Remotely-sensed measurements which are sensitive to hydrologically relevant variables (e.g. water vapor, soil moisture, etc.) • Mathematical models used for data assimilation: • Models of the physical system of interest • Models of the measurement process • Probabilistic descriptions of uncertain model inputs and measurement errors A description based on combined information should be better than one obtained from either measurements or model alone.

5. Key Features of Environmental Data Assimilation Problems State estimation -- System is described in terms of state variables(random vectors), which are characterized from available information Multiple data sources -- Estimates are often derived from different types of measurements (ground-based, remote sensing, etc.) measured over a range of time and space scales Spatially distributed dynamic systems -- Systems are often modeled with partial differential equations, usually nonlinear. Through discretization the resulting number of degrees of freedom (unknowns) can be very large. Uncertainty -- The models used in data assimilation applications are inevitably imperfect approximations to reality, model inputs may be uncertain, and measurement errors may be important. All of these sources of uncertainty need to be considered systematically in the data assimilation process.

6. State-space Framework for Data Assimilation • State-space concepts provide a convenient way to formulate data assimilation problems. Key idea is to describe system of interest in terms of following variables: • Input variables -- variables which account for forcing from outside the system or system properties which do not depend on the system state. • State variables -- dependent variables of differential equations used to describe the physical system of interest, also called prognostic variables. • Measurements -- variables that are observed (with measurement error) and are either directly or indirectly related to states. • Output variables -- variables that depend on state and input variables, also called diagnostic variables.

7. Basic Probabilistic Concepts in Data Assimilation • Uncertain forcing (u) and parameter (a) inputs: • Postulated unconditional PDFs:fu ( u)andfa(a ) • Uncertain States (y): • Derived (from state eq.) unconditional PDF: fy ( y) • Uncertain measurements (z): • Measurement PDF (error structure): fz ( z) • Knowledge of state after measurements included: • Characterized by conditional PDF: fy| z (y| z) • (Bayes Theorem)

8. Components of a Typical Hydrologic Data Assimilation Problem Time-varying input u(t) (e.g. precip.) State y (t) (e.g. soil moist.) Hydrologic system Measurement system Specified (mean) True True True Random fluctuations Output zi (e.g. radiobrightness) Random error, n Random fluctuations Specified (mean) Measured Time-invariant input  (e.g. sat. hydr. cond.) Data assimilation algorithm Means and covariances of true inputs and output measurement errors Estimated states and outputs State Eq: Measurement Eq: The data assimilation algorithm uses specified information about input uncertainty and measurement errors to combine model predictions and measurements. Resulting estimates are extensive in time and space and make best use of available information.

9. Characterizing Uncertain Systems p[y(t)| Zi] Std. Dev. y: p(y) p(y | Zi) Prior Conditional (Posterior) Mode Mean y(t) Zi What is a “good characterization” of the system state y(t), given the vector Zi = [z1, ..., zi] of all measurements taken through ti? Theposterior probability densityp(y| Zi) is the ideal estimate since it contains everything we know about the state y given Ziand other model inputsu and a. In practice, we must settle for partial information about this density • Variational DA: Derive mode of p[y(t)| Zi] by solving batch least-squares problem • Sequential DA: Derive recursive approximation of conditional mean (and covariance?) of p[y(t)| Zi]

10. Monte Carlo Approach: Ensemble Filtering p[y(ti)|Zi] Update with new measurement (zi+1) p[y(ti+1)|Zi+1] Update with new measurement ( zi+1 ) y j(ti+1| Zi+1) y j(ti| Zi) y j(ti+1| Zi) Propagate forward in time Propagate forward in time Time p[y(ti+1)|Zi] ti ti+1 Time ti+1 ti Divide filtering problem into two steps – propagation and update. Characterize random states with an ensemble (j = 1, … , J) of random replicates: Evolution of posterior probability density Evolution of random replicates in ensemble It is not practical to construct and update complete multivariate probability density. Ensemble filtering propagates only replicates (no statistics). But how should update be performed?

11. The Ensemble Kalman Filter (EnKF) Propagation stepfor each replicate (y j): Update step for each replicate: meas. residual K = Kalman gain derived from propagated ensemble sample covariance. K=Cyz [Czz+Cn]-1 After each replicate is updated it is propagated to next measurement time. No need to update covariance (i.e. traditional Kalman filter)—results in large computational savings.

12. Application: Microwave Measurement of Soil Moisture 1 sand silt 0.9 clay 0.8 microwave emissivity [-] 0.7 0.6 0.5 0 0.2 0.4 0.6 0.8 1 saturation [-] Land surface microwave emission (at L-band) is sensitive to surfacesoil moisture (~ 5 cm). • Measurement Limitations: • indirect measurement of soil moisture – inversion? • sparse in time (~ 1 measurement per day) – interpolation? • spatially coarse (~10s of kilometers) – downscaling? • contains information about surface moisture only (want rootzone soil moisture) – extrapolation?

13. Test Case: Application to SGP97 Experiment Site • Month-long experiment in central OK in summer 1997 (~10,000 km2 area) • Daily airborne L-band microwave observations (17 out of 30 days) to test feasibility of soil moisture estimation from space • Ground-truth soil moisture sampling performed daily at validation sites Can we use EnKF to map rootzone soil moisture fields and associated surface fluxes from microwave measurements? Margulis et al., 2002; 2005

14. Key Features of Problem • Off-the-shelf models • Hydrologic: NOAH LSM • Radiative Transfer: Jackson et al. (1999) • Spatially-distributed states and parameters • Dealing with model nonlinearities and input uncertainties • Real-time (sequential) estimation • Next generation satellite observations (L-band passive microwave)

15. Spatially variable model inputs Sand fraction 50 km 0 0.2 0.4 0.6 0.8 NOAH soil class NOAH vegetation class Meteor. Stations RTM Inputs Clay fraction El Reno 0 2 4 6 8 0 2 4 6 8 10 12 NOAH Model Inputs 0 0.05 0.1 Estimation region ~ 40 by 280 km (11 by 70 pixels--4 km resolution)

16. Illustrative Results: Sequential Updating 0.06 Estimated Vol. Soil Moisture Estimated Vol. Soil Moisture Error 0.04 (Before update) (After update) (Before update) (After update) 0.02 0.3 Day 169 0.2 0.02 0.1 0.015 0.3 0.01 Day 179 0.2 0.03 0.1 0.02 0.3 Day 184 0.01 0.2 0.1 • Left columns show estimated soil moisture fields before and after assimilating Tb • Right columns show estimated error in soil moisture fields from ensemble • Information in observations used to not only update mean fields, but reduces uncertainty

17. Illustrative Results: Downscaling 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 Observing System Simulation Experiments (OSSEs) Used to Investigate Impact of Coarse Measurements Microwave Observations (Tb in ºK): 4 km 12 km 40 km True Vol. Soil Moisture Field Day 178 Generate obs. at different meas. resolutions Estimated Vol. Soil Moisture Fields: Space-time averaged results

18. Illustrative Results: Interpolation 0.5 CF08 0.4 0.3 0.2 0.1 Ground truth gravimetric meas. 0 Individual replicates 170 175 180 185 190 195 Estimate (Cond. mean) Open Loop (Uncond. Mean) Volumetric Soil Moisture Day of Year ER Precip. Comparison of Estimates to Real Ground-truth Time series Microwave obs. times

19. Illustrative Results: Extrapolation/Flux Estimation Surface evaporation flux (latent heat) is a function of entire rootzone moisture, not just surface. Is information in radiobrightness propagating to sub surface? Note “spin-up” effect of filter during first 10 days Over time, information from Tb about surface conditions propagates downward through rootzone

20. Summary of Results • Data assimilation (in this case using the EnKF) allows for merging of model and data. Key benefits of this framework: • inversion of electromagnetic measurement into estimates of hydrologic states of interest (soil moisture) • downscaling of coarse microwave radiobrightness measurement resolution to estimation scale (similar potential for upscaling?) • value added data products which are essentially continuous in time/space (interpolation between sparse measurements) • extrapolation/propagation of information to unobserved portions of domain (subsurface states) via incorporation of model physics • estimates of additional outputs of interest (e.g. fluxes) which are difficult to measure directly • estimates of uncertainty about mean estimate (via error propagation through system)

21. CENS Example: Wastewater Reuse in Mojave Desert Where does the County Sanitation District (CSD) of Los Angeles put 4 million gallons per day of treated wastewater in a landlocked region? Can we use embedded sensors to track infiltration plume, assess nitrate concentrations, apply feedback control? Reclaimed wastewater irrigation pivot plots Palmdale, CA wastewater treatment plant (slide courtesy of Prof. Tom Harmon)

22. Distributed Monitoring and Adaptive Management Approach Monitoring network design: How many sensors can we get away with? How do we optimally place them? Interpolating between sensors/extrapolating to depth: Distributed parameter models Stochastic approaches mote image by Jason Fisher (Cal-CLEANER) (slide courtesy of Prof. Tom Harmon)

23. Site characterization At the field scale: rigorous characterization sampling being done geostatistical parameterization techniques indicator kriging (probability Ks exceeds...) ordinary kriging (Ks) (slide courtesy of Prof. Tom Harmon)

24. Proposed Research/Experiments Task 1: Model and EnKF Interface Design/Implementation Implementation of stochastic version of hydrologic flow/transport model Input error model analysis using site characterization studies EnKF “wrapper” design Task 2: Network Design with Observing System Simulation Experiments Model used to generate different measurement scenarios Evaluation of scenarios using OSSEs to determine optimal sensor locations, sensor numbers, etc. (via minimization of state estimation error) Data Assimilation (specifically the EnKF) proposed as a potential tool for investigating these research and operational implementation questions

25. Proposed Research/Experiments (cont.) Task 3: Real-time State and Parameter Estimation After network deployment, use as real-time state estimation tool Take advantage of early-life of sensors (accurate/stable error structure) to calibrate model parameters Use real-time state estimates for feedback control Task 4: Real-time Network Monitoring and Maintenance What about degradation of sensor network over time? Once model parameters are estimated, can measurement error be parameterized to detect changes in measurement error structure?

26. Summary • Data assimilation provides a very general framework for merging measurements and models • inversion, interpolation/extrapolation, uncertainty propagation, etc. • In hydrology, these techniques have primarily been used in the context of remote sensing due to limited availability of in-situ measurements • Problems where embedded sensor networks can be deployed are ideal candidates for application of these techniques where the ultimate goal is to maximize extraction of information content from measurements.

27. Acknowledgments Funding for Research: NSF Water Cycle Research Grant Collaborators: Dara Entekhabi (MIT) Dennis McLaughlin (MIT)

28. Some Helpful Data Assimilation References McLaughlin, D., 1995: Recent developments in hydrologic data assimilation, U.S. Natl. Rep. Int. Union Geod. Geophys. 1991-1994, Reviews in Geophysics, 33, 977-984. Margulis, S.A., D. McLaughlin, D. Entekhabi, and S. Dunne, 2002: Land data assimilation and soil moisture estimation using measurements from the Southern Great Plains 1997 field experiment, Water Resources Research, 38(12), 1299, doi:10.1029/2001WR001114. Evensen, G., 2003: The Ensemble Kalman Filter: theoretical formulation and practical implementation, Ocean Dynamics, 53, 343-367.