Development of a Regularized Dynamic System Response Curve for Real-Time Flood Forecasting Correction

Abstract

The dynamic system response curve (DSRC) is commonly applied as a real-time flood forecasting error correction method to improve the accuracy of real-time flood forecasting. It has been widely recognized that the least squares (OLS/LS) method, employed by DSRC, breaks down ill-posed problems, and therefore, the DSRC method may lead to deterioration in performance caused by meaningless solutions. To address this problem, a diagnostically theoretical analysis was conducted to investigate the relationship between the numerical solution of the Fredholm equation of the first kind and the DSRC method. The analysis clearly demonstrates the derivation of the problem and has implications for an improved approach. To overcome the unstable problem, a new method using regularization techniques (Tikhonov regularization and L-Curve criterion) is proposed. Moreover, in this study, to improve the performance of hydrological models, the new method is used as an error correction method to correct a variable from a hydrological model. The proposed method incorporates the information from a hydrological model structure. Based on the analysis of the hydrological model, the free water storage of the Xinanjiang rainfall-runoff (XAJ) model is corrected to improve the model’s performance. A numerical example and a real case study are presented to compare the two methods. Results from the numerical example indicate that the mean Nash–Sutcliffe efficiency value (NSE) of the regularized DSRC method (RDSRC) decreased from 0.99 to 0.55, while the mean NSE of DSRC decreased from 0.98 to −1.84 when the noise level was increased. The overall performance measured by four different criteria clearly demonstrates the robustness of the RDSRC method. Similar results were obtained for the real case study. The mean NSE of 35 flood events obtained by RDSRC method was 0.92, which is significantly higher than the mean NSE of DSRC (0.7). The results demonstrate that the RDSRC method is much more robust than the DSRC method. The applicability and usefulness of the RDSRC approach for real-time flood forecasting is demonstrated via the numerical example and the real case study.

1. Introduction

Flood forecasting is critical in many applications, such as flood warning [1], operation of a multi-reservoir system [2] and water resource management [3]. Conceptual rainfall-runoff models simplify and conceptualize these complex processes using a set of simple mathematical equations [4,5]. However, it is generally accepted that the satisfactory application of such a model can be hampered by many factors, including errors in the inputs [6,7], inadequacy of the model which generally includes errors in the model structure [8], errors in the model parameters [9] and errors in the model initial conditions [10,11]. The error correction (updating) methods have been made to handle the problem [12]. Several authors have proposed correction methods to update the river flow (discharge). Madsen and Skotner [13] incorporated an efficient data assimilation procedure into the update of river flow. Bao et al. [14] proposed a correction method for the roughness of the hydrodynamic model. Shamseldin and O’Connor [15] proposed an approach for updating river flow based on the neural network. Auto Regressive (AR) time series models have been used to improve the performance of predictions by predicting future discrepancies between model simulations and observations [16,17,18]. Sequential data assimilation (DA), which incorporates information from observations into model simulations, has demonstrated its applicability in improving model performance [19,20,21]. A few sequential DA techniques have been introduced in the hydrologic literature which mainly include the dynamic identifiability analysis (DYNIA) [22], the parameter estimation approach based on information localization [23], the Bayesian recursive estimation technique (BaRE) [24], etc. The standard Kalman filter is proposed to provide a recursive solution to the discrete-data linear filter problem [25]. Several nonlinear extensions of Kalman filter, such as extended Kalman filter (EKF) [26,27] and ensemble Kalman filter (EnKF) [28,29,30] have been utilized due to the nonlinearity of the hydrologic system. Particle filtering (PF), which is also known as sequential Monto Carlo sampling (SMC), is another alternative data assimilation technique [31] used in the hydrological literature [29].

The dynamic system response curve (DSRC) method was proposed by Bao et al. [32] to update the simulated runoff. Further, the method was used to update the areal mean precipitation [33]. In the DSRC method, the hydrologic model is approximated by first-order Taylor linearization. Then, the correction is achieved by solving the corresponding equations using the least squares technique. While the DSRC method may be a simple and efficient approach, there has been some discussion about the flaws of the method [34]. One of the main concerns is its ill-posed property, which can result in instability of the least squares solution and hence, poor performance. For such situations, the so-called regularization method is used to stabilize the least squares solutions in many branches of science and engineering [35].

In this paper, a new diagnostically theoretical analysis of the standard DSRC method is presented. The analysis is illustrated with a proof in the first part of this paper, where the relationship between the DSRC method and the discretized form of the first kind Fredholm equation is demonstrated. Based on the proof, a new method using regularization techniques (Tikhonov regularization and L-Curve criterion) is proposed to overcome the unstable problem of DSRC. To improve the performance of hydrological models, the new method is used in this study as an error correction method to correct a variable of a hydrological model. The applicability and usefulness of RDSRC for error correction of flood forecasting is demonstrated with a numerical example and a real case study.

2. Materials and Methods

2.1. DSRC Method and Regularized DSRC Method (RDSRC)

2.1.1. DSRC Method and Its Flaws

In presenting the new approach based on regularization techniques, a brief sketch of the DSRC method is presented. More specifications of derivation, proofs, and applications can be found in the following references [32,33]. A mathematical form of the nonlinear hydrological system is given as

$$\mathbf{Q}=f(\mathbf{P},\mathbf{E},\mathbf{X},\mathsf{\theta })$$

(1)

where $\mathbf{P}={[p_1,p_2,p_3,\cdots ,p_n]}^T$ and $\mathbf{E}={[e_1,e_2,e_3,\cdots ,e_n]}^T$ represent time-ordered vectors of areal mean rainfall and areal mean evapotranspiration, respectively. $\mathbf{Q}={[q_1,q_2,q_3,\cdots ,q_n]}^T$ is defined as the vector of simulated discharge (output of the model). $\mathbf{X}={[x_1,x_2,x_3,\cdots ,x_m]}^T$ denotes the vector of model variables at time 1, 2, …, m, and $\mathsf{\theta }={[{\theta }_1,{\theta }_2,\dots ,{\theta }_k]}^T$ represents the vector of parameters of the hydrological model.

Since our main purpose is to update the model’s variables, it is assumed that the simulations are briefly determined by the model’s variables with the equation

$$\mathbf{Q}(\mathbf{X})={(q_1(\mathbf{X}),q_2(\mathbf{X}),\cdots ,q_n(\mathbf{X}))}^T$$

(2)

Equation (2) can be linearized using first-order Taylor expansion for non-linear relationships and the relationship can be rewritten as Equation (3) [32,33]

$$\mathbf{Q}(\mathbf{X}+\mathsf{\Delta }\mathbf{X})\approx \mathbf{Q}(\mathbf{X})+J_Q(\mathbf{X})\mathsf{\Delta }\mathbf{X}$$

(3)

with

$$\mathsf{\Delta }\mathbf{X}=\left[\mathsf{\Delta }{x}_1,\mathsf{\Delta }{x}_2,\dots ,\mathsf{\Delta }{x}_m\right]$$

(4)

where $J_Q(\mathbf{X})$ is a Jacobian matrix and $\mathsf{\Delta }{x}_i$ is the increment (correction) of $x_i$. The Jacobian matrix consists of elements which are also known as sensitivity coefficients in Yeh [36]:

$$J_Q(\mathbf{X})=\left(\begin{array}{l}\frac{\partial q_1(X)}{\partial x_1}\frac{\partial q_1(X)}{\partial x_2}\cdots \frac{\partial q_1(X)}{\partial x_m}\\ \frac{\partial q_2(X)}{\partial x_1}\frac{\partial q_2(X)}{\partial x_2}\cdots \frac{\partial q_2(X)}{\partial x_m}\\ ⋮⋮\ddots ⋮\\ \frac{\partial q_n(X)}{\partial x_1}\frac{\partial q_n(X)}{\partial x_2}\cdots \frac{\partial q_n(X)}{\partial x_m}\end{array}\right)$$

(5)

where $\frac{\partial q_i(\mathbf{X})}{\partial x_j}$ represents the partial derivative of $q_i(\mathbf{X})$ with respect to $x_j$. This matrix enables us to get access to the “hidden” information in the hydrological model (the relationship between the model variable to be updated and the model output of the hydrological model) by calculating the partial derivative values. After that, this relationship (information) is utilized to correct the model’s variable. Specifically, each element in the matrix indicates the relationship between a change in $x_i$ ($\mathsf{\Delta }{x}_i$) and the corresponding change in outputs $q_j(\mathbf{X})$ (j = 1, 2, …, n). In other words, the column “i” represents the influence of updating $x_i$. In practice, the elements (sensitivity coefficients) are obtained by a simple approach proposed by Becker and Yeh [37].

The response of the perturbation to model variable X is defined as [32]

$$\mathsf{\Delta }\mathbf{Q}=\mathbf{Q}(\mathbf{X}+\mathsf{\Delta }\mathbf{X})-\mathbf{Q}(\mathbf{X})={\left(\mathsf{\Delta }{q}_1,\dots ,\mathsf{\Delta }{q}_n\right)}^T$$

(6)

and Equation (3) is rewritten in the form

$$J_Q(\mathbf{X})\mathsf{\Delta }\mathbf{X}\approx \mathsf{\Delta }\mathbf{Q}$$

(7)

It is assumed that $\mathsf{\Delta }\mathbf{Q}$ can be substituted by $Q_{obs}-\mathbf{Q}$, approximately, and Equation (7) turns into

$$J_Q(\mathbf{X})\mathsf{\Delta }\mathbf{X}\approx Q_{obs}-\mathbf{Q}$$

(8)

Use the ordinary least squares method to give the solution of formula [32,33,38,39]

$$\mathsf{\Delta }\tilde{\mathbf{X}}={(J_Q{(X)}^T{J}_Q(\mathbf{X}))}^{-1}{J}_Q{(\mathbf{X})}^T(Q_{obs}-\mathbf{Q})$$

(9)

The solution vector is $\mathsf{\Delta }\tilde{\mathbf{X}}={\left(\mathsf{\Delta }{\tilde{x}}_1,\dots ,\mathsf{\Delta }{\tilde{x}}_m\right)}^T$. Thus, the variable can be updated by $\mathsf{\Delta }\tilde{\mathbf{X}}$ if both $\mathsf{\Delta }\mathbf{Q}$ and $J_Q(\mathbf{X})$ are known. Note that both approximations in Equations (3) and (8) can introduce errors to the right-hand side. However, when using the DSRC method, even small perturbations to the right-hand side may cause instability of the solution and poor performance [34], and hence an analysis of the DSRC approach is motivated by these problems.

The basis for our analysis in this section is a Fredholm integral equation of the first kind with the kernel function, $K(x,y)$ [40,41]

$$f(x)={\displaystyle {\int }_a^{b}K(x,y)\phi (y)dy}c\le x\le d$$

(10)

where $\phi (y)$ is the unknown function, and $f(x)$ is known. To investigate the relationship between the DSRC method and the integral equation, the authors assume that the kernel, $\phi (y)$ and $f(x)$ satisfy the following relationship

$$K(p,t)=\frac{\partial q_p(\mathbf{X})}{\partial x_t}$$

(11)

$$\phi (t)=\mathsf{\Delta }{x}_t$$

(12)

and

$$f(p)=\mathsf{\Delta }{q}_p$$

(13)

where $\frac{\partial q_p(\mathbf{X})}{\partial x_t}$ is defined in Equation (5). By discretizing the equation at a set of points, the numerical solution can be expressed as

$$\mathbf{F}\approx \mathbf{K}\mathsf{\Phi }$$

(14)

with

$$\mathbf{F}=\left[\begin{array}{c}\mathsf{\Delta }{q}_1\\ \mathsf{\Delta }{q}_2\\ ⋮\\ \mathsf{\Delta }{q}_n\end{array}\right],\mathbf{K}=\left[\begin{array}{cccc}\frac{\partial q_1(X)}{\partial x_1}& \frac{\partial q_1(X)}{\partial x_2}& \cdots & \frac{\partial q_1}{\partial x_m}\\ \frac{\partial q_2}{\partial x_1}& \frac{\partial q_2}{\partial x_2}& \cdots & \frac{\partial q_2}{\partial x_m}\\ ⋮& ⋮& \cdots & ⋮\\ \frac{\partial q_n}{\partial x_1}& \frac{\partial q_n}{\partial x_2}& \cdots & \frac{\partial q_n}{\partial x_m}\end{array}\right]\mathrm{and}\mathsf{\Phi }=\left[\begin{array}{c}\mathsf{\Delta }{x}_1\\ \mathsf{\Delta }{x}_2\\ ⋮\\ \mathsf{\Delta }{x}_m\end{array}\right]$$

(15)

The detailed deviation is provided in the Supplementary Materials. Thus, it has been demonstrated that under certain assumptions, the basic equation (Equation (8)) of the DSRC method has the same form as the numerical solution of the Fredholm equation of the first kind.

It is generally accepted that using direct schemes to compute the integral equation is not adequate due to the noise of $f(p)$ in [40,41]. Under the perturbation in $f(p)$, the solution to an integral solution is unstable. Therefore, even if there is a solution to Equation (10), a small change to $f(x)$ can lead to an arbitrarily large change in $\phi (y)$. In the view of Tikhonov [42,43], a problem is ill-posed if the stability of the solution is violated, and any numerical methods for the solution will carry the inherent ill-posed nature to the corresponding discrete problem [44,45]. Unfortunately, even if the round-off errors in the solution process are ignored, $f(p)$ is usually contaminated with unavoidable perturbation, as discussed above. Hence, the solution of DSRC may be unstable and useless. For DSRC, this phenomenon typically appears in applications when solving the problem (Equation (8)) with the least squares method. The solution of the DSRC method often makes no sense when it is put into the context of flood updating. Hence, it is necessary to choose a regularization method to dampen down the instabilities.

2.1.2. Regularization Techniques and RDSRC Method

In this section, consider the linear system, which is rewritten from Equation (7) in a more simple and general form

$$\mathbf{A}x\approx \mathbf{b}\mathbf{A}\in R^{n\times m},\mathbf{b}\in R^n$$

(16)

where $A$ is a coefficient matrix with $n\ge m$, and $b$ is the vector of measurement (also called the right-hand side). Due to its ill-posedness, the direct solution to the problem is unstable, and a small perturbation of the data can cause large perturbations in the solution [46]. The approach for solving such problems should take cognizance of the difficulties; regularization techniques are required to stabilize and solve the discrete ill-posed problem. A well-known method for dealing with ill-posed problems is the Tikhonov regularization approach [42]. The regularized solution, ${\mathbf{x}}_{\lambda }$, is the solution to the following least squares problem as a function of the regularization parameter

$$\min ({\Vert \mathbf{A}\mathbf{x}-\mathbf{b}\Vert }_2^2+{\lambda }^2{\Vert \mathbf{L}\mathbf{x}\Vert }_2^2)$$

(17)

where the regularization parameter, $\lambda$, is applied to handle the weight of minimization of the regularization part with respect to the minimization of the residual norm [45]. The regularization error will decrease as $\lambda$ decreases, but the regularized solutions are more sensitive to the perturbations. That is, the choice of the regularization parameter can be interpreted as a compromise between the accuracy and the sensitivity. In the context of updating simulations obtained by the XAJ model, the rationale under this compromise is that the solution for correcting the model variable to the updating problem typically has a finite norm which is restricted by the hydrologic model. By using Equation (19), the proposed method merges the information about the model structure with the information from model simulations and independent observations for error correction of the hydrological model. Specifically, for the free water storage of the XAJ model, the updated value should range from 0 to the free water storage capacity (upper bound).

L is the regularization matrix which is used to define a norm or seminorm for controlling the size and smoothness of the solution. Specifically, in the standard form of the Tikhonov regularization, the matrix, $\mathbf{L}$, is considered to be I (identity matrix). In many applications, choosing a suitable matrix, $\mathbf{L}\ne \mathbf{I}$, is sometimes necessary to obtain better solutions [47]. However, adequate handling of the choice for the regularization matrix has hardly been discussed, and it is beyond the goal of this paper to present a comparison of these forms. Therefore, throughout the paper, the regularization matrix is assumed to be equal to the identity matrix. Hence, the Tikhonov regularization in standard form is

$$\min ({\Vert \mathbf{A}\mathbf{x}-\mathbf{b}\Vert }_2^2+{\lambda }^2{\Vert \mathbf{x}\Vert }_2^2)$$

(18)

and the corresponding solution is

$$\widehat{\mathbf{X}}={({\mathbf{A}}^T\mathbf{A}+{\lambda }^2\mathbf{I})}^{-1}{\mathbf{A}}^T\mathbf{b}$$

(19)

Using this formula, Equation (9) is replaced by

$$\mathsf{\Delta }\tilde{\mathbf{X}}={({\mathbf{J}}_Q{(\mathbf{X})}^T{\mathbf{J}}_Q(\mathbf{X})+{\lambda }^2\mathbf{I})}^{-1}{\mathbf{J}}_Q{(\mathbf{X})}^T({\mathbf{Q}}_{obs}-\mathbf{Q})$$

(20)

An important issue to consider is having the proper criterion for choosing the regularization parameter which determines the number of the smallest singular values to be neglected [48]. Unfortunately, there is no general criterion for this problem even though several approaches, such as the Generalized Cross-Validation (GCV) [49]; the Discrepancy Principle [50] and the L-Curve Criterion [51] have been proposed. The numerical experiment presented in [51] shows that the L-Curve Criterion has an overall more stable performance than the GCV method, and hence, in this study we used a combination of the L-Curve Criterion. Let ${\mathbf{x}}_{\lambda }$ donate the regularized solution, the L-curve is a plot for all valid regularization parameters of the residual norm, $\Vert \mathbf{A}{\mathbf{x}}_{\lambda }-\mathbf{b}\Vert$, versus the norm, $\Vert \mathbf{L}{\mathbf{x}}_{\lambda }\Vert$ [51,52,53]. The underlying idea for the L-curve method is to seek a balance between minimizing both the solution size and the residual size. The plot of the L-curve is typically a combination of two regions: the uppermost part, which is vertical, and the rightmost part, which is horizontal [45]. When the regularization parameter is small, the solution is sensitive to perturbations to the right-hand side and this results in the uppermost part. The horizontal part, where the regularization parameter is relatively big, corresponds to the solutions which are dominated by regularization errors. Then, the regularization parameter is chosen by locating the corner (the point of maximum curvature) of the curve. For detailed information about the L-Curve Criterion, refer to these references [51,52]. Together with Equation (20), we have presented the proposed method. For lack of a better name, we will distinguish the proposed method from the DSRC method by calling it regularized DSRC (RDSRC).

2.2. Hydrological Model

The Xinanjiang rainfall-runoff (XAJ) model is a conceptual rainfall-runoff model proposed by Zhao [54]. The XAJ model is one of the most widely used hydrologic models in humid and semi-humid regions in the south and east of China [55]. For the XAJ model, the basin is divided into several sub-basins and flood routing from the sub-basin outlets to the total basin outlet is implemented by applying the Muskingum flood routing method to successive sub-reaches [54]. The basic structure of the XAJ model consists of four parts: the three-layer evapotranspiration layer (TLE); the saturated runoff production layer (SRP); the runoff separation layer (SOR); and the flow concentration layer (FC). A block diagram illustrating the conceptual structure of the XAJ model is shown in Figure 1. Detailed information about the model’s parameters has been summarized in

Table 1. Parameters of the XAJ model.

Module	Function	Methods	Parameter	Meaning	Unit
TLE	Evapotranspiration	Three-layer soil moisture model	K	Ratio of potential evapotranspiration to pan evaporation	-
			WUM	Areal mean tension water capacity of the upper layer	mm
			WLM	Areal mean tension water capacity of the lower layer	mm
			WDM	Areal mean tension water capacity of the deeper layer	mm
			C	Coefficient of deep evapotranspiration	-
			IM	Ratio of impervious area	-
SRP	Runoff production	Runoff formation on storage repletion	WM	Areal mean tension water capacity	mm
			B	Exponent of the tension water capacity distribution curve	-
SOR	Runoff separation	Free water storage model	SM	Areal mean free water capacity of the surface soil layer	mm
			EX	Exponent of the free water capacity curve	-
			KI	Outflow coefficients of the free water storage to interflow	-
			KG	Outflow coefficients of the free water storage to groundwater	-
FC	Runoff concentration	Linear reservoir	CS	Recession constant of the surface water storage	-
			CI	Recession constant of the interflow storage	-
			CG	Recession constant of the groundwater storage	-
MSSR	Flood routing	Muskingum method	KE	Storage-time constant	-
			XE	Weight factor	-

.

As mentioned above, this paper focuses on updating the free water storage of the XAJ model to improve the accuracy of the model’s simulation. Therefore, an overview of the SOR layer is first presented. In the runoff separation layer, runoff (R) is separated into surface runoff (RS), sub-surface runoff (RI) and groundwater (RG) using the free water storage capacity curve. Implementation of the separation process is achieved with the following equations [54]:

$$RS(m)=f_1(S(m-1))$$

(21)

$$S(m)=f_2(S(m-1))$$

(22)

$$RI(m)=f_3(S(m))$$

(23)

$$RG(m)=f_4(S(m))$$

(24)

The meanings of these variables are presented in

Table 2. Descriptions of the variables.

Variable	Units	Description
RS	mm	Surface runoff
RI	mm	Interflow runoff
RG	mm	Groundwater runoff
S	mm	Free water storage

. We have replaced the mathematical equations of the model with the notations f₁, f₂, …, f₄ for simplicity. These equations specify that RS, RI and RG are actually functions of the free water storage (S). In other words, the update of the runoff components is achieved by modifying the free water storage and accordingly, changing the total flow. Note that for all iterations, the free water storage must satisfy [54]

$$0\le S\le SM$$

(25)

where SM represents the free water storage capacity. For detailed information about the free water storage and the runoff separation module, references [54,55,56,57,58] are recommended.

To examine the implications of the theoretical considerations and the new method, both the DSRC method and the RDSRC method were tested with a numerical example (synthetic case study) and a real case study. Before entering into the experimental section, it should be pointed out that the main purpose of this study was to test the methods with a numerical example.

2.3. Synthetic Case

Data Basis

The exact (ideal) discharge (the observed discharge is contaminated by the measurement noise), the exact values of the model variables and errors for the real cases cannot be accessed. Instead, a numerical example based on the XAJ model, where all these values are precisely known, is considered here. In this numerical example, the inputs, outputs, parameters, the model variables and the error vectors (the measurement error and the free water storage error) are artificially generated and hence, precisely known. A brief introduction of the numerical example is presented here for completeness; for detailed information about the numerical example, we recommend reading research by Bao, Si and Qu [32], Si, Bao and Gupta [33].

Here, ${\mathbf{Q}}_{exact}$ represents the output vector of exact (ideal) values which is directly generated by the inputs via a free run of the XAJ model using a predefined set of conditions (input data, model parameters and initial values), and the vector of ideal values of the free water storage ${\mathbf{S}}_{exact}$ can be obtained simultaneously. ${\mathbf{Q}}_{cal}$ represents the outputs of the XAJ model with the additional errors, ${\mathbf{S}}_{error}$, added in free water storage (free of the measurement noise ${\delta }_{\mathbf{Q}}$) under the same conditions, ${\mathbf{Q}}_{obs}$ represents the output vector of the observed discharge (with the noise ${\delta }_{\mathbf{Q}}$), ${\mathbf{Q}}_{DSRC}$ represents the output vector of the estimated discharge obtained by the DSRC method, and ${\mathbf{Q}}_{RDSRC}$ represents the output vector of the estimated discharge obtained by the RDSRC method. Hence, in this test, the right-hand side in Equation (7) is given by

$$\mathsf{\Delta }\mathbf{Q}=({\mathbf{Q}}_{exact}+{\delta }_{\mathbf{Q}})-{\mathbf{Q}}_{cal}$$

(26)

$({\mathbf{Q}}_{exact}+{\delta }_{\mathbf{Q}})$ is the vector of the measurements (i.e., ${\mathbf{Q}}_{obs}$) that consists of the vector of the exact values, ${\mathbf{Q}}_{exact}$, and white noise ${\delta }_{\mathbf{Q}}$ (with zero mean). ${\delta }_{\mathbf{Q}}$ is assumed to be additional noise which means that the noise signal can be directly added to the original signal. The variance in ${\delta }_{\mathbf{Q}}$ is obtained by defining a new notion-error level. This notion, which is used to quantify the noise in this paper, is analogous to the term used by Hansen and O’Leary [51]

$$level=\frac{{\Vert \mathbf{e}\Vert }_2}{{\Vert \mathbf{b}\Vert }_2}$$

(27)

where ${\Vert ·\Vert }_2$ is the 2-norm, $\mathbf{e}$ is the vector of measurement (or model variable) errors with length n, and $\mathbf{b}$ is the given vector of ideal discharge (or model variable). Obviously, when provided with the value of Equation (27), the variance of the additional noise $\frac{leve{l}^2{\Vert \mathbf{b}\Vert }_2^2}{n}$ can be obtained, and white noise with zero mean can be generated. As is evident from Equations (9) and (20), the difference between the DSRC method and the RDSRC method lies in the term ${\lambda }^{2}I$, which serves to stabilize the solution obtained by the DSRC method. Therefore, only a new flowchart (Figure 2) is needed to illustrate the specific procedures of RDSRC method (using the regularization techniques) which is included in the box “RDSRC” in Figure 3. For general procedures of DSRC method, we recommend reading research by Bao, Si and Qu [32], Si, Bao and Gupta [33].

The synthetic study area consists of five precipitation stations and has an area of 24,000 km². It is assumed that the synthetic basin has uniform parameter fields. The model parameters used in this study are presented in

Table 3. Parameters of the XAJ model used in the numerical example.

Module
TLE and SRP		SOR		FC		MSSR
Parameter	Value	Parameter	Value	Parameter	Value	Parameter	Value
K	0.8	SM	30	CS	0.875	KE	1
WUM	20	EX	1.5	CI	0.925	XE	0.49
WLM	80	KI	0.35	CG	0.995
WDM	30	KG	0.35
C	0.16
IM	0.01
WM	130
B	0.4

. The parameters of TLE and SRP are presented in the same column since these modules are closely related.

To evaluate the robustness of the origin method and the method using regularization tools (Tikhonov regularization and L-curve method), a wide range of error levels ($level\in \{0,0.01,0.02,\dots ,0.7\}$) were considered. For each noise level, 100 perturbation error vectors (flood events) were generated to compare the DSRC method and the RDSRC method (7100 flood events in total). For the present study, the main purpose was to apply the regularization method to stabilize the solutions perturbed by the right-hand side; thus, once fixed, the error vector of free water storage should remain unchanged for all 7100 runs. The error vector of the free water storage was also assumed to be additional white noise with zero mean and the corresponding noise variance was obtained by setting Equation (27) equal to 0.7. To compare the stability of the two methods in detail, four criteria were used. Firstly, the relative error is shown by

$$RE=\frac{{\Vert {\mathbf{S}}_{exact}-\mathbf{S}\Vert }_2}{{\Vert {\mathbf{S}}_{exact}\Vert }_2}$$

(28)

where ${\mathbf{S}}_{exact}$ represents the vector of the exact values of the free water storage, $\mathbf{S}$ represents the vector of the calculated values (obtained by DSRC, RDSRC or free of updating). The Nash–Sutcliffe efficiency value (NSE) [59] is shown by

$$NSE=1-\frac{{\displaystyle \sum _{t=1}^N{[Q_{sim}(t)-Q_{obs}(t)]}^2}}{{\displaystyle \sum _{t=1}^N{[Q_{obs}(t)-\overline{Q_{obs}}]}^2}}$$

(29)

where $Q_{sim}(t)$ and $Q_{obs}(t)$ represent the simulated discharge and the observed discharge at time step t (with time length N) respectively, $\overline{Q_{obs}}$ is the mean of the measured discharge. The Nash–Sutcliffe efficiency value is a widely used measure that serves to compare the performance of a given model to the mean value ($\overline{Q_{obs}}$) of the observed values [60]. A NSE value of 1 indicates a perfect match between the output calculated by the model and the observed data. A NSE value of 0 means that the performance of the model is only as accurate as the mean value. The following criteria are used to evaluate the free water storage, as well as the discharge. Firstly, the mean bias error (MBE) is calculated by

$$MBE=\frac{1}{N}{\displaystyle {\sum }_{i=1}^N(x_{obs,i}-}{x}_{sim,i})$$

(30)

where $x_{obs,i}$ is the observed value at time i and $x_{sim,i}$ is the simulated value at time i. and the root mean square error (RMSE)

$$RMSE=\sqrt{\frac{1}{N}{\displaystyle {\sum }_{i=1}^N{\left(x_{obs,i}-x_{sim,i}\right)}^2}}$$

(31)

2.4. Real Case

2.4.1. Study Area

The Qilijie Basin is located in Fujian Province in southwest China and has an area of 14,787 km². The basin is dominated by a subtropical monsoon climate. The annual precipitation in the basins ranges from 1289 mm to 2451 mm, and the corresponding measured pan evaporation ranges from 735 mm to 895 mm.

2.4.2. Data Basis

Observed data from the outlet station, including daily data and hourly data for the period 1988–1999, were used for parameter calibration and simulation. Precipitation inputs were based on daily data and hourly data from 42 stations. Measured pan evaporation inputs were based on daily data from Songxi station in the middle part of the basin, since Songxi station is the only available station in the study area. Together with the precipitation, the daily measured pan evaporation of Songxi station was disaggregated to hourly time-steps using linear interpolation to serve as the input to the XAJ model. Observed outlet discharge data were based on hourly data from Qilijie station. Figure 4 depicts the location of Qilijie Basin and the stations.

The parameter optimization of the XAJ model has been described by many authors and we do not intend to include it here. For a more detailed description, the references [54,55,56,57,58] are recommended. The initial values of the parameters are given empirically based on the feasible range proposed in Zhao [54]. Daily data were mainly used for parameter calibration in the evapotranspiration module and the runoff production module. By comparing the simulated and observed hydrographs, the parameters were calibrated manually. Thirty-six flood events were selected from the period 1988–1999. The performance of the XAJ model was presented together with the results of the two methods. For simplicity, the parameter fields were assumed to be uniform.

Table 4. Parameters of the XAJ model for Qilijie basin.

Module
TLE and SRP		SOR		FC		MSSR
Parameter	Value	Parameter	Value	Parameter	Value	Parameter	Value
K	1.18	SM	34	CS	0.798	KE	1
WUM	20	EX	1.5	CI	0.9	XE	0.38
WLM	80	KI	0.379	CG	0.995
WDM	50	KG	0.321
C	0.16
IM	0.001
WM	150
B	0.4

lists the values of the main parameters.

3. Results and Discussion

3.1. Synthetic Case

A total of 7100 runs were made under the conditions (parameters, initial conditions and the error vectors) mentioned above. To eliminate the effects of the units and ranges of different variables, for the synthetic case, the criteria RMSE and MBE were replaced by RMSE/abs(ORMSE) (R_RMSE) and MBE/abs(OMBE) (R_MBE), where ORMSE and OMBE represent the RMSE value and the MBE value of the original performance (free of updating), respectively. With this replacement, the “benchmark” information (original performance) is also included. Figure 5 displays both the results of the free water storage and the discharge. The means and standard deviations for RMSE/abs(ORMSE) of the two methods (DSRC and RDSRC) are displayed this figure.

For the free water storage, with an increase in the level of noise, the mean R_RMSE of DSRC increased from 0.89 to 3.17 with the mean being equal to 2.38, while for RDSRC, the mean increased from 0.54 to 2.04, with the mean equal to 1.42. Similar results were obtained for the discharge. There was a significant increase in the R_RMSE of DSRC (from 0.28 to 3.01) with an increase in noise level. However, the R_RMSE of RDSRC has only increased from 0.23 to 1.23. The performance of the XAJ model using the RDSRC method (0.99) is still slightly better than the original performance of the XAJ model when the noise level increases to 0.54, as measured by R_RMSE_Q. To evaluate the performance of the RDSRC method and DSRC method in detail, a similar plot (Figure 6) for the criterion R_MBE is presented.

As expected, satisfactory performances, as measured by R_MBE, were obtained. For the RDSRC method, the R_MBE of the free water storage decreased from −0.22 to −1.34 (with the mean equal to −0.78) and the R_MBE of the discharge decreased from −0.19 to −1.18 (with the mean equal to −0.67). Clearly, there is considerable deterioration in performance when using the DSRC method. The R_MBE of the free water storage ranged from −3.18 to −0.27 with the mean equal to −1.91, and the R_MBE of the discharge ranged from −0.24 to −3.18 with the mean equal to −1.84. An interesting observation from this figure is that for both methods, the R_MBE (and hence, MBE) monotonously decreases when increasing the noise level. The reason for this phenomenon lies in the mathematical equations for recursively calculating the free water storage (S)

$$S(t)=(1-KI-KG)\cdot S(t-1)+(R(t)-RS(t))/FR(t)$$

(32)

with

$$FR(t)=R(t)/PE(t)$$

(33)

where PE is the net rainfall. Descriptions of other variables are presented in Table 2.

The sum of the parameters, KI and KG, is generally set as equal to a small positive value which ranges from 0.7 to 0.8 [54]. Throughout the paper, the sum is set to 0.7, and therefore, Equation (32) is rewritten as

$$S(t)=0.3\times S(t-1)+PE(t)\times (1-RS(t)/R(t))$$

(34)

Obviously, this equation indicates that the value of free water storage is very likely to be relatively small due to the small scaling factor (0.3), even if there is a increment caused by the remainder of R. Specifically, if there is no precipitation for two or more time steps (hours), the value of the free water storage will surely approximate 0 since SM (the upper limit of S) generally ranges from 5 mm to 50 mm. Together with the limits imposed on the free water storage (Equation (25)), it is likely that the idea of updating the free water storage is more applicable for “overestimated” situations, where the simulated discharge is greater than the observed discharge. However, for “underestimated” cases, where the observed discharge is greater than the simulated discharge, the inapplicability of this idea is obvious. For example, assume that the increment obtained by the RDSRC method is −5 and the free water storage at the current time step is 1. The only solution is to project the unconstrained estimation (−5 + 1 = −4) onto the boundary (0) since the estimated free water storage should remain positive. Obviously, this can definitely lead to underestimation of the exact value of the discharge. Therefore, there is a clear tendency for the MBE to decrease when the noise level is increased from 0 to 0.7 due to its intrinsic inability to update the free water storage. The mean and the standard deviation for the relative error (RE) of the free water storage are displayed in Figure 7a. Similar results were obtained for NSE (Figure 7b).

Figure 7a shows that an increase in the noise level has a significant effect of increasing the relative errors for both methods. With an increase in the level of noise, the mean RE of DSRC increased from 1.68 to 5.95, and the mean RE of RDSRC increased from 1.02 to 3.84. In addition, the mean of RDSRC (2.66) was smaller than the value of DSRC (4.46). Figure 7b clearly indicates both methods can obtain an improvement in simulations when the noise level is less than 0.17 ($level<0.17$). For the cases where the noise level ranges from 0.17 to 0.56, only the RDSRC method can obtain better mean NSE values than the original case (without updating). However, for the cases where the level is higher than 0.56, both methods result in complete deterioration. When the noise level increased, the mean NSE of RDSRC decreased from 0.99 to 0.55 and the mean NSE of DSRC decreased from 0.98 to −1.84. Note that similar results are also shown in Figure 5 and Figure 6.

Figure 7a,b indicate that with an increase in observation noise, there was a general tendency for the performance of both methods to deteriorate, due to the ill-posed problem. As expected from the theoretical analysis, the RDSRC method was much more robust than the DSRC method. Note that even without artificial perturbations to the right-hand side (i.e., level = 0), the performance of the RDSRC method was still slightly better than the performance of the DSRC method (smaller RE, smaller RMSE, bigger NSE, MBE closer to 0). The small shift in the results may be caused by the errors of approximation in the DSRC method. These errors will also contaminate the right-hand side and thus, cause deterioration in performance.

An interesting observation from Figure 7b is that there is a clear increment in the standard deviation (STD) with an increase in noise level for both methods. However, the STD obtained by RDSRC only slightly increases when compared to the STD obtained by DSRC (Figure 7b). This may suggest that the updated series obtained by RDSRC are closer to the ideal series than the measurements with additional noise. If the results obtained by the updating methods are closer to the measurements than the ideal flood event, the STD values of NSE will tend to be higher.

An important feature of Tikhonov regularization is that it filters (suppress) the contributions from the right-hand side (b) (contaminated by the noise) corresponding to the small singular values (hence filters out the perturbations), and hence, the solution is less sensitive to the perturbations [51,52,61]. Thus, it is of great importance that the updated series should be closer to the unique ideal flood event than the measurements, since this phenomenon indicates the method has filtered out the noise. This motivates an examination of which series (measured or ideal) is closer to the updated series. Note that the criterion NSE is only used to compare the updated value and the exact value (free of noise) in Figure 7b. Hence, in this study, the NSE is generalized to compare the simulations to the observations (with noise). Figure 8 displays the updated results of both methods using the measured values and the exact values as $Q_{obs}(t)$ in Equation (29), respectively. This examination resulted in four sets of results, listed in Figure 8.

Figure 8 shows that the results obtained by RDSRC were closer to the exact values compared with the measurements, and with an increase in noise level, there was an obvious tendency for the difference between the two lines to increase monotonically. This indicates that the RDSRC is rarely mistaken due to noise (insensitive to the noise). However, the NSE values obtained by DSRC with the exact series were only slightly better than those obtained with the measurements when the noise level was less than 0.32. Furthermore, for cases where the noise level is higher than 0.32, there is a tendency for the NSE values with the ideal series to decrease at a much faster rate than that of the measurements. This indicates that the DSRC method is completely “misled” by noise. Thus, the RDSRC method efficiently filters out the perturbations to the exact right-hand side by filtering out the components of the solution corresponding to small singular values, whereas the DSRC can only weakly filter out the noise, leading to a deterioration in performance when dealing with high noise levels.

3.2. Real Case

An efficient method for comparing the model’s performance obtained by two different methods (i.e., RDSRC and XAJ, DSRC and XAJ) is to present the scatter plot of a criterion for the two methods. Figure 9 shows such a plot of RMSE.

Figure 9 clearly indicates that the RDSRC method has significantly improved the model’s performance (for almost all of the 35 flood events), while for the DSRC method, the “corrected” performance was almost identical to the original performance (free of correction). The RMSE of the RDSRC method only ranged from 106.21 m³/s to 1139.60 m³/s with the mean equal to 398.81 m³/s, which is superior to the original performance (ranging from 117.26 m³/s to 1988.40 m³/s with the mean equal to 755.18 m³/s). However, the RMSE of DSRC method ranged from 218.76 m³/s to 1810.47 m³/s with the mean equal to 753.97 m³/s. To compare the applicability of DSRC and RDSRC in more detail, the results measured by MBE are presented.

Figure 10 illustrates that both RDSRC method and DSRC method tend to overestimate the discharge since the MBE values were negative. As discussed in the numerical example, this phenomenon may be caused by the fact that the calculated free water storage tends to decrease rapidly and approximates the lower limit (0). In this case, the negative correction (increment) obtained by the RDSRC (DSRC) method has no effect on the model’s output. Therefore, the idea of improving the model’s performance by updating the free water storage may break down when the calculated discharge is greater than the observed discharge. These results support the theoretical findings of the numerical example. Nevertheless, there is a significant improvement in the performance of RDSRC when compared to the results of the DSRC method and XAJ model.

Finally, the results measured by the criterion NSE are briefly reported. As expected, Figure 11b clearly shows that the RDSRC method has achieved considerably satisfactory results. However, for DSRC, the performance was almost identical to the XAJ model (without correction) and in some cases, DSRC resulted in deterioration (Figure 11a). RDSRC yielded satisfactory results for most of the flood events, with all the flood events having a mean NSE of 0.92, while for DSRC, the mean NSE was only 0.70, which is even less than the mean NSE in the XAJ model (0.74).

Together with the results of the numerical example, it was shown that the RDSRC method obtained an overall satisfactory (stable) performance. However, the DSRC method typically breaks down when the noise level is high and can even lead to deterioration in the model’s performance in real cases. This may be caused by the fact that for real cases, perturbations to the right-hand side are much more complicated than in the numerical example. In the previous example, some errors (i.e., parameter uncertainty, model inadequacy, etc.) were assumed to be perfect. However, in real cases, these errors are actually unavoidable and therefore, contribute to perturbations to the right-hand side. Therefore, the “worse” right-hand side can result in deterioration in performance when using DSRC.

Unfortunately, our capability to quantify model error still shows considerable deficits [62]. With these considerations, we note that the main reason for the significant deterioration in the real case performance of DSRC is the “terrible” right-hand side with additional perturbations. The robust performance of the RDSRC method (compared to the performance of DSRC) suggests that the regularization technique efficiently filters out the components of the solution corresponding to small singular values and hence, the solution is less sensitive to the perturbations. Therefore, the applicability and robustness of RDSRC method were demonstrated by both the numerical example and the real case study using a conceptual rainfall-runoff model.

Before entering into the conclusion, it should be pointed out that despite the robust and satisfactory performance of this model, it has been demonstrated that the inability to update the free water storage can deteriorate its performance. The fundamental reason for this inability is the constraints on the model variables due to their physical limitations. This problem has also been reported in the context of nonlinear state estimation [63,64,65]. Therefore, future work is needed to incorporate this important information into the proposed method (constraints on model variables). Besides the intrinsic problem of the free water storage of XAJ model, it should be pointed out that both the DSRC and RDSRC methods can only update a single variable; this can indeed result in significant deviation between the simulations and the observations in some cases [34]. The reason for this flaw is the assumption n ≥ m in Equation (5) of the mathematic algorithms (the least squares method). This assumption mainly serves to avoid the rank-deficient problem and therefore, poses a constraint on the number of updated variables. As a consequence, an unavoidable problem will arise when updating the free water storage; for example, the errors in the soil moisture can never be eliminated by updating the free water storage but contribute to perturbations on the right-hand side. Hence, for real cases, the use of DSRC possibly leads to a more unstable performance than in the numerical examples. Interestingly, even though the methods can only update a single variable, the effect of updating a specified variable can be transferred to other variables which are calculated based on the value of the specified variable. For example, a shift in S will surely cause changes in RS, RI, RG, etc., but it has no effect on the soil moisture (W). In considering this, it is likely that choosing variables from the module (box) with smaller numbers (Figure 1) will have an effect on more variables. However, as discussed above, the main idea of the DSRC method is to approximate the state-to-output transition equations through first-order linearization, and therefore, the nonlinearity of different variables is of great importance. In general, variables which are closer to the output of a model tend to have weaker nonlinearity. Unfortunately, there is little literature regarding an overall evaluation of the nonlinearity of variables in the XAJ model. An investigation based on correlation analysis was proposed by Lü, Hou, Horton, Zhu, Chen, Jia, Wang and Fu [58]. The results showed that the free water storage has strong linear correlations with the output discharge. Consequently, it is likely that a good method for choosing the updated variable for RDSRC (DSRC) must seek a good balance between nonlinearity and the number of variables to be corrected. These considerations, such as the multivariable updating capability, the principles for choosing a variable, and incorporating the constraints of model variables, should be addressed in future work.

4. Conclusions

In this study, a new version of the DSRC method using regularization techniques was proposed and applied to improve the performance of the XAJ model. A theoretical analysis was conducted to investigate the instability of the DSRC method. Both the DSRC method and the RDSC method were tested with a numerical example and a real case study.

The analysis showed that under certain assumptions, the DSRC method was almost identical to the numerical solution of the Fredholm integral equation of the first kind, indicating that the DSRC method has an intrinsic ill-posed characteristic which may lead to meaningless solutions. The new method incorporates regularization techniques to overcome the instability of DSRC. The results of numerical example showed that the RDSRC method can efficiently filter out the perturbations to the right-hand side while the DSRC method can hardly distinguish and filter out the noise. Theoretical considerations show that the perturbations to the right-hand side tend to be more complicated when moving from numerical examples to real cases, thereby causing a poor updating performance. These theoretical considerations were supported by the case study. The regularization technique can greatly help enhance the stability and performance of updating in real cases.

An analogous analysis of other studies showed that the same problems occur when updating the soil moisture [34]. This study presented the reason for the unstable performance and proposed a new method which has been demonstrated to be able to handle the unstable problem. Ultimately, further work is needed to achieve multi-variable correction. A method for incorporating the additional information (constraints on model variables) and a general criterion for choosing the variable(s) to be updated are also needed.