# Daily Data Assimilation of a Hydrologic Model Using the Ensemble Kalman
Filter

Sami A. Malek
UC Berkeley, Department of Electrical Engineering and Computer Sciences
<EMAIL_ADDRESS>
&Alexandre M. Bayen
UC Berkeley, Department of Electrical Engineering and Computer Sciences
UC Berkeley, Department of Civil Engineering
UC Berkeley, Institute for Transportation Studies
<EMAIL_ADDRESS>
&Steven D. Glaser
UC Berkeley, Department of Civil Engineering
<EMAIL_ADDRESS>

###### Abstract

Accurate runoff forecasting is crucial for reservoir operators as it allows
optimized water management, flood control and hydropower generation. Land
surface models in mountainous regions depend on climatic inputs such as
precipitation, temperature and solar radiation to model the water and energy
dynamics and produce runoff as output. With the rapid development of cheap
electronics applied in various systems, such as Wireless Sensor Networks
(WSNs), satellite and airborne technologies, the prospect of practically
measuring spatial Snow Water Equivalent in a dense temporal scale is
increasing. We present a framework for updating the Precipitation Runoff
Modeling System (PRMS) with Snow Water Equivalent (SWE) maps and runoff
measurements on a daily timescale based on the Ensemble Kalman Filter (ENKF).
Results show that by assimilating SWE daily, the modeled SWE gets updated
accordingly, however no improvement is observed at the runoff model output.
Instead, a deterioration consistently occurs. Augmenting the state space with
model parameters and runoff model output allows for filter update with
previous day measured runoff using the joint state-parameter method, and
showed a considerable improvement in the daily runoff output of up to 60%
reduction in RMSE for the wet water year 2011 relative to the no assimilation
scenario, and improvement of up to 28% compared to a naive autoregressive
AR(1) filter. Additional simulation years showed consistent improvement
compared to no assimilation, but varied relative to the previous day
autoregressive forecast during the dry year 2014.

_K_ eywords Filtering $\cdot$ EnKF $\cdot$ Data Assimilation $\cdot$ Runoff
Model $\cdot$ Sensors $\cdot$ Hydrology $\cdot$ Forecasting $\cdot$ Estimation
$\cdot$ PRMS

##  1 Introduction

This research work fits into a larger objective of attaining consistently more
accurate snowmelt runoff forecast resilient to climate change and extreme
weather conditions. Presently, reservoir and dam operators rely on forecasts
largely based on historical methods. Numerous studies show that climate change
effects are challenging the stationarity assumption of this statistical
approach. The past can no longer predict the future as well as it did before.
Other studies show a systematic receding date for snow melt during the year.
[1, 2, 3, 4]

A mandatory evacuation order was issued in February 2017 after the Oroville
dam crisis because of hazardous water overflow at the Oroville reservoir. The
Department of Water Resources (DWR) published a report showing that reservoir
water levels rose earlier and higher than previous 16 years. Oroville is the
pouring point of the North Fork Feather River where 3 days earlier, recently
deployed sensors showed an onset of rapid snow melt up to 1 cm per hour most
probably due to a rain on snow event. Unlike statistical models, a physically
based forecast model that incorporates these new data should be able to
capture such inflow surges.

In addition to flood control, runoff forecasting allows for better reservoir
management and potentially more efficient power generation. Stakeholders, such
as Pacific Gas & Electric (PG&E) have already started shifting towards
physically based forecast models. These models though less vulnerable to
climate change, require dense and accurate data as well as a methodology to
seamlessly integrate sensor measurements of different temporal and spatial
scales.

Recent efforts to increase hydrologic instrumentation in the East Branch of
the North Fork Feather River basin by a partnership between the University of
California at Berkeley, the Department of Water Resources and the California
Energy Commission (CEC) led to the deployment of four 1 km2 scale state-of-
the-art remote autonomous ground measurement systems based on Wireless Sensor
Networks technology to collect temperature, humidity, solar radiation,
snowdepth, soil temperature and soil moisture measurements at a 15 minute
scale. These new deployments would supplement the existing measurement system
consisting of 8 snow pillows reporting daily SWE measurements and the monthly
manual snow course locations. Moreover the NASA airborne snow observatory, in
partnership with DWR are using airborne LiDAR system to produce unprecedented
bi-weekly basin-wide SWE maps of the Tuolumne and Merced river basins at 50
meter spatial resolution [5].

The objective of this study is to design a practical data assimilation
framework for real-time runoff forecasting using PRMS that can leverage sensor
data from the existing, state-of-the-art and future hydrologic measurement
systems at different temporal and spatial scales.

After describing the PRMS model dynamics, we present existing filtering
techniques, then we develop a practical data assimilation framework to
integrate in real-time simulated measurements of SWE and real measurements of
streamflow in the model using the Ensemble Kalman Filter.

##  2 Precipitation Runoff Modeling System

The Precipitation Runoff Modeling System is a hydrologic model originally
described by Leavesley (1983). It simulates snow accumulation and ablation
processes as well as soil-zone and streamflow water dynamics. The watershed or
basin is discretized into Hydrologic Response Units (HRUs) based on
topography, land use, climate, soil properties, and geologic units information
[6]. Sample HRUs for the East Branch of the Feather River basin are shown in
Fig. 1.

Figure 1: Hydrologic Response Units (HRUs) of the East Branch of the Feather
River overlaid on a SWE map from 1998.

The model runs on a daily time discretization, with some processes having
half-day resolution. The model takes daily precipitation, minimum and maximum
air temperature and optionally short-wave solar radiation as inputs. PRMS is a
module based system: it can handle other data inputs as well as different
modeling modules. We restrict this study to what is currently being
operationally used in our specific application. Each HRU consists of a system
of abstract water containers (or storages) interconnected by water flows as
shown in Fig. 2. Variables with subscripts “h" in this report indicate that
those variables are unique for each HRU.

Figure 2: Schematic diagram of PRMS and its inputs [7]

The components are Interception ($Sint_{h}$), Snowpack ($SWE_{h}$ and
$H_{h}$), Impervious Zone Reservoir ($Simp_{h}$), Soil Zone Reservoir
($Ssz_{h}$, subdivided into Recharge $Ssre_{h}$ and Lower zone $Ssl_{h}$),
Subsurface Reservoir ($Sss_{h}$) and the Ground-Water Reservoir ($Sgw_{h}$).
For each water storage, the water mass balance equation is satisfied.

$\frac{\Delta\left(S\right)}{\Delta t}=\sum{}{}{F}$ (1)

where $S$ and $F$ denote water storage and flow respectively.

$\frac{\Delta\left(Sint\right)}{\Delta
t}=\sum{}{}{Fint_{in}}-\sum{}{}{Fint_{out}}$ (2)

where $Fint_{in}$ consists of precipitation, and $Fint_{out}$ consists of
evaporation, sublimation and througfall.

$\frac{\Delta\left(SWE\right)}{\Delta
t}=\sum{}{}{Fswe_{in}}-\sum{}{}{Fswe_{out}}$ (3)

where $Fswe_{in}$ consists of precipitation on non-vegetated areas,
throughfall from canopy covered regions, and $Fswe_{out}$ consists of
sublimation and snowmelt.

$\frac{\Delta\left(Ssz\right)}{\Delta
t}=\sum{}{}{Fsz_{in}}-\sum{}{}{Fsz_{out}}$ (4)

where $Fsz_{in}$ consists of snowmelt, rain, and flow from interception on
pervious regions, and $Fsz_{out}$ consists of surface runoff,
evapotranspiration and excess $Sz$ water.

$\frac{\Delta\left(Sss\right)}{\Delta
t}=\sum{}{}{Fss_{in}}-\sum{}{}{Fss_{out}}$ (5)

where $Fss_{in}$ consists of recharge from $Ssz$, subsurface flow to
streamflow, and percolation to the ground water reservoirs ($Sgw$).

$\frac{\Delta\left(Sgw\right)}{\Delta
t}=\sum{}{}{Fsgw_{in}}-\sum{}{}{Fsgw_{out}}$ (6)

where $Fsgw_{in}$ consists of recharge from Soil Zone and from subsurface
reservoirs and $Fsgw_{out}$ consists of ground water sink and groundwater flow
to sreamflow output.

For the Snowpack storage($SWE$), melt occurs when the snowpack reaches
isothermal conditions ($H$ = 0). An additional energy balance equation is
satisfied that deals with the atmospheric and temperature gradient and energy
fluxes ($f$) interactions and can be written as:

$\frac{\Delta\left(H\right)}{\Delta t}=\sum{}{}{f}$ (7)

Symbol | PRMS symbol | Type | Quantity | Units | Description
---|---|---|---|---|---
P | hru_ppt | daily | per HRU | inches | total precipitation
Tmax, i | tmax | daily | one |  oF | max temperature at HRU i
Tmin, i | tmin | daily | one |  oF | min temperature at HRU i
Table 1: Model Inputs Symbol | PRMS symbol | Quantity | Units | Description
---|---|---|---|---
Sint | intcp_stor | per HRU | inches | interception storage
Sliq | freeh2o | per HRU | inches | liquid water snowpack storage
$\rho$ | pk_den | per HRU | g/cm3 | snowpack average density
SWE | pkwater_equiv | per HRU | inches | total water snowpack storage
Sice | pk_ice | per HRU | inches | frozen water snowpack storage
D | pk_depth | per HRU | inches | snowpack depth
H | pk_def | per HRU | Langleys | snowpack heat deficit
Tpk | pk_temp | per HRU |  oF | snowpack temperature
fsca | snowcov_area | per HRU | - | snow cover area
Sgw | gwres_stor | per GWR | inches | ground water storage
Simp | imperv_stor | per HRU | inches | impervious storage
Ssre | soil_rechr | per HRU | inches | soil recharge storage
Ssz | soil_moist | per HRU | inches | total soil-zone storage (soil-moisture)
Sss | ssres_stor | per SSR | inches | sub-surface storage
SWEmax | pst | per HRU | inches | tracks max SWE of pack for fsca albedo
Table 2: Model States

###  2.1 Model Inputs

Model inputs consist of daily time-series of precipitation for each HRU and
minimum and maximum air temperature for 1 HRU that contains a temperature
station. The daily short-wave radiation ($Rsw$) is estimated internally by the
model. They are detailed below:

Short-Wave Radiation: A constant solar table consisting of daily estimates of
the potential (clear sky) short-wave solar radiation ($Rpsw$) for each HRU
(they are derived from calculations of duration of sun exposure and latitude
for each HRU and each day [7]). Such values are then adjusted daily based on a
modification of the “degree-day" method to get the short-wave radiation
described in more details in [7]. The method consists of deriving a degree-day
coefficient ($dd$) and then a ratio of actual-to-potential radiation for a
horizontal surface ($r$). The following parametric linear curve is used to
find dd:

$dd=\bm{\phi}_{m}\cdot Tmax_{f}+\bm{\beta}_{m}$ (8)

where $\bm{\phi}_{m}$ and $\bm{\beta}_{m}$ are parameters and $Tmax_{f}$ is an
input. From dd, $r_{h}$ is retrieved via a non-linear relationship of $dd$ in
function of $r_{h}$. Finally the actual solar radiation is adjusted according
to each HRU slope.

$Rsw_{h}=\frac{r_{h}\times\bm{\gamma}}{cos\left(arctan\left(\bm{\delta}_{h}\right)\right)}\cdot
Rpsw_{h}$ (9)

where

$\bm{\gamma}=\begin{cases}\bm{{\gamma}s}&\text{during summer days \& $P$
greater than $\bm{Pmin}$}\\\ \bm{{\gamma}w}&\text{during winter days \& $P$
greater than $\bm{Pmin}$}\\\ 1&\text{$P$ less than $\bm{Pmin}$}\end{cases}$

because $Rpsw_{h}$ is calculated for days without precipitation.

Air Temperature

The single temperature input at HRU i with elevation $e_{i}$ is spatially
distributed for the remaining HRUs according to the following equation (and
similarly for $T_{min}$) in degrees Fahrenheit.

$Tmax_{h}=\bm{Tmax}_{h_{i}}-\bm{{\lambda}1}_{month}\times\left(\frac{\bm{e}_{h}-\bm{e}_{h_{i}}}{1000}\right)-\bm{{\beta}1}_{h}$
(10)

$Tmin_{h}=\bm{Tmin}_{h{i}}-\bm{{\lambda}2}_{month}\times\left(\frac{\bm{e}_{h}-\bm{e}_{h_{i}}}{1000}\right)-\bm{{\beta}2}_{h}$
(11)

Precipitation Precipitation $\mathit{P}_{h}$ for each HRU are inputs to the
model. The precipitation phase is determined solely based on temperature. If
the HRU maximum air temperature ($Tmax_{h}$) is less than or equal to the
parameter $\bm{Ts}$, all precipitation is snow. If $Tmin_{h}$ and $Tmax_{h}$
are greater than or equal to $\bm{Ts}$ and $\bm{Tr}$ respectively,
precipitation is assumed to be all rain. Otherwise, $P$ is considered a
mixture of rain and snow. The rain fraction of $P$ is computed as:

$fr_{h}=\left(\frac{Tmax_{h}-\bm{Ts}}{Tmax_{h}-Tmin_{h}}\right)\cdot\bm{\zeta_{m}}$
(12)

A value greater than 1 is considered an all rain event.

###  2.2 Interception

Portions of the HRU with canopy cover intercept precipitation. Different
canopy cover density values are used for the summer ($\bm{ds}_{h}$) and winter
($\bm{dw}_{h}$). Only when the precipitation amount exceeds the canopy storage
capacity ($Sc$) does it reach the ground as throughfall ($Pt_{h}$). For
example the net precipitation reaching the ground during summer is:

$Pn_{h}=Pt_{h}\cdot\bm{ds}_{h}+\left(1.0-\bm{ds}_{h}\right)\times\mathit{P}_{h}$
(13)

where the throughfall $Pt_{h}$ is computed as:

$Pt_{h}=\begin{cases}P_{h}-\frac{Sc_{h}}{\bm{A}_{h}\cdot\bm{ds}_{h}}&:P_{h}>\frac{Sc_{h}}{\bm{A}_{h}\cdot\bm{ds}_{h}}\\\
0.0&:\text{otherwise}\end{cases}$ (14)

where the canopy storage capacity $Sc_{h}$ is:

$Sc_{h}=\left(\bm{Crs}_{h}-Sint_{h}\right)\cdot\bm{A}_{h}\cdot\bm{ds}_{h}$
(15)

Parameters $\bm{ds}_{h}$ and $\bm{Crw}_{h}$ are used during winter. A similar
process is modeled for snow. $Sint_{h}$ is a model state holding intercepted
precipitation. $\bm{Crs}_{h}$ and $\bm{Crw}_{h}$ are the summer and winter
rain interception storage capacity for the major vegetation type for each HRU.
$\bm{Cs}_{h}$ is used for snow. Intercepted rain, (or snow) evaporates (or
sublimates) at a free-water surface rate.

###  2.3 Evapotranspiration

The Jensen-Haise formulation [7] is used to compute the potential
evapotranspiration:

$ET_{h}=\bm{jc}\cdot\left(T_{h}-jc_{h}\right)\cdot\frac{Rsw_{h}}{2.54\cdot\left(597.3-\left(0.5653\cdot
T_{h}\right)\right)}$ (16)

where the coefficient $jc_{h}$ is approximated to be:

$jc_{h}=22-\frac{\bm{e_{h}}}{1000}$ (17)

As figure 2 shows, evapotranspiration occurs in multiple reservoirs.

###  2.4 Snowpack

Snowpack dynamics are modeled using water and energy balances for each HRU.
The snowpack is abstracted as 2 layers. The energy exchanges that occur
between the snowpack and the atmosphere are radiative, conductive and
convective. The energy balance at the snow-atmosphere interface is computed
twice per day, for both the day and night periods as:

$\frac{dE}{dt}=I_{h}+Qv_{h}+R_{h}$ (18)

$E$

is the surface energy of the snowpack

The net short wave radiation $R_{h}$ in (18) is computed from $Rsw_{h}$ (9)
after accounting for the reflected portion of the radiation using the surface
albedo $\alpha$ and that limited by the vegetative transmission coefficient
parameter $\bm{\psi_{h}}$.

$R_{h}=\left(1-\alpha_{h}\right)\cdot\bm{\psi_{h}}\cdot Rsw_{h}$ (19)

Where $\alpha$ is obtained from a non-linear curve that is a function of the
snow dynamics phase and time since last new snow in days.

The net incoming longwave radiation $I_{h}$ originates from the atmosphere and
land cover:

$I_{h}=\bm{cw_{h}}\cdot
Ip_{h}+\left(1-\bm{cw_{h}}\right)\cdot\bm{\epsilon}\cdot Ip_{h}$ (20)

where $\bm{\epsilon}$ is the emissivity, a precipitation situation-dependent
parameter (between 0.0 and 1.0 dependent on storm type), and $Ip_{h}$ is the
perfect black-body emission:

$Ip_{h}=5.85\cdot 10^{-8}\cdot\left(T_{h}+273.16\right)^{4}$ (21)

The combined convection and latent heat flux from condensation is modeled as
linear function of temperature.

$Qv_{h}=\bm{\omega}\cdot T_{h}$ (22)

When (18) is negative, heat flow occurs by conduction only (ex: potentially
refreeze snow at isothermal) and depends on the snowpack density and layer
temperature gradient:

$Qc_{h}=2\cdot\left[0.5\cdot\sigma\sqrt{\frac{0.0077\cdot\rho_{h}^{2}\cdot\Delta
t}{\sigma\cdot\rho_{h}}}\right]\cdot\left(T_{h}-Tpk_{h}\right)$ (23)

where $\sigma$ is the specific heat of ice (in Celsius
$Cal.m^{-3}.^{o}C^{-1}$). A positive energy exchange translates into snowmelt
at the surface of the snowpack, that transports heat $\Delta E$ into the lower
snowpack by mass transfer. The potential melt from such heat is:

$M_{h}=\frac{\Delta E}{203.2}\cdot sca_{h}$ (24)

where $sca_{h}=fsca_{h}\cdot\bm{A}_{h}$ is the HRU snow covered area.

The energy state of the lower layer is represented as a heat deficit. The
snowpack heat deficit ($H$) for each HRU represents the amount of heat
necessary to bring the snowpack to isothermal 0 oC. Snowmelt occurs only when
the heat deficit reaches zero and the snowpack free-water storage is exceeded.
Otherwise, potential melt is “refrozen" by the decrease in the heat deficit.
Storage within the snowpack is tracked in two states: ice (solid) and free
water (liquid), the sum of which is termed Snow Water Equivalent ($SWE$), and
quantifies the volume of water obtained if all snow is melted. Melt decreases
the ice storage of the snowpack, and once the amount of free water surpasses
the maximum capacity, water exits the snowpack. Sublimation loss from the
snowpack occurs according to the following equation:

$B_{h}=\bm{\xi}\cdot ET_{h}\cdot fsca_{h}$ (25)

where $ET_{h}$ is the potential evapotranspiration described in 2.3 The
process lowers the non-isothermal heat deficit by: $Tpk_{h}\cdot ET_{h}\cdot
1.27$ calories.

The evolution of the snow depth ($D$) is approximated from the ordinary
differential equation [7]:

$\frac{dD_{h}}{dt}+\bm{\tau}\cdot
D_{h}=\frac{Pns_{h}}{\bm{\rho_{init}}}+\left(\bm{\tau}\cdot\frac{SWE_{h}+Pns_{h}}{\bm{\rho_{max}}}\right)$
(26)

and is calculated as:

$\Delta D_{h}=\Delta
t\cdot\left(\frac{Pns_{h}}{\bm{\rho_{init}}}+\bm{\tau}\cdot\left(\frac{SWE_{h}+Pns_{h}}{\bm{\rho_{max}}}-D_{h}\right)\right)$
(27)

Where $\bm{\rho_{init}}$ and $\bm{\rho_{max}}$ are parameters for initial and
maximum density of new snowfall, respectively. $Pns_{h}$ is the net new
snowfall after interception. The snowpack density is then:

$\rho_{h}=\frac{SWE_{h}}{D_{h}}$ (28)

The snow-covered area is determined by a multi-modal depletion curve [7] that
describes the evolution of the fractional snow-covered area ($fsca$) in
function of the fraction of maximum $SWE$. Parameter $\bm{Acurve}$ specifies
the 11 values of the curve. The depletion curve is used when HRU $SWE$ is less
than the maximum $SWE$ corresponding to total snow cover specified by
parameter $\bm{SWEmax}$.

Water available at the soil surface (from melt and rain) proceed to both
infiltrating the soil reservoir and filling the impervious zone reservoir.

###  2.5 Impervious Zone Reservoir

The impervious zone reservoir ($Simp$) constitutes the fraction of the HRU
specified by the parameter $\bm{fi}$ with no soil-infiltration capacity. It
receives $\bm{f_{i}}$ fraction of the total water from snow melt ($M$) and
rain throughfall ($Pnr$) for each HRU and has maximum retention capacity of
$\bm{Simax}$, above which water flows directly as surface runoff $Fsri$
described in 2.7. The process is described in the following discretized
dynamic equation:

$Simp=Simp_{t-1}+(Pnr_{h}+M_{h})\cdot\bm{fi_{h}}-Fsri_{h}-ei_{h}$ (29)

where ei is the water depleted by evaporation:

$ei_{h}=ETa_{h}\cdot(1-fsca_{h})\cdot\bm{fi_{h}}$ (30)

where $ETa_{h}$ is the available potential for evapotranspiration remaining in
the system after accounting for evapotranspiration and sublimation losses that
already occured (interception,..).

###  2.6 Soil Zone Reservoir

The soil zone reservoir ($Ssz$) represents the active soil profile and its
depth is considered to be the average rooting depth of the predominant
vegetation on the HRU. It receives part of the water from snow melt and rain.
It can hold up to the parameter $\bm{Sszmax}$. The upper layer is termed the
recharge zone ($Ssre$) and the lower layer is termed the lower zone ($Ssl$).
While water is lost through both evaporation and transpiration from the upper
layer, it is only lost through transpiration from the lower layer. Evaporation
is limited by both how much potential of $ET$ is left in $ETa_{h}$ for the day
and by the availability of water. Water needs to fill the recharge layer
before it can proceed to the lower layer. The soil water content represented
by $Ssz$ is also defined as the “soil moisture". A maximum of
$\bm{Fzgwmax_{h}}$ percolates from $Ssz$ to $Sgw$, and the remainder to $Sss$
as excess.

###  2.7 Surface Runoff

Both the pervious soil zone and the impervious reservoirs contribute to
surface runoff $Fsr$ as $Fsrp$ and $Fsri$ respectively.

$Fsr_{h}=(1-\bm{fi_{h}})\cdot Fsrp_{h}+\bm{fi_{h}}\cdot Fsri_{h}$ (31)

The pervious runoff is modeled using the following non-linear parametric
equation in function of soil moisture, net rain precipitation and snowmelt
[7]:

$Fsrp_{h}=max\left[\bm{{\alpha}1_{h}}\cdot 10^{\bm{{\theta}1_{h}}\cdot
Ssz\cdot(0.5\cdot Pn_{h})},\bm{Asr}\right]\cdot\left(Pnr_{h}+M_{h}\right)$
(32)

Where $Pn_{h}$ is net precipitation (13), $\bm{{\alpha}1_{h}}$,
$\bm{{\theta}1_{h}}$ and $\bm{Asr}$ are model parameters. The impervious
runoff is equal to the amount of water exceeding $\bm{Simax}$ in (29).

###  2.8 Subsurface Reservoir

The subsurface reservoir models the relatively rapid water movement of the
unsaturated soil zones towards a streamflow channel. Such behavior is
highlighted typically during and after a rainfall or snowmelt event. A
reservoir routing system is used to model the flow from the subsurface
reservoir. The flow $Fss$ is the solution that satisfies both the mass
continuity equation:

$Fss_{s,t}=Fss_{s,t-1}-\frac{d\left(Sss\right)}{dt}$ (33)

and the empirical quadratic relation:

$Fss_{s,t}=\bm{\alpha}3_{s}\cdot Sss_{s}+\bm{\beta}3_{s}\cdot Sss^{2}_{s}$
(34)

where $\bm{\alpha 2}_{s}$ and $\bm{\bm{\beta}2_{s}}$ are routing coefficients.

Another $Sss$ discharge route transfers water to the $Sgw$ reservoir through
the parametric power equation:

$Fss\\_gw_{s}=\bm{\alpha
2_{s}}\cdot\left[\frac{Sss_{s}}{\bm{smax_{s}}}\right]^{\bm{\theta 2_{s}}}$
(35)

where $\bm{\alpha 2_{s}}$, $\bm{smax}$ and $\bm{\theta 2_{s}}$ are model
parameters.

###  2.9 Ground-Water Reservoir

Once the soil zone reservoir is saturated, excess water starts filling the
ground-water ($Sgw$) reservoir at a maximum daily recharge rate
$\bm{Fzgwmax_{h}}$, a parameter defined for each HRU.

$Fsz\\_gw=min(\bm{Fzgwmax_{h}},Ssz\\_excess)$ (36)

Water leaving the ground-water reservoir either flows laterally ($Fgw$)
contributing to the baseflow portion of the streamflow or is lost from the
system via the groundwater sink, both at a linear rate:

$Fgw_{h}=\bm{\alpha 4_{g}}\cdot Sgw_{g}$ (37) $Fgsnk_{h}=\bm{\alpha
5_{g}}\cdot Sgw_{g}$ (38)

###  2.10 Model Output: Streamflow

Basin streamflow, which consititues the model output, is computed as an area
weighted average of all the reservoir flow outputs described above:

$F_{basin}=\sum_{h=1}^{N_{h}}\frac{Fsr_{h}}{A_{h}}+\sum_{s=1}^{N_{s}}\frac{Fss_{s}}{A_{s}}+\sum_{g=1}^{N_{g}}\frac{Fgw_{g}}{A_{g}}$
(39)

Improving the daily forecast accuracy of $F_{basin}$ is the ultimate objective
of this report.

##  3 Filtering Techniques

Models are simplifications of physical processes and are inherently imperfect
representations of reality. Combining related measurements with modeled
estimates through data assimilation would produce more accurate outcomes.
Assimilation approaches were introduced in oceanography and meteorology [8, 9]
with different schemes used such as the variational methods [10] and the
filtering techniques including nudging [11], optimal interpolation [12] and
Kalman filtering with its variants that are presented below.

Bayesian Filtering

Bayesian inference [13] is the underlining principle for all the data
assimilation approaches we will describe in this section. They stem from the
two following probability rules: marginalization (sum) rule: Given
$f_{xy}(x,y)$, we define:

$f_{x}(x):=\sum_{y\epsilon Y}f_{xy}(x,y)$ (40)

and the conditioning, or product rule defined as:

$f_{x|y}(x,y):=\frac{f_{xy}(x,y)}{f_{y}(y)}$ (41)

for $f_{y}(y)\neq 0$

We can formulate the following forecast and update system in the Bayesian
framework, where $x$ represents the model state, $z$ represent the
measurements and $k$ is the discrete time step:

Prior update (or forecast):

$f(x(k)|z(1:k-1))=\int f(x(k)|x(k-1))\cdot f(x(k-1)|z(1:k-1))dx(k-1)$ (42)

where the first item in the integral is the process model and the second item
is the previous time step measurement update result.

Measurement update (or update):

$f(x(k)|z(1:k))=\frac{f((z(k)|x(k))\cdot f(x(k)|z(1:k-1)}{\int
f(z(k)|x(k))f(x(k)|z(1:k-1))dx(k)}$ (43)

where the first and second items of the numerator are the measurement model
and the prior (42), respectively, and the denominator is the normalization
$f(z(k)|z(1:k-1))$.

Kalman Filter It can be shown [14, 15, 16] that given the system of equations
with a linear process model A and independent Gaussian distributions for
errors $v(k-1)$ and $w(k)$ and prior: assuming initialization:
$\hat{x}_{a}(0)=x_{0},P_{a}(0)=P_{0}$

$x(k):=A(k)\cdot x(k-1)+u(k-1)+v(k-1)\newline $ (44) $y(k):=H(k)x(k)+w(k)$
(45)

the analytical solution to the Bayesian state estimation problem in (42) and
(43) is the Kalman Filter (KF) and can be written in a recursive form while
keeping track of the distributions’ mean and covariance as:
Prior update/Forecast step:

$\hat{x}_{f}(k)=A(k-1)\hat{x}_{a}(k-1)+u(k-1)$ (46)
$P_{f}(k)=A(k-1)P_{a}(k-1)A^{T}(k-1)+Q(k-1)$ (47)

A Posteriori update/Analysis step:

$\hat{x}_{a}(k)=\hat{x}_{f}(k)+P_{a}(k)H^{T}(k)R^{-1}(k)\left(z(k)-H(k)\hat{x}_{f}(k)\right)$
(48) $P_{a}(k)=\left(P_{f}^{-1}(k)+H^{T}(k)R^{-1}(k)H(k)\right)^{-1}$ (49)

without the need to store all previous observations and states. It is commonly
referred to $P_{a}(k)H^{T}(k)R^{-1}(k)$ in (48) as the Kalman gain $K(k)$. The
maximum a posteriori estimate (MAP) of the normally distributed analysis
coincides with the mean: $x^{MAP}_{a}(k)=\hat{x}_{a}(k)$ and thus is chosen as
the best estimate of the system state at time k.

The basic Kalman Filter (KF) algorithm that is restricted for only linear
models consists of two steps, a predict step followed by an update step. The
KF and its variants have numerous applications in the physical sciences,
primarily in the guidance [17], navigation and control of vehicles [18], time-
series analysis [19] and robotics [20].

Extended Kalman Filter

The Extended Kalman Filter (EKF) is an extension to the KF for non-linear
systems. The approach is to linearize the non-linear process model around the
mean using first order Taylor expansion. It thus becomes an approximation of
the exact solution. The Jacobian is then used to advance the error covariance
matrix in time. It requires the model to be continuously differentiable and is
shown to have bad performance when the system exhibits strong non-linearities.
For instance, the EKF was successfully used to integrate GPS-derived flow
velocity measurements from drifters to improve the non-linear hydrodynamics
state estimates in a controlled channel pilot experiment [21].

Unscented Kalman Filter

The Unscented Kalman Filter (UKF), also referred to as Sigma-Point KF, tries
to solve the linearity problem using another approach. It works by first
deterministically sampling using an algorithm sigma points from the
distribution and assigning to each point a weight. At every time step the
points are individually propagated through the non-linear model, and finally
the resulting mean and covariance is computed from the forecasted sigma
points. Unlike the EKF, not only the mean is advanced in time, but also many
points sampled from the distribution. The method requires N = 2r + 1 sigma
points where r is the system dimension. [22] shows that the UKF, with the same
complexity, is a better filter than the EKF for non-linear models capturing
sigma points’ mean and covariance accurately up to the third order (in Taylor
series expansion terms) while the EKF only to the first. Using UKF instead of
EKF halfed the RMSE in estimating vehicle and wheel angular velocity of the
anti-lock braking system (ABS) in [23].

Ensemble Kalman Filter

The EnKF represents the distribution as a random sample of points called
ensemble members. In the forecast step, each ensemble member is integrated
forward in time with an additional process noise. The forecasted error
covariance matrix is approximated at any time step by the ensemble sample
covariance. The error covariance matrix is thus implicitly propagated saving
the high computational requirements associated with its storage and forward
integration such as in KF and EKF. The update stage consists of updating every
ensemble member using sample covariance as an approximation of the true
covariance. The ensemble members are never resampled from the distribution
thus preserving the ensemble skewness, kurtosis, clustering, etc. In the EnKF
update, the updated ensemble is obtained by shifting and re-scaling the
forecast ensemble.

Particle Filter

The Sequential Importance Resampling (SIR) Particle Filter (PF), samples
points called particles from the distribution to generate system states
realizations. These particles are assigned weights that represent their
likelihood. The sampling is done randomly like the EnKF, but not necessarily
from normal distributions. On update, the particle weights are adjusted.
Higher weights are assigned to those that are supported by sensor data and
vice versa. The Particle Filter suffers from a problem with sample degeneracy
and impoverishment, in that some scenarios end up with very few high-weight
particles while the other majority of particles have almost zero weights.
Another disadvantage is that the number of particles required for good
performance scales exponentially with the model dimension. The method does not
require Gaussian distributions nor model linearity [24]. It thus works well
with distributions that have complex shape and multiple peaks. Example
applications include vehicle localization [25], indoor occupant positioning
using a Radio Frequency Identity RFID and receivers [26], etc.

Filter Selection

The Kalman Filter in its basic form cannot be applied to our modeling system
because the latter is non-linear. The model is also not continuously
differentiable, which eliminates the option of using the Extended Kalman
Filter. For example, different physical conditions can lead to the use of
totally different process equations. For our application, using either of UKF,
EnKF or PF is theoretically valid. According to [27], UKF would require $2337$
sigma points to fully represent the mean and variance of the system. In [28],
the UKF was used in a similar application as ours but using algorithms to
reduce the number of sigma points to $n+2$. Multiple studies in the literature
pertaining to geophysics, oceanography and hydrology have shown that the EnKF
performs well even when the ensemble size is orders of magnitude smaller than
the system dimension [29, 30]. EnKF outperformed PF in the assimilation of a
conceptual rainfall-runoff in [31] and the justification was that it is “less
sensitive to misspecification of the model and uncertainties", and in [32], it
had similar performance. Given the above specified reasons and the limitation
in resources available that will restrict the number of model realizations, we
choose the EnKF.

###  3.1 EnKF & State-space

The Ensemble Kalman Filter [33] was introduced in [34] as an alternative to
EKF to overcome specific difficulties with nonlinear state evolution models,
including non-differentiability of the model and closure problems. As
described above, EnKF uses Monte Carlo (or ensemble integrations) and has the
same standard update equations of the KF, with state mean and covariance
replaced by an estimated ensemble sample mean and covariance respectively
[35]. The EnKF and more generally statistical inversion schemes were greatly
advanced in [36], mainly in the context of Tomography. A framework for
sampling a smooth prior is presented, especially when structural information
about the ensemble (correlations between pixels) is known apriori. Such
techniques can also be used in the perturbations of spatial data such as
climatologic inputs, model error etc. For instance these techniques were used
to generate smooth initial spatial model states ensemble in [37] without
shocks. [38] shows that using a non-trivial state covariance matrix in the
state noise yields much better performance than a stationary one.

With the Kalman Smoother, state estimation was applied in Electrical Impedance
Tomography (EIT) to reconstruct rapidly time-varying cross-sectional body
scans inferred from boundary voltage and current measurements [39], [40].
Similarly, Electrical Resistance Tomography (ERT) measurements were used to
infer water saturation levels in soil using EKF assimilation of a weakly non-
linear model [41]. Time-dependent noise was used in [42] infer gas temperature
from resistance measurements.

[37] provided a framework for assimilating readily available GPS measurements
of vehicle speed and velocity into a traffic model [43] using the EnKF as part
of the the Mobile Millennium traffic-monitoring system[44]. Results of an
unprecedented scale experiment [37], where GPS speed and position data from
100-vehicles were assimilated into the traffic model, showed improved
estimates of the traffic state even with very low percentage of GPS-equipped
vehicles participating.

We now present the conceptual algorithm of the EnKF: Let s be the system
dimension. Hydrologic models typically take inputs. Let i be the number of
inputs and o be the number of observations. Let $\bm{x}$ be the state vector
and $\mathcal{M}$ the non-linear model. The discrete time k forward
integration equation can be represented as:

$\bm{x}^{k}=\mathcal{M}(\bm{x}^{k-1},\bm{b}^{k})$ (50)

Where the model observed states are:

$\bm{y}^{k}=\bm{H}\bm{x}^{k}$ (51)

The localized EnKF assimilation scheme can be summarized by the algorithm
below:

  1. 1.

Initialization: Draw N ensemble realizations $\bm{x}^{0}_{a}\left(e\right)$
(with $e\in[1,...,N]$) from a process with a mean $\bar{\bm{x}}^{0}_{a}$ and
covariance $\bm{C}^{0}_{a}$

  2. 2.

Forecast: At each time step n, update each of the N ensemble members according
to PRMS forward simulation model. Then update the ensemble mean and covariance
according to:

$\bm{x}^{k}_{f}(e)=\mathcal{M}[\bm{x}^{k-1}_{f}(e),\bm{b}^{k}(e)]+\bm{\eta}^{k}(e)$
(52) $\bar{\bm{x}}^{k}_{f}=\frac{1}{N}\sum_{e=1}^{N}\bm{x}^{k}_{f}(e)$
$\bm{C}^{k}_{ens,f}=\frac{1}{N-1}\sum_{e=1}^{N}\left(\bm{x}^{k}_{f}(e)-\bar{\bm{x}}^{k}_{f}\right)\left(\bm{x}^{k}_{f}(e)-\bar{\bm{x}}^{k}_{f}\right)^{T}$

  3. 3.

Analysis: Assuming a measurement is obtained at time n, localize sample
covariance, compute the Kalman gain, and update the network forecast:

$\bm{C}^{k}_{ens,f,L}=\mathcal{L}\circ\bm{C}^{k}_{ens,f}$
$\bm{K}_{ens}^{k}=\bm{C}^{k}_{ens,f,L}\left(\bm{H}^{k}\right)^{T}\left(\bm{H}^{k}\bm{C}^{k}_{ens,f,L}\left(\bm{H}^{k}\right)^{T}+\bm{C}^{k}_{obs}\right)^{-1}$
$\bm{x}^{k}_{a}(e)=\bm{x}^{k}_{f}(e)+\bm{K}_{ens}^{k}\left(\bm{y}_{obs}^{k}-\bm{H}^{k}\bm{x}^{k}_{f}(e)+\bm{\chi}^{k}(e)\right)$
(53)

  4. 4.

Return to 2.

where:

  1. $s^{k}_{f}$

is the ensemble of states forecast at time k, [s x N]

  2. $\mathcal{M}$

is the model operator (PRMS)

  3. $\bm{b}^{k}$

is the perturbed ensemble input, [i x N]

  4. $\bm{\eta}^{k}$

is the white process noise with covariance Q, [s x N]

  5. $\bar{\bm{s}}^{k}_{f}$

is the ensemble sample mean of the states forecast at time k, [s]

  6. $\bm{C}^{k}_{ens,f}$

is the ensemble states forecast sample covariance at time k, [s x s]

  7. $\bm{H}^{k}$

is the observation matrix at time k, [o x s]

  8. $\bm{K}_{ens}^{k}$

is the Kalman gain, [s x o]

  9. $\bm{R}^{k}_{obs}$

is the measurement error covariance, [o x o]

  10. $x^{k}_{a}$

is the ensemble states after analysis at time k, [s x N]

  11. $\bm{y}_{obs}^{k}$

is the vector of observations, [o x 1]

  12. $\bm{\chi}^{k}$

is the measurement perturbations with 0 mean & covariance $\bm{R}^{n}_{obs}$,
[o, N]

One can refer to [45] for more details and for practical implementation.

Model Constraints

The model described above has physical constraints for inputs, states,
observations and parameters as described in section 4. During the perturbation
of such quantities with normally distributed noise and after the update step
of the EnKF, invalid scenarios are likely to arise. To circumvent this issue,
a hard boundary check is implemented after each analysis step. Any values
outside the range are set to the corresponding maximum or minimum of that
state.

Major determinants of EnKF performance are the ensemble size and the error
statistics.

Ensemble Size
With an infinite ensemble size, the EnKF approximation will reach the true KF
solution. We choose an ensemble size of 100 for practical reasons detailed in
Section 4.

Additionally due to the small ensemble size compared to the state space size,
it might be necessary to use localization when using the EnKF when spurious
correlations between physically uncorrelated states (ex: geographically
distant states) affect the filter performance [46, 30]. Distance based
localization weights using linear, exponential or Gaussian decorrelation
functions are typically used to attenuate the correlations between distant
states and observations. For this study, we do not use any localization as the
impact on the outflow is assumed to be negligible compared to the scheme
introduced in section 5.2.

Inputs and Model Error Statistics
The additive noise $\bm{\eta}^{k}$ in the integration equation at time $k$
(52) is modeled as a fraction of the mean of the ensemble. This scheme was
chosen because the majority of states are storage states and have zero lower
bound. Simultaneously forcing such boundaries as described above and adding
daily zero mean normal noise would introduce bias to the model. The downside
of such error modeling is that the model is assumed to have perfect estimates
for states with zero means.

$\bm{\eta}_{k}\sim\mathcal{N}(\mu=0,\,\sigma=\alpha\bar{\bm{x}}^{k}_{f})$ (54)

Inputs (and potentially parameters) perturbation propagates model errors that
are consistent with the physics of the model, compared to only a simple daily
addition of white noise to the ensemble. Inputs and observations error
perturbations are modeled as independent normal distributions with standard
deviations computed as the following:

For input precipitation:

$\sigma_{P}=0.4\cdot P$ (55) $Pens\sim\mathcal{N}(\mu=P,\,\sigma^{2}=(0.4\cdot
P)^{2})$ (56)

For simulated SWE measurement:

$\sigma_{SWE}=0.1\cdot SWE$ (57)

The standard deviation for temperature is assumed to be a constant:

$\sigma_{T}=2^{o}C$ (58) $Tens\sim T+\mathcal{N}(\mu=0,\,\sigma^{2}=2^{2})$
(59)

Variance inflation is required when the model forecast error is known to be
underestimated. Without variance inflation, the filter tends to diverge and
the observation fails to influence the model forecast because of the
underestimated spread in the ensemble state forecast due to previous update
events. Multiple studies implement variance inflation differently [46] and
some use adaptive inflation that is pre-computed at analysis step as a
function of the error between the forecasted state and the observation [47].
In this study we will use a post-analysis inflation procedure with constant
$\alpha_{i}=0.9$ similar to [48].

$\bm{x}^{\prime}_{a,infl}=\left(1-\alpha_{i}\right)\bm{x}^{\prime}_{a}+\alpha_{i}\cdot\bm{x}^{\prime}_{f}$
(60)

where the prime $\bm{x}^{\prime}$ denotes the ensemble anomaly
$(\bm{x}-\bm{\bar{x}})$.

##  4 Model Features

As can be seen by equations in section 1, the model is extremely nonlinear
(ex: equations 21, 32, and 35…) and there is frequent use of thresholds to
determine what situational equations are appropriate to use (ex: equation 12).
Those threshold functions also make the resulting system not continuously
differentiable.

There are some explicit constraints for model inputs, states and parameters
shown in tables 1, 2 and S.1. These constraints must be satisfied for the
model to run and to get a physically realistic behavior.

The model is spatially distributed and runs independently on each sub-basin.
Each sub-basin is divided into HRUs (for surface), SSRs and GWRs (for sub-
surface) spatial units. The basin-level model dimension depends on the number
of these sub-units.

To quantify the model dimension in this study, we focus on the region of
interest, which is the East Branch of the North Fork of the Feather River sub-
basin, highlighted in Fig. 1. 111 HRUs have been apriori selected and they are
spatially co-located with the ground-layer subsurface reservoirs (SSRs) and
groundwater reservoirs (GWRs). We will thus call the 111 regions HRUs in the
remainder of this study. This implies that $A_{h}=A_{s}=A_{g}$ in (39). Each
HRU has a total of 12 states, shown in Table 2. In our case, model parameters
can be basin-wide (i.e. fixed and constant value for all HRUs) or different
for each HRU. Some parameters are constant and others change monthly or
seasonally.

The dimension of the model in the region of interest is then 1332.

Data
In this study, model inputs are generated by various cooperators and are
obtained from Pacific Gas & Electric (PG&E). They were estimated using the
Parameter-Elevation Relationships on Independent Slopes Model (PRISM)’s
historical monthly maps [49] and daily precipitation gauges. More details are
available in [6].

SWE measurements are simulated by averaging for each HRU from the 90 meter
resolution daily historical product in [50].

Other measurements available are the streamflow $F\\_basin$ in equation 39
downstream of the basin to check the output performance. The model parameters
used are those obtained from the previous calibration study of the model and
can be found in [6].

The average water year 2006 was chosen for simulation. Additional water years
(2011, 2014) are simulated and their results can be found in the Supplement.
2014 was a dry year while 2011 was a wet water year [51]. We first run the
model with the calibrated parameters [6] which are presently used
operationally without any assimilation to serve as a reference. Next, we
proceed to assimilate the SWE state of the model in Section 5.1. The
experiment results motivate a new joint-state parameter assimilation scheme
with feedback that is explained and simulated in Section 5.2.

Performance Metrics

The Root Mean Squared Error is used to compare the stream flow model output
from the assimilation experiments with the measured stream flow.

$RMSE=\frac{1}{T}\sum_{t=1}^{T}{\sqrt{\left(Fbasin_{meas,t}-Fbasin_{sim,t}\right)^{2}}}$
(61)

The RMSE of all experiments conducted are summarized in tables 4 and 5.

For ensemble size, we use the maximum possible ensemble number of 100 that is
practical to run on a typical user accessible computer. The limiting factor
for increasing the ensemble size in our case is the time needed to read and
write state values to an ensemble of PRMS files. The PRMS used is available as
a pre-compiled binaries of a Fortran code and is setup to read/write state
data from/to files. To improve the latency of these procedures, the files were
stored on RAM instead of disk and the ensemble model forecasts were generated
in parallel in a multi-threaded wrapper. Simulation time for one water year
with daily assimilation is around one hour.

##  5 Results

###  5.1 SWE State Assimilation

Figure 4 shows that the analysis SWE follows closely the “measured" SWE. This
is expected since inflation was used to increase the uncertainty of the SWE
state forecast, and thus daily observations of SWE have a relatively smaller
uncertainty during the majority of the simulation days. The remaining water
year graphs are available in the Supplement. Since the SWE observation
uncertainty was modeled as a percentage of the observation value, it peaks
during the peak snow season and is illustrated by the grey ensemble spread in
the plot.

Figure 3: Water year 2006 basin streamflow output. Figure 4: Water year 2006
basin average SWE.

Figure 3 shows a deterioration in the performance in terms of runoff output.
This is consistent for all simulation years, as indicated in Table 4 and in
Figures S.1 and S.4 in the Supplement. The main explanation for such outcome
is that the model parameters used are those that were calibrated based on
inputs and output only. The former are temperature and precipitation while the
latter is the measured runoff. No SWE data were involved in the calibration,
which implies that updating the model SWE to a more accurate estimate will not
necessarily improve the model output. We should also note that the SWE used as
observation is not directly measured, but a modeled product. However, we
assume this is not the reason of deterioration. Furthermore, the runoff output
ensemble spread is not as large as desired, since the measured runoff is not
within the ensemble for the majority of the days.

Given that the SWE ensemble spread is relatively large, that the parameters
used were calibrated without SWE knowledge, we conclude that it would be
advantageous to perturb the model parameters and estimate them as well in a
joint state-parameters assimilation scheme. Moreover a feedback consisting of
the previous-day measured runoff would be necessary to estimate those
parameters mostly involved in runoff generation dynamics.

###  5.2 SWE and Runoff Joint State-Parameter Assimilation

####  5.2.1 New State Space

Results of the previous experiment imply that parameters need to be altered.
Model parameters can be appended to the state vector and updated online on
assimilation events in what is termed a “Joint state-parameter estimation". In
such method, there is no need to apriori train the model to find the optimal
parameters which might not be time-independent. Instead, the range of each
parameter must be known. Parameters are indirectly updated through the cross-
correlation between the parameter and the “state" being observed in the Kalman
gain matrix.

Joint parameter-state assimilation has been used in many applications with
complex models and spans different fields to name a few: assimilation of
measured blood concentration into a metabolic model to update the numerous
model parameters as well as other unobserved states [52], assimilation of
measured velocity fields from GPS-equipped drifters into a shallow water model
to update model parameters [53], and assimilation of measured displacements
into gas storage geomechanical model to reduce uncertainty of some model
parameters [54]. Results in [41] showed that simultaneous parameter estimation
is necessary when the model perturbed by daily additive white noise is far
from the truth.

Parameters involved in model dynamics from model states to streamflow output
require streamflow measurements to be updated. Parameters that are included in
the state space were chosen based on the sensitivity results obtained from the
study in [55]. We choose the parameters that score high on daily streamflow
statistics of mean, CV and AR1 indicating sensitivity of parameters to each
process. Processes chosen were snowmelt (assimilating SWE) and runoff
(assimilating runoff), with emphasis on runoff. Chosen paramters include 9
parameters per HRU and 4 global parameters and are shown in Table 3. The state
vector is augmented with these parameters. Exponential parameters are excluded
(ex: smidx_exp) to maintain approximate Gaussian ensemble distributions, a
requirement of the EnKF.

Symbol | PRMS symbol | Quantity | Type | Range | Units | Description
---|---|---|---|---|---|---
Tms | tmax_allsnow | 1 | constant | -10, 40 |  oF | temperature below which P is all snow
Tmr | tmax_allrain | 1 | monthly | -8, 60 |  oF | temperature above which P is all rain
$\alpha$1 | smidx_coef | per HRU | constant | 0.001, 0.06 | - | linear surface runoff (Fsr) coeficient
Asr | carea_max | per HRU | constant | 0.0, 1.0 | - | max areal fraction contributing to surface runoff Fsr
$\beta$ | dday_intcp | 1 | monthly | -60.0, 10.0 |  oday | intercept of degree-day equation
jc | jh_coef | 1 | monthly | 0.005, 0.06 |  /oF | coef. used in ETp
Sszmax | soilmoist_max | per HRU | constant | 0.001, 60.0 | inches | maximum soilzone water holding capacity
$\alpha$4 | gwflow_coef | per GWR | constant | 0.001, 0.5 | -/day | linear coefficient routing Sgw to streamflow
Fzgwmax | soil2gw_max | per HRU | constant | 0.0, 5.0 | inches | max soil excess water routed to gwStr
$\alpha$5 | gwsink_coef | per GWR | constant | 0.0, 1.0 | -/day | linear coefficient for groundwater sink
$\alpha$2 | ssr2gw_rate | per SSR | constant | 0.05, 0.8 | -/day | linear coefficient routing Sss to Sgw
$\alpha$3 | ssrcoef_sq | per SSR | constant | 0.0, 1.0 | - | coeficient routing Sss to streamflow
$\beta$3 | ssrcoef_lin | per SSR | constant | 0.0, 1.0 | -/day | linear coeficient routing Sss to streamflow

Table 3: Perturbed and updated parameters. Figure 5: Water year 2006 basin
average streamflow.

Initial parameter perturbations are modeled proportionally to their range,
such that:

$\sigma_{p}=0.25\cdot|p_{max}-p_{min}|$ (62)

Parameters are inflated post-analysis similarly to states as described in
(60). Parameters do not change during model integration, thus when the
parameter ensemble nears collapse - ie. when the standard deviation becomes
less than $\sigma_{target}$, inflation is performed by re-scaling the
parameter ensemble perturbations to $\sigma_{target}$:

$\sigma_{target}=0.25\cdot\sigma_{p}$ (63)

where $\sigma_{p}$ is the standard deviation of the initial parameter
perturbation.

Next, we additionally augment the state-space with the previous-day updated
runoff, so that it can be updated on the next analysis day by the previous-day
measured runoff. The current day simulated streamflow is also appended.
Finally, the dimension of the augmented state-space becomes 2337.

The observation vector $y_{obs}$ is augmented by the previous-day measured
runoff and correspondingly, a new row is appended to the observation matrix
$H$ in the analysis equation (52). The previous-day measured runoff is
available in practical applications, because we are updating the SWE state
daily. The measured streamflow is assumed near perfect and perturbed with an
independent normally distributed error with a scaling standard deviation:

$\sigma_{F}=0.005\cdot Fbasin_{meas}$ (64)

###  5.3 Results

As reference, we compare the results to both the estimated streamflow with no
assimilation, which is the green plot shown in Figure 5, and an additional
estimate which is simply the previous-day measured streamflow (not
visualized), also termed autoregressive AR(1) filter with equation:

$Fbasin_{sim,t}=Fbasin_{meas,t-1}$ (65)

Results of both experiments are summarized in Tables 4 and 5 and show a
consistent improvement in simulated streamflow compared to the no assimilation
parameter-calibrated case, with up to 60% reduction in RMSE for the wet water
year 2011. Results were occasionally better than the naive previous-day AR(1)
forecast such as the average water year 2006 (28%), but worse during the dry
water year 2014. For all assimilation experiments, the RMSE of the basin-mean
SWE state was substantially reduced ($>90\%$) as expected, since it is the
state being “measured" with $10\%$ uncertainty. SWE results constitute a
validation of the EnKF analysis procedure.

Water Year | no assim. | SWE update |  $\Delta$(%) | Joint SWE/Fbasin |  $\Delta$(%) | AR(1) |  $\Delta$(%)
---|---|---|---|---|---|---|---
2006 | 1228 | 1845 | 50 | 889 | -28 | 1235 | 1
2011 | 870 | 1891 | 117 | 360 | -59 | 430 | -51
2014 | 198 | 996 | 404 | 177 | -10 | 126 | -36

Table 4: Experiments performance in terms of streamflow output RMSE

Water Year | no assim. | SWE update |  $\Delta$(%) | Joint SWE & Stream |  $\Delta$(%)
---|---|---|---|---|---
2006 | 1.08 | 0.02 | -98 | 0.04 | -96
2011 | 1.97 | 0.05 | -97 | 0.05 | -98
2014 | 0.20 | 0.01 | -93 | 0.01 | -95

Table 5: Experiments validation in terms of basin-mean SWE RMSE.

##  6 Conclusion

The EnKF data assimilation framework presented succeeded in updating the
modeled SWE during analysis step. SWE state assimilation alone does not
improve runoff in a heavily parameterized model, where the parameters have
been calibrated based on inputs/output without taking into consideration the
SWE state being updated. Furthermore, results show that deterioration in
accuracy can occur. We postulate SWE assimilation in a non-parametric model is
likely to improve streamflow forecast accuracy.

Joint state-parameter assimilation using the previous day measured stream flow
as feedback shows a substantial improvement in the accuracy of the daily
estimate of streamflow (30% reduction in RMSE for water year 2006) compared to
the calibrated no-assimlation scenario, and compared to the simple previous-
day AR(1) estimator (30% reduction in RMSE for water year 2006) during average
water years.

In reality, streamflow measurements are not as accurate as assumed in this
study. It is typically computed from the river stage measurements assuming
constant cross section area and other approximations that should be accounted
for. Streamflow measurement uncertainty, if well known, can be modeled and
accounted for in 64. Nevertheless, what this study shows is that with a near-
perfect measurement of streamflow, the EnKF framework presented improves the
accuracy of the streamflow estimate. More accurate measurement methods do
exist such as current meters or the “Acoustic Doppler Current Profiler", but
are more expensive.

Although the previous day estimator (AR(1)) was better for the dry year 2014,
it has disadvantages over using the physical model-based assimilation
framework in that it is not a reliable method, especially when the previous
streamflow measurement is not available or severely inaccurate due to sensor
failure or environmental events, whereas the EnKF framework would provide a
more robust and physically sound estimate when measurements are missing. In
fact, the only change required when previous day measurement of streamflow is
not available is to delete one row from the observation matrix. Moreover, the
EnKF with PRMS framework is preferred to the naive AR(1) method given the
forecast time window the former can potentially achieve. We suggest that a
more sophisticated fully data-driven model that considers at least input data
is required to match its potential long-term forecast accuracy such as an
artificial recurrent neural network (ex: Long short-term memory, LSTMs).

We suggest as future work to replicate the study for all years where data is
available. A longer forecast window study is also of interest to stakeholders
and would be a natural continuation of this report, where multiple days in the
future are first simulated with model inputs from weather predictions
(forecast step), after which their streamflow outputs are updated with
measurements of SWE and streamflow available at the current day (update step).
This would increase the state dimensions only by the number of days in the
prediction window, keeping the assimilation tractable. With parameters updated
by measurements of streamflow and states such as SWE, current-day measurements
of SWE (and potentially other states) should have a strong impact on the
streamflow prediction accuracy for a long prediction window.

## Acknowledgement

This work was partially supported by the Civil Systems department at UC
Berkeley, California Energy Commission through the grant “Improving
Hydrological Snow pack Forecasting for Hydropower Generation Using Intelligent
Information Systems” (EPC-14-067), Pacific Gas and Electric Co (PG&E), the
California Department of Water Resources and UC Water.

I would like to thank Prof. Steven Glaser for recommending me to pursue this
research work and to thank Prof. Alexandre M. Bayen for his advice and words
of encouragement.

Finally, special thanks to Kevin Richards (PG&E) for sharing streamflow
measurements and Feather river data pertaining to PRMS such as HRUs and
calibrated parameters.

## References

  * [1] California Department of Water Resources, “Climate Change,” 2017. [Online]. Available: http://www.water.ca.gov/climatechange/
  * [2] R. Bales, N. P. Molotch, T. H. Painter, M. D. Dettinger, R. Rice, and J. Dozier, “Mountain hydrology of the western United States,” _Water Resources Research_ , vol. 42, p. W08432, 2006.
  * [3] I. T. Stewart, D. R. Cayan, and M. D. Dettinger, “Changes in snowmelt runoff timing in western north america under a ‘business as usual’ climate change scenario,” _Climatic Change_ , vol. 62, no. 1, pp. 217–232, Jan 2004. [Online]. Available: https://doi.org/10.1023/B:CLIM.0000013702.22656.e8
  * [4] ——, “Changes toward earlier streamflow timing across western north america,” _Journal of Climate_ , vol. 18, no. 8, pp. 1136–1155, 2005. [Online]. Available: https://doi.org/10.1175/JCLI3321.1
  * [5] T. H. Painter, D. F. Berisford, J. W. Boardman, K. J. Bormann, J. S. Deems, F. Gehrke, A. Hedrick, M. Joyce, R. Laidlaw, D. Marks, C. Mattmann, B. McGurk, P. Ramirez, M. Richardson, S. M. K. Skiles, F. C. Seidel, and A. Winstral, “The Airborne Snow Observatory: Fusion of scanning lidar, imaging spectrometer, and physically-based modeling for mapping snow water equivalent and snow albedo,” _Remote Sensing of Environment_ , vol. 184, pp. 139–152, 2016. [Online]. Available: http://dx.doi.org/10.1016/j.rse.2016.06.018
  * [6] K. M. Koczot, A. E. Jeton, Bruce, McGurk, and M. D. Dettinger, “Precipitation-runoff processes in the feather river basin, northeastern california, and streamflow predictability, water years 1971-97.” [Online]. Available: https://pubs.usgs.gov/sir/2004/5202/
  * [7] S. L. Markstrom, R. S. Regan, L. E. Hay, and R. J. Viger, “Prms-iv, the precipitation-runoff modeling system, version 4.” [Online]. Available: https://pubs.usgs.gov/tm/6b7/
  * [8] R. A. Anthes, “Data assimilation and initialization of hurricane prediction models,” _Journal of the Atmospheric Sciences_ , vol. 31, no. 3, pp. 702–719, 1974. [Online]. Available: https://doi.org/10.1175/1520-0469(1974)031<0702:DAAIOH>2.0.CO;2
  * [9] F. COUNILLON and L. BERTINO, “Ensemble optimal interpolation: multivariate properties in the gulf of mexico,” _Tellus A_ , vol. 61, no. 2, pp. 296–308, 2009. [Online]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1600-0870.2008.00383.x
  * [10] F.-X. L. DIMET and O. TALAGRAND, “Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects,” _Tellus A_ , vol. 38A, no. 2, pp. 97–110, 1986. [Online]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1600-0870.1986.tb00459.x
  * [11] C. Paniconi, M. Marrocu, M. Putti, and M. Verbunt, “Newtonian nudging for a richards equation-based distributed hydrological model,” _Advances in Water Resources_ , vol. 26, no. 2, pp. 161 – 178, 2003. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0309170802000994
  * [12] A. T. Lorenzo, M. Morzfeld, W. F. Holmgren, and A. D. Cronin, “Optimal interpolation of satellite and ground data for irradiance nowcasting at city scales,” _Solar Energy_ , vol. 144, pp. 466 – 474, 2017. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0038092X17300555
  * [13] M. Bayes and M. Price, “An essay towards solving a problem in the doctrine of chances. by the late rev. mr. bayes, f. r. s. communicated by mr. price, in a letter to john canton, a. m. f. r. s.” _Philosophical Transactions (1683-1775)_ , vol. 53, pp. 370–418, 1763. [Online]. Available: http://www.jstor.org/stable/105741
  * [14] R. Kalman, “A new approach to linear filtering and prediction problems,” _Journal of Basic Engineering (ASME)_ , vol. 82D, pp. 35–45, 01 1960.
  * [15] N/A, “The kalman filter as best linear estimator,” University Lecture, 2018. [Online]. Available: https://berkeley-me233.github.io/
  * [16] J. Kaipio and E. Somersalo, _Statistical and computational inverse problems_. Springer, 2010.
  * [17] J. Fonseca, R. Oliveira, J. Abreu, E. Ferreira, and M. Machado, “Kalman filter embedded in fpga to improve tracking performance in ballistic rockets,” 04 2013, pp. 606–610.
  * [18] R. Kasper and S. Schmidt, “Sensor-data-fusion for an autonomous vehicle using a kalman-filter,” in _2008 6th International Symposium on Intelligent Systems and Informatics_ , Sep. 2008, pp. 1–5.
  * [19] T. Choudhry and H. Wu, “Forecasting the weekly time-varying beta of uk firms: Garch models vs. kalman filter method,” _The European Journal of Finance_ , vol. 15, no. 4, pp. 437–444, 2009. [Online]. Available: https://doi.org/10.1080/13518470802604499
  * [20] P. Zarchan and H. Musoff, _Fundamentals of Kalman filtering: a practical approach_. American Inst. of Aeronautics and Astronautics, 2005.
  * [21] A. Tinka, M. Rafiee, and A. M. Bayen, “Floating sensor networks for river studies,” _IEEE Systems Journal_ , vol. 7, no. 1, pp. 36–49, March 2013.
  * [22] E. Wan and R. Van Der Merwe, “The unscented kalman filter for nonlinear estimation,” vol. 153-158, 02 2000, pp. 153 – 158.
  * [23] D. Kun, L. Kaijun, and X. Qunsheng, “Application of unscented kalman filter for the state estimation of anti-lock braking system,” in _2006 IEEE International Conference on Vehicular Electronics and Safety_ , Dec 2006, pp. 130–133.
  * [24] M. S. Arulampalam, S. Maskell, N. Gordon, and T. Clapp, “A tutorial on particle filters for online nonlinear/non-gaussian bayesian tracking,” _IEEE Transactions on Signal Processing_ , vol. 50, no. 2, pp. 174–188, Feb 2002.
  * [25] A. U. Peker, O. Tosun, and T. Acarman, “Particle filter vehicle localization and map-matching using map topology,” in _2011 IEEE Intelligent Vehicles Symposium (IV)_ , June 2011, pp. 248–253.
  * [26] K. Weekly, H. Zou, L. Xie, Q. Jia, and A. M. Bayen, “Indoor occupant positioning system using active rfid deployment and particle filters,” in _2014 IEEE International Conference on Distributed Computing in Sensor Systems_ , May 2014, pp. 35–42.
  * [27] Y. Cheng and Z. Liu, “Optimized selection of sigma points in the unscented kalman filter,” in _2011 International Conference on Electrical and Control Engineering_ , Sept 2011, pp. 3073–3075.
  * [28] P. Jiang, Y. Sun, and W. Bao, “State estimation of conceptual hydrological models using unscented Kalman filter,” _Hydrology Research_ , vol. 50, no. 2, pp. 479–497, 12 2018. [Online]. Available: https://doi.org/10.2166/nh.2018.038
  * [29] Zhang, Hongjuan, H. Franssen, Harrie-Jan, J. A., and Harry, “State and parameter estimation of two land surface models using the ensemble kalman filter and the particle filter,” Sep 2017. [Online]. Available: https://www.hydrol-earth-syst-sci.net/21/4927/2017
  * [30] J. Rasmussen, H. Madsen, K. H. Jensen, and J. C. Refsgaard, “Data assimilation in integrated hydrological modeling using ensemble kalman filtering: evaluating the effect of ensemble size and localization on filter performance,” _Hydrology and Earth System Sciences_ , vol. 19, no. 7, pp. 2999–3013, 2015. [Online]. Available: https://www.hydrol-earth-syst-sci.net/19/2999/2015/
  * [31] A. H. Weerts and G. Y. H. El Serafy, “Particle filtering and ensemble kalman filtering for state updating with hydrological conceptual rainfall-runoff models,” _Water Resources Research_ , vol. 42, no. 9, 2006. [Online]. Available: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2005WR004093
  * [32] D. Plaza, K. R. De, and V. Pauwels, “From the ensemble kalman filter to the particle filter: a comparative study in rainfall-runoff models*,” _IFAC Proceedings Volumes_ , vol. 44, no. 1, pp. 14 201 – 14 206, 2011, 18th IFAC World Congress. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S1474667016459082
  * [33] G. Evensen, “The ensemble kalman filter: theoretical formulation and practical implementation,” _Ocean Dynamics_ , vol. 53, no. 4, pp. 343–367, Nov 2003\. [Online]. Available: https://doi.org/10.1007/s10236-003-0036-9
  * [34] ——, “Sequential data assimilation with a nonlinear quasi-geostrophic model using monte carlo methods to forecast error statistics,” _Journal of Geophysical Research: Oceans_ , vol. 99, no. C5, pp. 10 143–10 162, 1994. [Online]. Available: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/94JC00572
  * [35] G. Burgers, P. Jan van Leeuwen, and G. Evensen, “Analysis scheme in the ensemble kalman filter,” _Monthly Weather Review_ , vol. 126, no. 6, pp. 1719–1724, 1998. [Online]. Available: https://doi.org/10.1175/1520-0493(1998)126<1719:ASITEK>2.0.CO;2
  * [36] J. Kaipio and E. Somersalo, _Statistical and computational inverse problems_. Springer, 2010.
  * [37] D. Work, O.-P. Tossavainen, S. Blandin, A. Bayen, T. Iwuchukwu, and K. Tracton, “An ensemble kalman filtering approach to highway traffic estimation using gps enabled mobile devices,” 01 2008, pp. 5062–5068.
  * [38] A. Seppänen, M. Vauhkonen, E. Somersalo, and J. P. Kaipio, “State space models in process tomography — approximation of state noise covariance,” _Inverse Problems in Engineering_ , vol. 9, no. 5, pp. 561–585, 2001. [Online]. Available: https://doi.org/10.1080/174159701088027781
  * [39] J. P. KAIPIO, P. A. KARJALAINEN, E. SOMERSALO, and M. VAUHKONEN, “State estimation in time-varying electrical impedance tomography,” _Annals of the New York Academy of Sciences_ , vol. 873, no. 1, pp. 430–439, 1999. [Online]. Available: https://nyaspubs.onlinelibrary.wiley.com/doi/abs/10.1111/j.1749-6632.1999.tb09492.x
  * [40] A. Seppänen, M. Vauhkonen, P. J Vauhkonen, E. Somersalo, and J. Kaipio, “State estimation with fluid dynamical evolution models in process tomography: Eit application,” _Inverse Problems - INVERSE PROBL_ , vol. 17, 08 2000.
  * [41] A. Lehikoinen, S. Finsterle, A. Voutilainen, M. Kowalsky, and J. Kaipio, “Dynamical inversion of geophysical ert data: state estimation in the vadose zone,” _Inverse Problems in Science and Engineering_ , vol. 17, no. 6, pp. 715–736, 2009. [Online]. Available: https://doi.org/10.1080/17415970802475951
  * [42] D. Baroudi, J. Kaipio, and E. Somersalo, “Dynamical electric wire tomography: a time series approach,” _Inverse Problems_ , vol. 14, no. 4, pp. 799–813, aug 1998.
  * [43] D. B. Work, S. Blandin, O.-P. Tossavainen, B. Piccoli, and A. M. Bayen, “A Traffic Model for Velocity Data Assimilation,” _Applied Mathematics Research eXpress_ , vol. 2010, no. 1, pp. 1–35, 04 2010. [Online]. Available: https://doi.org/10.1093/amrx/abq002
  * [44] “Mobile millennium.” [Online]. Available: http://traffic.berkeley.edu/
  * [45] G. Evensen, _Data Assimilation The Ensemble Kalman Filter_. Springer Berlin, 2014.
  * [46] A. J. Majda and X. T. Tong, “Performance of ensemble kalman filters in large dimensions,” _Communications on Pure and Applied Mathematics_ , vol. 71, no. 5, pp. 892–937. [Online]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1002/cpa.21722
  * [47] H. H. Bauser, D. Berg, O. Klein, and K. Roth, “Inflation method for ensemble kalman filter in soil hydrology,” _Hydrology and Earth System Sciences Discussions_ , pp. 1–18, 03 2018.
  * [48] Y. Shi, K. J. Davis, F. Zhang, C. J. Duffy, and X. Yu, “Parameter estimation of a physically-based land surface hydrologic model using an ensemble kalman filter: A multivariate real-data experiment,” _Advances in Water Resources_ , vol. 83, pp. 421 – 427, 2015. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0309170815001359
  * [49] C. Daly, M. Halbleib, J. I. Smith, W. P. Gibson, M. K. Doggett, G. H. Taylor, J. Curtis, and P. P. Pasteris, “Physiographically sensitive mapping of climatological temperature and precipitation across the conterminous united states,” _International Journal of Climatology_ , vol. 28, no. 15, pp. 2031–2064, 2008. [Online]. Available: https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/joc.1688
  * [50] S. A. Margulis, G. Cortés, M. Girotto, and M. Durand, “A landsat-era sierra nevada snow reanalysis (1985–2015),” _Journal of Hydrometeorology_ , vol. 17, no. 4, pp. 1203–1221, 2016. [Online]. Available: https://doi.org/10.1175/JHM-D-15-0177.1
  * [51] L. S. Huning and S. A. Margulis, “Climatology of seasonal snowfall accumulation across the sierra nevada (usa): Accumulation rates, distributions, and variability,” _Water Resources Research_ , vol. 53, no. 7, pp. 6033–6049. [Online]. Available: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2017WR020915
  * [52] A. Arnold, D. Calvetti, and E. Somersalo, “Parameter estimation for stiff deterministic dynamical systems via ensemble kalman filter,” _Inverse Problems_ , vol. 30, no. 10, p. 105008, sep 2014.
  * [53] M. Rafiee, A. Tinka, J. Thai, and A. M. Bayen, “Combined state-parameter estimation for shallow water equations,” in _Proceedings of the 2011 American Control Conference_ , June 2011, pp. 1333–1339.
  * [54] C. Zoccarato, D. Baù, M. Ferronato, G. Gambolati, A. Alzraiee, and P. Teatini, “Data assimilation of surface displacements to improve geomechanical parameters of gas storage reservoirs,” _Journal of Geophysical Research: Solid Earth_ , vol. 121, no. 3, pp. 1441–1461, 2016. [Online]. Available: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2015JB012090
  * [55] H. Steven L., L. E., Clark, and M. P., “Towards simplification of hydrologic modeling: identification of dominant processes,” Nov 2016. [Online]. Available: https://www.hydrol-earth-syst-sci.net/20/4655/2016/