A computationally efficient subgrid model for coupled surface and groundwater flows

Coastal Engineering 157 (2020) 103665

Contents lists available at ScienceDirect

Coastal Engineering journal homepage: www.elsevier.com/locate/coastaleng

A computationally efficient subgrid model for coupled surface and groundwater flows Yujie Chen a,b ,∗, Fengyan Shi b , James T. Kirby b , Guoxiang Wu a , Bingchen Liang a a b

College of Engineering, Ocean University of China, Qingdao 266100, China Center for Applied Coastal Research, University of Delaware, Newark, DE 19716, USA

ARTICLE

INFO

Keywords: Subgrid model Groundwater flow Surface water and groundwater interaction

ABSTRACT We present the development of a high-efficiency and high-resolution subgrid model for simulating surface water and groundwater interaction in coastal regions. The governing equations are derived based on an assumption of the fully-saturated subsurface flow and a phase-averaging method, in which two phase functions are defined respectively for the surface phase and subsurface phase at the subgrid level. Both the surface water and groundwater processes are modeled by solving the phase-averaged equations on a coarse grid, taking into account the subgrid effects. Model validation with laboratory data indicates that the subgrid model is capable of simulating flow processes within permeable structures. The model was applied to simulating flooding and draining processes and the interaction of surface flows and groundwater flows in a meso-tidal salt marsh. The model results demonstrate that the subgrid model can be used in a complex salt marsh system with numerous narrow channels, creeks and artificial ditches, with accuracy similar to that of the full model running directly on the high-resolution grid, but computational costs can be two orders of magnitude smaller than for the full grid model. Tidal influence on the phreatic surface of the groundwater in the marsh system was discussed.

1. Introduction Recently, increasing attention has been paid to the study of interactions between groundwater and surface water in the research fields of hydrology, water resources management, and watershed ecosystems (Sophocleous, 2002; Danielopol et al., 2003; Liang et al., 2007; Yuan et al., 2011; Kong et al., 2010; Casulli, 2017). The importance of integrating surface water and groundwater into a single system has been enabled by different model coupling approaches (Sophocleous, 2002; Winter, 1998). In the modeling community, the surface water flow and groundwater flow are traditionally treated as two separate processes, partly due to the different timescales of surface/groundwater movements, and partly due to the difficulties in measuring and modeling their interactions (Liang et al., 2007). The inherent deficiencies of the separation have been recognized by previous studies (e.g., Freeze and Harlan, 1969; Freeze, 1972; Yen and Riggins, 1991) , in modeling hydrological processes in lakes, streams and wetlands, where the surface and subsurface water are hydraulically interconnected. In the recent study of large-scale watershed floods, it has been found that the water motion in the near-surface soil layer (i.e., 0–5 m) can have a similar timescale to that of the surface water flow (Liang et al., 2007; VanderKwaak and Loague, 2001; Gunduz and Aral, 2005). Salt marshes are also characterized by strong dynamic interactions between

surface water and groundwater, which form the basis for wetland functionality. Xin et al. (2011) pointed out that the cycles of flooding and draining processes induced by the tide on the marsh platform lead to dynamic, complex pore water flows, which are an important driving mechanism to export nutrients to the coastal sea. There are various types of numerical models for simulating surface flow and subsurface flow interactions. Integrated surface flow and subsurface flow models can be developed based on different governing equations with various levels of complexity. The most sophisticated models solve the 3-D Reynolds-averaged hydrostatic Navier–Stokes equations for surface water flow and the 3-D Richards’ equation for subsurface flow, such as in Spanoudaki et al. (2009), Yuan et al. (2011), and others. The two sets of governing equations were coupled through the kinematic boundary condition at the interface between the surface water and subsurface domains, either in an iterative or noniterative fashion. However, this 3-D approach may not be applicable to real-world simulations due to excessive computational demand. Compared to the 3-D approach, more computationally efficient models have been developed based on simplified governing equations for the surface and subsurface flows. For surface flows, simplified diffusion approximation equation or depth-integrated shallow water equations can be used. The diffusion wave approximation-based model

∗ Corresponding author at: College of Engineering, Ocean University of China, Qingdao 266100, China. E-mail address: [email protected] (Y. Chen).

https://doi.org/10.1016/j.coastaleng.2020.103665 Received 11 September 2019; Received in revised form 3 December 2019; Accepted 13 February 2020 Available online 17 February 2020 0378-3839/© 2020 Elsevier B.V. All rights reserved.

Coastal Engineering 157 (2020) 103665

Y. Chen et al.

The phase-averaging method introduced by Defina (2000) is generalized so that it can deal with two or multiple model components, such as surface flow, subsurface flow, as well as impermeable structures, in the subgrid scale. Noting that the initial idea of Defina’s subgrid model was to deal with the small-scale irregular bottom topography, we propose a similar subgrid concept that can treat both irregular surface water bottom and subsurface geometric features in small-scales. To obtain a theoretical formulation for the subsurface flows which is consistent with Defina’s equations, we re-derive the subgrid model equations which govern both surface and subsurface flows in a single system. In Section 2, we extend the phase-averaging method to include two phase functions for surface and groundwater flows. The numerical scheme is introduced in Section 3. In Section 4, a numerical test of a one-dimensional horizontal flow is reproduced to validate the subgrid model, and the convergence rates of our model are discussed. The model is also used to simulate flooding/draining processes and the interaction between surface and groundwater flows in a small meso-tidal marsh on Delaware Bay. Conclusions are presented in Section 5.

is the most efficient model among all types of surface flow models because the diffusion equation does not have the gravity wave characteristics which limit the numerical time stepping. It has been widely applied to modeling overland flows in inland watersheds (e.g., Panday and Huyakorn, 2004; Gunduz and Aral, 2005). However, the oversimplified diffusion approximation may not be applied to coastal areas, which are typically subjected to strong dynamic oceanic forces such as tides and waves (Kong et al., 2010). In coastal areas such as salt marshes, the surface water models based on the depth-integrated shallow water equations (or the Saint Venant equations in 1-D) are usually used (e.g., Liang et al., 2007). For subsurface flows, the 3-D Richards’ equation can be simplified to the 2-D Boussinesq equations according to the Dupuit–Forchheimer assumption, assumptions of saturated conditions, and saturation-independent hydraulic conductivity. The 2-D approximation can reduce the computational cost greatly and thus has been widely used in coupled surface water and subsurface flow models (e.g., Liang et al., 2007; Yuan et al., 2008; Casulli and Zanolli, 2010; Casulli, 2015, 2017). Modeling unsaturated subsurface flow is another aspect of simulating surface and subsurface flow interactions. When the unsaturated subsurface flow is considered, the 3-D Richards’ equation should be solved, resulting in a significant increase in computational cost. Recently, Casulli (2017) proposed an interesting strategy for solving the coupled 3D surface water and variably saturated subsurface equations more efficiently. He integrated the continuity equation over the entire water column in terms of both the surface and subsurface vertical integral fluxes, resulting in an exact mass conservation equation with respect to the piezometric head. Numerically, the vertical integrated continuity equation forms a well-posed mildly nonlinear system where the piezometric head is the only unknown variable. The equation with a single dependent variable can be solved more efficiently than the traditional coupling approach. However, the 3-D Richards’ equation still cannot be avoided due to the 3-D nature of the saturation distribution. Therefore, for 2-D integrated surface flow and subsurface flow models, it is assumed that the underground zone below surface water is always fully saturated (e.g., Liang et al., 2007; Yuan et al., 2008; Casulli, 2015), and the water surface and velocity components are continuous across the surface and subsurface interface. The computational efficiency of a coupled surface flow and subsurface flow model is the key issue for real-world applications. As more and more LiDAR-derived high-resolution Digital Elevation Model (DEM) data become available, it is increasingly desirable to run a highresolution surface flow model in a large domain and for a long period of time to get the small-scale geometrical impact on surface flows and sedimentation trends. There is the same computational efficiency problem for a subsurface model component when it is coupled with a high-resolution surface water model. In the previous modeling studies, semi-implicit numerical methods are known to be effective approaches to stable solutions at a minimal computational cost (e.g., Casulli, 1990; Kong et al., 2010). The subgrid techniques, which incorporate smallscale geometrical details into a regular grid model, were also developed to increase computational efficiency. In particular, Casulli (2009) developed a subgrid-based algorithm which can deterministically represent the mass balance in wetting and drying regions in the subgrid scale. Stelling (2012) and Volp et al. (2013) proposed a subgrid model based on a finite-volume form of the shallow water equations which are discretized on a relatively coarse grid, but take into account highresolution effects in the coarse grid computation. Wu et al. (2016) developed a subgrid model which combines the concept of Defina’s phased-averaged subgrid model (Defina, 2000) and the idea of Volp et al. (2013) in treatment of the fine grid effects into the coarse grid calculations. In applications of the interaction between surface and subsurface flow, Casulli (2017) incorporated subgrid geometrical details into the coupled surface and subsurface flow model. In this study, we extend the phased-averaged subgrid surface flow model developed by Wu et al. (2016) to include the groundwater component for modeling the surface water and groundwater interaction.

2. Derivation of governing equations The aim of developing a phase-averaged subgrid model is to formulate a set of governing equations which can be discretized on a coarse grid but include the feedback of small-scale geometrical effects. In Defina (2000), a phase function, representing ‘dry’ or ‘wet’ area, was introduced. By averaging the wet and dry areas over a spatial elementary area, and depth-integrating over the 3-D shallow water equations, the 2-D shallow water equations can be obtained with the phase-averaged quantities. The resulting equations can deal with partially dry/wet area in the coarse grid scale with the wet fraction at the water surface derived from the phase-averaging. Compared to the regular 2D depth-integrated shallow water equations without phaseaveraging, the phase-averaged equations contain extra terms which represent the subgrid effects or the feedback from the subgrid-scale processes. In the mass conservation equation, the wet fraction value appears in front of the time derivative of surface elevation as a ‘porous’ effect induced by partially dry area within an elementary region. In the momentum equations, there are several new terms derived from the phase-averaging and depth-integrating processes. The new diffusion terms resulted from the phase-averaging and depth-integrating of the Reynolds stress terms can be re-parameterized using the so-called effective viscosity (Defina, 2000). The new advective terms contain the momentum correction factor analogous to the shear-dispersion correction in Svendsen and Putrevu (1994) who did the similar averaging process over the momentum equations. The most significant term in the new momentum equations is the bottom stress term which represents the phase-averaged bottom friction. In Wu et al. (2016), the bottom friction term was calculated using the integration of bottom stresses in the subgrid scale, assuming the momentum balance within a coarse grid cell is essentially based on the steady-state solution. The steady-state solution gives that the surface gradient force is balanced by the bottom friction, or the friction slope 𝑆 can be evaluated by 𝑆 = 𝐶𝑑 𝑢2𝑠 ∕𝑔ℎ𝑠 , where 𝐶𝑑 is the friction coefficient, 𝑔 is the gravity acceleration, 𝑢𝑠 and ℎ𝑠 are flow velocity and the water depth, respectively, at each subgrid cell. The bottom friction integrated in the subgrid scale is essentially the feedback from small-scale geometrical effects to the coarse grid model. Following the same procedure, we introduce both the surface flow layer and the subsurface flow layer in a computational domain. Fig. 1 shows a vertical (𝑥, 𝑧) transect of a representative elementary area 𝐴 in the Cartesian coordinate system 𝐱 = (𝑥, 𝑦, 𝑧). The surface water (shaded blue) lies in the interval −𝑠 ≤ 𝑧 ≤ 𝜂 when 𝜂 ≥ −𝑠, and is absent where 𝜂 ≤ −ℎ. The subsurface flow (shaded tan) is confined within −ℎ ≤ 𝑧 ≤ min(𝜂, −𝑠), where the elevation 𝑧 = −ℎ represents an impermeable substrate. Here, 𝜂 has two definitions, one is the water surface elevation in the surface flow region, the other is the phreatic surface elevation/piezometric head in the subsurface region. 2

Coastal Engineering 157 (2020) 103665

Y. Chen et al.

directions. By phase-averaging over the mass conservation equation (6), we get ⟨𝜑1 ∇ ⋅ 𝐮𝑠 ⟩ = ⟨∇ ⋅ (𝜑1 𝐮𝑠 )⟩ − ⟨𝐮𝑠 ⋅ ∇𝜑1 ⟩ = 0,

(7)

Using (5), (7) can be rewritten as 𝜕𝑊 𝑠 − ⟨𝐮𝑠 ⋅ ∇𝜑1 ⟩ = 0, (8) 𝜕𝑧 𝑠 𝑠 where 𝐔𝐻 and 𝑊 are the phase-averaged surface water velocities in the horizontal and vertical directions, respectively. ∇𝐻 is the horizontal gradient operator. Integrating (8) over the depth from 𝑧 = −∞ to 𝑧 = 𝜂 ) 𝜂 ( 𝜂 𝜕𝑊 𝑠 𝑑𝑧 − ⟨𝐮𝑠 ⋅ ∇𝜑1 ⟩𝑑𝑧 = 0. (9) ∇𝐻 ⋅ 𝐔𝑠𝐻 + ∫−∞ ∫−∞ 𝜕𝑧 ∇𝐻 ⋅ 𝐔𝑠𝐻 +

Note that ∇𝜑1 is zero except at the soil surface. ⟨𝐮𝐬 ⋅∇𝜑1 ⟩ represents the infiltration/exfiltration at the soil surface. According to Leibniz’s rule, (9) can be rewritten as Fig. 1. Section view of the representative elementary area (Blue area: surface water; Tan area: groundwater).

𝜂

∇𝐻 ⋅

∫−∞

𝜂

𝐔𝑠𝐻 𝑑𝑧 − 𝐔𝑠𝐻 |𝜂 ⋅ ∇𝐻 𝜂 + 𝑊 𝑠 |𝜂 −

∫−∞

⟨𝐮𝑠 ⋅ ∇𝜑1 ⟩𝑑𝑧 = 0.

(10)

The kinematic boundary condition at the water surface is Following Casulli (2017), 𝜂 is assumed to be continuous across the interface between surface flow and subsurface flow. This assumption makes it easy to calculate the pressure gradient at the interface, where Dupuit’s assumption can be used on the subsurface flow side. In contrast to Casulli (2017), we assume that the subsurface layer is saturated below the local phreatic surface, so that the 3-D model equations for subsurface flows can be simplified as 2-D equations as described below. According to Defina (2000), a further assumption is made that the phreatic surface elevation, 𝜂, is constant within the elementary area, meaning no horizontal variation of pressure, based on the hydrostatic assumption. We introduce two phase functions which can be used for the surface phase and the subsurface phase separately. For the surface phase, the phase function is defined by { 1 𝑧 > −𝑠 𝜑1 (𝐱) = (1) 0 else For the subsurface phase, the phase function is { 1 𝑧 ≤ −𝑠 𝜑2 (𝐱) = 0 else

𝜕𝜂 + 𝐮𝑠𝐻 |𝜂 ⋅ ∇𝐻 𝜂 − 𝑤𝑠 |𝜂 = 𝑆, (11) 𝜕𝑡 where 𝑆 is the source term at the surface, may represent the precipitation/vapor rate. Phase-averaging (11) and making use of (4) and (5) yields 𝜕𝜂 + 𝐔𝑠𝐻 |𝜂 ⋅ ∇𝐻 𝜂 − 𝑊 𝑠 |𝜂 = 𝜗1 |𝜂 𝑆, 𝜕𝑡 where 𝜗1 is the surface water fraction at any vertical level: 𝜗1 |𝜂

𝜗1 (𝐱) = ⟨𝜑1 (𝐱)⟩ =

1 𝜑 (𝐱)𝑓 (𝐱, 𝑡)𝑑𝐴 𝐴 ∬𝐴 𝑖

𝜕𝜂 + 𝐔𝑠𝐻 |𝜂 ⋅ ∇𝐻 𝜂 − 𝑊 𝑠 |𝜂 = 𝛩1 𝑆. 𝜕𝑡 Substituting (14) into (10) gives 𝛩1

𝛩1

(2)

𝜂 𝜂 𝜕𝜂 + ∇𝐻 ⋅ 𝐔𝑠𝐻 𝑑𝑧 − ⟨𝐮𝑠 ⋅ ∇𝜑1 ⟩𝑑𝑧 = 𝛩1 𝑆. ∫ ∫ 𝜕𝑡 −∞ −∞

(3)

(5)

Note that ∇𝜑2 is zero except at the soil surface and the impervious bed. At the impervious bed, ∇𝜑2 is normal to the flow direction, so (∇𝜑2 )𝑧=−ℎ = 0. As the same as (9), ⟨𝐮𝐠 ⋅ ∇𝜑2 ⟩ represents the infiltration/exfiltration at the soil surface. According to Leibniz’s rule, (17) can be rewritten as

2.1.1. Surface flow The mass conservation equation in a 3-D domain can be expressed by

𝜂

∇𝐻 ⋅

∇ ⋅ 𝐮 = 0, 𝐮𝑠

(15)

𝜕𝜃 + ∇ ⋅ 𝐮𝑔 = 0, (16) 𝜕𝑡 where the superscript ()𝑔 represents the groundwater component, hereafter. 𝜃 = 𝜖𝐸 is the moisture content, in which 𝜖 is the porosity of soil and 𝐸 is the water saturation. For the case that the subsurface flow is saturated, the first term in (16) vanishes. By phase-averaging using (2) and depth-integrating over the depth from 𝑧 = −∞ to 𝑧 = 𝜂 ) 𝜂 ( 𝜂 𝜕𝑊 𝑔 𝑑𝑧 − ⟨𝐮𝑔 ⋅ ∇𝜑2 ⟩𝑑𝑧 = 0. (17) ∇𝐻 ⋅ 𝐔𝑔𝐻 + ∫−∞ ∫−∞ 𝜕𝑧

2.1. Mass conservation

𝑠

(14)

2.1.2. Subsurface flow Based on the Richards’ (1931) equation for porous media, the continuity equation governing the subsurface flow can be written as

where 𝐴 is the representative elementary area, and 𝑖 = (1, 2). The averaging procedure is assumed to satisfy both Leibniz’ rule and Gauss’ rule, that is ⟨ ⟩ 𝜕𝜑𝑖 (𝐱)𝑓 (𝐱, 𝑡) 𝜕 (4) = ⟨𝜑𝑖 (𝐱)𝑓 (𝐱, 𝑡)⟩, 𝜕𝑡 𝜕𝑡 ⟨∇𝜑𝑖 (𝐱)𝑓 (𝐱, 𝑡)⟩ = ∇⟨𝜑𝑖 (𝐱)𝑓 (𝐱, 𝑡)⟩.

(13)

Note that we do not use the terminology ‘porosity’ which was used in Wu et al. (2016) in order to avoid the soil porosity defined in the following section. Following Wu et al. (2016), we denote the surface water fraction at the water surface as 𝛩1 = 𝜗1 |𝜂 . Eq. (12) is then rewritten as

The phase-averaging process is conducted in each flow field. For a generic variable 𝑓 (𝐱, 𝑡), the phase-averaging can be expressed as ⟨𝜑𝑖 (𝐱)𝑓 (𝐱, 𝐭)⟩ =

1 𝜑 (𝐱)𝑑𝐴. 𝐴 ∬𝐴 1

(12)

(6)

∫−∞

𝐔𝑔𝐻 𝑑𝑧 − 𝐔𝑔𝐻 |𝜂 ⋅ ∇𝐻 𝜂 + 𝑊 𝑔 |𝜂

𝜂

−

()𝑠

where is the velocity vector of the surface flow. The superscript represents the surface flow component, hereafter. For convenience, we express the velocity vector as 𝐮𝑠 = 𝐮𝑠𝐻 +𝑤𝑠 𝐤, where 𝐮𝑠𝐻 is the horizontal velocity vector 𝐮𝑠𝐻 = 𝑢𝑠𝑥 𝐢 + 𝑢𝑠𝑦 𝐣, and ( 𝐢, 𝐣, 𝐤) are unit vectors in (𝑥, 𝑦, 𝑧)

∫−∞

⟨𝐮𝑔 ⋅ ∇𝜑2 ⟩𝑑𝑧 = 0.

(18)

The kinematic boundary condition at 𝑧 = 𝜂 requires 𝜖𝑠 3

𝜕𝜂 + 𝐮𝑔𝐻 |𝜂 ⋅ ∇𝐻 𝜂 − 𝑤𝑔 |𝜂 = 𝑆 𝜕𝑡

(19)

Coastal Engineering 157 (2020) 103665

Y. Chen et al.

depth-integrating used in Defina (2000), the phase-averaged momentum conservation equations can be written as ( ) ( ) 𝑃 𝑠 𝑄𝑠 𝜕𝑃 𝑠 𝜕 𝑃 𝑠𝑃 𝑠 𝜕 𝜖𝑥𝑦 + 𝜖𝑥𝑥 + 𝜕𝑡 𝜕𝑥 𝑌 𝜕𝑦 𝑌 𝜕𝜂 𝜕𝑅𝑥𝑥 𝜕𝑅𝑥𝑦 𝜏𝑠𝑥 𝜏𝑏𝑥 + 𝑔𝑌 − − − + =0 (30) 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜌 𝜌 ( ) ( ) 𝜕𝑄𝑠 𝑃 𝑠 𝑄𝑠 𝑄𝑠 𝑄𝑠 𝜕 𝜕 + 𝜖𝑥𝑦 + 𝜖𝑦𝑦 𝜕𝑡 𝜕𝑥 𝑌 𝜕𝑦 𝑌 𝜕𝜂 𝜕𝑅𝑥𝑦 𝜕𝑅𝑦𝑦 𝜏𝑠𝑦 𝜏𝑏𝑦 + 𝑔𝑌 − − − + =0 (31) 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜌 𝜌

where 𝜖𝑠 = 𝜖|𝜂 , the soil porosity at 𝑧 = 𝜂. 𝑆 is the surface source term representing the source/sink rate. Here, 𝑆 is the same term as used for the surface flow. The phase-averaged kinematic boundary condition at 𝑧 = 𝜂 can be written as 𝜕𝜂 + 𝐔𝑔𝐻 |𝜂 ⋅ ∇𝐻 𝜂 − 𝑊 𝑔 |𝜂 = 𝛩2 𝑆, (20) 𝜖 𝑠 𝛩2 𝜕𝑡 where 𝛩2 is the subsurface fraction at 𝑧 = 𝜂. Substituting (20) into (18) gives 𝜖 𝑠 𝛩2

𝜂 𝜂 𝜕𝜂 + ∇𝐻 ⋅ 𝐔𝑔 𝑑𝑧 − ⟨𝐮𝑔 ⋅ ∇𝜑2 ⟩𝑑𝑧 = 𝛩2 𝑆. ∫−∞ 𝐻 ∫−∞ 𝜕𝑡

(21)

where 𝑌 is the effective water depth defined similarly as in Defina (2000),

2.1.3. Integrated mass conservation equation The integrated mass conservation equation can be obtained by the sum of (15) and (21). As mentioned in the section above, ∇𝜑2 is zero except at the soil surface and impervious bed. At the impervious bed (𝐮𝑔 ⋅ ∇𝜑2 )|−ℎ = 0

𝜂

𝑌 =

(32)

(𝜏𝑠𝑥 , 𝜏𝑠𝑦 ) and (𝜏𝑏𝑥 , 𝜏𝑏𝑦 ) represent the phase-averaged wind stress on the water surface and the bottom friction, respectively. 𝑅𝑖𝑗 , where 𝑖, 𝑗 = 𝑥, 𝑦, represents shear stress tensor which can be calculated in terms of an effective viscosity 𝜈𝑒 ) ( 𝑠) ( 𝑠 ( 𝑠) 𝜕𝑄𝑠 𝜕𝑄 𝜕𝑃 𝜕𝑃 + , 𝑅𝑦𝑦 = 2𝜈𝑒 , 𝑅𝑥𝑥 = 2𝜈𝑒 , 𝑅𝑥𝑦 = 𝜈𝑒 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑦

(22)

∇𝜑1 is zero except at the soil surface. Using the flow continuity for the infiltration/exfiltration at the soil surface, i.e., (𝐮𝑠 ⋅ ∇𝜑1 + 𝐮𝑔 ⋅ ∇𝜑2 )|−𝑠 = 0,

1 𝜑 𝑑𝐴𝑑𝑧 𝐴 ∫−∞ ∬𝐴 1

(23)

(33) we get 𝜂

∫−∞

In (30) and (31), 𝜖𝑖𝑗 are the momentum correction factors which are from the depth-averaging of the phase-averaged advection terms. The correction factors cannot be computed directly unless a 3D velocity distribution is known. In Defina, 𝜖𝑖𝑗 was estimated to be slightly larger than 1. It was also pointed out that, since the convective inertia is usually negligible at small water depth or partially wet area, 𝜖𝑖𝑗 =1 can be conveniently assumed. The bottom friction 𝝉 𝑏 is estimated from the subgrid representations following Wu et al. (2016). Here, the general form of bottom friction proposed in Wu et al. (2016) is used. √ (𝜏𝑥 , 𝜏𝑦 ) = 𝜌𝐶𝑑 (𝑃 𝑠 )2 + (𝑄𝑠 )2 (𝑃 𝑠 , 𝑄𝑠 ), (34)

𝜂

⟨𝐮𝑠 ⋅ ∇𝜑1 ⟩𝑑𝑧 +

∫−∞

⟨𝐮𝑔 ⋅ ∇𝜑2 ⟩𝑑𝑧 = 0.

(24)

The sum of (15) and (21) gives (𝛩1 + 𝜖𝑠 𝛩2 )

𝜂 𝜕𝜂 + ∇𝐻 ⋅ (𝐔𝑠 + 𝐔𝑔𝐻 )𝑑𝑧 = 𝑆 ∫−∞ 𝐻 𝜕𝑡

(25)

In the derivation of (25), (𝛩1 + 𝛩2 ) = 1 is used to get the total 𝜂 precipitation/vapor 𝑆. ∫−∞ (𝐔𝑠𝐻 + 𝐔𝑔𝐻 )𝑑𝑧 is the total phase-averaged momentum flux contributed from both the surface and subsurface flows. We use (𝑃 𝑠 , 𝑄𝑠 ) and (𝑃 𝑔 , 𝑄𝑔 ) to denote the (𝑥, 𝑦) components of the surface water flux and the subsurface water flux, respectively. The mass conservation equation (25) can be rewritten as (𝛩1 + 𝜖𝑠 𝛩2 )

𝜕𝜂 𝜕(𝑃 𝑠 + 𝑃 𝑔 ) 𝜕(𝑄𝑠 + 𝑄𝑔 ) + + =𝑆 𝜕𝑡 𝜕𝑥 𝜕𝑦

where 𝜌 is the water density. 𝐶𝑑 is the so-called equivalent bottom friction coefficient (Volp et al., 2013) in the coarse grid and can be calculated based on quantities at the subgrid level:

(26) 𝐶𝑑 =

2.2. Momentum equation

in which √ √ 𝐶𝑑𝑠 ∕ℎ𝑠 𝛼𝑠 = ℎ ℎ𝑠 ∕𝐶𝑑𝑠 𝑑𝐴, ∬𝐴 𝑠 𝐴𝑌

2.2.1. Surface flow Based on the hydrostatic approximation, the momentum equation in the domain of the surface water can be expressed as 𝜕𝐮𝑠𝐻 𝜕𝑡

+ 𝐮𝑠 ⋅ ∇𝐮𝑠𝐻 = −𝑔∇𝐻 𝜂 + ∇ ⋅ 𝝉 𝐻

(27)

(36)

2.2.2. Subsurface flow The governing equation for the horizontal water velocities in the subsurface layer is based on Darcy’s theory (Darcy, 1856) and can be written as

(28)

The advection term (the second term in (28)) can be rewritten as

𝐮𝑔𝐻 = −𝜅∇𝐻 𝜂

⟨𝜑1 𝐮𝑠 ⋅ ∇𝐮𝑠𝐻 ⟩ = ⟨𝜑1 ∇ ⋅ (𝐮𝑠 𝐮𝑠𝐻 )⟩

(37)

where 𝜅 is the hydraulic conductivity. For a simple 2D subsurface model, 𝜅 can be treated as a depth-averaged value so that

= ∇ ⋅ ⟨𝜑1 𝐮𝑠 𝐮𝑠𝐻 ⟩ − ⟨∇𝜑1 ⋅ (𝐮𝑠 𝐮𝑠𝐻 )⟩ = ∇ ⋅ ⟨𝜑1 𝐮𝑠 𝐮𝑠𝐻 ⟩

(35)

ℎ𝑠 is the local water depth at subgrid point, and 𝐶𝑑𝑠 is the friction coefficient used at the subgrid points. Note that the expression of 𝛼𝑠 (36) is slightly different from that in Wu et al. (2016) due to the use of the flux form (𝑃 𝑠 , 𝑄𝑠 ), instead of the velocity form, in (34).

where 𝑔 is the gravitational acceleration. 𝝉 𝐻 = (𝝉 𝑥 , 𝝉 𝑦 ), the turbulent stresses in (𝑥, 𝑦) directions. Phase-averaging (27) gives 𝜕 ⟨𝜑 𝐮𝑠 ⟩ + ⟨𝜑1 𝐮𝑠 ⋅ ∇𝐮𝑠𝐻 ⟩ = −⟨𝜑1 𝑔∇𝐻 𝜂⟩ + ⟨𝜑1 ∇ ⋅ 𝝉 𝐻 ⟩ 𝜕𝑡 1 𝐻

𝐶𝑑𝑠 1 𝑑𝑧, 𝐴 ∬𝐴 𝛼𝑠2

(29)

In the derivation of (29), the incompressible flow condition is used. In addition, we assume ⟨∇𝜑1 ⋅ (𝐮𝑠 𝐮𝑠𝐻 )⟩ = 0, implying the momentum exchange between the surface layer and the subsurface layer is neglected. It is a reasonable assumption which has been adopted by many researchers such as Casulli (2015, 2017). With this assumption, the momentum equation is the same as Defina’s equation (Defina, 2000). Following exactly the procedure of the phase-averaging and

𝑃 𝑔 = −𝜅𝐷𝑔

𝜕𝜂 𝜕𝑥

(38)

𝑄𝑔 = −𝜅𝐷𝑔

𝜕𝜂 𝜕𝑦

(39)

where 𝐷𝑔 is the effective thickness of the subsurface layer defined similarly as the effective water depth, 𝜂

𝐷𝑔 = 4

1 𝜑 𝑑𝐴𝑑𝑧 𝐴 ∫−∞ ∬𝐴 2

(40)

Coastal Engineering 157 (2020) 103665

Y. Chen et al.

Fig. 2. Initial setup of the idealized 1-D flow model (Black solid line: initial water level (0.274 m); Blue trapezoid: sand embankment; Green dotted line: the waterlevel-collected location; Red dotted line: the velocity-collected location; Gray dotted line: the schematic diagram of tidal fluctuation).

Because there is no subgrid friction involved for the subsurface flow, the momentum equation is only solved at the coarse grid level.

The values of 𝛩1 , 𝛩2 , 𝑌 and 𝐶𝑑 are calculated at the subgrid level. The pre-storage method proposed in Wu et al. (2016) can be used to pre-calculate those values to save computational time. We assume that the surface is fully covered by either water or soil in the entire computational domain, i.e., 𝛩1 + 𝛩2 = 1. Therefore, there is no wetting and drying process because a soil cell is treated as a wet point with a porosity. The wetting–drying masks are not needed in the model.

3. Numerical method A mixed difference–differential parabolic equation with respect to the surface elevation can be derived following Wu et al. (2016). Eqs. (26), (30), (31), (38), and (39) are governing equations for the coupled surface and subsurface system and solved on the coarse grid. Integrated quantities such as 𝛩1 , 𝛩2 , 𝑌 and 𝐶𝑑 are calculated at the subgrid level. Eqs. (30) and (31) are first discretized in time from 𝑛 to 𝑛 + 1 to get (𝑃 𝑠 )𝑛+1

𝜕𝜂 = −𝐴𝑛 + 𝐵𝑛 𝜕𝑥

(41)

𝜕𝜂 + 𝐶𝑛 𝜕𝑦

(42)

(𝑄𝑠 )𝑛+1 = −𝐴𝑛 where ( 𝑛

𝐴 = 1 + 𝐶𝑑 (

√

𝑔𝑌 𝛥𝑡

4. Model tests In this section, the subgrid model with the coupled surface and subsurface flows was first validated against laboratory measurements and simulation results from Casulli’s (2015) model. The model was used to reproduce a one-dimensional horizontal flow and water exchange between an idealized tidal basin and an adjacent lagoon (Ebrahimi et al., 2007; Casulli, 2015). Then, the model was applied to a meso-tidal marsh, Brockonbridge Marsh in Delaware, U.S.A.

)𝑛 (43)

4.1. Experiment of surface/groundwater flows in a tidal lagoon

(𝑃 𝑠 )2 + (𝑄𝑠 )2 𝛥𝑡

)𝑛 𝛥𝑡 ⋅ √ 1 + 𝐶𝑑 (𝑃 𝑠 )2 + (𝑄𝑠 )2 𝛥𝑡 ( 𝑠 𝑠) ]𝑛 [ ( ) 𝜕𝑅𝑥𝑥 𝜕𝑅𝑥𝑦 𝜏𝑠𝑥 𝑃 𝑠 𝑃 𝑄 𝜕 𝜕 𝑃 𝑠𝑃 𝑠 − − + + + + 𝜕𝑥 𝑌 𝜕𝑦 𝑌 𝜕𝑥 𝜕𝑦 𝜌 𝛥𝑡 ( )𝑛 𝛥𝑡 𝐶𝑛 = ⋅ √ 1 + 𝐶𝑑 (𝑃 𝑠 )2 + (𝑄𝑠 )2 𝛥𝑡 [ ( 𝑠 𝑠) ( 𝑠 𝑠) ] 𝜕𝑅𝑥𝑦 𝜕𝑅𝑦𝑦 𝜏𝑠𝑦 𝑄𝑠 𝑛 𝑃 𝑄 𝑄𝑄 𝜕 𝜕 − − + + + + 𝜕𝑥 𝑌 𝜕𝑦 𝑌 𝜕𝑥 𝜕𝑦 𝜌 𝛥𝑡

The laboratory experiment of surface and groundwater flows and water exchange between an idealized tidal basin and an adjacent lagoon was conducted in the Hydraulics Laboratory, at Cardiff University (Ebrahimi et al., 2007). The experiment has become a standard benchmark test to verify the numerical models of surface water and groundwater interaction and flows in porous media. As shown in Fig. 2, a trapezoidal sandy embankment was built in the 6.48 m-long flume with the flat bottom. It can be expressed as 𝑠(𝑥) = 𝑚𝑖𝑛[0, 𝑚𝑎𝑥(−0.33, 21 |𝑥 − 1.40| − 0.35)], where 𝑠 is the height of the embankment. The model embankment was characterized by an average grain diameter of the sands 𝐷50 = 1 mm, a constant porosity 𝜖 = 0.3, and hydraulic conductivity 𝜅 = 0.01𝑚∕𝑠. The harmonic tides with constant amplitude of 0.06 m and period of 355 s were produced at the left end of the flume by a vertically oscillating weir. The initial water level on both sides of the embankment was initially 0.274 m, which was set as the High Water Level (HWL), while the Mean Water Level (MWL) was 0.214 m and the Low Water Level (LWL) is 0.154 m. Water elevations were collected at the middle points of the physical model cross-section on both sides of the embankment (Green dotted lines in Fig. 2) and for a period of 10 tidal cycles. Velocities were recorded over three tidal periods at the location denoted by the red dotted line, which is 0.12 m in font of the dike edge. The velocity data were collected at the depth of 0.128 mm under the MWL. The horizontal velocities which were time-averaged over 10 s were used to compare the model. Similar to Casulli (2015) who simulated the 1-D process of surface and groundwater flows using the model with a grid resolution of 0.02 m, we set up our model in a 1-D horizontal region of 𝑥 ∈ [−2.56, 2.56] as shown in Fig. 2. The model started with the initial water elevation of 0.274 m in the entire domain. The tidal elevation was specified

𝐵𝑛 =

(44)

(45)

𝜕𝜂 Note that the time level for 𝜕𝑥 is not specified in (41) and (42) and will be fixed in the following parabolic equation. An implicit scheme is used to treat the bottom friction terms in (34) as reflected in the √ denominator (1 + 𝐶𝑑 (𝑃 𝑠 )2 + (𝑄𝑠 )2 ) in (43)–(45). Substituting (41), (42), (38) and (39) into the mass conservation equation (26) yields a mixed difference–differential equation [ ] [ ] 𝜕𝜂 𝜕𝜂 𝜕𝜂 𝜕 𝜕 𝜕𝐵 𝑛 𝜕𝐶 𝑛 (𝛩1 + 𝜖𝑠 𝛩2 ) = (𝐴𝑛 + 𝜅𝐷𝑔 ) + (𝐴𝑛 + 𝜅𝐷𝑔 ) − − 𝜕𝑡 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕𝑥 𝜕𝑦

(46) Eq. (46) is the parabolic-type equation of surface elevation, integrating both the surface flow and subsurface flow contributions. The contribu𝜕𝜂 tion of the subsurface flow is reflected by the terms 𝜅𝐷𝑔 𝜕𝑥 and 𝜅𝐷𝑔 𝜕𝜂 , 𝜕𝑦 which are the flux components of the subsurface flow formulated in (38) and (39). The equation can be solved using the Odd-Even Hopscotch (OEH) numerical scheme, the same scheme used in Wu et al. (2016). After 𝜂 is obtained, the momentum fluxes can be calculated by (30) and (31) for the surface layer, and (38) and (39) for the subsurface layer. 5

Coastal Engineering 157 (2020) 103665

Y. Chen et al.

Fig. 5. Convergence rates with the coarse grid refinement. The subgrid ratio is 4.

Fig. 3 compares the modeled results with the laboratory data. The top panel shows that the modeled surface elevations are nearly identical to the data in the entire time series at the measurement location in the tidal basin, indicating that the model predicted the correct tidal oscillation in the basin using the elevation-nudging boundary condition. The middle panel shows the model/data comparison of surface elevation in the tidal lagoon region. The surface elevation initially drops from the still water level (𝑧=0.274 m), and then oscillates around the water level (𝑧 = 0.224) higher than MWL of the tidal basin. The water level fluctuation in the lagoon has more than 75% of amplitude reduction, and approximately 90◦ of phase lag comparing to those in the basin, consistent with the measured data as demonstrated in the figure. The model slightly over-predicted the troughs of the low tide. The bottom panel shows the model/data comparison of the current velocity at the measurement location denoted in Fig. 2 (red dotted line). The model predicted correct amplitude and phase of the current velocity, and the time series of modeled velocity look smoother compared to the measured data. To compare the results between the subgrid model and the full grid model running directly on the high-resolution grid, we carried out the same simulation using the full grid model with the grid resolution of 0.02 m. The results of surface elevation and velocity from the full grid model appeared to be identical to the subgrid model and were not shown here. Because there is no measured data of the phreatic surface available inside the embankment, we compared our results of the phreatic surface with the results from Casulli’s (2015) model. Fig. 4 shows the snapshots of free surface and phreatic surface at low tide (t=3375 s) and high tide (t=3550 s), respectively, from the present model (black solid lines) and Casulli’s model (red dashed lines). The two results look nearly identical, suggesting the present subgrid model performed similarly to other full grid model in this case. Beside the baseline case, more test cases with either different coarse grid sizes or subgrid ratios were conducted to test the grid convergence. Two sets of tests were performed, one is for testing the convergence of the coarse grid with a constant subgrid ratio, the other is for the convergence of the subgrid with a constant coarse grid size. For the former case, we adopted a sequence of different coarse grid spacing, 𝛥𝑥 × 2𝑝 , where 𝛥𝑥 = 0.01𝑚 and 𝑝 = 0, 1, 2, 3, 4. A constant subgrid ratio, 4, was used. Following (Shi et al., 2003), the convergence rates can be demonstrated by the overall root-mean-square (RMS) differences of modeled surface elevations at each grid point in the entire computational time and defined by √ ∑𝑛 ∑𝑚 2 𝑖=1 (𝜂𝑝 (𝑖, 𝑘) − 𝜂𝑠 (𝑖, 𝑘)) 𝑘=1 𝑅𝑀𝑆𝑝 = (47) 𝑚×𝑛

Fig. 3. Model/data comparisons of surface elevation at the measurement location in the tidal basin (top) and in the tidal lagoon (middle), and current velocity at the tidal basin (bottom).

Fig. 4. The snapshots of surface/phreatic surface elevation from the present subgrid model and Casulli’s (2015) model.

at the left boundary, while a vertical wall boundary condition was applied at the right boundary. The initial and boundary conditions can be expressed as: 𝜂(0, 𝑥) = 0.214, 𝜂(𝑡, 0) = 𝑎 cos( 2𝜋 𝑡) + 0.214, 𝑢(𝑡, 𝐿) = 0, 𝑇 where 𝑎 = 0.06 m, 𝑇 = 355 s. We set up a baseline case with the ratio of the subgrid size to the coarse grid size (subgrid ratio, hereafter) as 4, which means 0.08 m of grid resolution for the coarse grid, and 0.02 m of grid resolution for the subgrid. The results from the baseline case were compared with the laboratory measurement data.

where 𝜂𝑝 denotes the computed surface elevation in the case 𝑝, 𝑚 is the number of grid points (for 1-D) in the case of p = 4, 𝑛 is the number of time steps and the output time spacing is 10 s. 6

Coastal Engineering 157 (2020) 103665

Y. Chen et al.

For the latter set of the convergence tests, we kept the coarse grid size constant as 0.16 m and performed six cases with different subgrid ratios as 2, 4, 8, 16, 32, and 64, respectively. The RMS differences were calculated at the coarse grid using (47). Fig. 6 shows the RMS differences versus the subgrid ratios. The RMS difference decreases with the increasing of the subgrid ratio, indicating the convergence of the subgrid refinement. However, giving that the plot is in the semi-logarithmic coordinates (log for subgrid ratios), the convergence rate is much lower than that of the coarse grid refinement. The low convergence rate for the subgrid refinement suggest that the subgrid method may not affect the accuracy of the numerical schemes used in the coarse grid. 4.2. Application to Brockonbridge Marsh, Delaware Brockonbridge Marsh is a meso-tidal sault marsh located about 38 km upstream from the month of Delaware Bay (Fig. 7). A main channel connects the marsh system with Delaware Bay, and its width decreases with increasing distance inland, from about 25–30 m at the month to less than 10 m about 4 km inland. The marsh platforms are cut through by numerous second-order channels, creeks and the man-made ‘‘mosquito" ditches, which can be as narrow as 1 m, or less. Typical fronting berms at the channel edges, where the elevations are higher than the adjacent marsh platforms, serve as a levee system. The recent observation in Bombay Hook Marsh, which is close to Brockonbridge Marsh, revealed that the narrow cuts distributed either on the marsh platforms or around channels and ditches have significant effects on the flooding and draining processes on marsh platforms (Deb et al. 2019). A high resolution model which can resolve the small-scale feature of

Fig. 6. Convergence rates with subgrid refinement. The coarse grid size is 0.16 cm.

Fig. 5 shows the RMS differences versus the coarse grid spacing. A linear trend can be found in the logarithmic coordinates, indicating the nearly constant convergence rate for the refinement of the coarse grid with the constant subgrid ratio. It was expected because the subgrid only affects the bottom friction term by providing the equivalent bottom friction coefficient (eq. (35)) for the coarse grid without involving any grid discretization.

Fig. 7. The location (b) and topography of Brockonbridge Marsh, Delaware. Observation sites are marked with white dots, Sites A–F are located in the main channel, and Site G is in the tertiary channel. S1–S3 are dry land locations selected for demonstrating the groundwater processes. Detailed topography in the rectangle in (b) is shown in (c). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 7

Coastal Engineering 157 (2020) 103665

Y. Chen et al.

Fig. 8. Comparisons of modeled and measured water level time-series. The black solid lines represent the measurements, blue dot dashed lines represent the 2 m-resolution full grid model, and red dot dashed lines represent 8 m-resolution model with 2 m-subgrid model, and dark green dot dashed lines represent 8 m full grid model.

the groundwater interaction. Three locations, S1–S3, on uplands were selected adjacent to Site F as shown in Fig. 7(b) and the zoom-in view in Fig. 7(c). The model setup followed Wu et al. (2016) and vegetation effects were taken into account using the equivalent drag coefficient defined in Wu et al. (2017). Three cases were carried out with different model configurations, including the original full grid models with 2 m and 8 m grid resolutions, respectively, and the subgrid model with 8 m coarse grid and 2 m subgrid. The open water boundaries to the northeast of the domain are forced with water levels recorded at Site A. The parameters used in the groundwater model are relatively random due to the lack of measured data. We used a marsh soil porosity of 0.425, which was used in Gardner and Wilson (2006) and a saturated hydraulic conductivity of 10−3 based on the typical value for marsh deposit (Harvey and Odum, 1990). The initial phreatic surface elevations were set to the mean tidal level (MTL) in the entire domain though they are expected to be higher than MTL, considering inland hydraulic conditions. The influence of the initial condition is discussed later in this section. Fig. 8 shows the comparison of model results and the observed water level time-series at Site A–G. The subgrid model with 8m/2m subgrid configuration performed similarly to the 2 m full grid model, and both models were in good agreement with the data. In particular, the models accurately simulated the large flooding event during the storm on March 25–26, which persisted much longer with increasing distance into Brockonbridge Marsh, with the gauge signal at Site F taking several days to return to normal water levels. On the contrary, the 8 m full grid model failed to predict the draining processes after the storm due to the coarse grid which cannot resolve the narrow creeks

the marsh topography is necessary for modeling the hydrodynamic processes in the marsh system. To increase the computational efficiency, Wu et al. (2016, 2017) applied the subgrid model to simulation of the flooding and draining processes in Brockonbridge Marsh. Their model tests showed that up to two-order of magnitude speed-up can be achieved using the subgrid model with 8 m-coarse grid resolution and 2 m-subgrid resolution (subgrid ratio of 4) versus the full grid model running directly on the 2 m-grid. The accuracy of the subgrid model was similar to the high resolution full grid model. In this study, we carried out the similar tests using the coupled surface and subsurface subgrid model in the marsh system where the hydrodynamic data are available. Although there is no measured data for groundwater processes, the tests may demonstrate the model capability of simulating the interaction between the surface flows and groundwater flows as well as the computational efficiency of the subgrid approach. The high-resolution bathymetric/topographic data are based on LiDAR scanning on the platforms and field surveys inside the channels (Mieras et al., 2014) as shown in Fig. 7(b). A field campaign was conducted between 20 March and 4 April 2013 (Mieras et al., 2014; Pieterse et al., 2016), including a storm event from 25 to 27 March. Water-level records were available from pressure gauges at seven sites labeled as Sites A–G in Fig. 7(b), among which, Site A–F were deployed along the main channel, and Site G was located in a small tertiary channel. In Wu et al. (2016), a special attention was paid to the results at Site G because the water elevation was not modeled correctly using the 8-m resolution full grid model, while was successfully predicted by the subgrid model with 8 m coarse grid and 2 m subgrid. In the present study, we focused on the simulation of the surface water and 8

Coastal Engineering 157 (2020) 103665

Y. Chen et al.

Fig. 9. Surface elevation (left) and phreatic surface elevation at high tide (03/30/2013 12:00).

Fig. 10. Surface elevation (left) and phreatic surface elevation at low tide (03/30/2013 18:00).

and cuts and small-scale feature of the marsh platforms. These results are basically identical to the previous results from Wu et al. (2016, 2017) when the same model parameters were used, indicating that the groundwater processes do not alter the surface flow processes in the short time period. Figs. 9 and 10 illustrate the subgrid model results (8 m/2 m configuration) of surface elevation for the surface water and phreatic

surface elevation for the groundwater at the high tide (03/30/2013 12:00) and the low tide (03/30/2013 18:00), respectively. The flooding and draining processes are indicated by the surface elevation changes and wetting/drying on the marsh platforms between the high tide and low tide. Compared to the intra-tidal variation of the surface water, the phreatic surface elevations for the groundwater are relatively static between the high-tide and low-tide. The tidal influence on the 9

Coastal Engineering 157 (2020) 103665

Y. Chen et al.

still keep expanding slowly and could not get the equilibrium stage in the present simulation period. In addition, extra sources such as precipitation, evaporation, as well as river discharges may be important but were neglected in the simulation. To better present the influence of the tide on the phreatic surface, Fig. 11 shows the spatial distribution of the envelopes of the phreatic surface variation. The envelopes at the channel edges get the maximum, close to the range of water level fluctuation inside the channel. They decrease with increasing distances from the main channel. This behavior is similar to the tidal influence on the water table in beaches (e.g. Nielsen, 1990). To demonstrate the surface flow effects on the groundwater on dry lands, Fig. 12 shows the time series of the phreatic surface elevation at Sites S1–S3, which are about 100 m from the main channel and stay dry in the entire period of the simulation. The time series of surface elevation at the closest channel point, Site, F, are also illustrated in the figure for a comparison. The results show that the groundwater phreatic surfaces at the three dry points initially increase in the three-day spinup period and then fluctuate around a certain higher water level. The fluctuation magnitudes at the three sites are different, depending on their locations. The farther the site is from the main channel, the smaller the fluctuation is at the site. The phase lags of the phreatic surface fluctuation relative to the water level at Site F can also be observed. The phase lag at S3 is the largest among S1–S3 because it is the farthest site from the channel. The computational efficiency of the present subgrid model is consistent with the original subgrid model in Wu et al. (2016) and will not discussed in detail here. We observed that the groundwater calculation only accounts for about 4% of the total computational time, suggesting that a similar high speed-up rate to that in Wu et al. (2016) can be achieved for the subgrid model versus the full grid model. In the application of Brockonbridge Marsh, a speed-up rate of 36.80 was observed for the present subgrid model with the subgrid ratio of 4 and without the pre-storage method, similar to the speed-up (36.87) in Wu et al. (2016).

Fig. 11. Envelope of phreatic surface elevation obtained from the simulation in the period of 03/25–04/04, 2013.

groundwater gets maximum close to the channel edge and decreases with increasing distance from the channel. It should be mentioned that the phreatic surface distribution shown here may not be realistic due to unrealistic boundary and initial conditions of the groundwater. The time scale of a large groundwater system is much longer than the tidal scale, especially in the areas far from the channel system. For a short-term simulation, the initial distribution of water table may be crucial for modeling the evolution of the phreatic surface distribution in the large marsh system. Our tests showed that, with the initial phreatic surface condition set to MTL, the phreatic surface close to the channel can reach a certain high level in about three-day spin-up period. However, the phreatic surfaces

5. Conclusions In this study, we extended the phased-averaged subgrid surface flow model developed by Wu et al. (2016) to include the groundwater component for modeling the surface water and groundwater interaction. The subgrid model equations were derived based on the phase-averaging method initially introduced by Defina (2000). To include the groundwater component, two phase functions are defined respectively for the surface phase and subsurface phase. The resulting equations include the mass conservation equation integrating the surface flow and subsurface flow, and momentum equations governing

Fig. 12. The time series of the phreatic surface elevation at Sites F and S1–S3. 10

Coastal Engineering 157 (2020) 103665

Y. Chen et al.

References

the surface flow and subsurface flow, respectively. The momentum equations for the surface flow are consistent with the ones in Wu et al. (2016), and thus the same approaches to subgrid configuration and numerical method can be applied in the present subgrid model. For the groundwater flow, the simple Boussinesq equations (Darcy formulas) were used only on the coarse grid. The subgrid model was validated against the data from the laboratory experiment of surface and groundwater flows in a tidal lagoon. The model/data comparison indicated the subgrid model is capable of simulating flow processes involving permeable structures with the similar accuracy to the full grid model. The convergence tests suggested the nearly constant convergence rate for the refinement of the coarse grid and a constant subgrid ratio. The convergence rates for the subgrid refinement are slow. The model was applied to simulating the flooding/draining processes and the interaction of surface flows and groundwater flows in a meso-tidal salt marsh, Brockonbridge Marsh in Delaware. The results indicated that the surface elevations at the field measurement gauges were very well reproduced by the subgrid model. The groundwater processes do not alter the water surface level in the channel. Tidal influences on the phreatic surface of the groundwater are maximum at the channel edges and decrease with increasing distance from the channel. We pointed out the phreatic surface distribution from the present application may not be realistic due the unrealistic boundary and initial conditions. However, the model test in Brockonbridge Marsh demonstrated the potential use of the subgrid model in simulating the interaction of the surface water and groundwater, and the high model efficiency for modeling a large-scale marsh system with narrow channels, creeks and artificial ditches. For general coastal applications, the present subgrid model offers a promising approach to making use of high-resolution bathymetry measurements, like LiDAR-derived digital elevation model data, in simulating the flooding and draining processes on highly-complex topography, including the interaction between surface and groundwater flows. The phase-averaging method used in this study provides a general approach to deriving subgrid models with multiple model components, such as surface flow, groundwater flow, unsaturated subsurface flow, as well as impermeable structures.

Casulli, V., 1990. Semi-implicit finite difference methods for the two-dimensional shallow water equations. J. Comput. Phys. 86 (1), 56–74. Casulli, V., 2009. A high-resolution wetting and drying algorithm for free-surface hydrodynamics. Internat. J. Numer. Methods Fluids 60 (4), 391–408. Casulli, V., 2015. A conservative semi-implicit method for coupled surface–subsurface flows in regional scale. Internat. J. Numer. Methods Fluids 79 (4), 199–214. Casulli, V., 2017. A coupled surface-subsurface model for hydrostatic flows under saturated and variably saturated conditions. Internat. J. Numer. Methods Fluids 85 (8), 449–464. Casulli, V., Zanolli, P., 2010. A nested Newton-type algorithm for finite volume methods solving Richards’ equation in mixed form. SIAM J. Sci. Comput. 32 (4), 2255–2273. Danielopol, D.L., Griebler, C., Gunatilaka, A., Notenboom, J., 2003. Present state and future prospects for groundwater ecosystems. Environ. Conserv. 30 (2), 104–130. Darcy, H., 1856. Les fontaines publiques de la ville de Dijon: exposition et application.... Victor Dalmont. Defina, A., 2000. Two-dimensional shallow flow equations for partially dry areas. Water Resour. Res. 36 (11), 3251–3264. Ebrahimi, K., Falconer, R.A., Lin, B., 2007. Flow and solute fluxes in integrated wetland and coastal systems. Environ. Model. Softw. 22 (9), 1337–1348. Freeze, R., 1972. Role of subsurface flow in generating surface runoff: 1. Base flow contributions to channel flow. Water Resour. Res. 8 (3), 609–623. Freeze, R., Harlan, R., 1969. Blueprint for a physically-based, digitally-simulated hydrologic response model. J. Hydrol. 9 (3), 237–258. Gardner, L.R., Wilson, A.M., 2006. Comparison of four numerical models for simulating seepage from salt marsh sediments. Estuar. Coast. Shelf Sci. 69 (3–4), 427–437. Gunduz, O., Aral, M.M., 2005. River networks and groundwater flow: a simultaneous solution of a coupled system. J. Hydrol. 301 (1–4), 216–234. Harvey, J.W., Odum, W.E., 1990. The influence of tidal marshes on upland groundwater discharge to estuaries. Biogeochemistry 10 (3), 217–236. Kong, J., Xin, P., Song, Z.-y., Li, L., 2010. A new model for coupling surface and subsurface water flows: with an application to a lagoon. J. Hydrol. 390 (1–2), 116–120. Liang, D., Falconer, R.A., Lin, B., 2007. Coupling surface and subsurface flows in a depth averaged flood wave model. J. Hydrol. 337 (1–2), 147–158. Mieras, R., Shi, F., Kirby, J.T., 2014. A high-resolution numerical model investigation into the response of a channelized salt marsh to a storm surge event. Research report CACR-14-07, Center for Applied Coastal Research, University of Delaware. Nielsen, P., 1990. Tidal dynamics of the water table in beaches. Water Resour. Res. 26 (9), 2127–2134. Panday, S., Huyakorn, P.S., 2004. A fully coupled physically-based spatially-distributed model for evaluating surface/subsurface flow. Adv. water Resour. 27 (4), 361–382. Pieterse, C.M., de Jonge, R., Berendsen, R.L., 2016. The soil-borne supremacy. Trends Plant Sci. 21 (3), 171–173. Shi, F., Svendsen, I.A., Kirby, J.T., Smith, J.M., 2003. A curvilinear version of a quasi-3D nearshore circulation model. Coast. Eng. 49 (1–2), 99–124. Sophocleous, M., 2002. Interactions between groundwater and surface water: the state of the science. Hydrogeol. J. 10 (1), 52–67. Spanoudaki, K., Stamou, A.I., Nanou-Giannarou, A., 2009. Development and verification of a 3-D integrated surface water–groundwater model. J. Hydrol. 375 (3–4), 410–427. Stelling, G.S., 2012. Quadtree flood simulations with sub-grid digital elevation models. Proc. Inst. Civ. Eng. 165 (10), 567. Svendsen, I.A., Putrevu, U., 1994. Nearshore mixing and dispersion. Proc. R. Soc. Lond. Ser. A Math. Phys. Sci. 445 (1925), 561–576. VanderKwaak, J.E., Loague, K., 2001. Hydrologic-response simulations for the R-5 catchment with a comprehensive physics-based model. Water Resour. Res. 37 (4), 999–1013. Volp, N., Prooijen, B., Stelling, G., 2013. A finite volume approach for shallow water flow accounting for high-resolution bathymetry and roughness data. Water Resour. Res. 49 (7), 4126–4135. Winter, T.C., 1998. Ground Water and Surface Water: A Single Resource, vol. 1139. DIANE Publishing Inc.. Wu, G., Li, H., Liang, B., Shi, F., Kirby, J.T., Mieras, R., 2017. Subgrid modeling of salt marsh hydrodynamics with effects of vegetation and vegetation zonation. Earth Surface Process. Landforms 42 (12), 1755–1768. Wu, G., Shi, F., Kirby, J.T., Mieras, R., Liang, B., Li, H., Shi, J., 2016. A pre-storage, subgrid model for simulating flooding and draining processes in salt marshes. Coast. Eng. 108, 65–78. Xin, P., Yuan, L.-R., Li, L., Barry, D.A., 2011. Tidally driven multiscale pore water flow in a creek-marsh system. Water Resour. Res. 47 (7). Yen, B.C., Riggins, R., 1991. Time scales for surface-subsurface flow modeling. In: Irrigation and Drainage. ASCE, pp. 351–358. Yuan, D., Lin, B., Falconer, R., 2008. Simulating moving boundary using a linked groundwater and surface water flow model. J. Hydrol. 349 (3–4), 524–535. Yuan, L.-R., Xin, P., Kong, J., Li, L., Lockington, D., 2011. A coupled model for simulating surface water and groundwater interactions in coastal wetlands. Hydrol. Process. 25 (23), 3533–3546.

Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. CRediT authorship contribution statement Yujie Chen: Software, Validation, Formal analysis, Data curation, Writing - original draft, Visualization. Fengyan Shi: Conceptualization, Methodology, Data curation, Writing - review & editing, Funding acquisition. James T. Kirby: Writing - review & editing, Investigation, Project administration. Guoxiang Wu: Data curation, Writing - review & editing. Bingchen Liang: Data curation, Writing - review & editing. Acknowledgments This study was supported by Delaware Sea Grant, USA project (R/HRCC-2 and R/RCE-2). Yujie Chen was supported by the Chinese Scholarship Council and the University of Delaware during her Ph.D. study. Guoxiang Wu was supported by the National Natural Science Foundation of China (No. 51709243). This research was supported in part through the use of Information Technologies (IT) resources at the University of Delaware. 11