{'fraction_non_alphanumeric': 0.04894883176430418, 'fraction_numerical': 0.016974781717964977, 'mean_word_length': 4.030425714636241, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A two-layer approach for Coupling 1D/2D Shallow Water Flow Models

In this paper, we propose a novel approach for coupling 2D/1D shallow water flow models. Efficiently coupling these models is vital for simulating the flow and flooding of open channels. Currently, existing methods couple the models either at the channel lateral boundaries (lateral methods) or at the location, along the channel flow direction, where the two sub-domains intersect (frontal methods). We classify these methods as horizontal methods since the coupling points are on the horizontal plane. The limitations of these methods include their inability to recover the 2D channel flow structure during flooding without losing efficiency, and their inability to switch between lateral and frontal methods, based on the problem. Here, we propose a new paradigm which is to think of the channel flow as a two layer flow. This leads to the vertical coupling method (VCM) which is able to recover 2D channel flow structure without losing efficiency, switch types depending on the flow and is a superset of some existing methods. The VCM is based on a user chosen elevation above which the user considers the channel to be full. From this elevation, all other flow quantities are defined including the two layers. The flows in the lower and upper layers are assumed to be 1D and 2D respectively, and the appropriate flow models with exchange terms are derived. The 2D shallow water flow model is retained for the floodplains. Finite volume methods (FVM) are formulated for the flood model; the FVM with the operator splitting approach are also formulated to solve and couple the two layer channel flow models. We show that the resulting method (i) is well-balanced (ii) preserves the"no-numerical"flooding property (iii) preserves conservation properties and (iv) adapts to the flow situation.


Introduction
One dimensional Saint Venant models are often used to simulate open channel flows but they become inadequate once the channel overflows. However, multi-dimensional (even 2D shallow water) simulations are computationally expensive. This has led to the development of methods to couple 1D channel simulations with 2D floodplain simulations. River/floodplain coupling simulations started as quasi 2D models in which floodplains are represented by storage cells, then coupled with an existing 1D river model [9,4]. These approaches would not allow to simulate the fluid dynamics in the floodplains [10].
In [5], the coupling of a 1D channel model with a full 2D floodplain model was achieved by including the 2D numerical fluxes into the finite volume scheme for the 1D model, while [8] utilized the theory of characteristics to couple 1D/2D models through suitable matching conditions defined at the 2D/1D interface. In [20], the 1D river model and a 2D non-inertia model were also coupled where the water level differences between the flows in the two domains are used to calculate the interacting discharges in the sub-domains.
Methods based on post-processing the separately computed solutions are presented in [15], see also [14]. At 2D/1D interface, each model computes its own solution from which the total water volume in a 1D cell and all its adjacent 2D cells is computed. Then, the water height for 2D cells and the wetted Area for the associated 1D cell are found. In [17], the methods have been applied to Tiber River, Rome. In [13], the 1D model including the coupling terms were classically derived from the full 3D inviscid Euler's equations and an optimal control process applied to couple the models. These models have been numerically treated with the finite volume method [10] where the discrete exchange term, which lead to globally well-balanced scheme, is proposed. This approach superposes a 2D grid over the 1D channel grid and convergence is achieved using a Schwartz-like iterative algorithm.
A major difficulty in coupling 2D/1D shallow water models is the computation of the channel flow lateral discharge. In [11], this channel lateral discharge was set to zero, while [12] adopted an iterative technique that uses the solution of successive Riemann problems to estimate the transverse velocity. This difficulty in computing the channel lateral discharge remains challenging to compute accurate fluxed between the channel and the floodplain.
A further issue is the assumption that the flow will remain one dimensional in the channel even in the region, where a flooding event is occurring. In [19] a method to compute different lateral discharges at each channel boundary was purposed. However, the free-surface and xvelocity component were still assumed to be laterally constant across the channel and although the approach does improve the accuracy of the method, it still does not recover the complete 2D flow structure during flooding.
Efficient methods that recover the 2D flow structure during flooding but revert back to 1D simulation if no flooding is occurring could solve both the issues mentioned above. Frontal coupling methods in which the floodplain extends into the channel in parts of the domain, recover 2D flow structure but they loose efficiency because they compute the 2D solutions at all times. Moreover, most of the existing methods need to know the location of possible flooding a-priori and can not take into account that flooding locations may vary with time. So methods that can adapt to flooding regions a-posteriori are desirable.
The goal of this paper is to propose a method, the vertical coupling method (VCM), that (i) recovers the 2D flow structure during flooding while reverting back to 1D simulation in the channel regions where no flooding occurs. (ii) automatically detects flooding regions. (iii) can be easily added to well established 1D channel and 2D floodplain flow solvers.
The VCM is based on partitioning the flow in an overflowing channel into two vertical layers where the flows in the lower and upper layers are simulated using 1D and 2D models, respectively. We then derive a coupled models for the two layers which represents the channel flow. The upper layer model can be easily coupled to a 2D floodplain model. The method is parameterized by prescribing a height function z w b (x) along the channel which determines at which water height the channel is in danger of overflowing. Only when this water height is reached will the upper layer model be activated, below that layer the channel flow is simulated using a standard 1D model. Different choices for z w b (x) will lead to different method and special choices will lead to some standard lateral or frontal coupling methods found in the literature. In this sense VCM is a superset of some existing methods.
The rest of the paper is organised as follows: The two layer channel model is presented in section 2. We start by presenting the notation used to derive the VCM in section 2.1, derive the lower and upper layer flow models in sections 2.2 and 2.3 respectively, the complete two layer channel flow model is summarised in section 2.4. For the floodplain flow we use a standard shallow water model which is summarized in section 2.5. A numerical algorithm for the coupled channel flow models is presented in detail in section 3 and the properties of the method are considered in section 4 where we prove that the method is well-balanced, preserves no-numerical flooding and is mass conservative. Numerical experiments are presented in section 5 to evaluate the performance of the method compared to simpler approaches found in the literature. Finally, we give a summary in section 6.  Let Ω H ⊂ R 2 be a fixed 2D horizontal domain with fixed bottom elevation, z b ( X), X = (x, y) ∈ Ω H . Let H( X, t) ≥ 0 denote the depth of water at point, X ∈ Ω H at time, t ≥ 0, so that is the water free-surface elevation at point, X, at time t. Then, at time, t the flow occupies the 3D domain, Ω t defined by A cross section of the flow domain, Ω t at a fixed x is shown in figure 1(a).
2.1.1 Channel Geometry: Figure 1(b) shows the flow cross section indicating the floodplain part and the channel cross section. To simplify the derivation of the models we assume its length lies along the x-axis (frontal direction) and the width along the y-axis (lateral direction). We define The functions y l (x, z) and y r (x, z) are the y-coordinates of the left and right lateral wall boundaries respectively at the elevation, z. The function, B(x, z) gives the channel cross sectional lateral width, i.e., and for convenience we set see figure 1(b). c by the authors, 1677 2.1.2 The VCM Background Given the channel geometry above, the formulation of the VCM starts by choosing an elevation, z w b (x) ≥ Z b (x) above which the channel is considered full, see figure 1(b). This elevation is completely decided by the user; the only constraint is for it not to be less than Z b (x). So, z w b (x) becomes the channel top and is referred to as the maximum wall elevation defined below. Definition 2.1 (Maximum channel wall elevation, z w b (x)) The maximum channel wall elevation at cross section x, is the minimum elevation of the channel banks above which flooding is said to have occurred [18].
Once z w b (x) has been chosen, other quantities are derived. The lateral width and the y-coordinates of the lateral wall boundaries at the top (z = z w b (x)) becomes B(x) := B(x, z w b (x)) and y w l,r (x) := y l,r (x, z w b (x)) respectively. For the VCM we then partition the flow into the channel flow within the region y w l (x) < y < y w r (x), and the flow in the floodplain the remaining region, namely −∞ < y < y w l (x) and y w r (x) < y < ∞, see figure 1(b). As noted above, the floodplain flows will be simulated using a standard 2D shallow water models (given in section 2.5), hence we concentrate on deriving the models for the flow in the channel region.

Remark 2.1
The idea of the VCM method is to use a 1D channel model for regions in the domain where the channel is not full while using a two layer model otherwise, where only the lower layer is evolved using a 1D model while the upper layer uses a 2D model. As said above the choice of when to switch from a one layer to a two layer model depends on the choice of z w b (x) and thus is up to the user. The closer the chosen z w b (x) is to Z b (x), the smaller the channel flow region y w l (x) < y < y w r (x), hence the larger the floodplain flow region (see figure 1(b)). In particular, choosing z w b (x) = Z b (x) in some part of the channel, will result in an approach where the 2D floodplain model is used in parts of the channel (a frontal type coupling approach), while taking z w b (x) = ∞ will result in a one layer 1D model being used in that region leading to a horizontal type coupling [18,19]. Thus, different choices of z w b (x) lead to different coupling approaches.
Having identified the channel and floodplain flow regions, we now concentrate on the channel region. First, we note that the following condition holds for fixed x: We also extend the definition of the width functions, B(x, z), y l,r (x, z) to the region above the top (that is where z ≥ z w b (x)) as follows: see figure 1(b). This results in a channel with straight vertical walls above the height where the channel is assumed to be full.
For the scheme later on we need to assume that the channel geometry is known such that we can reconstruct the water height given the wetted cross section: Definition 2.2 (Water height for given wetted cross section) Let A = A(x) be a given function describing the wetted area in the channel defined by the bottom topography z b (x, ·). We define the height function H( X;

Remark 2.2
With the above definition H( X) − z b ( X) is independent of y and represent a 1D flow with the wetter area given by A. Due to the way we extended y l,r (x, z) for z ≥ z w b (x) it is easy to see that the above definition provides a unique height function H even for large value of A, i.e., even in the case where the given wetted cross section leads to a overfull channel.

Important Quantities for the Channel Flow
We now proceed to define other important quantities for the channel flow.

Definition 2.3 (Channel Depth)
The channel depth, β( X) is the laterally varying height between the channel bed and the chosen elevation, z w b (x), that is see figure 1(b).
We also introduce the critical area defined as follows.

Definition 2.4 (Critical Area)
The critical area, A c (x) is the wetted cross sectional area of an exactly filled cross section. That is, the wetted cross sectional area when the water level is exactly at the chosen elevation, z w b (x). It is defined by

Remark 2.3
If the channel is exactly full (A = A c ), then the water depth is exactly equal to the channel depth, β, and consequently H( X; A) = β( X).

Definition 2.5 (Channel Flow Lateral Boundaries)
Let y * l (x, t) and y * r (x, t) denote the y-coordinates of the left and right channel flow lateral boundaries, i.e., y * l (x, t) := min{y : η( X, t) > z b ( X), y ≥ y l (x, z w b (x))}, y * r (x, t) := max{y : η( X, t) > z b ( X), y ≤ y r (x, z w b (x))}. (11) Note that if there is water everywhere in the channel, then y * l,r (x, t) = y w l,r (x).

Remark 2.4 (Channel Assumption)
We assume that the channel never goes dry anywhere between the lateral wall boundaries, i.e., there are no islands in the channel: Definition 2.6 (Average Free Surface) With the flow lateral boundaries known, we can define the laterally averaged free surface elevation, namelȳ Next, we define precisely what we mean by a full channel.

Definition 2.7 (Full Channel)
We say a channel cross section at x is full if H( X, t) ≥ β( X) ∀y ∈ y w l (x), y w r (x)) , This means that if the channel is full the following statements hold: • given the total wetted area A(x, t) in the cross section and the critical area A c as given in definition the following inequality holds (recall definition 2.4): • alsoη(x, t) ≥ z w b (x) since by definition

Vertically Partitioned Channel Flow
The main idea of the VCM is to partition the channel flow into two layers, based on the chosen channel top elevation, z w b (x). Since we want to approximate the flow with a 1D model if the channel is not full, we make the following assumption to allow consistency with a 1D formulations.

Definition 2.8 (1D Consistency Assumption)
If the channel is not full, then the lateral variation in free surface elevation, η( X, t) is negligible and the free surface can be taken to be its lateral average,η(x, t).

Remark 2.5
This simply means that we assume the free surface to be laterally flat. This is exactly the 1D assumption and allows to apply 1D modelling whenever the channel is not full (a similar assumption is made for example for the 1D models in [15,5,9]).
We can now identify the channel flow domain, Ω ct at time, t as: Definition 2.10 (Flow Partitions) Given the total water depth, H( X, t), we partition the flow into the two layers see figure 1(b). The water depths in the lower and upper layers are given by h 1 ( X, t) = min(H( X, t), β( X)), y l (x, η 1 (x, t)) ≤ y ≤ y r (x, η 1 (x, t)), h 2 ( X, t) = H( X, t) − h 1 ( X, t), y l (x, η 1 (x, t)) ≤ y ≤ y r (x, η 1 (x, t)), Note that h 1 ( X, t) + z b ( X) = η 1 (x, t) defines the top of the lower layer and so by our 1D consistency assumptions h 1 is independent of y.
(a) Full cross section; η 1 (x, t) equal to the channel wall elevation z w b (x) and lower layer water depth, h 1 (x, y, t) equal to the the channel depth β(x, y). Free surface η(x, y, t) remains non flat laterally. Channel flow boundaries y l,r (x, η 1 ) equal to the channel top lateral walls, y w l,r (x).
(b) Non full cross section; laterally flat free surfacē η(x, t), η 1 (x, t) equal to the laterally flat free surfacē η(x, t), which is less than the channel wall elevation z w b (x), h 1 (x, y, t) equal to the total water depth, H(x, y, t) which is less than the channel depth β(x, y). The channel flow lateral interval, |y r (x, η 1 ) − y l (x, η 1 )| is less than the total channel lateral interval |y w r (x) − y w l (x)|.
where (u, v, w) T is the fluid velocity vector at point ( X, z) ∈ Ω ct at time t. Furthermore, the following boundary conditions hold

The Lower Layer Flow Model
The lower layer flow is assumed to be always one dimensional and thus the lateral variations in the free-surface elevation do not have any impact on the lower layer flow. Therefore, using the averaged free-surface,η(x, t), instead of η( X, t). We thus replace the horizontal moment equation (22) and the kinematic boundary condition (25) by We define the area of wetted cross section, A 1 (x, t), the section averaged volumetric discharge, Q 1 (x, t) and section averaged velocity, u 1 (x, t), for the lower layer, as follows: Before we proceed, let us state the following important relations which are easily proven (see [18] for details): (9) and (14) respectively, and

Lemma 2.1 Let A c and A be defined in
Theorem 2.1 Let A c , A, A 1 and A * c be defined in (9), (14), (29) and (31) respectively, then

Mass Conservation Equation for
Lower Layer: First we derive the model equations for the lower layer. Integrating the equation (21) over the lower layer cross section, and applying the Leibniz rule using the kinematic boundary condition (24), we obtain where Remark 2.6 Note that´y r (x,η 1 (x,t)) y l (x,η 1 (x,t)) S( X, t)dy is the averaged mass exchange term between the two layers. If the channel is not full, i.e., η 1 (x, t) =η(x, t) then S( X, t) = 0 using the boundary condition (27)

Momentum Equation for
Lower Layer: To derive the momentum equation for the lower layer, we integrate equation (26) over the lower layer cross section and using standard arguments arrive at where the last integral on the right is the momentum exchange term between the two layers. c by the authors, 1677

Upper Layer Flow Model
The upper layer flow is allowed to remain fully two-dimensional, so the Free-Surface Euler Equations, (21) -(23) and (25) are applicable. However, the kinematic boundary condition on the bottom does not apply here because the bottom of the upper layer is η 1 (x, t) which is not a physical boundary that fluid particles cannot cross. Define the following quantities: So that the velocities are where q 2 ( X, t) = (q 2x ( X, t), q 2y ( X, t)) T is the upper layer 2D discharge vector and u 2 ( X, t) = (u 2 ( X, t), v 2 ( X, t)) T is the velocity vector in the upper layer.
In the following, we derive the equations for the 2D quantities. Integrating equation (21) vertically over the upper layer, we have Also, integrating equation (22) vertically over η 1 (x, t) ≤ z ≤ η( X, t), applying the kinematic boundary condition (equation (25)) and simplifying, we have Similarly, integrating equation (23) vertically over η 1 (x, t) ≤ z ≤ η( X, t), applying the kinematic boundary condition (equation (25)) and simplifying, we have The mass exchange S( X, t) between the layers is defined in equation (35).

Summary of Coupled Channel Flow Models
The two layer channel model derived so far consisting of the 1D lower layer model Sdy, and the 2D upper layer model, where the fluxes are given by Note that the models (42) and (43) are not closed since the exchange term S and interface velocity u η 1 are not known. Following a similar idea as used in [3] we will solve the system using a two step approach where in the first step we solve the equations without the exchange term and and then use this intermediate solution to approximate the mass and momentum exchange between the layers. But even then there is a difficulty in solving the 1D lower layer model since the free-surface elevation term,η(x, t), appearing in (44) is different from the actual top level, η 1 (x, t) of the lower layer flow, which is represented by A 1 . This makes directly deriving a wellbalance scheme for the lower layer model challenging and more importantly we cannot directly reuse an existing 1D channel solver for the solution. Since this is one of our major goals, we base our numerical scheme on a reformulation of the above coupled channel model: Integrating (43) over the cross section y l , y r we arrive at where assuming that the average horizontal velocities in both layers are very similar, i.e., ( A 2 ) 2 ≈ 0 we can add the equations for A 1 , Q 1 and A 2 , Q 2 together arriving at a model for the 1D full channel So with this further approximation we end up with the following systems to model the flow within the channel: where Remark 2. 7 The system (49) is a standard 1D approximation to the channel flow with Φ providing the exchange terms between the channel and the floodplain along the horizontal boundary. Thus a standard model can be used to approximate (49). The advantage of our approach is that the second set of equations (50) provides additional 2D information in the channel that can be used to accurately compute the exchange terms Φ.

Floodplain Flow Model
We describe the flow in the floodplains using the 2D Shallow water equations, namely where (52)

Numerical Algorithm for the Coupled Channel Flow Models
In this section, we describe in detail the numerical schemes and algorithms to solve the two-layer models, (49) and (50) for the channel flow. The 1D model (49) is solved on a 1D (quasi 2D) channel mesh (see figure 3(a)), while the 2D upper layer model (50) is solved on a 2D channel mesh (see figure 3(b)). The results are effectively combined to derive the update values for the channel flow. Let Ω 1 h be a 1D (quasi 2D) channel mesh with cells K i = (x i−1/2 , x i+1/2 ), see figure 3(a). To consider the upper layer channel model, (43), we further discretize each 1D channel cell K i into 2D cells, T ij , j = 1, ..., N 1 y , see figure 3(b). This forms the 2D channel mesh Ω 2 h , so that

Background
For the 1D cell, K i , we define the following discrete channel quantities: as the channel wall elevation, channel top width, critical area, 1D bed elevation, and channel depth, respectively, evaluated at the center x i of each 1D cell, K i (see figure 4(a)). We also define the discrete quantities denoting approximations to the average total wetted cross section and section-averaged discharge in K i .
For the associated 2D cells, T i j in K i , we define two similar sets of discrete quantities being the bed elevation and channel depth at the cell center, X i j of T i j , see figure 4(b). Note the identity, where q 2 = (q 2x , q 2y ) T and |T i j | is the size(area) of T i j . Finally we also define full 2D cell averages given by which are cell averages for the full 2D data (sum of lower and upper layer) in 2D cells, T i j ∈ Ω 2 h .  Given a lateral distribution η n i, j of the free-surface elevation in the 2D cells T i j ⊂ K i for j = 1, . . . , N 1 y we assume that we can compute the wetted cross sectional area, A n i in K i by This requires information about the channel geometry, for instance, if the channel has a rectangular cross-section, with a laterally constant bottom elevation, z bi j = Z bi for all j, then A is given by Next, we formulate a discrete 1D consistency assumption similar to definition 2.8: c by the authors, 1677

Definition 3.1 (Discrete Consistency requirement)
The solution, (H, q x , q y ) n i, j at t n satisfies a discrete consistency requirement, if the following condition holds. If ∃ j * ∈ {1, 2, ..., N 1 y } such that η n i, j * < η β i , then We will later show that if the initial data satisfies this discrete consistency requirement, then the requirement is satisfied for all time by our scheme.

Solution of the two-layer Channel Models
We now solve the two layer channel models, (49), (50) to evolve the full 2D data, (H, q) i j from the initial to the final time. To this end we first initialize the value (A 0 i , Q 0 i ) i for the full 1D model and the values (h 2,i, j i 0 , q 0 2,i, j ) i j for the 2D upper layer model (in the region where the channel is full). To evolve these values to the next time level we first update the 2D upper layer values ignoring the exchange term between the two layers but taking the horizontal fluxes with the floodplain into account. Next we update the values of the full 1D model, using the channel-floodplain lateral fluxes computed by the 2D upper layer solver to completely solve (49). Finally, we correct the solution in the upper layer by estimating the exchange term between the two layers. Note that the first step only leads to computational cost in regions where the channel is full (i.e. h n 2,i, j > 0). Otherwise this step can be skipped and only the 1D step has to be carried out. Exchange terms with the floodplain are also correctly taken into account in the first step.

3.2.1
Initializing the two-layer model from given full 2D data Given full 2D data, (H, q) 0 i, j at the initial time and z bi j in each 2D channel cell T i j ∈ Ω 2 h , we initialize the 1D model, (49) and the 2D upper layer model, (50) as follows: and where u and ∆y i, j is the average lateral width of T i, j .

Step 1: Solution of Upper Layer Model without Exchange Terms
We now propose a method to evolve the solution for the model, (45). Let us emphasis that any appropriate method for a 2D shallow water model can be used here since (45) can be interpreted as a standard 2D shallow water equation with topography, z b ( X) + h 1 ( X, t) referred to as the apparent topography in [6,7]. In our simulation we use an apparent topography hydrostatic reconstruction scheme [1,6,7]: Let T j be a 2D cell in a 2D channel mesh and T k be its neighbour. Let e jk be the edge between T j and T k , and n jk be the unit vector normal outwards to T j . Furthermore, let |T j | and |e jk | be the area of T j and length of e jk respectively and let E j be the set of all edges of T j . Let h n 1, j , β j , z b, j and Π n 2, j = (h n 2,j , q n 2, j ) T be the lower layer water depth, channel depth, bottom elevation and upper layer cell average vector in 2D channel cell, T j ; while h n 1,k , β k , z b,k and Π n 2,k = (h n 2,k , q n 2,k ) T are those in 2D cell, T k . Note that the neighbour cell, T k could actually be in the floodplain; the only requirement is that the current cell, in which we compute the upper layer solution, must be in the channel.
Define the apparent topography: η n 1,p = h n 1,p + z b,p , for p = j, k.

Step 2: Complete Solution of the Full 1D Model
We now solve the full 1D model, (49). This involves two sub-steps, namely

Step 2.1: Black-Box Stage
Here, we use any available 1D solver to solve the full 1D model without the lateral flux term, Φ, namely with the data, (A, Q) n i . In this paper, we adopt the scheme from [15,16], see also [18,19]. This gives the approximation, (A n+1 i , Q n+1 i ).

3.4.2
Step 2.1: Add the floodplain coupling term The discrete coupling term, Φ i is the vector consisting of the first two components of the 2D lateral fluxes already computed by the upper layer solver in section 3.3. Therefore, the final update value for the full 1D model, (49) is then given by where ∆t is the time step.

Step 3: Final Update of the Upper Layer Model
With the intermediate solution, (h 2 , q 2 ) n+1 * i, j known, we now completely solve the upper layer model including the exchange term, (50). The approximate solution of (50) is the approximate solution of the system with the initial data, (h 2 , q 2 ) n+1 * i,j , where Then using forward Euler time discretization, the approximate solution of (71) is To complete the description of the scheme, we need to define the exchange terms and interface velocities S i,j and u η 1 ,i, j . But there are no equations for these terms. However they can be determined by requiring that the following conditions be satisfied: These conditions allow us to first obtain the heights h n+1 2,i, j , then the exchange terms, and finally the interface velocities u η 1 ,i, j are calculated so that we can compute q n+1 2,i, j .

Step 3.1: Approximating the Lower Layer Flow By (33) we have
given A n+1 i computed above. Note that we can not directly apply A n+1 1,i = A n+1 i − j h n+1 2,i j ∆y since we do not yet know h n+1 2,i j . However, by mass/momentum conservation, the following must hold Hence the intermediate lower layer wetted area is given by Using this lower layer update, A n 1,i and the intermediate solutions, A n+1 * 1,i and h n+1 * 2,i j we can now compute the upper layer heights in the next step.

Step 3.2: Upper Layer Heights Using the definition of A n+1
1,i from (73) we have to distinguish two cases A n+1 and consider two further cases.
(c) Remove area of reduced height from total gap area: Equations (76), (78) and (79) compute h n+1 2,i,j in T i j ⊂ K i for all cells.

Step 3.3: Upper Layer Discharge and Lower Layer Discharge
Having computed h n+1 2,i, j , we compute the exchange terms, S i, j using the first equation in (72), hence and compute the interface velocity, u η 1 ,i, j following [3] namely Thus the interface velocity is equal to the upper layer velocity if mass is flowing from the upper to the lower layer and otherwise equal to the lower layer velocity.
Using the above definitions of the exchange terms and interface velocity, we define the upper layer discharge q Implicitly we can also compute the lower layer solution (A n+1 1,i , Q n+1 1,i ) by using (73) and local momentum conservation which leads to the lower layer discharge given by 3.6 Obtaining the full 2D data H n+1 i, j , q n+1 x,i, j and q n+1 y,i, j for T i, j ⊂ K i This final step is only required for postprocessing and for studying the properties of the scheme in the following section.
The full water depth, H n+1 i,j is given by where The full 2D x-discharge is computed as Finally, the y-discharge is computed as follows : where This completes the numerical algorithm for the channel flow, for the floodplain flow we use a hydrostatic reconstruction scheme as described for the upper layer model, see [18] for details. At the channel/floodplain interface, the lateral fluxes are effortlessly computed since in the 2D channel cell we have complete 2D value (H i, j , q i, j ) n at least if the channel is not full. In this important case, we have an approximation of the vertical discharge which is not the case for most existing methods where only a 1D model is used in the channel which does not provide values for the lateral discharge. If the channel is not full but the surrounding floodplain is not dry, the channel 1D flow data (A n i , Q n i ) is used to compute the coupling fluxes.

Properties of the Scheme
We now consider some important properties of the scheme proposed above.

Definition 4.1 (Consistency of Distribution Operation)
We require a distribution operation such as the one in section 3.2.1 to be consistent in the sense that whenever the channel is not full (H i, j ≤ β i j for all T ij ⊂ K i ), then the following conditions must hold: h 2i, j , q 2x,i, j and q 2y,i, j each equal zero ∀T i j ⊂ K i ;.
This ensure that we do not solve a two-layer problem when the channel is not full.

Definition 4.2 (No-Numerical Flooding Properties)
In order to satisfy the no-numerical flooding property, the full 2D data (H i,j , q xi, j , q yi, j ) n must satisfy the following condition for all n: i Either the lower layer is full or the upper layer is empty, i.e., ii If h n 2,i,j = 0, then (a) q n 2x,i, (which is constant laterally) and (c) q n y,i,j = 0; iii If ∃ j * ∈ {1, 2, ..., N 1 y } such that η n i, j * < η β i , then η n i, j = η n i, j * ∀j.
Condition (i) means that the lower layer is either full (h n+1 1,i, j = β i, j ) or upper layer is dry h n+1 2,i, j = 0. In other words, there should not be gap between the two layers. A consequence of this property is also that no overflowing of the channel can occure unless it is full. Condition (ii) states that if the upper layer is dry, the upper layer velocities and discharges must vanish, and the full layer flow velocity, u must be laterally uniform. Condition (iii) means that if the channel is not full, then the free-surface must be flat, i.e., not vary in the lateral direction. Thus the discrete consistency requirement (definition 3.1) is satisfied at all time steps.
In the following we will show that if the solution at time t n satisfies the above no numerical flooding properties then this is also true at time t n+1 . In addition, we also prove that the scheme is well-balanced and conserves mass under suitable conditions on the intermediate solutions. Proof: We need to prove that the ascertion in definition 4.1 is true. Let the channel not be full, that is, H i, j ≤ β i j . Then, we have h 2,i, j = H i, j − min(H i, j , β i j ) = H i, j − H i, j = 0, and q 2x,i, j = q 2y,i, j = 0 (see (63)). Therefore, we have h n+1 2,i,j = 0 or h n+1 1,i, j = β i j in either case. Hence, ii Let h n+1 2,i, j = 0, then (a) (see (72)).
That is η n+1 i, j is independent of j, hence also equal to η n+1 i, j * . Therefore, Hence the free surface is laterally flat. Therefore, the discrete consistency requirement (definition 3.1) holds at t n+1 , and the scheme satisfies the no-numerical flooding property.
Secondly, we show that A n+1 2,i = A n 2,i using (74) and that A n+1 * 1,i = A n 1,i : Hence, A n+1 i = A n i . Thirdly, we use the above results to show that h n+1 2,i, j = h n 2,i, j .
Hence by (78) and (79), we have h n+1 Fourtly, we show that q n+1 2,i, j = q n 2,i, j and Q n+1 Hence, we have shown that (H, q) n+1 i, j = (H, q) n i, j , so the method is well-balanced. Next, we prove that the proposed schemes for the intermediate solutions are mass conservative and also well-balanced for lake at rest. Since the hydrostatic reconstruction method is mass conservative, [2,1] and the 2D model, (45) does not introduce any source term to the height equation, then the scheme (67) is mass conservative. We state and prove a theorem below to show that this scheme, like the standard hydrostatic reconstruction [1] method, preserves well-balance of lake at rest.

Theorem 4.4
The upper layer scheme, (67) is well-balanced with respect to lake at rest.
Proof: The proof follows the same lines given in [2,1], for details see [18]. Since the 2D solver, (67), is also well-balanced, so (h n+1 * 2,i, j , q n+1 * 2,i, j ) T = (h n 2,i, j , q n 2,i, j ) T . Therefore, the lower layer scheme, (75) becomes By (74), we have (90) ). This completes the theoretical aspect of the VCM. In the next section, we present some numerical experiments to evaluate its performance compared to other coupling methods.

Numerical Results
We consider two simple test cases to verify the performance of the proposed method. We use the full 2D simulation results as the reference solution, and compare the results of the VCM with those of the Horizontal Coupling Method (HCM) [19] and the Flux-Based Method (FBM) [5].
Recall that in addition to the usual data required for the simulation, e..g, domain size and intitial conditions, the VCM method also requires a choice of the function z w b (x) used to determin when a channel is considered full and flodding might occure. To apply the VCM to this problem, we consider the following z w b (x): The free surface elevation at the last time step is shown in figure 5 for the full 2D, VCM and HCM. The figure shows that the coupling methods, VCM and HCM, approximates the full flow field with good accuracy, however, one can also see that the VCM computes better approximation than the HCM. Especially the 2D flow structure in the right part of the channel where the upper layer is active, is correctly captured by the VCM, unlike the HCM. The accuracy of the vertical coupling method is further illustrated in figure 6 which displays the time evolution of the free surface elevation, at the selected probe points. It can be seen that the VCM captures the full 2D results better than both the HCM and FBM, at all the probe points. Again, from figure 5, one can see that the VCM recovers 2D flow structure within the channel unlike the HCM and FBM. For the this test case, the vertical coupling method results in about 48% gain in computational time over the full 2D simulation. The fastest method (FMB) leads to a gain of about 55%. Hence in this test case, the VCM truly has shown good accuracy improving on the simple flux based coupling method while retaining most of the gain in efficiency.
To apply the vertical coupling method to this problem, we consider the following smoother c by the authors, 1677 version of the wall elevation: elsewhere. (94)

Result of test 2:
For all four methods we simulated this problem with 55 × 90 grid cells in the floodplain, while the grid resolutions in the channel are exactly the same as in the first test case. The simulation was run for ten seconds and CFL of 0.95.
The free surface elevation and velocity magnitude are shown in figures 7 and 8 respectively. We can see that that the VCM computes a better approximation of the full 2D solution than the HCM which in turn, is more accurate than the FBM. Moreover, the VCM reproduces a non laterally constant free surface elevation and velocity within the channel, which can not be achieved by either HCM or FBM. To further understand the results, the free surface elevation, the x-component and y-component velocity are plotted for selected probe points in figures 9. It can be seen that the vertical coupling method is more accurate than the other methods for all three flow quantities at all probe points and almost throughout the duration of simulation. Again, the vertical coupling method really captures the flow structure of the full 2D simulation.
Finally, figures 5 and 8 show that the VCM, unlike the other methods, recovers the 2D flow structure within the channel at the flooding regions. Again, the VCM continues to compute 1D solutions at non flooding regions; this demonstrates the self-adaptive nature of the method.
In this example VCM was about 75% slower then the computationally much simpler FBM method and about 35% more efficient then the full 2D simulation.
For this problem the VCM clearly improved the accuracy of the simulation both within the channel and in the floodplain compared to simpler coupling method but at a increased computational cost. But note that in this examples the percentage of the domain where a 1D assumption for the flow is valid is quite small so that it is not surprising that the gain in computational efficiency between full 2D and VCM is not so large. For a simulation of a very large network of rivers where the 2D region might be very small compared to the entire computational domain, the difference in efficiency of one coupling method over another will be far less significant, so that accuracy becomes the deciding factor. In this regard, the VCM is the best method of the three considered here.

Conclusion
In this paper we investigated a new method for coupling 2D and 1D shallow water models. We focused on the need to efficiently recover 2D channel flow structure during a flooding event to achieve high accuracy while maintaining the efficiency of the 1D channel model as much as possible. To this end we presented a vertical coupling approach which adds a second 2D layer inside the channel in regions where overflow can be expected to occur. This makes coupling the channel flow to the floodplain straightforward since 2D information is always available also within the channel when required. Any standard 1D channel flow and 2D floodplain solver can be used as building blocks for the new VCM method and only a slight modification of standard 2D solvers is required for the evolution of the second layer. We proved that the resulting method retains many properties of the 1D and 2D solvers used, e.g., mass conservation and well balancing. In addition we also studied a no-numerical flooding property. Our numerical results show that the VCM, in most cases, outperforms the other methods studied in this paper and does accurately recover 2D flow structures also within the channel when flooding occurs.
More detailed investigations with more complex channel geometries will be the focus of our further research. As mention in this paper, VCM actually consists of a very large family of methods. It depends on the choice of z w b (x) which is used to determine when a channel is considered full and the 2D layer model is used. For example, taking z w b (x) = ∞ everywhere leads to a horizontal flux coupling method similar to the FBM method. On the other hand taking z w b (x) = 0 in regions in danger of flooding and z w b (x) very large away from these regions using a smooth transition, results in a frontal type coupling method where a standard 2D solver is used in the region where z w b (x) = 0 and on the boundary of this region the VCM will lead to a blending type approach between the 2D and the 1D regions. The width of the blending region will depend on how z w b (x) changes from 0 to ∞. Further tests are required to understand the influence of different choices of z w b (x) and are the focus of ongoing work.