*Article* **Improvement in Ridge Coefficient Optimization Criterion for Ridge Estimation-Based Dynamic System Response Curve Method in Flood Forecasting**

**Kexin Liu 1,2,\*, Weimin Bao 1, Yufeng Hu 2, Yiqun Sun 1, Dongjing Li 2, Kuang Li <sup>2</sup> and Lili Liang <sup>3</sup>**


**Abstract:** The ridge estimation-based dynamic system response curve (DSRC-R) method, which is an improvement of the dynamic system response curve (DSRC) method via the ridge estimation method, has illustrated its good robustness. However, the optimization criterion for the ridge coefficient in the DSRC-R method still needs further study. In view of this, a new optimization criterion called the balance and random degree criterion considering the sum of squares of flow errors (BSR) is proposed in this paper according to the properties of model-simulated residuals. In this criterion, two indexes, namely, the random degree of simulated residuals and the balance degree of simulated residuals, are introduced to describe the independence and the zero mean property of simulated residuals, respectively. Therefore, the BSR criterion is constructed by combining the sum of squares of flow errors with the two indexes. The BSR criterion, L-curve criterion and the minimum sum of squares of flow errors (MSSFE) criterion are tested on both synthetic cases and real-data cases. The results show that the BSR criterion is better than the L-curve criterion in minimizing the sum of squares of flow residuals and increasing the ridge coefficient optimization speed. Moreover, the BSR criterion has an advantage over the MSSFE criterion in making the estimated rainfall error more stable.

**Keywords:** flood forecasting; error correction; residual property; ridge coefficient criterion

#### **1. Introduction**

Flood forecasting, an important non-structural measure, plays an important role in regional flood control, flood warning, risk decision making, etc. [1–3]. The hydrological model simplifies and conceptualizes the flood process with a set of equations aiming at obtaining the outlet flow. However, not all problems can be solved with such a model as the flood forecasting accuracy is often hampered and influenced by many error factors existing in the hydrologic system, including the errors in the model inputs, the errors in the model initial condition, the errors in the model simplification and the errors in the model parameters. Therefore, many scholars have devoted themselves to the research of error correction methods. For example, the autoregressive (AR) model estimates the flow error existing in a certain forecasting period by using the correlation of error series, and it was later developed into improved methods such as the recursive autoregressive model and the forgetting factor recursive autoregressive model [4–6]; Kalman filtering (KF) technology is widely used to update hydrological element time series in flood forecasting, and many improved types have been gradually formed, including the extended Kalman filter (EKF) [7] and the ensemble Kalman filter (EnKF) [8]. Data assimilation technology has also shown satisfying results in improving the prediction accuracy of models, including dynamic identifiability analysis (DYNIA) [9], and the Bayesian recursive estimation technique

**Citation:** Liu, K.; Bao, W.; Hu, Y.; Sun, Y.; Li, D.; Li, K.; Liang, L. Improvement in Ridge Coefficient Optimization Criterion for Ridge Estimation-Based Dynamic System Response Curve Method in Flood Forecasting. *Water* **2021**, *13*, 3483. https://doi.org/10.3390/w13243483

Academic Editor: Renato Morbidelli

Received: 5 November 2021 Accepted: 1 December 2021 Published: 7 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

(BaRE) [10]. Robust theory [11] and the comprehensive correction method [12] are also used in flood forecasting error correction. However, error correction technology still needs to be improved; for example, the autoregressive model assumes that the error series has a linear correlation, but its performance near the flood peak is often not satisfying [12]. Additionally, some studies have also shown that the EnKF does not perform very well when the error structure is far away from a Gaussian distribution [13]. Scholars are still trying to solve the above problems.

A new error correction method called the dynamic system response curve (DSRC) method has been proposed by Bao et al. [14]. This method constructs a feedback model which is conceptualized on updating the hydrologic element series by tracing back to the source of the error. With the help of the first-order Taylor linearization to approximate the hydrologic model, error correction is achieved by solving the corresponding equations using the least square method. The DSRC method was initially applied to correct single hydrological elements, including runoff [15], rainfall [16] and model state variables [17]. Then, it was used to correct several hydrological elements comprehensively [18]. However, some studies [17–20] found that the correction results are not always stable, reflected in the excessive correction of hydrological element series and none-smooth simulated flow hydrographs. To solve the above problems, the DSRC-R method was developed by Si et al. [19] from the point of the regularization known as ridge estimation, and this method has improved the stability of correction results to some degree. Nevertheless, the selection criterion of the ridge coefficient still needs further study and improvement. Previous studies [17–19] often chose the ridge coefficient based on the L-curve criterion [21]. However, it has been found in practice that two aspects still need attention. One is that the L-curve seems insufficient to reflect the properties of model-simulated residual errors including the independence and the zero mean property of simulated residuals, which hinders the performance of DSRC-R in some cases; the other is that the application of the L-curve criterion takes a long time, which is not conducive to the real-time performance of flood forecasting. The L-curve criterion involves derivative calculation, and the difference method makes the operation efficiency lower.

Therefore, in this paper, we analyze the L-curve criterion and introduce the concepts of the random degree of simulated residuals and the balance degree of simulated residuals to describe the properties of the model-simulated residuals, and then a new criterion called the balance and random degree criterion considering the sum of squares of flow errors (BSR) is proposed. The new criterion takes the independence and the zero mean property of simulated residuals into account, which is conducive to obtaining a ridge coefficient in line with the statistical characteristics of residuals. The new criterion does not involve derivative calculation; thus, it can greatly shorten the search time of the ridge coefficient in optimization, improve the operational efficiency and enhance the real-time flood forecasting performance.

#### **2. Methodology**

#### *2.1. DSRC Method*

The main idea behind the DSRC method is that it firstly retrieves the rainfall errors from the outlet flow errors, then updates the rainfall series and finally reruns the model with the updated rainfall series. In this method, given a hydrological model *Q* = *Q*(*P*) that generates outlet flow *Q* as a function of rainfall *P*, the variation process of the outlet flow change by the unit perturbation in rainfall is called the system response curve. Based on this, multi-time system response curves form the system response matrix *S*, and this matrix sets up the relation between rainfall errors Δ*P* and flow errors Δ*Q*; then, the estimation of rainfall errors Δˆ *PLS* can be computed via the least square method. In this round, the input rainfall series is updated with Δˆ *PLS*, and then the model is rerun with the updated rainfall series to correct the forecasting results. In this study, the DSRC method was combined with the Xinanjiang (XAJ) model, a hydrological model constructed by Professor Zhao Renjun of Hohai University which is widely used in flood forecasting in humid areas of China [22]. Additionally, here, we mainly talk about the calculation process of the DSRC method shown in Figure 1 rather than the XAJ model, as this paper will introduce it in Section 3. Additionally, the theoretical derivation of the DSRC method can be gained from [16].

**Figure 1.** Schematic diagram of the DSRC method. *QC*, *QO* and *Q <sup>C</sup>* are model-simulated flow, observed flow and updated model-simulated flow, respectively; Δ*Q* is the model-simulated deviation series of the outlet flow; Δˆ *P* is the estimated rainfall error series; *PO* is the initial rainfall series; *P* is the updated rainfall series; and *E* is pan evaporation.

According to [16], the rainfall error estimation Δˆ *PLS* is expressed as the following Equation (1):

$$
\Delta \hat{P}\_{LS} = \left(\mathbf{S}^T \mathbf{S}\right)^{-1} \mathbf{S}^T \Delta \mathbf{Q} \tag{1}
$$

where Δ*Q* is the model-simulated deviation series of the outlet flow, Δ*Q* = *QO* − *QC*; *QO* = [*QO*,1, *QO*,2, ··· *QO*,*M*] *<sup>T</sup>* is the observed flow series; and *QC* = [*QC*,1, *QC*,2, ··· *QC*,M] *T* is the simulated flow series computed from the observed rainfall series.

In Equation (1), *S* is the system response matrix defined as

$$S = \begin{bmatrix} \frac{\partial Q\_1(P)}{\partial p\_1} & \frac{\partial Q\_1(P)}{\partial p\_2} & \cdots & \frac{\partial Q\_1(P)}{\partial p\_N} \\ \frac{\partial Q\_2(P)}{\partial p\_1} & \frac{\partial Q\_2(P)}{\partial p\_2} & \cdots & \frac{\partial Q\_2(P)}{\partial p\_N} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial Q\_M(P)}{\partial p\_1} & \frac{\partial Q\_M(P)}{\partial p\_2} & \cdots & \frac{\partial Q\_M(P)}{\partial p\_N} \end{bmatrix} \tag{2}$$

where *p*1, *p*<sup>2</sup> ··· *pn* are the initial rainfall values; and the indices *M* and *N* represent the lengths of observed flow and rainfall, respectively (*<sup>M</sup>* <sup>≥</sup> *<sup>N</sup>*). *<sup>∂</sup>Qi*(*P*) *<sup>∂</sup>pj* represents the influence of the *j*-th rainfall on the *i*-th outlet flow. When *i* < *j*, it is obvious that *<sup>∂</sup>Qi*(*P*) *<sup>∂</sup>pj* = 0 because rainfall does not affect the outlet flow that occurs before it. *<sup>∂</sup>Qi*(*P*) *<sup>∂</sup>pj* is generally obtained by the difference, which is *<sup>∂</sup>Qi*(*P*) *<sup>∂</sup>pj* <sup>=</sup> *<sup>Q</sup>*(*p*1,...,*pj*+Δ*p*,...,*pN*)−*Q*(*p*1,...,*pj*,...,*pN*) <sup>Δ</sup>*<sup>p</sup>* .

Accordingly, the updated rainfall series *P* can be expressed as Equation (3):

$$P' = P\_O + \Lambda \hat{P} \tag{3}$$

where *P* is the updated rainfall series; *PO* is the initial rainfall series; and Δˆ *P* is the rainfall error estimation series. In the DSRC method, Δˆ *P* is replaced by Δˆ *PLS* mentioned above.

In order to improve the forecasting accuracy, the updated rainfall series *P* is introduced into the hydrological model for recalculation, and then the updated model-simulated flow series *Q <sup>C</sup>* is obtained by Equation (4):

$$\mathcal{Q}'\_{\mathbb{C}} = \mathcal{Q}(P') \tag{4}$$

where *Q <sup>C</sup>* = [*Q <sup>C</sup>*,1, *Q <sup>C</sup>*,2, ··· , *Q <sup>C</sup>*,*N*] *<sup>T</sup>* is the updated model-simulated flow series.

#### *2.2. DSRC-R Method*

Correction results from the DSRC method are sometimes unstable. Relevant studies [19,20] have pointed out that the DSRC method is prone to be ill-conditioned, which generates unstable results when this method is applied to small basins, or when the length of flow information is short. Therefore, Si et al. [19] combined the DSRC method with the ridge estimation method and proposed the DSRC-R method, which is more robust than the DSRC method. In previous studies, the ridge coefficient *β* was often selected via the L-curve criterion [17–19]; however, the correction results were not always stable. Thus, the criterion for obtaining the appropriate ridge coefficient *β* in the DSRC-R method needs further study. Here, we directly provide the formula of rainfall error estimation Δˆ *PRE* as the following Equation (5). For more details about the derivation process of DSRC-R, please refer to [19].

$$
\Delta \hat{P}\_{RE} = \left(\mathbf{S}^T \mathbf{S} + \beta I\right)^{-1} \mathbf{S}^T \Delta \mathbf{Q} \tag{5}
$$

where *β* is the ridge coefficient; *I* is the identity matrix; and Δˆ *PRE* is the rainfall error estimation series of the DSRC-R method.

Δˆ *PRE* can be introduced into Equation (3) to update the rainfall series, and then the model can be rerun with *P* to correct the forecasted flow. The flow chart of the DSRC-R method is shown in Figure 2.

**Figure 2.** Flow chart of the DSRC-R method. *Si* is the system response curve of the *i*-th rainfall, that is, the *i*-th column of matrix S.

#### *2.3. L-Curve Criterion*

The ridge coefficient *β* is significantly important in the DSRC-R method. In some previous studies [17–19], the L-curve criterion was used to determine the value of *β*. According to [23], the main idea behind the L-curve criterion in selecting the appropriate ridge estimation coefficient *β* in the DSRC-R method can be summarized as: the balance of lg*<sup>Q</sup>*(*PO* <sup>+</sup> <sup>Δ</sup><sup>ˆ</sup> *PRE*) <sup>−</sup> *QO* and log <sup>Δ</sup><sup>ˆ</sup> *PRE*, and the appropriate ridge coefficient lies at the corner of the curve, usually with the largest curvature. That is, the technique to find the appropriate ridge coefficient can be expressed as Equation (6).

$$\max\_{\boldsymbol{\mathcal{S}}} \left( \frac{|\boldsymbol{f}' \mathbf{g}'' - \boldsymbol{f}'' \mathbf{g}'|}{\left[ \left( \boldsymbol{f}' \right)^2 + \left( \mathbf{g}' \right)^2 \right]^{3/2}} \right) \tag{6}$$

where *<sup>f</sup>*(*β*) = log <sup>Δ</sup><sup>ˆ</sup> *PRE* 2 2; *<sup>g</sup>*(*β*) = lg*<sup>Q</sup>*(*PO* <sup>+</sup> <sup>Δ</sup><sup>ˆ</sup> *PRE*) − *QO* 2 2; and •<sup>2</sup> <sup>2</sup> is the modular square.

It can be proved (see Appendices A and B) that *f*(*β*) and *g*(*β*) can be expressed as Equations (7) and (8).

$$f(\beta) = \log[\sum\_{i=1}^{n} \left(\frac{k\_i}{\lambda\_i + \beta}\right)^2] \tag{7}$$

$$\lg(\beta) = \lg \left\| Q(P\_O + \sum\_{i=1}^n \frac{k\_i}{\lambda\_i + \beta} v\_i) - Q\_O \right\|\_2^2 \tag{8}$$

where *<sup>λ</sup>i*(*<sup>i</sup>* = 1, ··· , *<sup>N</sup>*) is the eigenvalue of the matrix *<sup>S</sup>TS*; *vi*(*<sup>i</sup>* = 1, ··· , *<sup>N</sup>*), orthogonal to each other, is the unit eigenvector corresponding to *λi*(*i* = 1, ··· , *N*); and *ki*(*i* = 1, ··· , *n*) is a group of coefficients that enable *<sup>S</sup>T*Δ*<sup>Q</sup>* to be linearly expressed by *vi*(*<sup>i</sup>* = 1, ··· , *<sup>N</sup>*), that is, *<sup>S</sup>T*Δ*<sup>Q</sup>* <sup>=</sup> *<sup>n</sup>* ∑ *i*=1 *kivi*.

The L-curve criterion has a good effect on selecting the ridge coefficient in the DSRC-R method, but there are still some problems that are worthy of attention. First, the result is sometimes unsatisfactory. The reason is that the L-curve criterion seems to insufficiently reflect the properties of the model-simulated residuals, although it pays attention to the balance of lg*<sup>Q</sup>*(*PO* <sup>+</sup> <sup>Δ</sup><sup>ˆ</sup> *PRE*) <sup>−</sup> *QO* and log <sup>Δ</sup><sup>ˆ</sup> *PRE*. Second, the optimization of the ridge coefficient consumes too much time. The L-curve criterion involves first-order and secondorder derivatives, as shown in Equation (6). Additionally, the explicit expression of *g*(*β*) cannot be obtained at present, meaning its derivative can only be obtained by a difference method; thus, it will take a long time and is not conducive to the real-time performance of flood forecasting.

#### *2.4. New Optimization Criterion (BSR)*

In previous studies [17–19], the L-curve criterion was generally adopted to find the suitable ridge estimation coefficient *β* in the DSRC-R method. Nevertheless, the L-curve criterion has some shortcomings such as insufficient consideration of model-simulated residuals, imperfect utilization of information and huge consumption of time. Therefore, this study takes the properties of the model-simulated residuals into consideration and then explores a new optimization criterion which is more suitable for the DSRC-R method.

For any model, it is always expected that the simulated residual series satisfies the zero mean property and non-correlative statistical property; in other words, let the mean of the residual series and correlation coefficient be as small as possible. This shows us that the criterion for determining the ridge coefficient should consider Equations (9) and (10). We use Equation (9) to express the zero mean property of the residuals, and this indicator is called the balance degree of simulated residuals (BDSR). We use Equation (10) to express the correlation of the residual series, and its reciprocal is called the random degree of simulated residuals (RDSR), which is shown in Equation (11). Additionally, RDSR indicates

the independence of the residuals, and we generally expect a smaller value of the BDSR indicator and a bigger value of the RDSR indicator.

$$BDSR = abs \left[ \sum\_{i=1}^{m} \left( Q\_{O,i} - Q\_{C,i}' \right) \right] \tag{9}$$

$$r\_{\varepsilon} = abs[\frac{\sum\_{i=1}^{m-1} (\Delta Q\_i' - \overline{\Delta Q\_X})(\Delta Q\_{i+1}' - \overline{\Delta Q\_Y})}{\sqrt{\sum\_{i=1}^{m-1} (\Delta Q\_i' - \overline{\Delta Q\_X})^2 \sum\_{i=1}^{m-1} (\Delta Q\_{i+1}' - \overline{\Delta Q\_Y})^2}}] \tag{10}$$

$$RDSR = \frac{1}{r\_c} \tag{11}$$

where *abs* represents the absolute value sign; *re* is the absolute value of the correlation coefficient of adjacent residuals;Δ*Q <sup>i</sup>* = *QO*,*<sup>i</sup>* − *Q C*,*i* ; Δ*QX* = <sup>1</sup> *m*−1 *m*−1 ∑ *i*=1 Δ*Q i* ; and Δ*QY* =

$$\underset{\stackrel{m-1}{\underline{m}-1}}{\overset{1}{\underline{m}-1}} \underset{\stackrel{\underline{n}-1}{\underline{\qquad}}}{\Delta Q} \mathbf{Q}'\_{i+1}.$$

Traditional methods often take the minimum sum of squares of flow errors (MSSFE) as the objective function for calibration parameters, but the effective information contained in this method is not sufficient to obtain the ridge coefficient. The reason is that the derivation process of the DSRC-R method utilizes the least square method, and when the system is linear, the least sum of squares of flow errors equals the least square method; thus, the value of *β* should be zero in this circumstance. Although the DSRC-R method belongs to non-linear system inversion methods, the value of *β* still has a decreasing tendency, and this is not conducive to the stability of the method. Therefore, we need to further excavate more useful information in the simulated errors.

In this paper, we take the independence and the zero mean property of simulated residuals into account, combine these two points with the sum of squares of flow errors (SSFE) and lastly explore a new criterion called the balance and random degree criterion considering the sum of squares of flow errors (BSR criterion). Generally, we hope to find a large value of RDSR which is more consistent with the property of the residuals and is conducive to avoiding system errors; moreover, we hope to find a small value of BDSR which can satisfy the zero mean property of residuals and can decrease the flood volume errors. Based on this, considering the SSFE indicator, we propose the BSR criterion, the mathematical form of which is provided in Equation (12). Overall, the new BSR criterion considers the traditional sum of squares of residuals; furthermore, it is possible to find the ridge coefficient which satisfies the properties of flow residuals.

$$\min\_{\beta} \frac{(BDSR + 1)SSFE}{RDSR} \tag{12}$$

where "+1" is used to avoid the value of BDSR being zero; and *SSFE* <sup>=</sup> *<sup>m</sup>* ∑ *i*=1 *QO*,*<sup>i</sup>* − *Q C*,*i* 2 .

Since the BSR criterion pays more attention to the independence and the zero mean property of simulated residuals than the L-curve criterion, it is more likely to select a ridge coefficient *β* which satisfies the properties of residuals, and thus a better performance can be achieved with the DSRC-R method; moreover, the BSR criterion does not involve derivation calculation, meaning it can improve operational efficiency and save much more time.

In real-time flood forecasting, in order to quickly obtain an appropriate *β*, we need to utilize an automatic optimization method. Additionally, in this paper, we adopt the particle swarm optimization algorithm, which was first proposed by Kennedy and Eberhart and constructed on the concept of mimicking the social behavior of birds [24–27]. This algorithm has been widely used in many types of optimization problems. For more details about particle swarm algorithms, please refer to [24].

#### *2.5. The Entire Research Process*

In order to make the whole research process more clear, we created a flow chart, as shown in Figure 3. The figure shows the entire research process including the proposal of the BSR criterion, the research of a synthetic case and a real case and the comparison of three criteria (BSR, L-curve and MSSFE).

**Figure 3.** The research flow chart showing the entire research process.

#### **3. Case Study**

This research projecting synthetic and real-data studies aimed at comparing the performance of the DSRC-R method under three criteria which include the BSR criterion, the L-curve criterion and the MSSFE criterion.

#### *3.1. Model Description*

The selected hydrological model in this research is the Xinanjiang (XAJ) model by Professor Zhao Renjun of Hohai University [22], which is one of the most widely used conceptual hydrological models in China. The XAJ model, including the inputs of observed precipitation as well as pan evaporation, and the outputs of forecasted flow as well as evaporation, can be used in different spatial and temporal scales and be divided into four layers: the first layer utilizing the three-layer evapotranspiration (TLE) model to realize basin evaporation; the second layer utilizing the saturated runoff production (SRP) model to realize the basin runoff production; the third layer utilizing the free water storage model to realize runoff separation; and the fourth layer utilizing the linear reservoir method and the Muskingum method to realize the basin flow concentration. When applying the XAJ model, firstly, divide the basin into several sub-basins and then compute the runoff and outlet flow in every sub-basin; lastly, gather the flow of each sub-basin at the outlet of the basin. For

more details about the XAJ model, please refer to [22]. The structure of the XAJ model is shown in Figure 4, and the meaning of parameters in each layer is shown in Table 1.

**Figure 4.** Schematic diagram of the Xinanjiang (XAJ) model. TLE represents the three-layer evapotranspiration model; SRP represents the saturated runoff production layer; SOR represents the runoff separation layer; and FC represents the flow concentration layer.



#### *3.2. Synthetic Case*

The synthetic case was to design a typical artificial basin whose catchment area, station distribution, model parameters, hydrological features of the basin condition, "observed" precipitation, evaporation, "observed" outlet flow, error factors and any other information about the basin are all known in order to compare different schemes expediently. The synthetic case included 10,000 synthetic precipitations and corresponding floods in order to compare the performance of the DSRC-R method under three criteria (the BSR criterion, the L-curve criterion and the MSSFE criterion).

A major point of the synthetic case was to obtain the "observed" flow. Here, we referred to [16] and utilized Equation (13) to obtain it. By using different *PO* and Δ*P*, we can obtain a different *QO*.

$$Q\_{\mathcal{O}} = Q(P\_{\mathcal{O}} + \Delta P) + \varepsilon \tag{13}$$

where *QO* is the "observed" flow series; *PO* is the initial precipitation series; Δ*P* is the given error series, and each value in the series Δ*P* does not exceed 30% of the corresponding value in *PO*; and *e* is Gaussian white noise which cannot exceed 5% of the initial value in this study.

This synthetic case assumed that the basin area is 1000 km2 and there are 8 precipitation stations in the basin. The value of parameters from each layer is shown in the following Table 2.



#### 3.2.1. Data

A major point of the synthetic case was to generate different initial series of precipitations *PO*. In order to increase the diversity of *PO*, we applied the following method. Firstly, we chose 55 typical areal precipitation processes from a real basin, then transformed the position of the rainfall peak in each precipitation and eventually formed 500 synthetic typical precipitation processes. When generating synthetic rainfall, we selected a synthetic typical precipitation, randomly adjusted each rainfall period ranging less than 30% of the typical rainfall and then obtained the proportion of each time interval of the rainfall series; then, we randomly generated the total rainfall and allocated it to each time interval according to the above proportion, thus forming one synthetic precipitation, that is, one initial series of precipitation *PO*. Then, we introduced *PO* and the given Δ*P* into Equation (13) to obtain the "observed" flow series. In this case, we constructed 10,000 synthetic precipitations and corresponding floods, and the total rainfall of each flood ranged from 10 to 200 mm.

#### 3.2.2. Statistical Indicators

In the synthetic case, what we consider most is the performance of the DSRC-R method under different criteria rather than the contrast of the results between the DSRC and DSRC-R methods, as was accomplished in [19]. The criteria include the BSR criterion, the L-curve criterion and the MSSFE criterion. The relevant statistical indicators include the relative error of flood peak (RPF), relative error of runoff depth (RRD), Nash–Sutcliffe effiency coefficent (NSE), time needed to update a flood (TU) and root mean square error (RMSE). RMSE can be utilized to evaluate the robustness of the DSRC-R method under different criteria, and this index was applied in [17]. The smaller the value of the RMSE indicator, the more robust the DSRC-R method will be.

The statistical indicators between synthetic cases and real-data cases are different. In synthetic cases, RMSE is one of the indicators; however, it is not covered in real-data cases as there is no way to obtain the "true" precipitation. The mathematical definitions of each statistical indicator are expressed as follows:

Relative error of flood peak (RPF):

$$\text{RPF} = (Q\_{\text{OP}} - Q\_{\text{CP}}) / Q\_{\text{OP}} \times 100\% \tag{14}$$

Relative error of runoff depth (RRD):

$$\text{RRD} = (R\_O - R\_C) / R\_O \times 100\% \tag{15}$$

Nash–Sutcliffe effiency coefficent (NSE):

$$\text{NSE} = 1 - \frac{\sum\_{i=1}^{N} \left(Q\_{C,i} - Q\_{O,i}\right)^2}{\sum\_{i=1}^{N} \left(Q\_{O,i} - \overline{Q}\_{O}\right)^2} \tag{16}$$

Root mean square error (RMSE):

$$\text{RMSE} = \sqrt{\frac{1}{n} \sum\_{i=1}^{n} \left( p\_{T,i} - p\_i' \right)^2} \tag{17}$$

where *QCP* is the forecasted value of the flood peak; *QOP* is the observed value of the flood peak; *RC* is the forecasted depth of runoff; *RO* is the observed depth of runoff; *QO* is the average value of flow; *QC*,*<sup>i</sup>* and *QO*,*<sup>i</sup>* are the forecasted flow and observed flow in the *i*-th time interval; and *pT*,*<sup>i</sup>* and *p <sup>i</sup>* are the original precipitation and the updated precipitation in the *i*-th time interval. In the synthetic case, *pT*,*<sup>i</sup>* = *pO*,*<sup>i</sup>* + Δ*pi*, and *p <sup>i</sup>* <sup>=</sup> *pO*,*<sup>i</sup>* <sup>+</sup> <sup>Δ</sup>ˆ*pi* . The larger the NSE (NSE ≤ 1), the higher the forecasting accuracy, the smaller the RMSE and the more robust the DSRC-R method will be.

#### 3.2.3. Computational Process of DSRC Method and DSRC-R Method

The mechanism of the DSRC method is that it firstly utilizes the error information of the outlet flow to invert and estimate the rainfall error, then updates the original rainfall series and lastly reruns the model with the updated rainfall series to correct the forecasting result. The specific steps of the DSRC method are as follows:


According to relevant research [19,20], when the flow data are insufficient, the DSRC method will tend to be unstable, characterized by wide fluctuations in the precipitation error estimated series and potentially oscillation, which will influence the flow correction effect. Therefore, a more robust method, DSRC-R, is proposed through combination with the ridge estimation method. Although the DSRC-R method ensures the stability of the error estimate, the ridge coefficient varies with floods. Thus, the ridge coefficient selection should be in accordance with the instant data. The steps of ridge coefficient optimization are as follows:


In order to evaluate the correction effect of the three criteria, the synthetic case not only includes some flow indicators (RPF, RRD, NSE) but also includes some precipitation indicators (RMSE). Additionally, RMSE is proposed to quantitatively describe the robustness of the method and the oscillation phenomena. The smaller the RMSE, the more robust the precipitation error estimation will be. RMSE has been applied in previous studies and has been successful. However, in the actual case, due to the inability to acquire the "true" value of rainfall, the performance of the three criteria cannot be evaluated by the RMSE indicator.

#### 3.2.4. Results and Discussion

The performance of the DSRC-R method on 10,000 synthetic floods under three criteria is shown in Table 3, and one typical synthetic flood is shown in Figure 5. As can be seen in Table 3, although the three criteria (BSR, L-curve and MSSFE) have a certain effect, some differences still exist. In terms of the indicators RPF and RRD, the BSR criterion has the best performance (1.90% and 1.01%). The MSSFE criterion takes second place in this regard (1.97% and 1.27%), and the L-curve criterion has the worst performance (2.19% and 1.30%). In terms of the indicator NSE, the BSR criterion and MSSFE criterion have the same result, with a value of 0.999, which outnumbers that of the L-curve criterion, with 0.001. In terms of operational efficiency, the L-curve criterion consumes much more time, where the average TU of a flood is 12.21 s. The value of the ridge coefficient *β* under the three criteria has a big difference, and *β* tends to be smaller (average value of 64.26) when it applies the MSSFE criterion. In terms of the indicator RMSE, the value under the BSR criterion is 0.759, which is significantly less than the value (1.040) under the MSSFE criterion, indicating that the BSR criterion is more conducive to improving the robustness of the DSRC-R method than the MSSFE criterion.


**Table 3.** The results of the synthetic case.

<sup>1</sup> The values of indicators (RPF, RRD, NSE, RDSR, BDSR, TU, *β*, RMSE) in the table are the average values of 10,000 synthetic floods.

**Figure 5.** The performance of the DSRC-R method in a typical flood under three criteria. (**a**) The contrast between the real value of the rainfall error and the estimation value under the MSSFE criterion; (**b**) the contrast between the real value of the rainfall error and the estimation value under the L-curve criterion; (**c**) the contrast between the real value of the rainfall error and the estimation value under the BSR criterion; (**d**) the contrast of a typical flood forecasted flow, where *QO* is the "observed" flow, *QC* is the forecasted flow, *QL* is the updated forecasted flow under the L-curve criterion, *QSS* is the updated forecasted flow under the MSSFE criterion and *QBSR* is the updated forecasted flow under the BSR criterion.

As it is depicted in Figure 5a–c, the estimated value of the rainfall error under the BSR criterion is the closest to the "true" value, the total difference is 0.8 mm (1.2 mm errors under the L-curve criterion and 0.9 mm errors under the MSSFE criterion), the correlation coefficient reaches 0.962 (0.875 under the L-curve criterion and 0.896 under the MSSFE criterion) and the spots are evenly distributed on both sides of the 1:1 line. All of the above contribute to the DSRC-R method achieving the best performance under the BSR criterion. As it is shown in Figure 5d, the indicators RPT and RRD have values of 0.9% and 0.6%, respectively, under the BSR criterion. RPT has a value of 2.7% and RRD a value of 1.4% under the MSSFE criterion. Lastly, RPT has a value of 4.5% and RRD a value of 3.8% under the L-curve criterion. In a typical flood, the optimal value of *β* under the MSSFE criterion is 2.47, which is significantly less than 651.33 under the BSR criterion. This result shows that the instability of the method is not sufficiently alleviated. Therefore, the points in Figure 5a are scattered. The L-curve criterion does not fully consider the properties of the simulated residuals, meaning the flow correction result in Figure 5d is not satisfactory.

Compared with the L-curve criterion, the BSR criterion improves the performance of the DSRC-R method. This is because the BSR criterion takes more account of the properties

of model-simulated residuals, including the mutual independence between residuals and the zero mean property of residuals. As is shown in Table 3, the average BDSR under the BSR criterion is 66.9 m3/s, which is significantly less than 80.8 m3/s under the Lcurve criterion, indicating that the former can better reflect the zero mean property of the simulated residuals; moreover, the average RDSR under the BSR criterion is 12.88, which is significantly greater than 5.57 under the L-curve criterion, indicating that the former is more conducive to the mutual independence of simulated residuals. The above analysis shows that the BSR criterion is more conducive to the optimization of the ridge coefficient, in line with the properties of simulated residuals, than the L-curve criterion. Therefore, the BSR criterion improves the performance of the DSRC-R method. In terms of operation efficiency, the BSR criterion takes less time than the L-curve criterion because the former does not involve derivative calculation, while the latter involves derivative difference calculation.

Compared with the BSR criterion, the MSSFE criterion tends to make the value of *β* smaller (average value of 64.26), which is not conducive to improving the robustness of the DSRC-R method. The average RSME under the MSSFE criterion is 1.040, which is significantly greater than 0.759 under the BSR criterion, indicating that the BSR criterion makes the DSRC-R method more stable. This is because the MSSFE method is equivalent to the least square method in linear systems, meaning the value of *β* should be 0 if the DSRC-R method is applied in a linear system. Although the XAJ model is a non-linear system, the value of *β* will still tend to be small, which is not conducive to improving ill-conditioned problems. However, the BSR criterion introduces BDSR and RDSR, and this makes more use of the effective information contained in the simulated errors, which is conducive to avoiding a value of *β* that is too small, making the DSRC-R method more stable.

To sum up, the BSR criterion has greater advantages among the three criteria. It pays more attention to extracting effective information from simulated errors of the outlet flow (in fact, any type of error will eventually be reflected here). It takes more account of the mutual independence and zero mean property of residuals, which is conducive to selecting a more reasonable value of *β* and improving the performance of the DSRC-R method.

#### *3.3. Real Case*

The research basin in this study is Tankeng basin, with a total area of 3330 km2 and 15 precipitation stations, which is located at the tributary of the Ou River in Zhejiang Province. Tankeng basin is in the subtropical monsoon climate zone, and it enjoys a temperate climate with well-marked seasons and plenty of rainfall and sunshine. Runoff from Tankeng basin is mainly supplied by precipitation. The annual average precipitation of Tankeng basin is between 1500 and 2100 mm, and multi-annual average evaporation is 969.9 mm. The precipitation spatial distribution is uneven, as is the annual precipitation temporal distribution; based on this, the entire year can be divided into three parts with the first part "spring rain" ranging from March to April, the second part "plum rain" ranging from May to June and the third part "thunderstorm" ranging from July to September. In Tankeng basin, multi-annual average flow is 120 m3/s, and the maximal runoff happens in June, with the proportion of the entire annual runoff reaching up to 19.7%. More details about Tankeng basin and its rainfall station distribution are illustrated in Figure 6.

**Figure 6.** Map showing the location of Tankeng basin and depicting the stations. The map describes the location, longitude and latitude range and shape of Tankeng basin and shows the location of the rainfall stations and the flow station.

#### 3.3.1. Data

There are two time scales for observed data including the day scale and the hour scale. Evaporation data are day-scale information, while precipitation and flow data are both day-scale information and hour-scale information. This research collected observed data such as evaporation, observed flow and precipitation from 1980 to 2005. The rainfall data were from 15 rainfall stations, evaporation data came from one evaporation station and observed flow data came from the outlet flow of Tankeng basin.

Tankeng basin has been used in previous studies using the XAJ model [20,28], and thus no calibration was required here. The parameters of the XAJ model at Tankeng basin are shown in Table 4.


**Table 4.** Parameters of the XAJ model for Tankeng basin.

#### 3.3.2. Results and Discussion

The statistical indicators of the real case are different from those of the synthetic case. In the real case, the "true" value of rainfall cannot be obtained, meaning RMSE cannot be applied. Here, the indicators include RPF, RRD, NSE and TU.

In this study, 31 historical floods in the basin were selected to compare the results of three criteria including the L-curve criterion, MSSFE criterion and BSR criterion. In order to show the performance of the DSRC-R method under each criterion, the statistical indicators of each flood are listed in Table A1 (see Appendix C). In order to more clearly compare the results of the BSR criterion and the other two criteria, scatter diagrams of each indicator of 31 floods were constructed, as shown in Figure 7. As it is illustrated in Table A1, in terms of the flood peak, the average RPF is 5.42% under the L-curve criterion, 3.95% under the MSSFE criterion and 3.49% under the BSR criterion; therefore, it is apparent that the BSR criterion has a better performance. In terms of the runoff depth, the average RRD is 2.77% under the BSR criterion, which is 3.91% and 1.31% lower, respectively, than that under the L-curve and MSSFE criteria. In terms of NSE, the BSR criterion has the maximal NSE of 0.940, while the NSE under the L-curve criterion is 0.933 and 0.938 under the MSSFE criterion. In terms of operational efficiency, the average time consumed under the BSR

criterion is 12.7 s, which is 25.7 s lower than that under the L-curve criterion, and this illustrates that the operation efficiency is dramatically improved when the BSR criterion is adopted.

**Figure 7.** The contrast of application effects of 31 floods under three criteria. (**a**) The contrast of RPF of each flood between the BSR criterion and L-curve criterion; (**b**) the contrast of RRD of each flood between the BSR criterion and L-curve criterion; (**c**) the contrast of NSE of each flood between the BSR criterion and L-curve criterion; (**d**) the contrast of RPF of each flood between the BSR criterion and MSSFE criterion; (**e**) the contrast of RRD of each flood between the BSR criterion and MSSFE criterion; (**f**) the contrast of NSE of each flood between the BSR criterion and MSSFE criterion.

Figure 7a–c show the comparison of the statistical indicators of 31 floods between the BSR criterion and L-curve criterion. Figure 7a,b are inclined to the lower side of the 1:1 line, and Figure 7c is inclined to the upper side of the 1:1 line, indicating that the BSR criterion is more effective than the L-curve criterion for most floods. In terms of the flood peak and runoff depth, the BSR criterion is more effective in 24 floods and 29 floods, respectively; additionally, the NSE of 23 floods is greater under the BSR criterion. Figure 7d–f show the comparison of flood indicators between the BSR criterion and MSSFE criterion. Similar to Figure 7a–c, under the BSR criterion, the forecasting accuracy of 21 flood peaks and 25 flood runoff depths is higher, and NSE of 18 floods is larger, indicating that the BSR criterion has more advantages over the MSSFE criterion for most floods.

Compared with the L-curve criterion, the BSR criterion improves the performance of the DSRC-R method. This is because the BSR criterion introduces BDSR and RDSR, thus more fully considering the properties of the model-simulated residuals including the zero mean property and the mutual independence. The results in Table A1 show that the average BDSR corresponding to the BSR criterion is 648.1 m3/s, which is far less than 1269.9 m3/s of the L-curve criterion, indicating that the BSR criterion is more conducive to the zero mean property of the simulated residuals. At the same time, the average RDSR corresponding to the BSR criterion is 1.69, which is greater than 1.28 of the L-curve criterion, indicating that the BSR criterion tends to meet the mutual independence of the simulated residuals. Therefore, under the BSR criterion, the DSRC-R method has a better performance. In the real case, the mean value of *β* under the MSSFE criterion is 327.67, which is far less than 2048.97 under the BSR criterion. This is because the MSSFE criterion only considers the sum of squares of errors, which tends to make *β* small, while the BSR criterion focuses

on the nature of the simulated residuals and extracts more effective information from the outlet flow error information, which is conducive to avoiding a small value of *β*. Therefore, the application of the BSR criterion in the real case is better than that of the MSSFE criterion.

#### **4. Conclusions and Prospect**

A previous study [19] proposed the DSRC-R method and verified that it has stronger robustness. However, the selection criterion of the ridge coefficient, usually the L-curve criterion [17–19], has only received limited attention. This essay constructed the BSR criterion based on the properties of model-simulated residuals, utilizing the indicator RDSR to quantitatively describe the independence of the residual series while utilizing the indicator BDSR to quantitatively describe the zero mean property of residuals. Additionally, we then contrasted the performance of the DSRC-R method under three different criteria (BSR criterion, L-curve criterion and MSSFE criterion) through synthetic and real-data studies.

From the results, we found that among the three criteria, the BSR criterion is more suitable for the DSRC-R method. Compared with the L-curve criterion, the BSR criterion improves the performance of the DSRC-R method. This is because the BSR criterion introduces RDSR and BDSR, which quantitatively describe the mutual independence of model-simulated residuals and zero mean property of model-simulated residuals, respectively. Moreover, the BSR criterion saves more time than the L-curve criterion because the BSR criterion does not involve derivative calculation. In addition, compared with the traditional MSSFE criterion, the BSR criterion is more conducive to enhancing the robustness of the DSRC-R method. The MSSFE criterion tends to make the ridge coefficient *β* smaller, which is unfavorable to the performance of the DSRC-R method. Meanwhile, the BSR criterion is conducive to avoiding a small ridge coefficient by extracting more effective information contained in the simulated errors, and this makes the DSRC-R method more robust and improves its performance.

Further research is needed. The BSR criterion proposed in this paper improves the performance of the DSRC-R method, and this seems to benefit from the rational use of outlet flow information. In recent years, data assimilation technologies combined with radar information and remote sensing information have been continuously emerging and have received a significant amount of attention. However, making full use of outlet flow information to update hydrological elements deserves more attention. This is because hydrological elements such as rainfall, evaporation and soil moisture will eventually be reflected in the outlet flow. Therefore, much more effort needs to be made in this regard.

**Author Contributions:** Conceptualization, W.B. and K.L. (Kexin Liu); methodology, W.B. and K.L. (Kexin Liu); software, Y.H. and Y.S.; formal analysis, K.L. (Kexin Liu); data curation, L.L.; writing original draft preparation, K.L. (Kuang Li); writing—review and editing, D.L.; visualization, Y.S.; supervision, Y.H.; project administration, Y.H.; funding acquisition, L.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Fundamental Research Funds of IWHR (grant No. WH0145B032024), the Special Project of National Science and Technology Basic Resources Investigation (grant No. 2021xjkk04005) and the National Key R&D Program of China (grant No. 2018YFC0406400).

**Data Availability Statement:** Data used for this study can be made available upon request.

**Acknowledgments:** We gratefully acknowledge financial support from the Fundamental Research Funds of IWHR and the National Key R&D Program of China.

**Conflicts of Interest:** The authors declare no conflict of interest.

**Appendix A. Proof of Δˆ** *PRE* **<sup>=</sup>** *<sup>n</sup>* ∑ *i***=1** *ki <sup>λ</sup>i***+***<sup>β</sup> vi* **in the DSRC-R Method**

Let *<sup>λ</sup>i*(*<sup>i</sup>* = 1, ··· , *<sup>n</sup>*) be the eigenvalue of matrix *<sup>U</sup>TU*.

Next, *<sup>λ</sup><sup>i</sup>* + *<sup>β</sup>*(*<sup>i</sup>* = 1, ··· , *<sup>n</sup>*) is the eigenvalue of matrix *<sup>U</sup>TU* + *<sup>β</sup>I*. Observe that matrix *STS* and matrix *UTU* + *βI* have the same eigenvectors.

We already have *STS* as a symmetric matrix; therefore, the N mutually orthogonal unit eigenvector *vi*(*i* = 1, ··· , *n*) exists. Hence

$$S^T S v\_i = \lambda\_i v\_i \text{.} \ (S^T S + \beta I) v\_i = (\lambda\_i + \beta) v\_i$$

Let *A*= (*v*1, ··· , *vn*) be the orthogonal matrix.

Next, *A* is the full rank; thus, *Ak* = *ST*Δ*Q* must have a unique solution, ... Then, *<sup>S</sup>T*Δ*<sup>Q</sup>* has a linear expression with *vi*(*<sup>i</sup>* = 1, ··· , *<sup>n</sup>*) and can be solved with only one group *ki*(*i* = 1, ··· , *n*). Hence

$$\begin{array}{rcl} (\prescript{}{}{S}^{T}\mathbf{S} + \beta I) \Delta \hat{P}\_{RE} &=& \prescript{}{\*}{}^{T}\Delta Q \\ &=& \sum\_{i=1}^{n} k\_{i} v\_{i} \\ &=& \sum\_{i=1}^{n} \frac{k\_{i}}{\lambda\_{i} + \beta} [(\lambda\_{i} + \beta) v\_{i}] \\ &=& \sum\_{i=1}^{n} \frac{k\_{i}}{\lambda\_{i} + \beta} [(\prescript{}{}{S}^{T}\mathbf{S} + \beta I) v\_{i}] \\ &=& (\prescript{}{}{S}^{T}\mathbf{S} + \beta I) \sum\_{i=1}^{n} \frac{k\_{i}}{\lambda\_{i} + \beta} v\_{i} \\ &\ddots (\prescript{}{}{S}^{T}\mathbf{S} + \beta I) \hat{\Delta P}\_{RE} = (\prescript{}{}{S}^{T}\mathbf{S} + \beta I) \sum\_{i=1}^{n} \frac{k\_{i}}{\lambda\_{i} + \beta} v\_{i} \end{array}$$

Notice that *STS* + *βI* is reversible. Hence, Δˆ *PRE* <sup>=</sup> *<sup>n</sup>* ∑ *i*=1 *ki <sup>λ</sup>i*+*<sup>β</sup> vi*.

**Appendix B. Proof of** *f* **(***β***) = log[** *n* ∑ *i***=1 (** *ki <sup>λ</sup>i***+***<sup>β</sup>* **) 2 ]**

Notice that if we want to prove the equation *f*(*β*) = log[ *n* ∑ *i*=1 ( *ki <sup>λ</sup>i*+*<sup>β</sup>* ) 2 ], firstly, we need to prove <sup>Δ</sup><sup>ˆ</sup> *PRE* 2 <sup>2</sup> <sup>=</sup> *<sup>n</sup>* ∑ *i*=1 ( *ki <sup>λ</sup>i*+*<sup>β</sup>* ) 2 .

Furthermore,

$$\left\| \| \hat{\Delta} \mathcal{P}\_{RE} \| \right\|\_{2}^{2} = \left( \hat{\Delta} \mathcal{P}\_{RE}, \hat{\Delta} \mathcal{P}\_{RE} \right) = \left( \sum\_{i=1}^{n} \frac{k\_{i}}{\lambda\_{i} + \beta} \upsilon\_{i\prime}, \sum\_{i=1}^{n} \frac{k\_{i}}{\lambda\_{i} + \beta} \upsilon\_{i} \right) = \sum\_{i=1}^{n} \sum\_{j=1}^{n} \frac{k\_{i}}{\lambda\_{i} + \beta} \frac{k\_{j}}{\lambda\_{j} + \beta} (\upsilon\_{i\prime}, \upsilon\_{j}).$$

where (*vi*, *vj*) is the inner product of *vi*, *vj*.

Notice that *vi* (*i* = 1, ..., n) is the mutually orthogonal unit vector.

$$\begin{array}{c} \text{\(\,.\,(v\_i, v\_i) = 0\text{ when } i \neq j \text{ and } \left(v\_i, v\_i\right) = 1 \text{ (\$i = 1, 2 \cdots \$n\$)\)}\\ \text{\(\,.\,\|\,\|\hat{\Delta} P\_{RE}\|\|\_{2}^{2} = \sum\_{i=1}^{n} \left(\frac{k\_i}{\lambda\_i + \beta}\right)^{2}\\ \text{\(\,.\,}\, f(\beta) = \log[\sum\_{i=1}^{n} \left(\frac{k\_i}{\lambda\_i + \beta}\right)^{2}]. \end{array}$$




1 The average values of RPF and RRD actually represent the average of the absolute values of RPF and RRD, respectively.

#### **References**

