Estimation of the Interaction Between Groundwater and Surface Water Based on Flow Routing Using an Improved Nonlinear Muskingum-Cunge Method

The interaction between groundwater (GW) and surface water (SW) not only sustains runoff in dry seasons but also plays an important role in river floods. Lateral inflow is the recharge of groundwater to surface water during a river flood; this recharge is part of the GW-SW exchange. Hydrological engineers proposed the idea of modelling flood routing using the Muskingum-Cunge method, in which the GW-SW exchange is not fully considered. This study proposes an improved nonlinear Muskingum-Cunge flood routing model that considers lateral inflow; the new method is denoted as NMCL1 and NMCL2 and can simulate flood routing and calculate the GW-SW exchange. In addition, both the linear and nonlinear lateral inflows (with the channel inflows) are discussed, and the stable lateral inflows that occur due to the GW-SW exchange are considered for the first time. A sensitivity analysis shows that different parameters have different effects on the simulation results. Three different flood cases documented in the literature are selected to compare the four classical and two updated Muskingum-Cunge methods. Two different floods of the River Wye are selected to verify the accuracy of the calibrated model. The simulation results of the improved Muskingum-Cunge method are compared with the temperature inversion results measured from the Zhongtian River, China, to indicate the feasibility and reliability of the improved method. A comparison shows that, for several cases, the proposed method is capable of obtaining optimal simulation results. The proposed method inherits the ability of the Maskingum-Cunge method to simulate flood routing. Moreover, it can quantify the GW-SW exchange, and the reliability of the estimations is owed to the nonlinearity and sign flexibility of the calculated exchange process.

analysis and selects different groups of literature data to estimate the model parameters, the 115 channel routing, and the exchange volume of GW and SW. Section 4 applies the new model for a 116 field study conducted by the authors, and evaluates the lateral inflow and exchange volume 117 estimated by the traditional MC method and the proposed one, as well as the thermal method. 118 Section 5 discusses the application and deficiency of the proposed methods. Lastly, Section 6 119 reports the main conclusions. . This assumption is not valid for some reaches of the river, motivating the nonlinear 151 Muskingum model (Mohan, 1997), such as the nonlinear Maskingum method (named NMC1) 152 proposed by Gill (1978): where the new parameter m is the exponent of the power-law relationship between the 156 accumulated storage and the weighted flow. When m=1, equation (5) reduces to the specific 157 linear relation. K in the nonlinear models is the storage parameter for the river reach assumed to 158 be equal to the travel time of flood peaks that move through the reach. Also, x is the 159 dimensionless weight which indicates a weight between inflow and outflow on storage depends 160 on the channel type. The graphical method used for historical inflow and outflow hydrographs is 161 inappropriate for the nonlinear Muskingum models. Hence, the determination of K and x is 162 difficult, requiring advanced methods (Moghaddam et al., 2016). 163 (2) MC models with lateral flow (MCL1) 164 All the models mentioned above ignore the lateral flow existing in the river reach in the 165 actual flood events. O'Donnell (1985) assumed that the lateral flow entering the river reach was 166 directly proportional to the inflow with a proportionality factor  , and then proposed the first 167 MC model with the lateral flow: Integrating the nonlinear continuity and storage equations, which takes lateral flow into 172 consideration, one obtains: Tung (1985) and Geem (2006) 2.2 Improve MC method to calculate GW-SW interaction (NMCL1 and NMCL2) 181 The MC methods reviewed above assume that the lateral inflow is linearly proportional to 182 the inflow. The zero-value inflow leads to zero lateral inflow, which cannot explain the potential 183 lateral inflow due to GW-SW interaction. To solve this issue, it is hereby first assumed that a 184 stable GW-SW interaction process exists before the flood event, and add a constant representing 185 the stable exchange volume in the storage (named NMCL1): where e represents the stable lateral inflow due to GW-SW exchange. We can then estimate the 188 sum of GW-SW interaction and lateral inflow using: Another option is to consider the nonlinear relationship between the lateral inflow and the 191 channel inflow, leading to the following storage (named NMCL2): Both the linear and nonlinear relations between the inflow and the lateral inflow are 198 considered in this study and will be checked against real-world data in the next section. 199 When using the MC methods mentioned above to simulate flood routing, the primary objective is 200 to minimize the discrepancy between the simulated and measured outflow. For this purpose, the 201 sum of the squared errors (SSQ), the Nash efficiency coefficient (NSE), and the root mean 202 square error (RMSE) are used as the objective functions.  tests are needed before starting the actual search process to determine 1  and 2  .

212
The particle swarm optimization (PSO) approach, which is an intelligent cluster 213 optimization algorithm (Eberhart and Kennedy, 1995), is applied for parameter optimization.

Model Validation
We choose the improved NMCL1 and NMCL2 methods, which consider the nonlinear lateral 218 inflow (i.e., the portion changing nonlinearly with the river inflow) and the stable lateral inflow 219 due to GW-SW exchange. 220 The observed flow at the River Wyre in UK (O'Donnell, 1985) is selected here for 221 sensitivity analysis. This example involves inflow with multi peaks (Figure 1) Table 1, it can also be 253 seen that the sensitivity of parameter x is only next to that of parameter m , which is ten times 254 larger than that of parameter K . Here x is the weighting factor allocating the contribution of the 255 wedge to the total storage. A larger x (>0) means more storage for the wedge, so that more river 256 inflow (from the channel) can directly transfer to outflow, resulting in a stronger fluctuation of 257 the outflow, especially during the advance of a flood event. 258 The change of K also has big influence on the model results (Figure 2( shows that the smaller the parameter K is, the greater the sensitivity. K is the storage constant 265 expressing the ratio between storage and discharge, which may also be viewed as the lag or 266 travel time through the reach (Karahan, 2012). When K is smaller, the flood rise and fall 267 difference of the reach is larger, and the water storage capacity of the channel is also reduced. 268 When K is equal to 1, the channel water storage is very small and the flood fluctuation increases. 269 At this time, the inflow of the river is constant, and the outflow needs to be increased to balance The parameter  shown in Figure 2(e) and 2(g) and Table 1 is more sensitive to the outflow 291 than e . This sensitivity is stronger for the NMCL2 method than the NMCL1 method. This is expanded, the results obtained above need to be re-evaluated. Using the optimal parameters, the results simulated by different MC methods are compared.
312 Table 2 lists the error analysis of flood routing simulation results using different methods for the 313 three cases.  315 The first example (Wilson, 1974) was a smooth single-peak hydrograph with low lateral   phase of rain event, the GW-SW interaction gradually returns to the initial state before rainfall. 353 The flood data observed at the River Wyre, UK (O'Donnell, 1985) has considerable lateral 354 inflow. We use the six MC methods mentioned above to simulate the flood routing. 355    . We tried the 477 three popular methods and found that these three methods generated deviating results, which we 478 attribute to the sensitivity of the Ar method to non-ideal input data (see Hatch et al. (2006)). 479 Therefore, we chose to use the Ar method because of its robustness and consistent results.

497
The measured water level data are used to calibrate model parameters in the PSO optimal 498 algorithm, and the flood routing and exchange volume are simulated by different MC methods. According to Figures 8 and 9, only the GW-SW exchange estimated by the NMCL2 539 methods is slightly greater than 0 at the initial stage, indicating that the quantitative difference of 540 these methods in the interaction of surface water and groundwater is small, and the amount of 541 surface water gained from groundwater and the lateral inflow is slightly greater than that due to This feature reflects the observed fact that the river water level rises due to rainfall, which 553 simultaneously enhances the loss of surface water and strengthens the exchange between GW 554 and SW. The trend of the simulation results of the two methods is consistent. Compared with the 555 thermal based method, the flood routing method proposed by this study is more accurate and 556 easier to operate, which can simulate the overall change of GW-SW exchange volumes. When 557 the field investigation is limited, more reliable dynamic data of flood routing can be obtained by 558 numerical simulation. In addition, before the field survey, the numerical simulation method can 559 also be used to predict the flow routing, so as to make the field work more feasible.       Qout are the measured inflow and outflow, respectively. QMC1 is the estimated outflow using MCL1, QMC2 is the estimated outflow using MC2, QNMC1 is the estimated outflow using NMC1, QMCL1 is the estimated outflow using MCL1, QNMCL1 is the estimated outflow using NMCL1, QNMCL2 is the estimated outflow using NMCL2. QEX represents the sum of the transient/conventional lateral inflow and the vertical inflow due to the interaction between GW and SW. and Qout are the measured inflow and outflow, respectively. QMC1 is the estimated outflow using MCL1, QMC2 is the estimated outflow using MC2, QNMC1 is the estimated outflow using NMC1, QMCL1 is the estimated outflow using MCL1, QNMCL1 is the estimated outflow using NMCL1, QNMCL2 is the estimated outflow using NMCL2. QEX represents the sum of the transient/conventional lateral inflow and the vertical inflow due to the interaction between GW and SW.