On the Modeling of Bank Storage in a Groundwater Model : The April , 1983 , Flood Event in the Neuwieder Becken ( Middle Rhine )

For predictive numerical simulations of subsurface floods (groundwater head rise due to high water in a contiguous river), it is important to know how to represent the bank storage process in a numerical groundwater model. Whilst leakage approaches are frequently used for modeling bank storage, another option is the application of a head boundary condition. In order to get a better understanding of the bank storage process, we analyze the bank storage event in the Neuwieder Becken (Middle Rhine) in April 1983, which has been reported by Ubell (1987). We found the leakage function to be nonlinear and hysteretic. The evaluation of different model variants for Ubell’s bank storage event shows that both a head boundary condition and a leakage boundary condition are appropriate modeling approaches. For practical reasons, the leakage boundary condition is preferred. A linear leakage function represents the bank storage process for the analyzed event sufficiently. A hysteretic course of the leakage function can be achieved in a three-dimensional groundwater model by layering the hydraulic soil properties in the vicinity of the bank.


Introduction
The August 2002 flood was one of the most disastrous floods in Germany in recent years.Many cities and villages in the catchment of the River Elbe were affected.A special characteristic of this flood event was the severe groundwater head rise: at some locations, the groundwater rose up to six meters over a short period of time.This level had never been measured before [1].In particular, the inhabitants of the City of Dresden were impacted by this type of subsurface flood.We define a "subsurface flood" (unterirdisches Hochwasser or Grund-Hochwasser; see [2,3]) as a groundwater head rise that is caused by high water in a river and distinguish it from groundwater flooding, which is related to heavy rainfall events (see, e.g., [4,5]).For the Dresden flood event, Huber et al. [6] and Kreibich and Thieken [7] quote that a significant percentage of the total economical losses due to the flood was related to groundwater.
Subsurface floods are mainly a problem of urban areas that are located in the vicinity of large rivers.Damage arises from groundwater that enters cellars and basements and destroys walls and interiors by wetting.When critical infrastructure, like power grid facilities, is affected [8], a power breakdown can be a consequence of a subsurface flood.Groundwater that enters sewage pipes can cause an overload of the sewage system [1,2].Finally, rising groundwater levels can cause the failure of foundations due to uplift forces.This can affect apartment buildings and office buildings, but also road and railway tunnels.
A special problem of subsurface floods is that the groundwater rise is often unexpected.This is for two reasons: Firstly, the subsurface flood can propagate far away from the river.Therefore, subsurface floods can affect locations that are usually not affected by the visible surface inundations.Secondly, the subsurface flood follows the river flood with delay.Groundwater levels keep rising when the flood wave in the river is already decreasing.Beyer [8] reports the case of a low-lying sports hall in Dresden during the August flood 2002.The basement plate of this sports hall is located 3.6 m below the ground surface, and the building is located in an area that was not inundated.Not until the parquet floor bulged several decimeters due to groundwater entering the hall had the threat even been noticed.As an urgent measure, water tanks and sand bags were brought into the hall to prevent the ground plate from cracking.
Subsurface flood hazard maps help urban planners to prepare cities located in the vicinity of large rivers for protection against subsurface floods.Decision makers benefit from such maps during flood protection operations, and inhabitants can inform themselves if they live in areas that are at risk of subsurface floods.The essential information of such a subsurface flood hazard map is the minimal depth to the groundwater table that probably is reached during a flood.Groundwater models are a helpful tool to prepare such subsurface flood hazard maps [9].
Bank storage (see [10]) is one important mechanism of water exchange between a river and the neighboring aquifer and, consequently, an important mechanism to account for in studies related to subsurface floods.Bank storage events are reported in several studies [11][12][13][14].Within the frame of predictive subsurface flood modeling, it is important to understand the corresponding interchange processes and to know how they can be represented in a numerical groundwater model.In groundwater models, the river-aquifer interaction is usually represented with a head boundary condition or a leakage boundary condition (see, e.g., [15][16][17]).The question we address in this article is which one of these two approaches is more suitable for large rivers that strongly interact with the adjacent aquifer.In Section 2, we summarize the mathematical background of these two options and discuss selected leakage functions from the literature.
Ubell [12] reports a bank storage event that occurred during the Rhine flood in April 1983, in the "Neuwieder Becken" (Middle Rhine).We review this bank storage event in Section 3. In Section 4, we describe a groundwater model that we set up for the Ubell bank storage event with the flow simulation program Feflow (see [18,19]).The water exchange between river and aquifer is modeled with a head boundary condition and a leakage boundary condition, which leads to two model variants.The groundwater model should reproduce the groundwater levels observed by Ubell [12].Based on the simulation results of these two model variants, we discuss the two model approaches for modeling bank storage in a groundwater model.The article closes with a recommendation for modeling bank storage and an outlook on the further needs for research.
A second order boundary condition (q-boundary condition, Neumann-type) is usually not applicable, because the amount of water exchange per time step is unknown in most cases.

Head Boundary Condition
Modeling the river-aquifer interaction with a head boundary condition means setting h R for the unknown groundwater level h at the boundary condition node: The boundary condition causes a withdrawal or an allocation of water at the boundary node according to the applied boundary condition value.The exchange flux also depends on the groundwater flow situation in the vicinity of the boundary condition node.Freeze [20] limits the applicability of head-boundary conditions to shallow rivers with non-varying depths and where steep groundwater gradients and large leakage areas are existent.For subsurface floods, these conditions will not apply in general.Cunningham and Sinclair [21] recommend head-boundary conditions for rivers that are fully penetrating the aquifer.This condition will apply to many of the large rivers where subsurface floods occur.Ubell [22] uses a head boundary condition in his modeling study.Modeling bank storage with a head boundary condition does not necessarily rely on the assumption of an interface layer between river and aquifer.If there is no interface, h R is the river stage.If a sediment layer separates the river and aquifer, the sediment layer must be modeled directly, i.e., with elements that comprise the hydraulic conductivity properties of the interface layer, or h R represents the groundwater level on the aquifer side of the interface layer.Boundary condition data for the latter case can be obtained by pre-model calculation.

Introduction
A leakage boundary condition uses a leakage function with the river stage h R as the boundary condition value.Figure 1 shows a schematic view of an idealized stream and the neighboring aquifer.The river stage h R,1 is higher than the groundwater level h G .In this case, the water infiltrates towards the aquifer.If the river stage is lower than the groundwater level, the hydraulic gradient points towards the stream, and water exfiltrates from the aquifer.Mathematical leakage models or leakage approaches relate the head difference between the river stage h R (m) and the groundwater level h G (m) to a leakage flux q ( m /s) per unit of the interface area between the river and the aquifer (the leakage rate) with a so-called leakage function: The groundwater level h G and the leakage rate q are initially unknown and obtained as model results.A time-variant groundwater level in a transient groundwater model accounts for the impact of storage on the interchange process.In this article, positive values of q and ∆h are defined as infiltration (inflow of river water into the aquifer), and negative values of q describe exfiltration.
. Schematic view of two different flow conditions in a cross-section of an idealized stream and the corresponding leakage function (head difference ∆h over the exchange rate q, c l,in and c l,ex denote the leakage coefficient for infiltration and exfiltration, respectively).River stage h R,1 gives a positive head difference and leads to infiltration, and the river stage below groundwater level h R,2 leads to exfiltration (modified and extended after [23]).

Leakage Functions
Certainly, the leakage function that is most frequently used is given by the following equation: Within this article, this approach is called the Darcy model, because the conceptional model assumes saturated Darcian flow through the interface layer.In Equation ( 4), c l ( 1 /s) is a leakage coefficient.The leakage coefficient takes positive values and is defined as a ratio of hydraulic conductivity k S ( m /s) and the thickness d (m) of an interface layer that is assumed to separate the river and the aquifer: (see, e.g., [16,17,[24][25][26]). Often, data for d and k S are not available and highly spatially variable.Consequently, in many groundwater modeling studies, the leakage coefficient c l is treated as a model parameter and determined by calibration.
For many rivers in central Europe, water exfiltrates during normal or low flow conditions, while during high water conditions, water infiltrates into the aquifer.If Equation ( 4) is used both for infiltration and exfiltration, the flow resistance must be the same for both flow directions.Brunke and Gonser [27] point out that the stream bed can act like a mechanical filter.Depending on flow conditions and suspended load in the stream, colmation of the bed can reduce the infiltration rate [27][28][29].In order to account for different resistance properties, different leakage coefficients are applied for infiltration (c l,in ) and exfiltration (c l,ex ) (see, e.g., [18,30]): The diagram in Figure 1 shows the relationship between the leakage flux and the head difference according to the Darcy leakage approach.Different leakage coefficients are used for infiltration and exfiltration, respectively.
The total leakage flux Q ( m 3 /s) is the product of the exchange rate q and the interface area A (m 2 ): With Equation (7), the Darcy approach describes a linear relationship between Q and ∆h if the interface area is assumed to be constant.The leakage function becomes nonlinear if the interface area A is calculated in dependence of the wetted perimeter of the river bed p (m) (see Figure 1), as proposed by Pinder and Sauer [10], Cunningham and Sinclair [21], Panday and Huyakorn [31], Krom and Graham [32] and Rushton [26]: In this equation, l (m) stands for the length of the river section, which usually remains constant during the flow simulation.The wetted perimeter p depends on the river stage and the river morphology.
An empirical nonlinear relationship between exchange rate and head difference is presented by Rushton and Tomlinson [33] (see also [34]): for ∆h < 0 C 1 ( m 2 /s), C 2 ( 1 /m) and C 3 ( m 2 /s) are positive leakage parameters.For small changes in the head difference, a large increase of the exchange rate is achieved if the parameters C 1 and C 3 are adjusted appropriately.C 2 limits the exchange rate for large values of ∆h.
The leakage rate q ( m 2 /s) is related to the length of the river section.Thus, the total leakage flux is obtained as: for this case.The diagram in Figure 2 shows example infiltration curves for the Darcy approach, the wetted perimeter approach and Rushton and Tomlinson's empirical approach.For all curves, the groundwater level h G has been set constant to 4.5 m.The river stage h R varies between 2.5 m and 8.5 m.In order to show the different characters of the three approaches, the leakage parameters and the channel geometry (Table 1) have been chosen in such a way that all three leakage functions give flux values in the same order of magnitude for a given head difference.The approach presented by Rushton & Tomlinson [33] leads to disproportional small flux values in case of large head differences.The linear Darcy approach does not differentiate between large and small head differences, and the wetted perimeter approach gives disproportional high flux values for large head differences.For the wetted perimeter approach curve, a trapezoidal cross-section with a bed elevation of 0 m above the datum has been applied.4) and ( 6)), the approach by Rushton and Tomlinson [33] (Equation ( 9)) and the leakage approach that uses the wetted perimeter (Equation ( 8)), evaluated here with a trapezoidal cross-section.Subsurface flood events are characterized by high infiltration into the aquifer during a comparatively short period of time.High water in the river initially causes large head differences.Figure 2 shows that leakage flux values that correspond with high head differences can be very different, depending on which leakage approach is applied.Which leakage approach is suitable depends on the study area, the intensity of the surface water-groundwater interaction and the objectives of the study.For hydrological modeling studies that aim to estimate the runoff induced by a specific rainfall event, a simple approach can be sufficient, while for groundwater models that are bordered by large rivers and that are built in order to assess the impact of groundwater extractions on the groundwater balance, a more advanced approach can be necessary.The third order boundary condition feature of several groundwater flow simulation programs, like the MODFLOW river package [35], FESSim [36] or Feflow [23], is based on the Darcy approach, and several textbooks (e.g., [15][16][17]) recommend this approach for the stream-aquifer interaction; therefore, it is assumed that in many practical model applications, the Darcy approach is used both for small and large rivers and streams.The approach proposed by Rushton and Tomlinson [33] limits the infiltration for high water level differences, so this approach will apply for small rivers and streams with small variations in water levels.Krom and Graham [32] report an application case where they use the wetted-perimeter approach for a braided river with wide cross-sections, where a small increase of discharge is attended by a large increase of the wetted perimeter.In theory, the wetted perimeter approach is physically more accurate than the Darcy approach, because it takes the change of the interface area into account.

The Leakage Boundary Condition of the Groundwater Flow Simulation Program Feflow
The groundwater flow simulation program Feflow provides a leakage boundary condition based on Darcian flow; Equation (4).The leakage rate q = c l ∆h depending on the river stage h R and the unknown groundwater level h G is set as a source or sink on the boundary condition node.Different parameters for infiltration and exfiltration can be set (Equation ( 6)).The total interface area A is calculated as: with: for each boundary condition node i individually within the flow simulation procedure.Therefore, the Feflow leakage approach is a wetted perimeter approach like the one given with Equation ( 8).However, the wetted perimeter p * is not related to the river stage, but to the groundwater level h G , or, in other words, the wetted perimeter is related to the aquifer's side of the interface layer.In a Feflow model, the leakage parameters are a property of the finite elements.The boundary condition is set at a node [30] and uses the average value of the leakage coefficients from the adjacent elements.Boundary condition nodes are activated only if the water level is higher or equal than the elevation of the node.In the situation sketched in Figure 3a, only Nodes 1 and 2 contribute to the leakage flux.Figure 3b shows a different situation, where boundary condition Node 3 also contributes.For Node 3, the average value of leakage coefficients c 1 and c 2 applies.Because the leakage coefficient c 1 is greater than c 2 , the leakage rate q is disproportionally higher than in the case shown in Figure 3a.The larger value of the wetted perimeter p * furthermore increases the total leakage flux Q (not shown in the figure) for the second case.The examples in Figure 3 show how different values of the wetted perimeter p * and the total leakage coefficient c l can produce different leakage fluxes Q for the same head difference ∆h.

Study Area
With a length of 1232 km, the Rhine is one of the largest European rivers.At the gauge in Emmerich, the average discharge is 2330 m 3 /s, and the Rhine is 730 m wide here [37].The Neuwieder Becken is a broad valley at the Rhine located between the cities of Koblenz and Andernach (see Figure 4).Here, the Rhine is about 400 m wide.The Neuwieder Becken covers an area of 7 km in width and 15 km in length and is a result of tectonic settlement within the "Rheinisches Schiefergebirge".Sedimentary and erosion processes generated several terraces (the two large ones are the Ältere Niederterrasse and the Jüngere Niederterrasse).This results in a localized groundwater system with minimal interaction with other groundwater systems [38].Above the tertiary aquifer base of hard rock [12], there is a thick gravel aquifer with a high permeable material that is covered with layers of sand and sandy loam.Loose pumice stones originating from the Laacher See volcano [39] can also be found.In Figure 5, a cross-sectional hydrogeological profile is given.[12]).The observation wells U01, U03 and U05 are referred to in Figure 6.
Ubell [12] gives hydraulic conductivity values k f of between 2 • 10 −2 m /s and 4 • 10 −3 m /s for the Niederterrasse aquifer, the average saturated thickness being 10-15 m and the effective porosity n eff 0.20.The groundwater in the Niederterrasse aquifer is primarily influenced by the Rhine water [12].Groundwater recharge originates mainly from rainfall and base flow from the hillside.A water supply company extracts groundwater for water supply.However, due to the high permeability of the aquifer material, the influence of groundwater extractions on the groundwater table is limited locally.
Near the village of Urmitz at Rhine kilometer 602.4 (official chainage), a chain of seven observation points has been installed [12,40].The location of the observation points is shown in Figures 4 and 5.These observation points are not located in the area of influence of the groundwater extractions [40].4 (modified after Ubell [12]; elevation given in meters above official sea level reference m ü.NN).

Hydrographs during the Flood Event in April, 1983
Figure 6 shows hydrographs of the Rhine river stage and groundwater levels at three groundwater observation points located in the profile "Urmitz" with different distances to the river bank (see Table 2).The hydrographs are related to the flood event in April, 1983.The diagrams give evidence for a strong interaction between the river and neighboring aquifer: according to the distance between the observation point and the river bank, the groundwater level hydrographs follow the river stage, damped and with delay.River and aquifer are always connected, because the groundwater level never falls below the river bed.The fact that the impact of rising river stage reaches that far into the hinterland is facilitated by the high conductivity property of the aquifer and the long duration of the flood event.Rushton [26] points out that the river coefficient or river conductance, which is defined as: "depends primarily on the horizontal hydraulic conductivity of the aquifer system rather than the vertical hydraulic conductivity of the riverbed deposits".A rapid increase of the groundwater level is attended by an increase of saturated thickness that, again, facilitates the entrance and spreading of river water towards the hinterland. .Hydrograph of the Rhine gauge (in meters above official sea level reference m NN) and three groundwater observation points with different distances to the river bank, respectively (values taken from [12]); leakage flux for the 1 km broad profile "Urmitz" as day-averaged values (values taken from [12]) and the head difference for the observation point located at the bank.
Table 2. Distance of groundwater observation points to the river bank for the observation well chain near Urmitz (derived from [12]).The leakage flux has been determined by Ubell [12] with his volume balance method (Volumenbilanzmethode, see [11,12,41,42] and Appendix).The course of leakage flux over time does neither coincide exactly with the Rhine water stage nor with any groundwater observation hydrograph: the maximum exchange flux occurs on Day 12, while the maximum water level in the river Rhine occurs later, on Day 13.The maximum groundwater level for groundwater observation point U01 reaches its maximum value one day later again.

A Leakage Function for Ubell's Rhine Flood Event
In order to derive a leakage function for the flood event in April, 1983, head differences between the river stage and the groundwater level observed at observation point U01 have been calculated (Figure 6, bottom).Observation point U01 has been chosen, because this observation well is located directly at the bank and, thus, is the most feasible of all observation wells according to the leakage approach conceptual model, where the head difference ∆h is defined between both sides of the assumed interface layer.The head difference values mainly comprise two sources of inaccuracy: firstly, the river stage used for the head difference has been measured at a different location and has been transferred to the chain of observation points [12].Secondly, the groundwater level is related to an observation well that does not necessarily represent the edge of the interface layer.If there is a distance between the aquifer side of the interface layer and the location the head difference ∆h is computed, the groundwater level h G reacts to the leakage flux q with delay (see [43]) and not directly.However, as observation point U01 is located at the river bank, we assume this effect to be small.Leakage flux values were taken from Ubell [12].It is generally difficult to obtain the amount of water exchanged between river and aquifer, because it is nearly impossible to measure the leakage flux directly (see Appendix).Ubell's method requires that the flow direction must always be parallel to the chain of observation points.For the current example, this condition is violated for normal flow conditions when the groundwater flows parallel to the river [38].However, the leakage function covers a time period of high water conditions, where the flow direction is oriented along the line of observation points [38].For the calculation of the volume stored in the aquifer per period of time, the storage coefficient has been assumed to be spatially and temporarily constant.Finally, the method neither allows for the separation of groundwater recharge processes and backwater effects from base flow from the hillside, nor does it account for the influences of groundwater withdrawal and the effect of the gravel pits on the groundwater storage.However, the impact of those processes is assumed to be small in comparison to the impact of the infiltration of river water [40].
Figure 7 shows the pairs of values of Q and ∆h for infiltration and exfiltration and smoothing functions.For positive leakage flux, a linear leakage function seems to represent the infiltration process quite well.However, the curve does not hit the point of origin.This is not in accordance with the leakage model concept.Inaccuracies of the volume balance method that have been mentioned above can be one reason for this behavior.Another explanation is that the water exchange is governed by processes not covered by the leakage conceptual model.These processes are discussed at the end of this section.
For the exfiltration process, the slope of the regression line is slightly lower than the slope of the regression line for the infiltration process.The axis intercept is larger than in the case of infiltration, while the coefficient of determination R 2 calculated for exfiltration is smaller than it is in the case of infiltration.An exponential function appears to match better with the observed data.However, the interpretation that the leakage flux is independent of the head difference seems also to be suitable here.In this case, the leakage model would not be appropriate to describe the interchange process.As mentioned in Section 2.3.2, it is often assumed that the leakage flux for a certain head difference is smaller for infiltration than for exfiltration, due to a river bed colmation (see, e.g., [27]).The leakage function presented here does not support this assumption for the current case.
In the case of infiltration, the absolute head differences range between zero and 1.28 m, while for exfiltration, the maximum absolute head difference is only 0.45 m.Large absolute head differences coincide with large leakage fluxes.This matches well with the view in which a groundwater level rise occurs within a short time period due to infiltration, but then decays very slowly (see, e.g., [1,2]).An evaluation of the Q-∆h relation with respect to the chronological sequence of the flood event is given in Figure 8.The curve has a hysteretic course: up to Day 10, the curve increases with an approximately linear course.From Day 10 to Day 11, the head difference rises only a little, while the leakage flux increases disproportionally high.On Day 12, the head difference already begins to decrease, but the leakage flux continues to rise.Therefore, a limitation of the leakage flux for large head differences ∆h, as given with Equation (9), seems less appropriate than a leakage function that returns disproportionally high values for large head differences, like the wetted perimeter approach.
There are several possible physical explanations for a hysteretic leakage function: • Both the groundwater level and the river stage change over time.Consequently, the wetted perimeter of the river cross-sectional profile changes as well.Therefore, for the same head difference, two different leakage fluxes are possible.• The properties of the interface layer between river and aquifer change during the flood event (see, e.g., [44]).Sedimentation of a suspended load or biogenous stabilization processes create a colmation layer on the river bed during normal or low flow conditions.A flood wave can destroy this colmation layer because of increased shear stress, and reduce the sediment layer's thickness or even totally erode the sediment layer.In addition, the loosening of the river bed material may increase the hydraulic conductivity of the river bed, which increases the exchange of water between the river and the aquifer.With the decay of the flood event, the river bed is clogged again by sedimentation and stabilization processes and allows the leakage flux to decrease slowly [28,[45][46][47].Rushton [48] (Section 12.2) reports a hysteretic course of the river coefficient over river stage for the semi-arid River Yobe in northeastern Nigeria within a cycle of the hydrological year.• The temperature of infiltrating river water is higher than the temperature of the pre-event groundwater in the aquifer (see [25]).With an increasing temperature, the hydraulic conductivity increases, which again allows more river water enter the aquifer.• Preferential flow paths created by animals or ice crushes can be activated during the flood (see [49]).Within this section, we present a groundwater model for Ubell's bank storage case.The objective is to reproduce the observed groundwater levels and the derived leakage flux from Ubell [12] with the groundwater model.The river-aquifer interaction is modeled in two different ways: as a leakage boundary condition and as a head boundary condition.By comparing the simulated leakage function of these two model variants with the one derived from observed values, we assess these two ways of modeling the river-aquifer interaction in a groundwater model.

The Groundwater Model
Ubell [22] already presents a groundwater model of the vertical cross-section through the aquifer along the chain of observation points near Urmitz (see Figures 4 and 5).He uses this model to quantify the impact of river-aquifer exchange processes on the water balance of the River Rhine.
Our groundwater model represents the same model area, and the conceptual model of the hydrogeological system follows Ubell [22] using the equations for saturated flow in porous media: (see, e.g., [17,50]).In this equation, h denotes the groundwater head at time t; k i the hydraulic conductivity; and x 1 , x 2 and x 3 are Cartesian space coordinates.The specific storage S sp incorporates the specific yield or drainable porosity S y and the specific storage (compressibility).For the current case, the specific storage is neglectable against the drainable porosity, because of unconfined conditions.Equation ( 15) is solved with the finite element method using the groundwater flow simulation program Feflow [18,19].
For the current model, the flow equations are solved in three dimensions, but boundary conditions are set in such a way that there is no velocity component perpendicular to the cross-section, so the model is in principle a two-dimensional vertical model.The computational grid and the location of boundary condition nodes are shown in Figures 9 and 10 for the two model variants.The computational grid has been generated automatically, and it consists of one times 247 brick elements in the horizontal plain.Each element has a base area of 13.415 m by 13.415 m, which adds to a length of 3300 m, plus one element length for the river.The vertical dimension is discretized by four element layers.On the left model boundary, two element columns represent the river (the water), the river bed and aquifer material below the river bed and the bank.The element height of the topmost element layer is variable and changes with the free groundwater surface.The movement of the groundwater table is modeled with the BASD method (best adaption to stratigraphic data; see [51]).This method incorporates the free-surface location, which is a priori unknown, as two additional boundary conditions of the three-dimensional formulation of the flow equations.In the head boundary condition model variant, we set the river stage at the nodes that represent the river itself, while for the leakage boundary condition variant, we applied the boundary condition on those nodes that represent the surface of the interface area.Following Ubell [22], on the right boundary (the Niederterrassenrand), a no flow boundary condition is set, and neither groundwater recharge nor leakage from the aquifer base are considered in the model.

Calibration Data and Calibration Parameters
A model calibration aims for a good match between simulated values and observed values (the calibration data) by adjusting specific model parameters (the calibration parameters).For the current case, the hydrographs from the seven observation points [12] are used as calibration data.The locations of the observation points are given in Figures 4 and 5, and three hydrographs from the observation points are given in Figure 6.In order to assess the calibration quality, the Nash-Sutcliffe efficiency [52] is calculated for each observation point hydrograph.The hydrographs represent the groundwater level at a point in space over time.Ubell [12] also constructed groundwater tables from the observed data for the whole modeling area, i.e., the cross-sectional profile, for three points in time.These groundwater tables are used as calibration data for visual assessment.
The hydraulic conductivity k and the leakage parameters c l are used as calibration parameters and leakage parameters only in the case of the leakage boundary condition model variant.The finite elements have been divided into eight element zones, as shown in Figure 11.Zones 1, 2 and 3 represent the river bank, and Zones 4 to 8 have been defined based on a preliminary study [53].Zone 5 and 7 account for the change of the groundwater gradient in Ubell's groundwater tables (see Figure 12) that indicates discontinuities in the aquifer.Table 3 indicates the usage of calibration parameter zones according to the calibration method and the model variant.Figure 12.Groundwater tables simulated with the manually calibrated model with the head boundary condition (a) and leakage boundary condition (b); groundwater tables from observed measurements from Ubell [12], for three points in time, respectively (after [53]).
A drainable porosity value of S y = 17.5 % for elements representing the aquifer has been determined in the preliminary analysis (see [53]).In the vicinity of observation point U05, the storativity value has been set to S y = 80 % in order to account for a gravel pit.River elements are assigned a value of S y = 0.

Calibration Procedure
The model variants have been calibrated manually by carrying out the following steps repeatedly: (1) comparison of simulated and observed hydrographs, as well as simulated and constructed groundwater tables; (2) identification of the calibration parameters to be changed and a decision for how to change them; (3) adjustment of one or multiple calibration parameters according to the decision made in Step ( 2); (4) simulation run.
Within the first calibration runs, the modeler gathers experience for how sensitive the model reacts with respect to certain parameter values.For the decisions on how to change a certain calibration parameter (Step (2)), it has been assumed that conductivity values in the river bed are higher than at the bank according to the findings of Giebel and Hommes [40] and Ubell [22].Decisions were based on the match between observed and simulated hydrographs at the observation points and the match of simulation results to the groundwater tables from Ubell [22].
In order to verify the calibration result, an automatic calibration has been carried out as well.For the automatic calibration, a Feflow module [54] that contains the PEST algorithm [55] has been used.PEST regards the calibration as an optimization problem.With a gradient type optimization algorithm, PEST minimizes the least squares of difference between observed values (the calibration data) and simulated values.Since the optimization algorithm does not necessarily find a global minimum of the objective function, initial parameter values have an impact on the calibration result.According to Ubell [22], we used a homogenous hydraulic conductivity of k = 3.0 • 10 -2 m /s as the initial value; for the leakage boundary condition model variant, an initial leakage parameter value of 1.0 1 /s has been applied.In order to avoid calibration results that physically make no sense, a confidence interval has been set as a constraint for each calibration parameter: For the automatic calibration, the parameter zones 4 to 8 were treated as one zone, because for parameter Zones 5 and 7, no corresponding observation hydrograph is available.Test runs showed that the optimization algorithm found results that physically make no sense if there are zones where no calibration data is assigned to.
The Feflow calibration module cannot differentiate between c l,ex and c l,in .Therefore, for the leakage model variant, the automatic calibration has been carried out separately for the infiltration phase and the exfiltration phase of the bank storage event.

Calibration Results
In Table 4, the Nash-Sutcliffe efficiency for the seven observation points for each model variant is given, both for the manual calibration and the automatic calibration.Figure 13 shows a graphical representation.The Nash-Sutcliffe efficiency takes values between −∞ and one, where a Nash-Sutcliffe efficiency of one indicates an exact match between observed and simulated values.In the case of a negative Nash-Sutcliffe efficiency, the observations are better represented by their average value than by simulated values [52].4) for the observation points obtained with manual and automatic calibration for the two model variants.
For observation points U01 to U06, a very good calibration quality for all model variants and both calibration methods has been achieved.The Nash-Sutcliffe efficiencies for observation point U07 indicates a good or acceptable match for the manually calibrated models and the automatically calibrated head boundary condition model as well; however, for the automatically calibrated leakage model, the obtained Nash-Sutcliffe efficiency of 0.103 is not sufficient.Regarding the sum of the Nash-Sutcliffe efficiencies (Table 5), the manual calibration seems to perform better than the automatic calibration; this is not surprising, as in the manual calibration more parameter zones were applied than in the automatic calibration.The manual calibrated leakage model variant gives the best calibration result in terms of Nash-Sutcliffe efficiency.
Figure 12 shows a comparison of simulated groundwater tables along the cross-sectional profile from the model with the head boundary condition (Figure 12a) and the leakage boundary condition model variant (Figure 12b) with the groundwater tables Ubell [12] derived from observations.On 13 April, the river stage reaches its maximum.On this day, the groundwater table simulated with the head boundary condition model matches quite well with the groundwater table based on observed values for the whole profile.Two and a half days later, on Day 15, however, the simulated water table lays below the one based on observed values for positions between 0 m and 1400 m.For the point in the simulation time of 26 April, 24:00 h, the match between the simulated groundwater table and the one based on observed data is even worse.Changing the hydraulic conductivity for bank elements would improve the data match for 26 April, but would also affect the water exchange during the infiltration phase.This would worsen the match between simulation results and observations for Days 13 and 15 of April, 1983.Because the leakage approach used here allows one to address the leakage parameters c l for infiltration and exfiltration independently, with the leakage boundary condition model a better match between the simulated groundwater tables and the groundwater tables derived from observed data could be achieved (Figure 12b).As shown in Section 3.3, infiltration and exfiltration follow different flow patterns, so assigning different parameter values for infiltration and exfiltration seems more appropriate for the current case.The water table in the hinterland of the river bank (positions further than 2600 m from the river in Figure 12) for 26 April is generally not very well represented by simulation results.At this point in time, low flow conditions with groundwater currents parallel to the river occur (see [38]) that are not represented in the models, because the model represents a cross-section perpendicular to the river only.This also might explain the low Nash-Sutcliffe efficiency values obtained for observation point U07.
Figure 9 shows the parameter value distribution of the head boundary condition model variant obtained by manual and automatic calibration.Contour levels of the finite elements indicate the value of hydraulic conductivity in the pictures.The corresponding parameter values are given in Table 3.The parameter distribution differs with the calibration method: while in the manual calibration larger conductivity values are allocated to the river bed, the PEST calibration algorithm has allocated high conductivity values to the upper bank.A similar distribution pattern of calibration parameter values has been obtained with the leakage boundary condition model variant (Figure 10): for the infiltration process, high leakage parameter values are located at the river bed in the manually calibrated model, while in the automatically calibrated model, the high parameter values are located in the bank.
The parameter value distribution has an effect on the simulated Q-∆h-relationship: Figure 14a shows the leakage functions of the manually calibrated models, and Figure 14b shows the two leakage functions obtained with the automatically calibrated models.Only the simulation results of the automatically calibrated models show a distinct hysteretic course.In the beginning and at the end of the flood event, the top element layer contributes only marginally to the total leakage flux.With a rising groundwater level, the calculated interface area A (Equation ( 12)) increases, but also the top element layer has more impact on the leakage flux calculation (see Figure 3).The optimization algorithm addressed this element layer in order to achieve disproportionally high leakage fluxes during high water, which results in a hysteretic leakage function.In the manual calibration, however, the upper element layers were assigned low values of hydraulic conductivity, according to the findings of Giebel and Hommes [40] and Ubell [22], so the impact of a larger interface area during high groundwater levels is not as great.Figure 14a also shows that the head boundary condition model variant leads to a uniform slope for negative and positive leakage fluxes, while the leakage function obtained from the manually calibrated leakage boundary condition model gives a different slope for infiltration and exfiltration, respectively.Around the origin, there are outlier values in the curves from models with leakage functions.These outliers arise from the fact that the groundwater level is still unknown when the leakage parameter is chosen according to the flow direction.Feflow does not carry out iterations, even if it turns out that the parameter for the wrong flow direction has been applied.

Discussion
The Nash-Sutcliffe efficiency and the comparison with constructed groundwater tables indicate that the manually calibrated leakage model variant performs better for the current case than the corresponding model with a head boundary condition.The fact that the leakage boundary condition model variant allows one to address infiltration and exfiltration independently results in a leakage function that matches better with the one derived from observed data, but is also advantageous for the practical calibration work.However, the leakage approach contains the assumption that the conductivity properties of the river bed change with the direction of exchange water flow, which is not necessarily the case, as discussed in Section 3.3.
The manually calibrated models result in a leakage function with a linear and (nearly) non-hysteretic course.The high Nash-Sutcliffe efficiencies obtained with the manually calibrated models indicate that the hysteretic course of the leakage function is not essential to reproduce the measured groundwater levels at the observation points.In the automatically calibrated models, the hysteretic course of the leakage function has been achieved by addressing the elements located at the upper layers, which are only active during high water.The resulting leakage functions match better with the leakage function derived from observed values.
The groundwater model geometry does not represent the real cross-section of the river very well.Doble et al. [56] show that the slope of the river bank itself and the corresponding variation of the interface area have an impact on the exchange of water between river and aquifer.A disproportional high increase of the wetted perimeter p for high river stages, for example when the flow situation changes from bank full discharge towards an inundation of forelands, is not incorporated in the model, because the vertical mesh boundary is required for three-dimensional models in Feflow.High permeability values in the upper parts of the bank (as obtained with the automatic calibration) can compensate for this technical limitation.This allows one also to account for preferential flow paths created by animals or ice crush (see Section 3.3) in the upper part of the bank, which are activated during high water only.However, for the current case, there is no evidence for such a preferential flow path in the literature, and the parameter distribution obtained by automatic calibration does not match qualitatively with the available information on the conductivity properties of the bank and river bed [22,40], whereafter the higher conductivity values are located in the river bed.
A possible change of hydraulic conductivity properties during the flood event due to erosion and sedimentation processes and biogenous stabilization are not taken into account in the groundwater model.A coupled simulation that incorporates river morphology [47] would be necessary to represent this processes.Furthermore, temperature effects have not been considered in the model.We have no information on the relevance of these processes for the current case, but both processes might be responsible for the hysteretic behavior.Doble et al. [56] show that the capillary fringe above the groundwater table has a damping effect on the rate of flux that infiltrates the bank, because the soil must become saturated before it can convey substantial volumes of water.On the other hand, it takes longer to drain water that is stored in the unsaturated zone than in the case of saturated flow.The groundwater models used in this study assume fully saturated flow.Consequently, the effects of the unsaturated zone on bank storage processes are not taken into account.We have chosen to assume fully saturated flow for two reasons: Firstly, the aquifer characteristics with its comparatively high values of hydraulic conductivity legitimate this assumption.The calibration results show that a saturated flow model is able to reproduce the bank storage event in a good way.Secondly, modeling of variable saturated flow is not feasible for many practical applications, because the extent of the modeling area is too large and variable saturated flow modeling is computationally too expensive or because there is a of lack of data for unsaturated flow properties.
Regrettably, none of the model variants represents the exfiltration process of the analyzed case satisfactorily.Figure 6 shows that the leakage flux still changes after Day 16, while the head difference remains nearly constant.The effects of the bank slope [56] that are not taken into account due to the assumption of a vertical bank may play a role here.There are no model elements that can be activated only for the process of exfiltration in the way that it is possible for the infiltration process under high water levels.However, in the case of exfiltration, only small absolute values of leakage flux occur, so for the current case, the impact of the exfiltration process on the groundwater flow situation is small in comparison to the infiltration process.Thus, the exfiltration process does not affect the calibration quality that much.
Different calibration results lead to a good or acceptable match between observed and simulated values, although the parameter distribution is different.This shows once more the non-uniqueness of the solution of the calibration problem (see [57]).The validation of groundwater models with independent data and the post audit of the modeling concept [58] become important if models are applied for predictive analysis.

Summary, Conclusions and Outlook
A subsurface flood is a severe groundwater head rise due to a flood event in a river.Bank storage is one of the important water exchange processes between a river and the contiguous aquifer that governs the subsurface flood propagation.For predictive subsurface flood simulations, it is important to know how to represent the water exchange in a groundwater model.Against this background, we addressed the question of how bank storage can be represented in a numerical groundwater model.
Ubell [12] reports a bank storage event for the Neuwieder Becken in the Middle Rhine area (Germany).The leakage function is hysteretic, and during the infiltration phase, higher leakage fluxes and larger head differences occurred than during the exfiltration phase.Note that we analyzed one bank storage event of a specific site only, so the outcomes of this study are not necessarily transferable to other sites.For Ubell's bank storage event, we set up a groundwater model where we applied the river stage as: • head boundary condition; and • leakage boundary condition.
With both approaches, the groundwater model could be calibrated to a good match between observed and computed groundwater levels.This is not surprising, because both approaches are based on a similar conceptual model.For practical applications, however, the leakage boundary condition is recommended, because infiltration and exfiltration can be treated as two different processes.This matches better with the observations and makes the practical calibration work easier.Another advantage is that the interface layer is represented independently of the mesh, which will be relevant especially for models on a larger scale (see, e.g., [30]).
With our model, we obtain a hysteretic characteristic of the leakage function with a vertical variation of leakage parameters or hydraulic conductivity properties, respectively.This approach can represent a variable interface area that changes with the river stage and preferential flow paths that are only active during high water, but not a dynamic change of leakage properties due to erosion, sedimentation and biogenous stabilization.However, the simulations show that a parameter distribution that results in a linear leakage function is sufficient to represent the water interchange between river and aquifer for the analyzed case.
In many conceptual models for the water exchange processes between river and aquifer that are applied in practice (i.e., the program features available in groundwater models), the river is only represented with its water level (see Sections 2.3.2).Consequently, some relevant processes are neglected.These approaches should be further developed in such a way that not only the river stage, but also: • heterogeneous conductivity properties of river bed and river bank; • the relation between river stage and interface area and • a dynamic change of leakage properties (as time series or in dependence of flow properties) are considered.Some new approaches have been published recently: Affolter et al. [14] apply leakage parameters that vary with the river discharges in order to account for the temporal changes of the hydraulic conductivity of the river bed during high water.Simpson and Meixner [47] present a quantitative approach to model the dynamics of stream bed hydraulic conductivity with the help of sediment transport modeling.The hydraulic conductivity of the stream bed is derived from the particle size distribution, layer-specific bulk density values and a pedotransfer function.
More field observations are desirable to verify existing approaches on bank storage modeling and to extend the conceptual understanding of bank storage processes.The analysis of the bank storage event presented in this paper is based on the volume balance of a two-dimensional vertical cross-section only.As a next step, it is planned to extend the coherent analysis of river stage measurements and groundwater level observations to the horizontal plane.This allows one to take the effects of groundwater flow components parallel to the course of the river into account as well.limitations arise from the fact that the devices must be installed in great water depths.High current velocities during high water conditions make direct seepage measurements a challenging task.
Rushton and Tomlinson [33] derive the leakage flux by base flow separation for normal conditions.However, Mattheß and Ubell [41] point out that usage of base flow separation methods requires a hydraulic gradient from the aquifer towards the river, i.e., effluent conditions, during the observation period.This condition is usually not satisfied during a high water period.
The estimation of the leakage flux by comparison of incoming and outgoing discharge of a river segment has a limited accuracy [26], because in most cases, the leakage flux is small compared to the river discharge, even if the river cuts an aquifer with high values of hydraulic conductivity.Another problem is the exact determination of the river discharge, because in the case of high water, the relation between the river stage and the river discharge is not unique due to transient conditions [59].

A2. The Volume Balance Method
Ubell [12] derives the leakage flux from observation wells that are installed in a line orthogonal to the river (see Figure 15).From the change in groundwater level, he concludes the volume change in the aquifer.Figure 15 shows the initial situation at time t 1 , where the river stage is h 01 , and the groundwater level for the observation points is given by h 11 to h 41 .At t 2 , the increased water level in the river causes water infiltration into the aquifer.The groundwater level rises to h 12 , h 22 and h 32 , while for observation point 4, the water level does not change.The groundwater storage between t 1 and t 2 is given by the area A. The change in water volume V can be calculated by: where n P is the effective porosity of the aquifer.The volume change per time unit is the leakage flux: with the length of the river section l.Ubell's method does not allow one to separate precipitation and the backwater of base flow.Ideal conditions for the application of this method are a groundwater table that is initially horizontal and no precipitation.Another important prerequisite is a groundwater flow pattern that is oriented along the chain of observation points for normal conditions as well as for flood conditions.

Figure 2 .
Figure 2. Leakage flux Q over head difference ∆h for the Darcy leakage approach (Equations (4) and (6)), the approach by Rushton and Tomlinson[33] (Equation (9)) and the leakage approach that uses the wetted perimeter (Equation (8)), evaluated here with a trapezoidal cross-section.

Figure 3 .
Figure 3. Wetted perimeter p * and leakage rate q for different water levels h R and h G , but the same head difference ∆h.The interface layer between the aquifer and the river is modeled as heterogeneous (c 1 > c 2 ).(a) Low river stage and low groundwater level; (b) high river stage and high groundwater level.

Figure 4 .
Figure 4. Overview of the Neuwieder Becken and the chain of groundwater observation points near Urmitz (after[12]).The observation wells U01, U03 and U05 are referred to in Figure6.

Figure 5 .
Figure 5. Hydrogeological cross-sectional profile of the Rhine valley near Urmitz along the chain of groundwater observation points shown in Figure4(modified after Ubell[12]; elevation given in meters above official sea level reference m ü.NN).

Figure 6
Figure 6.Hydrograph of the Rhine gauge (in meters above official sea level reference m NN) and three groundwater observation points with different distances to the river bank, respectively (values taken from[12]); leakage flux for the 1 km broad profile "Urmitz" as day-averaged values (values taken from[12]) and the head difference for the observation point located at the bank.

Figure 7 .
Figure 7. Head difference ∆h over leakage flux Q for observation point U01 (data taken from [12]) and the smoothing function.(a) Infiltration towards the aquifer; (b) exfiltration towards the River Rhine.

Figure 8 .
Figure 8. Head difference ∆h over leakage flux Q for groundwater observation point U01.The points are connected according to the time line of the flood event.Numbers at the data points indicate the days in April, 1983.

Figure 9 .
Figure 9. Finite element discretisation of the groundwater model representing the vertical cross-section of the aquifer along the chain of observation points near Urmitz.The pictures show the model variant where the river water interchange is modeled with a head boundary condition.The coloring of the finite elements indicate values of hydraulic conductivity after automatic and manual calibration, respectively (see Section 4.5 and Table 3).(a) Parameter distribution results of the manual calibration; (b) parameter distribution results of the automatic calibration.

Figure 10 .
Figure 10.Detail of the 3D groundwater model "Urmitz", where the river is modeled with a leakage boundary condition with calibration results for the leakage parameters c l,in .The coloring of the finite elements indicates the values of leakage parameters obtained by manual and automatic calibration (see Section 4.5 and Table3).(a) Results of manual calibration; (b) results of automatic calibration.

Figure 11 .
Figure 11.Zones for the calibration of the hydraulic conductivity k and leakage parameters c l .

Figure 13 .
Figure 13.Graphical representation of the Nash-Sutcliffe efficiency (Table4) for the observation points obtained with manual and automatic calibration for the two model variants.

Figure 14 .
Figure 14.Q-∆h-diagrams from the automatically and the manually calibrated model variants with leakage and head boundary conditions for the River Rhine.Numbers at the data points indicate the day in April 1983.(a) Manually calibrated; (b) automatically calibrated.

Table 1 .
Values used to calculate the leakage fluxes for the diagrams in Figure2.

Table 3 .
Assignment of parameter zones for the model/calibration variants and calibration results.

Table 4 .
Nash-Sutcliffe efficiency of the observation points obtained by manual and automatic calibration for the two model variants.

Table 5 .
Sum of Nash-Sutcliffe efficiencies over observation points U01 to U07 from