3.2.1 LL Dual EnKF

Suppose a dynamical system can be described by a vector of states x_t and outputs y_t and a vector of associated model parameters θ_t at any given time t. The uncertain system states and parameters are represented by an ensemble of states ${\left\{{\mathit{x}}_t^i\right\}}_{i=\mathrm{1}:n}$ and parameters ${\left\{{\mathit{\theta }}_t^i\right\}}_{i=\mathrm{1}:n}$ each with n members. The prior state and parameter distributions ${\left\{{\mathit{x}}_t^{i-}\right\}}_{i=\mathrm{1}:n}$ and ${\left\{{\mathit{\theta }}_t^{i-}\right\}}_{i=\mathrm{1}:n}$, respectively, represent our prior knowledge of the system, usually derived as the output from a numerical model. Suppose also that the system outputs are observed $\left({\mathit{y}}_t^o\right)$ but that there is also some uncertainty associated with these observations. The purpose of the data assimilation algorithm (here the EnKF) is to combine the prior estimates with measurements, based on their respective uncertainties, to obtain an improved estimate of the system states and parameters. A single cycle of the LL Dual EnKF procedure for a given time t is undertaken as follows. Note that in the following, the overbar notation is used to indicate the ensemble mean.

1. Propose a prior parameter ensemble. This involves generating a parameter ensemble using prior knowledge. In this case, our prior knowledge comes from the updated parameter ensemble from the previous time (${\mathit{\theta }}_{t-\mathrm{1}}^{i+})$ and how it has changed over recent time steps. The assumed parameter dynamics is a Gaussian random walk with time-varying mean and variance, given by

   $$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }{\mathit{\theta }}_t^{i-}\sim N\left({\mathit{\theta }}_{t-\mathrm{1}}^{i+}+{\mathit{m}}_t.\mathrm{\Delta }t,s^{\mathrm{2}}{\mathbf{\Sigma }}_{t-\mathrm{1}}^{\mathit{\theta }}\right)\mathrm{for}i=\mathrm{1}:n\text{(3)}& {\displaystyle }& {\displaystyle }{\mathbf{\Sigma }}_{t-\mathrm{1}}^{\mathit{\theta }}={\displaystyle \frac{\mathrm{1}}{n-\mathrm{1}}}\sum _{i=\mathrm{1}}^n\left({\mathit{\theta }}_{t-\mathrm{1}}^{i+}-\stackrel{\mathrm{‾}}{{\mathit{\theta }}_{t-\mathrm{1}}^{+}}\right){\left({\mathit{\theta }}_{t-\mathrm{1}}^{i+}-\stackrel{\mathrm{‾}}{{\mathit{\theta }}_{t-\mathrm{1}}^{+}}\right)}^{\mathrm{T}},\end{array}$$

   where ${\mathbf{\Sigma }}_{t-\mathrm{1}}^{\mathit{\theta }}$ is the sample covariance matrix of the updated parameter ensemble at time t−1; $\stackrel{\mathrm{‾}}{{\mathit{\theta }}_{t-\mathrm{1}}^{+}}$ indicates the ensemble mean of the updated parameters at time t−1; ()^T represents the transpose operator; and s² is a tuning parameter. The prior ensemble mean is determined as the linear extrapolation of the updated ensemble means from the previous two time steps; i.e.,

   $$\begin{array}{}\text{(4)}& {\displaystyle }& {\displaystyle }{\mathit{m}}_t\left[k\right]=\left\{\begin{array}{l}{\mathit{m}}_{t-\mathrm{1}}\left[k\right],\left|{\mathit{m}}_{t-\mathrm{1}}\left[k\right]\right|\le m_{\max }\\ {\mathit{m}}_{t-\mathrm{2}}\left[k\right],\left|{\mathit{m}}_{t-\mathrm{1}}\left[k\right]\right|>m_{\max }\end{array}\right.\text{(5)}& {\displaystyle }& {\displaystyle }{\mathit{m}}_{t-\mathrm{1}}={\displaystyle \frac{\stackrel{\mathrm{‾}}{{\mathit{\theta }}_{t-\mathrm{1}}^{+}}-\stackrel{\mathrm{‾}}{{\mathit{\theta }}_{t-\mathrm{2}}^{+}}}{\mathrm{\Delta }t}}\text{(6)}& {\displaystyle }& {\displaystyle }{\mathit{m}}_{t-\mathrm{2}}={\displaystyle \frac{\stackrel{\mathrm{‾}}{{\mathit{\theta }}_{t-\mathrm{2}}^{+}}-\stackrel{\mathrm{‾}}{{\mathit{\theta }}_{t-\mathrm{3}}^{+}}}{\mathrm{\Delta }t}},\end{array}$$

   where m_t[k] indicates the kth component of the vector m_t, the estimated rate of change. Note that the extrapolation is forced to be less than a predefined maximum rate of change m_max to minimize overfitting and avoid parameter drift due to isolated large updates. The maximum rate of change is model-specific and will depend on the modeler's judgement regarding expected extreme changes.

2. Consider observation and forcing uncertainty. This is done by perturbing measurements of forcings and system outputs with random noise sampled from a distribution representing the uncertainty in those measurements. The result is an ensemble of forcings (${\mathit{u}}_t^i)$ and observations (${\mathit{y}}_t^i)$ each with n members. For example, if random errors in measurements of system outputs (herein also referred to as observations) are characterized by a zero mean Gaussian distribution, the ensemble of observations is given by

   $$\begin{array}{}\text{(7)}& {\displaystyle }{\mathit{y}}_t^i\sim N\left({\mathit{y}}_t^o,{\mathbf{\Sigma }}_t^{y^o{y}^o}\right)\mathrm{for}i=\mathrm{1}:n,\end{array}$$

   where ${\mathit{y}}_t^o$ is the recorded measurement at time t and ${\mathbf{\Sigma }}_t^{y^o{y}^o}$ is the error covariance matrix of the measurements.

3. Generate simulations using prior parameters. The prior parameters from Step 1, ${\mathit{\theta }}_t^{i-}$, and updated states from the previous time, ${\mathit{x}}_{t-\mathrm{1}}^{i+}$, are forced through the model equations to generate an ensemble of model simulations of states (${\widehat{\mathit{x}}}_t^i)$ and outputs (${\widehat{\mathit{y}}}_t^i)$:

   $$\begin{array}{}\text{(8)}& {\displaystyle }& {\displaystyle }{\widehat{\mathit{x}}}_t^i=f\left({\mathit{x}}_{t-\mathrm{1}}^{i+},{\mathit{\theta }}_t^{i-},{\mathit{u}}_t^i\right)\mathrm{for}i=\mathrm{1}:n\text{(9)}& {\displaystyle }& {\displaystyle }{\widehat{\mathit{y}}}_t^i=h\left({\widehat{\mathit{x}}}_t^i,{\mathit{\theta }}_t^{i-}\right)\mathrm{for}i=\mathrm{1}:n.\end{array}$$

4. Perform the Kalman update of parameters. Parameters are updated using the Kalman update equation and the prior parameter and simulated output ensemble from Step 1 and 3:

   $$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle }{\mathit{\theta }}_t^{i+}={\mathit{\theta }}_t^{i-}+{\mathit{K}}_t^{\mathit{\theta }}\left({\mathit{y}}_t^i-{\widehat{\mathit{y}}}_t^i\right)\mathrm{for}i=\mathrm{1}:n\text{(11)}& {\displaystyle }& {\displaystyle }{\mathit{K}}_t^{\mathit{\theta }}={\mathbf{\Sigma }}_t^{\mathit{\theta }\widehat{y}}{\left[{\mathbf{\Sigma }}_t^{\widehat{y}\widehat{y}}+{\mathbf{\Sigma }}_t^{y^o{y}^o}\right]}^{-\mathrm{1}},\end{array}$$

   where ${\mathbf{\Sigma }}_t^{\mathit{\theta }\widehat{y}}$is a matrix of the sample cross-covariance between errors in parameters ${\mathit{\theta }}_t^{i-}$ and simulated output ${\widehat{\mathit{y}}}_t^i$ ; and ${\mathbf{\Sigma }}_t^{\widehat{y}\widehat{y}}$is the sample error covariance matrix of the simulated output:

   $$\begin{array}{}\text{(12)}& {\displaystyle }& {\displaystyle }{\mathbf{\Sigma }}_t^{\mathit{\theta }\widehat{y}}={\displaystyle \frac{\mathrm{1}}{n-\mathrm{1}}}\sum _{i=\mathrm{1}}^n\left({\mathit{\theta }}_t^{i-}-\stackrel{\mathrm{‾}}{{\mathit{\theta }}_t^{-}}\right){\left({\widehat{\mathit{y}}}_t^i-\stackrel{\mathrm{‾}}{{\widehat{\mathit{y}}}_t}\right)}^{\mathrm{T}}\text{(13)}& {\displaystyle }& {\displaystyle }{\mathbf{\Sigma }}_t^{\widehat{y}\widehat{y}}={\displaystyle \frac{\mathrm{1}}{n-\mathrm{1}}}\sum _{i=\mathrm{1}}^n\left({\widehat{\mathit{y}}}_t^i-\stackrel{\mathrm{‾}}{{\widehat{\mathit{y}}}_t}\right){\left({\widehat{\mathit{y}}}_t^i-\stackrel{\mathrm{‾}}{{\widehat{\mathit{y}}}_t}\right)}^{\mathrm{T}}.\end{array}$$

5. Generate simulations using updated parameters. Step 3 is repeated with the updated parameter ensemble ${\mathit{\theta }}_t^{i+}$ to generate the prior ensemble of model simulations of states (${\mathit{x}}_t^{i-})$ and outputs (${\stackrel{\mathrm{̃}}{\mathit{y}}}_t^i)$:

   $$\begin{array}{}\text{(14)}& {\displaystyle }& {\displaystyle }{\mathit{x}}_t^{i-}=f\left({\mathit{x}}_{t-\mathrm{1}}^{i+},{\mathit{\theta }}_t^{i+},{\mathit{u}}_t^i\right)\mathrm{for}i=\mathrm{1}:n\text{(15)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{y}}}_t^i=h\left({\mathit{x}}_t^{i-},{\mathit{\theta }}_t^{i+}\right)\mathrm{for}i=\mathrm{1}:n.\end{array}$$

6. Perform the Kalman update of states and outputs. Use the Kalman update equation for correlated measurement and process noise (Eqs. 15 to 18) and the simulated state (${\mathit{x}}_t^{i-})$ and output (${\stackrel{\mathrm{̃}}{\mathit{y}}}_t^i)$ ensembles from Step 5 to update them. Since the measurements have already been used to generate ${\stackrel{\mathrm{̃}}{\mathit{y}}}_t^i$, the errors in model simulations and measurements are now correlated. The standard Kalman update equation (as in the form of Eqs. 10 and 11) can no longer be used as it relies on the assumption that errors in measurements and model simulations are independent.

   $$\begin{array}{ll}\text{(16)}& {\displaystyle }& {\displaystyle }{\mathit{x}}_t^{i+}={\mathit{x}}_t^{i-}+{\mathit{K}}_t^x\left({\mathit{y}}_t^i-{\stackrel{\mathrm{̃}}{\mathit{y}}}_t^i\right)\mathrm{for}i=\mathrm{1}:n{\displaystyle }& {\displaystyle }{\mathit{K}}_t^x=\left[{\mathbf{\Sigma }}_t^{x\stackrel{\mathrm{̃}}{y}}+{\mathbf{\Sigma }}_t^{{\mathit{\epsilon }}_x{y}^o}\right]\\ \text{(17)}& {\displaystyle }& {\displaystyle }{\left[{\mathbf{\Sigma }}_t^{\stackrel{\mathrm{̃}}{y}\stackrel{\mathrm{̃}}{y}}+{\mathbf{\Sigma }}_t^{{\mathit{\epsilon }}_{\stackrel{\mathrm{̃}}{y}}{y}^o}+{\left({\mathbf{\Sigma }}_t^{{\mathit{\epsilon }}_{\stackrel{\mathrm{̃}}{y}}{y}^o}\right)}^{\mathrm{T}}+{\mathbf{\Sigma }}_t^{y^o{y}^o}\right]}^{-\mathrm{1}}\text{(18)}& {\displaystyle }& {\displaystyle }{\mathit{\epsilon }}_{x_t}^i={\mathit{x}}_t^{i-}-{\widehat{\mathit{x}}}_t^i\text{(19)}& {\displaystyle }& {\displaystyle }{\mathit{\epsilon }}_{{\stackrel{\mathrm{̃}}{y}}_t}^i={\stackrel{\mathrm{̃}}{\mathit{y}}}_t^i-{\widehat{\mathit{y}}}_t^i,\end{array}$$

   where ${\mathbf{\Sigma }}_t^{x\stackrel{\mathrm{̃}}{y}}$ is a matrix of the sample cross-covariance between simulated states ${\left\{{\mathit{x}}_t^{i-}\right\}}_{i=\mathrm{1}:n}$ and outputs ${\left\{{\stackrel{\mathrm{̃}}{\mathit{y}}}_t^i\right\}}_{i=\mathrm{1}:n}$ from Step 5; ${\mathbf{\Sigma }}_t^{{\mathit{\epsilon }}_x{y}^o}$ represents the sample covariance between ${\left\{{\mathit{\epsilon }}_{x_t}^i\right\}}_{i=\mathrm{1}:n}$ and the observations; and ${\mathbf{\Sigma }}_t^{{\mathit{\epsilon }}_{\stackrel{\mathrm{̃}}{y}}{y}^o}$ represents the sample covariance between the ${\left\{{\mathit{\epsilon }}_{{\stackrel{\mathrm{̃}}{y}}_t}^i\right\}}_{i=\mathrm{1}:n}$ and the observations.

The above algorithm specifies the updating of states and parameters at any given time, based on available observations. This allows one to retrospectively estimate time variations in model parameters, as well as provide one-time-step-ahead forecasts of states and outputs (as per Eqs. 8 and 9). Forecasts at longer time horizons (i.e., longer than one time step ahead) would be made by generating prior parameters and states as detailed in Steps 1 to 3, although the local linear extrapolations are only valid close to the current time point.

3.2.2 Application to the Nam Muc catchment

Joint state and parameter estimation was undertaken for the Nam Muc catchment over the period 1980 to 2004 by assimilating streamflow observations into the HyMOD and HBV models at a daily time step. Additionally, simulations using the time-invariant parameters obtained from calibration over the period 1973–1979 were generated for 1980 to 2004, for comparison. Estimating a large number of parameters from limited data is problematic in that the system is highly under-determined, making it difficult to ensure the estimated parameters are meaningful. Given the fairly low parameter dimensionality of HyMOD, all model parameters were allowed to vary in time, while for HBV we applied the Sobol method to identify the most sensitive parameters to be included in the time-varying parameter estimation. The Sobol method is a global sensitivity analysis method based on variance decomposition. It identifies the partial variance contribution of each parameter to the total variance of the hydrological model output (see for example Saltelli et al., 2008; Nossent et al. 2011). The method, implemented through the SAFE toolbox (Pianosi et al., 2015), found the lp and Maxbas parameters to be the least sensitive and least important in defining variations in catchment hydrology (see Table 3). These were held fixed (lp = 1 and Maxbas = 1 day) in the following analysis. Note that although the hl1 parameter was found to have low sensitivity, it was retained as a time-varying parameter due to its conceptual importance in separating interflow and near-surface flow (refer Fig. 3).

Table 3Variance-based sensitivity analysis results for HBV parameters: first-order sensitivity index representing the contribution of varying a single parameter to the variance of the model output. Lower values indicate lower sensitivity.

Download Print Version | Download XLSX

Table 4LL EnKF inputs for the HBV model case.

Download Print Version | Download XLSX

Unbiased normally distributed ensembles of the parameters and states are required to initialize the LL Dual EnKF. Initial parameter ensembles were generated by sampling from a Gaussian distribution with a mean equal to the calibrated parameters over the pre-change period and variance estimated from parameter sets with similar objective function values. Parameter sets with similar objective function values were obtained when using different starting points to the optimization algorithm during the model calibration stage. Initial state ensembles were also sampled from normal distributions with a mean equal to the simulated state at the end of the calibration period. An ensemble size of 100 members was adopted and assumed sufficiently large based on the findings of Moradkhani et al. (2005) and Aksoy et al. (2006). Due to the stochastic–dynamic nature of the method, ensemble statistics were calculated over 20 separate realizations of the LL Dual EnKF. The prior parameter-generating method described in Step 1 of Sect. 3.2 requires specification of the tuning parameter s² to define the variance of the perturbations. This was tuned by selecting the s² value that optimized the quality of forecast streamflow over the calibration period. Forecast quality was assessed using the logarithmic score (LS) (Good, 1952) of background streamflow predictions (${\stackrel{\mathrm{̃}}{y}}_t^i)$ using updated parameters (Eq. 15), which was averaged over the calibration period of length T:

where f(y) is the probability density function of the background streamflow predictions (represented by the empirical probability density function of the sample points ${\left\{{\stackrel{\mathrm{̃}}{y}}_t^i\right\}}_{i=\mathrm{1}:n})$; and $y_t^o$ is the measurement of the system outputs. The s² value that gave the largest $\stackrel{\mathrm{‾}}{\mathrm{LS}}$ was adopted for the assimilation period. The maximum allowable daily rate of change in the ensemble mean was based on assuming a linear rate of change within the entire feasible parameter space over a 3-year period.

Figure 4Parameter trajectories using the HBV model. The dark grey shaded areas indicate the middle 90 % of the ensemble, bounded by the 5th and 95th percentiles. The light grey shaded areas indicate the middle 50 % of the ensemble, bounded by the 25th and 75th percentiles. The ensemble mean is indicated by the blue line. The vertical green panel indicates the assumed time period of rapid deforestation.

Download

Figure 5Parameter trajectories using the HyMOD model. The dark grey shaded areas indicate the middle 90 % of the ensemble, bounded by the 5th and 95th percentiles. The light grey shaded areas indicate the middle 50 % of the ensemble, bounded by the 25th and 75th percentiles. The ensemble mean is indicated by the blue line. The vertical green panel indicates the assumed time period of rapid deforestation.

Download

As detailed in Sect. 3.2, observation and forcing uncertainty is considered by perturbing measurements with random noise. Here streamflow errors were assumed to be zero-mean normally distributed (truncated to ensure positivity) and heteroscedastic. The variance is defined as a proportion of the observed streamflow, to reflect the fact that larger flows tend to have greater errors than low flows:

where TN indicates the truncated normal distribution to ensure positive flows and d= 0.1. A multiplier of 0.1 was chosen based on estimates adopted for similar gauges in hydrologic data assimilation (DA) studies (e.g., Clark et al., 2008; Weerts and El Serafy, 2006; Xie et al., 2014).

Table 5LL EnKF inputs for the HYMOD model case.

Download Print Version | Download XLSX

Several studies have noted that a major source of rainfall uncertainty arises from scaling point rainfall to the catchment scale (Villarini and Krajewski, 2008; McMillan et al., 2011) and that multiplicative error models are suited to describing such errors (e.g., Kavetski et al., 2006). Rainfall uncertainties were therefore described using unbiased, lognormally distributed multipliers:

where P_t is the measured rainfall at time t; m and v are the mean and variance of the lognormally distributed rainfall multipliers M, respectively; and μ and σ² are the mean and variance of the normally distributed logarithm of the rainfall multipliers M. For unbiased perturbations, we let m= 1. The variance of the rainfall multipliers (v) was estimated by considering upper and lower bound error estimates in the Thiessen weights assigned to the four rainfall stations (see Sect. 2.1 for calculation of catchment-averaged rainfall, P_t). The resulting upper and lower bound catchment-averaged rainfall data were then used to estimate error parameters due to spatial variation in rainfall:

where P_upper,10 indicates catchment-averaged rainfall data estimated using the upper bound Thiessen weights with daily depth greater than 10 mm (similar for P_lower,10). A 10 mm rainfall depth threshold was chosen to avoid large rainfall fractions due to small rainfall depths. $\widehat{{\mathit{\sigma }}^{\mathrm{2}}}$ was found to be 0.05 in this case study. Similarly, we assume the dominant source of uncertainty in temperature data arises from spatial variation. Differences in temperature records at Lai Chau and Quynh Nhai (only available gauges with temperature records) were analyzed and found to be approximately normally distributed with a sample mean of 0.2 ^∘C and variance of 1.4 ^∘C. A perturbed temperature ensemble was then generated according to Eq. (28):

where $T_t^{\mathrm{avg}}$ represents catchment-averaged temperature data (see Sect. 2.1). Note that perturbations were taken to be unbiased (zero mean) as the sample mean of the differences in the temperature records was close to zero. The same perturbed input and observation sequences were used for the HyMOD and HBV runs for the sake of comparison. A summary of the values adopted for the various components of the LL Dual EnKF for each model is provided in Tables 4 and 5.

Figure 6Representative hydrographs of background streamflow from the LL Dual EnKF (black line), time-varying parameter model with no state updating (blue line), time-invariant parameter model with no DA (green line) and observed streamflow (red line). Results for HBV are shown in the top row and HyMOD in the bottom row. A pre-change year (1974) is shown on the left and a post-change year (1998) on the right.

Download