Revisiting Surface-Subsurface Exchange at Intertidal Zone with a Coupled 2D Hydrodynamic and 3D Variably-Saturated Groundwater Model

Abstract

A new high-performance numerical model (Frehg) is developed to simulate water flow in shallow coastal wetlands. Frehg solves the 2D depth-integrated, hydrostatic, Navier–Stokes equations (i.e., shallow-water equations) in the surface domain and the 3D variably-saturated Richards equation in the subsurface domain. The two domains are asynchronously coupled to model surface-subsurface exchange. The Frehg model is applied to evaluate model sensitivity to a variety of simplifications that are commonly adopted for shallow wetland models, especially the use of the diffusive wave approximation in place of the traditional Saint-Venant equations for surface flow. The results suggest that a dynamic model for momentum is preferred over diffusive wave model for shallow coastal wetlands and marshes because the latter fails to capture flow unsteadiness. Under the combined effects of evaporation and wetting/drying, using diffusive wave model leads to discrepancies in modeled surface-subsurface exchange flux in the intertidal zone where strong exchange processes occur. It indicates shallow wetland models should be built with (i) dynamic surface flow equations that capture the timing of inundation, (ii) complex topographic features that render accurate spatial extent of inundation, and (iii) variably-saturated subsurface flow solver that is capable of modeling moisture change in the subsurface due to evaporation and infiltration.

1. Introduction

Numerical modeling simplifications in hydrology, hydraulics and coastal engineering typically depend on (i) spatial/temporal scales of interest, (ii) relative importance of various hydrological processes, and (iii) data availability. For example, estuary and wetland models often consider only surface flow, thereby neglecting dynamics of surface-subsurface exchange (e.g., [1,2,3]). Coastal aquifer models are often built in the two-dimensional (2D) vertical (x-z) plane, neglecting the influence of complex topography in the x-y plane (e.g., [4,5,6,7]). Coupled surface-subsurface models in three dimensions (3D) often apply the diffusive wave approximation (DW) to surface flow instead of the shallow water equations (SWE) to reduce computational costs (e.g., [8,9]). The focus of the present work is to develop our understanding on the effects of using the DW approximation and other relevant simplifications for coupled surface-subsurface flow modeling in a coastal marsh.

Model simplifications with regard to surface and groundwater coupling are known to affect the accuracy of salinity predictions for coastal systems. For example, Li and Hodges [10] simulate surface water salinity transport in the Nueces Delta (TX, USA) with a hydrodynamic model that neglects groundwater coupling. The salinity is underestimated in shallow marsh regions, which is at least partially attributable to ignorance of porewater fluxes and the surface-subsurface exchange. Langevin et al. [11] report increased model-data agreement when surface-subsurface exchange is included in a model of the Everglades. Yu et al. [8] find that 3D topographical features, especially surface depressions, can affect groundwater salinization after ocean surge events. The challenge for coupled surface-groundwater models is their additional computational complexity, which has led a number of researchers to experiment with the DW approximation for surface water flow as a way to reduce overall computational costs.

In the present study, the Fine Resolution Environmental Hydrodynamic and Groundwater (Frehg) model is developed to simulate the overland flow, surface-subsurface exchange and the transport processes in shallow coastal wetlands. In the surface domain, Frehg solves the 2D, hydrostatic, depth-integrated Navier–Stokes equations (i.e., the SWE). For comparison purposes the DW option implemented in Frehg turns off the unsteady and advective inertia terms. In the subsurface domain, Frehg solves the 3D variably-saturated Richards equation with a predictor-corrector method, which is more robust than the commonly used Newton’s method [12]. The Frehg model is validated against benchmark problems from the literature. More complex exercises with the model include a synthetic sloping plane and real-world topography of the Trinity River Delta (Texas, USA). The latter exercises are used to investigate the sensitivity of model results to a variety of simplification techniques from the literature. The overall goals are to (i) provide insight into the simplifications used on different types of simulation problems, (ii) identify their effects on the simulation results, and (iii) provide validation of an open-source numerical code that can be used for hydrological research.

Although this work is motivated by the need to accurately capture salinity fluxes between surface and subsurface in coastal zones, herein we limit our focus to volume fluxes (i.e., neglecting forcing due to salinity-driven density gradients). This approach allows us to focus on the momentum consequences of the DW approximation alone. We argue that inadequacy of the DW approximation for homogenous density (as demonstrated herein) cannot be recovered by adding density-driven flow.

The remainder of this manuscript is arranged as follows: Section 2 provides insight into the meaning of the DW approximation. Section 3 describes the governing equations and numerical techniques used in the Frehg model. Section 4 introduces the test problems and the results. Section 5 discusses and summarizes the implications of the test results to wetland simulation with the coupled surface-subsurface models.

2. Background

To reduce computational complexity, the DW approximation neglects the inertia terms in conservation of momentum, which is reasonable for flow varying over large spatial and long time scales [13,14]. However, topography of coastal marshes has variability over a wide range of spatial scales [15], which makes it difficult to prove that large scales are dominant. Similarly, tidal frequencies or time scales for rapidly-moving weather fronts might not be “long enough” for valid DW scaling of system forcing. Interestingly, although the DW equations are clearly simpler than the SWE and less computationally complex, it is not necessarily true that they are always more computationally efficient. Cea and Blade [16] noted that the DW method could be slow at fine grid resolution because of more restrictive stability criteria. Similar findings are reported in Neal et al. [17].

The DW approximation is commonly used in modeling catchment rainfall runoff, and was recently compared to the SWE in Caviedes-Voullieme et al. [18]. Although the two methods produced similar hydrographs, the SWE provided more accurate depth and velocity predictions. Their results help illustrate the differences between the DW approximation used in catchment runoff versus a coastal marsh. For the former, the difference between the spatial gradient of the water surface and the landscape is a secondary portion of the flow forcing, which is unidirectional and primarily driven by the landscape slope. In contrast, the spatial gradients of the water surface (without signficant landscape slope) are the primary flow forcing in a coastal marsh driven by winds and tide.

The surface-water DW approximation in surface-subsurface coupled modeling is typically used when the focus is on the subsurface domain, which has longer residence times and its behavior over time scales of years or decades is of interest (e.g., [9,19]). Indeed, under such conditions an even simpler approach than DW has been used: representing the surface domain as a spatially-uniform tidal elevation (e.g., [6]). However, evidence shows that daily tidal spatial and temporal variation has influence on the timing and magnitudes of surface-subsurface exchange [4,20,21]. In particular, spatio-temporal surface water gradients across a marsh lead to oscillatory wetting and drying of the landscape, which affects both the area over which subsurface exchanges occur and the driving pressure gradients. Although the spatio-temporal inundation in a marsh is primarily driven by wind and tide, evaporation and rainfall events also play a role and add more complexity to the processes. For example, Geng and Boufadel [22] find that when evaporation exists, the groundwater salinity near the surface of the intertidal zone varies at tidal time scales. Their results imply that the combined effect of evaporation and frequent wetting/drying leads to stronger infiltration/exfiltration in pore waters. A further challenge is that surface-subsurface exchange is inherently 3D. Xin et al. [6] find that fluctuation of tidal creeks affect marsh salt budget at different time scales depending on the distance to the creeks. As such, simpler 2D simulations of vertical cross-sections cannot capture the flow and transport processes.

Validity of the DW approximation can be examined by scaling of conservation of momentum, presented previously by Ponce et al. [14] and Hunter et al. [13]. The focus in Hunter et al. [13] is the nonlinear advection terms, which they argue should be conditional on a Froude number ($Fr$) and a shallow-water Reynolds number ($Rs$):

$$Fr=\frac{U}{\sqrt{gH}}<1,Rs=\frac{H^{\frac{4}{3}}}{g{n}^{2}L}<1$$

(1)

where U is a representative flow velocity, H is a representative depth and L is a characteristic length scale. Here $Rs$ is the ratio between advective inertia and bottom friction and is presented above as the inverse of that proposed in Hunter et al. [13].

In contrast, Ponce et al. [14] focused on the validity of neglecting the unsteady term in momentum. Their approach used a ratio of dimensionless wave period to steady flow Froude number, with a final recommendation that could be written as

$$Rt=\frac{30L}{T\sqrt{gH}}<1$$

(2)

where T is the wave period of surface disturbance, the value 30 is an empirical limiter on wave celerity error, and other terms are as denoted above.

An interesting observation is that Equation (1) requires small $H/L$ whereas Equation (2) favors large $H/L$, which limits the range over which we can neglect both unsteady momentum and nonlinear advection. Furthermore, the water depth in a coastal marsh could vary by orders of magnitudes, ranging from <10${}^{-2}$ m at the wetting/drying front to ∼10 m at the tidal inlet. The spatial scales of topographical features also exhibit a wide range of variation [15]. Different types of upstream and downstream boundary conditions introduce additional factors that affect the applicability of the DW approximation [23]. As a consequence, validity of the DW approximation in real marsh-scale simulations cannot be determined simply from Equations (1) and (2).

The studies noted above highlight the complex interactions among topography, tide, wetting/drying, evaporation, variably-saturated subsurface flow, and surface-subsurface exchange. To the authors’ knowledge, numerical studies comparing these processes have not been presented in the literature. In particular, unsaturated subsurface flow is often modeled in 2D with simplified surface flow equations (e.g., [19,22]) while 3D surface-subsurface models often assume fully saturated groundwater flow (e.g., [11,24,25]). The reason might be the extensive computational cost and simulating challenges incurred by solving the variably-saturated Richards Equation [26].

For the present work, we seeks to understand the consequences of using the DW approximation in coupled surface-subsurface simulations. Rather than using simple scaling arguments, the validity of the DW approximation is examined by numerically calculating and comparing each term in the momentum equations [27]. This approach allows examination of force balancing at different times and different locations in the marsh. Model sensitivity to the inclusion/exclusion of evaporation and wind is also investigated. It will be shown that discrepancies of using DW model could be exacerbated when evaporation and wind forces are considered.

3. Methods

The Frehg model (described herein) is a descendant of the Fine Resolution Environmental Hydrodynamic model (Frehd) used in Li and Hodges [10], Hodges [28], Hutschenreuter et al. [29]. The Frehd model solves the full nonhydrostatic 3D Navier–Stokes equations with options for model simplifications to the 3D SWE and 2D depth-averaged SWE. In the present study, a 3D variably-saturated groundwater flow model is coded and coupled with the 2D hydrostatic, depth-integrated approach of Frehd to form the Frehg model. Transport of passive scalars is provided by an optional module within Frehg, but is not used or described herein.

3.1. Surface Flow Module

For the default operation of Frehg, surface flow is modeled with the 2D depth-integrated SWE for continuity:

$$\frac{\partial \eta }{\partial t}+\frac{\partial h{u}_j}{\partial x_j}=0$$

(3)

and momentum:

$$\frac{\partial u_j}{\partial t}+u_k\frac{\partial u_j}{\partial x_k}=-\left|\overrightarrow{g}\right|\frac{\partial \eta }{\partial x_j}+\frac{\partial }{\partial x_k}\left(\nu \frac{\partial u_j}{\partial x_k}\right)-\frac{|\overrightarrow{g}|n^2{u}_j\left|\overrightarrow{u}\right|}{2h^{4/3}}+\frac{{\tau }_j}{\rho h}$$

(4)

where $j,k\in \left\{1,2\right\}$, u is velocity, $\eta$ is surface elevation, h is depth, $\nu$ is eddy viscosity, $\tau$ is wind stress, $\rho$ is water density, $\overrightarrow{g}={[0,0,g]}^T$ represents gravity and n is the Manning’s roughness coefficient. Note that the Einstein summation convention is reduced to implied sums over $\left\{1,2\right\}$ for repeating subscripts. From left to right, Equation (4) consists of the unsteady term, the nonlinear advection term, the hydrostatic pressure term, the horizontal turbulent diffusion term, the bottom friction term and the wind term. It has been shown that horizontal turbulent momentum diffusion is relatively insignificant (compared to bottom friction) for shallow marsh models [10,30] and will be neglected herein, but is retained above for completeness and because it is an optional term in the Frehg code.

The wind stress vector, $\overrightarrow{\tau }$, is estimated as:

$$\overrightarrow{\tau }=C_w{\rho }_a{\left(u_w-\left|\overrightarrow{u}\right|cos\beta \right)}^2{\left[cos\omega ,sin\omega \right]}^{\mathrm{T}}$$

(5)

where, $C_w$ is the wind stress coefficient, ${\rho }_a$ is air density, $u_w$ is wind speed, $\overrightarrow{u}$ is the water velocity vector in the model x-y coordinate system, $\beta$ is the angle between wind direction and water flow direction, $\omega$ is the angle between wind direction and the “north” of the model’s grid system. From Equation (4), the wind term becomes dominant when depth h is small, which could cause numerical instability. In Frehg, this issue is handled by applying an exponential decay function to wind stress as h approaches zero [10]. This approach prevents model instability when a strong wind stress is applied over shallow depths. Note that wind-driven surface flow is commonly modeled in shallow coastal marshes (e.g., [31,32]), but wind effects are excluded in the DW approximation.

In Frehd and Frehg, Equations (3) and (4) are solved with a semi-implicit, hybrid finite difference/finite volume method that is similar to Casulli [33], Casulli and Cattani [34]. In short, Equation (4) is discretized such that time-stepping of velocity $u_j$ is a function of all other variables and parameters, with $|\overrightarrow{u}|$ handled explicitly and $\eta$ in an implicit form. The $u_j$ time-marching relationships are substituted into Equation (3), which creates a linear implicit system for the surface elevation $\eta$. The system is solved to obtain $\eta$ at the new time level and back-substituted into the discrete equations for time-marching $u_j$ to update the velocities. The surface flow solution in Frehg uses an Arakawa C grid, with solution of $u_1$ on the x face, $u_2$ on the y face, and $\eta$ at cell centers. A detailed description of the numerical schemes can be found in Li and Hodges [10] and citations therein.

3.2. Diffusive Wave Approximation

The diffusive wave approximation to the 2D SWE is made by neglecting the inertia terms (unsteady and nonlinear advection) and the wind term in Equation (4). Here we also neglect the insignificant momentum diffusion term so that DW momentum is a balance relationship between the driving water surface gradient and the bottom friction:

$$\frac{\partial \eta }{\partial x_j}=-\frac{n^2{u}_j\left|\overrightarrow{u}\right|}{2h^{4/3}}$$

(6)

which implies that spatio-temporal gradients in u are instantaneously compensated by the free-surface gradients affecting continuity, Equation (3).

To implement the DW approximation within the semi-implicit time-marching schema of Frehg is straightforward. The $u_j$ in Equation (3) is replaced with

$$u_j=-\left(\frac{2h^{4/3}}{n^2\left|\overrightarrow{u}\right|}\right)\frac{\partial \eta }{\partial x_j}$$

(7)

where $\eta$ is treated implicitly and the parenthetical factor is treated explicitly in the semi-implicit time-marching algorithm. The time-linearization implied by treating h as lagged coefficient (hence removing an implicit nonlinearity with $\eta$) is consistent with the original SWE semi-implicit scheme of Casulli and Cattani [34] that neglects a similar 2nd-order unsteady nonlinear coupling [35].

For completeness, it is useful to present another simplification used in the literature: the Kinematic Wave (KW) approximation. The KW adds to the approximations used in DW by assuming the surface water depth gradient, $\partial h/\partial x$ is negligible. By definition,

$$\frac{\partial \eta }{\partial x}\equiv \frac{\partial h}{\partial x}-S_0$$

(8)

where $S_0$ is the bottom slope. So the KW approximation reduces momentum to a simple relationship between velocity, depth, roughness and bottom slope:

$$\frac{n^2{u}_j\left|\overrightarrow{u}\right|}{2h^{4/3}}=S_0$$

(9)

Herein, results from prior KW studies in the literature are used in Section 4.2 as part of model validation.

3.3. Subsurface Flow Module

In Frehg, variably-saturated groundwater flow is modeled with the 3D Richards equation, which can be written in either the pressure head form:

$$\left(\frac{S_s\theta }{\varphi }+S_c\left(\psi \right)\right)\frac{\partial \psi }{\partial t}=\frac{\partial q_j}{\partial x_j}$$

(10)

or in the mixed form:

$$\frac{S_s\theta }{\varphi }\frac{\partial \psi }{\partial t}+\frac{\partial \theta }{\partial t}=\frac{\partial q_j}{\partial x_j}$$

(11)

where $j\in \left\{1,2,3\right\}$, $\psi$ is pressure head, $\theta$ is water content, $\varphi$ is porosity, $S_s$ is specific storage, $S_c=\partial \psi /\partial \theta$ is specific capacity and $q_j$ is the flux between grid cells. Note that $q_j$ can be modeled using Darcy’s law:

$$q_j=K\left(\psi \right)\left(\frac{\partial \psi }{\partial x_j}-g_j\right)$$

(12)

where $K\left(\psi \right)$ is the hydraulic conductivity and $g_j$ is the vector gravity.

Frehg uses the Mualem–van Genuchten model [36,37] to represent the nonlinear relationship between K and $\psi$, $\theta$ and $\psi$:

$$\theta ={\theta }_r+\left({\theta }_s-{\theta }_r\right)S_a$$

(13)

$$K=K_s{S}_a^{\frac{1}{2}}{\left[1-{\left(1-S_a^{\frac{1}{m}}\right)}^m\right]}^2$$

(14)

where $S_a={\left(1+{\left|\alpha \psi \right|}^{\widehat{n}}\right)}^{-m}$ is soil water saturation, $m=1-1/\widehat{n}$, $\alpha$ and $\widehat{n}$ are soil parameters, ${\theta }_s$ and ${\theta }_r$ are saturated and residual water contents, $K_s$ is the saturated hydraulic conductivity.

The solution of discrete versions of Equations (10)–(12) uses the predictor-corrector-allocation (PCA) scheme developed in Li et al. [12]. In short, the predictor step discretizes Equation (10) with a linearized implicit scheme (i.e., $\psi$ is treated implicitly but $\theta$, K and $S_c$ are all explicit), which results a solution for $\psi$ that is non-conservative in unsaturated zones due to linearization errors. In the corrector step, $\psi$ obtained from Equation (10) is substituted into Equation (11) to explicitly update water content, $\theta$. Finally, the solution of the predictor step is used in the saturated zone and the solution of the corrector step is used in the unsaturated zone. Excess moisture at the saturated-unsaturated interface (due to mismatch between the predictor and the corrector solutions) is redistributed to adjacent grid cells to achieve strict mass conservation. The PCA method is based on predictor-corrector methods of Kirkland et al. [38], Lai and Ogden [39]. The new method is conservative and has proven more robust than commonly-used Picard or Newton type iterative schemes.

Although Li et al. [12] only introduce a 1D PCA scheme, it is easily extended into 3 dimensions because computation of the Darcy fluxes Equation (12) is independent in each direction. The 3D moisture allocation is performed by distributing excess water based on the pressure gradient in each direction, which is similar to the 1D split distribution method described in Li et al. [12].

3.4. Surface-Subsurface Coupling

In Frehg, both surface and subsurface domains use rectangular grids of same sizes (i.e., uniform $\mathsf{\Delta }x$ and $\mathsf{\Delta }y$ throughout the two domains, but $\mathsf{\Delta }x$ might not equal $\mathsf{\Delta }y$). The surface and subsurface domains are coupled asynchronously, meaning that flow in each domain is solved independently with surface-subsurface exchange acting as boundary conditions. This coupling approach maximizes the flexibility of the Frehg model, and allows it to be reconfigured to model either surface or subsurface flow alone.

For coupled flow, the solution procedure of one complete time step is illustrated in Figure 1. The superscripts “n”, “n + 1” and “*” represent time level. The surface flow equations are solved first (with evaporation/rainfall included) to obtain a temporary solution of surface depth, $h^{*}$. The depth (if non-zero) is used as a Dirichlet pressure head boundary condition of Equation (10) for the subsurface domain. The exchange flux (positive upward) at the interface of the two domains is modeled as:

$$q_{\mathrm{exchange}}=K_{sz}\left(\frac{{\psi }_{\mathrm{top}}-h^{*}}{0.5\mathsf{\Delta }{z}_{\mathrm{top}}}-g\right)$$

(15)

where $K_{sz}$ is the saturated hydraulic conductivity in the vertical direction, ${\psi }_{\mathrm{top}}$ is the pressure head of the topmost grid cell of the subsurface domain, and $\mathsf{\Delta }{z}_{\mathrm{top}}$ is the vertical grid resolution of that cell.

For the surface domain, Frehg uses a simple depth-increment coupling, where the exchange flux is converted into an equivalent exchange depth and is added to the free surface elevation:

$$h_{\mathrm{exchange}}=q_{\mathrm{exchange}}\mathsf{\Delta }t\varphi A_z$$

(16)

where $\mathsf{\Delta }t$ is the size of the time step, $A_z$ is the area of the surface grid cell. We use area instead of the product of the grid side lengths (e.g., $\mathsf{\Delta }x\mathsf{\Delta }y$) because Frehg has the capability of modeling subgrid-scale topographic effects [40], in which $A_z\ne \mathsf{\Delta }x\mathsf{\Delta }y$. An alternative approach (used in the original Frehd for externally-specified high inflows) uses the exchange flux directly as a boundary condition on the $\eta$ linear solver of the semi-implicit method. For Frehg the simpler depth-increment approach is preferred as it is consistent with the approach to rainfall and evaporation (Section 3.5) and has proven stable in all our tested simulations. With the exchange depth added, the final surface elevation (${\eta }^{n+1}={\eta }^{*}+h_{\mathrm{exchange}}$) and depth ($h^{n+1}=h^{*}+h_{\mathrm{exchange}}$) are substituted into Equation (4) or (6) to get the surface velocities and complete the current time step.

3.5. Rainfall and Evaporation

In Frehg, rainfall is added to the surface domain by directly increasing the free surface elevation. If the surface cell is dry, rainfall first leads to surface ponding, then enters the subsurface domain according to Equation (15). Similarly, evaporation from the wet surface domain is modeled by directly reducing surface elevation based on the evaporation flux, $q_e$. For wet surface cells, $q_e$ is often computed as functions of net radiation and wind speed [41]. Since time-varying radiation and wind data are not used in the present study, the $q_e$ for wet surface is simply set as a constant. Detailed discussion on surface evaporation models is beyond the scope of this manuscript. Evaporation from the subsurface domain (i.e., direct losses from a dry cell) is estimated with the bulk aerodynamic model [42]:

$$q_e=\frac{{\rho }_a\left[{\alpha }_1{q}_{\mathrm{sat}}-q_a\right]}{\rho R_a}$$

(17)

where ${\rho }_a$ is air density, $q_{\mathrm{sat}}$ is the saturated specific humidity, which is a function of soil surface temperature, ${\alpha }_1=min\left(1,\frac{1.8\theta }{\theta +0.3}\right)$, $q_a$ is the air specific humidity, $R_a$ represents the aerodynamic resistance, which is a function of wind velocity. Note that the expressions of ${\alpha }_1$ and $R_a$ described in Geng and Boufadel [42] are based on measurements at specific field sites reported in Barton [43] and Liu et al. [44], respectively. Since no measured soil properties, meteorological and thermodynamic conditions are available in the present tests, we estimate ${\alpha }_1$, $q_{\mathrm{sat}}$ and $R_a$ following the data-fitted equations in Geng and Boufadel [42]. The meteorological and thermodynamic conditions are assumed constant. The only variable in Equation (17) that affects the dry-surface evaporation rate in our test simulations is $\theta$, the soil water content. A similar approach was used in Xin et al. [6]. Figure 2 shows the evaporation flux used in the present study with a surface temperature of 20 deg. C and a wind speed of 2 m/s. A low-saturation cutoff is applied in Frehg that enforces zero evaporation when water content reaches the pre-defined residual water content, ${\theta }_r$, of the Mualem-van Genuchten model, Equation (13).

3.6. Thin-Layer Treatment

In surface-water modeling of shallow marshes with wetting/drying, a minimum allowable depth ($h_{\min }$) is typically set to avoid numerical instabilities. Such problems can occur when the hydrostatic pressure gradient ($\partial \eta /\partial x$) applied to an infinitesimal layer causes unphysically high velocities [10]. The use of a simple cutoff $h_{\min }$ presents few difficulties when modeling wind and tidally-driven systems as the associated mass loss is typically smaller than the convergence residual of the linear solver for $\eta$. However, coupling with the subsurface and rainfall can be problematic for simple cutoff algorithms. For example, if $h_{\min }={10}^{-5}$ m (as used herein) with a model time step of $\mathsf{\Delta }t=60$ s, then a rainfall rate below 0.6 mm/h can never develop any ponding depth, no matter what the rainfall duration. This issue is particularly important when daily cumulative rainfall data are simply spread evenly over a model day—e.g., with the above example values a significant rainfall rate of 1.4 cm/day would simply be lost. A similar effect can occur when seepage from subsurface to surface is at low flux rates.

To address such issues, Frehg introduces “ponding storage” that tracks surface water volume in nominally dry cells. These are cells that considered dry in the surface water transport solution until the cumulative ponded volume exceeds $h_{min}$. Because the exfiltration depths of less than ${10}^{-5}$ m have limited importance, for simplicity the small ponding storage is ignored in the subsurface solution.

A further mass conservation problem for surface-subsurface exchange occurs when the predicted $q_{\mathrm{exchange}}$ is negative (infiltrating) and the volume $q_{\mathrm{exchange}}\mathsf{\Delta }t$ is greater than the available volume (V) in the surface water. In such a case, Frehg sets $q_e=-V\mathsf{\Delta }{t}^{-1}$ and reclassifies the cell as dry.

3.7. Comparing SWE and DW Models

Similar to Arico and Nasello [27], a set of three numerical dimensionless numbers are defined to evaluate the proximity between the SWE and DW models. These dimensionless numbers allow us to track the spatial-temporal variation of the relative importance for each term in the momentum equation. Consider a simple 1D flow and introduce the standard discrete notational convention that subscripts represent the location in the 1D grid space with i and $i+1$ as cell centers and $i+1/2$ as cell faces. Discrete numerical dimensionless numbers can be defined on the cell faces (where velocity is computed) as:

$${\widehat{Fr}}_{i+\frac{1}{2}}^n=\left|\frac{u_{i+\frac{1}{2}}^n(u_{i+\frac{1}{2}}^n-u_{up}^n)}{g\left({\eta }_{i+1}^n-{\eta }_i^n\right)}\right|$$

(18)

$${\widehat{Rt}}_{i+\frac{1}{2}}^n=\left|\frac{\mathsf{\Delta }x}{g\mathsf{\Delta }t}\frac{(u_{i+\frac{1}{2}}^n-u_{i+\frac{1}{2}}^{n-1})}{\left({\eta }_{i+1}^n-{\eta }_i^n\right)}\right|$$

(19)

$${\widehat{Rf}}_{i+\frac{1}{2}}^n=\left|\frac{n^2\mathsf{\Delta }x{(u_{i+\frac{1}{2}}^n)}^2}{2{(h_{i+\frac{1}{2}}^n)}^{\frac{4}{3}}\left({\eta }_{i+1}^n-{\eta }_i^n\right)}\right|$$

(20)

where $\widehat{Fr}$ is the square of the discrete Froude number, which is the ratio between the advective inertia and the hydrostatic pressure gradient, $\widehat{Rt}$ is the ratio between the unsteady term and the hydrostatic pressure gradient, $\widehat{Rf}$ is the ratio between the friction term and the hydrostatic pressure. n as the superscript represents time level, which needs to be distinguished from Manning’s n. Note that $u_{up}^n$ represents the upwind face velocity (e.g., $u_{i-1/2}$ for flow in the $+x$ direction), consistent with the 1st-order upwind scheme in Frehg. To examine the validity of the DW approximation, $\widehat{Fr}$, $\widehat{Rt}$ and $\widehat{Rf}$ are calculated using the results simulated with the SWE. For DW to be an appropriate simplification, we expect SWE results to provide $\widehat{Rf}\sim 1$ with $\widehat{Fr},\widehat{Rt}\ll 1$; i.e., friction approximately balances the hydrostatic pressure gradient while the advective inertia and unsteady terms are negligible.

4. Tests and Results

4.1. Overview

Three test problems are used to validate Frehg and examine performance with both DW and SWE algorithms. Validation of Frehg is in Section 4.2, which presents a test case from the literature for rainfall on a sloping plane with uniform hydraulic conductivity. The behavior of a simplified intertidal zone is examined in Section 4.3 to provide insight into the difference between the DW and SWE approximations that can be readily quantified. Frehg is exercised over a coarse-resolution model of the Trinity River Delta (Texas, USA) in Section 4.4 to examine the combined effects of evaporation, tide-/wind-induced wetting/drying and surface-subsurface exchange in a large-scale simulation. For ease of exposition, each of the following sections presents the model setup, results, and focused discussion of the individual case. Discussion and synthesis across the test cases is provided in Section 5.

4.2. Rainfall on a Sloping Plane

An idealized system with rainfall on a sloping plane has been used by Sulis et al. [45], Maxwell et al. [46], Wu et al. [47], and provides validation of the surface-subsurface coupling approach in Frehg. The system uses a simple sloping plane in the x direction with uniform properties in the transverse y direction, which can be modeled as a 2D x-z vertical section. The setup of the model is provided in

Table 1. Common parameters for test of rainfall on a sloping plane.

Parameter	Value	Units
L	800	m
$S_0$	0.0005	–
n	0.01986	m${}^{-1/3}$ s
$\varphi$	0.4	–
${\theta }_r$	0.08	–
$\alpha$	1.0	m${}^{-1}$
$\widehat{n}$	2.0	–
$\mathsf{\Delta }x$	80	m
$\mathsf{\Delta }z$	0.1	m
$\mathsf{\Delta }t$	2.0	s
$q_{\mathrm{rain}}$	$5.5\times {10}^{-6}$	m s${}^{-1}$
$T_{\mathrm{rain}}$	200	min
$T_{\mathrm{recession}}$	100	min

and

Table 2. Parameters varied for three test cases of rainfall on sloping plane.

Parameter	(1)	(2)	(3)	Units
$K_{s\left(i\right)}$	$1.16\times {10}^{-4}$	$1.16\times {10}^{-5}$	$1.16\times {10}^{-6}$	m s${}^{-1}$
$d_{i\left(i\right)}$	0.5	1.0	1.0	m

. The physical system is characterized by the plane length (L), slope ($S_0$), and the model parameters discussed in Section 3, above. The vertical and horizontal saturated hydraulic conductivities are set to a uniform value ($K_s$) in each test, with three different values tested ($K_{s\left(1\right)}$, $K_{s\left(2\right)}$, $K_{s\left(3\right)}$). The discretization uses uniform values for the horizontal grid ($\mathsf{\Delta }x)$, vertical grid ($\mathsf{\Delta }z$), and time step ($\mathsf{\Delta }t$). The initial conditions are a uniform water table depth ($d_i$) and zero water depth above the sloping plane. Boundary conditions include a rainfall rate ($q_{\mathrm{rain}}$), rainfall duration ($T_{\mathrm{rain}}$), recession duration ($T_{\mathrm{recession}}$) and a fixed surface elevation at the downstream boundary. Surface runoff at $(x,z)=(400,0)$ m is used for model-model comparison, which is consistent with Sulis et al. [45], Maxwell et al. [46]. Note that to use the Dirichlet downstream boundary condition for the free surface, the total plane length L is doubled compared to that in Sulis et al. [45], but the subsurface domain beyond $x=400$ m is deactivated. The fixed free surface elevation at outlet is much lower than $z=0$ m, which guarantees the measured runoff at $(x,z)=(400,0)$ m is unaffected by the downstream boundary conditions. The selection of the rainfall rate and the three $K_s$ values implies that $K_{s\left(1\right)}$ will require the subsurface to saturate before runoff is generated, whereas $K_{s\left(2\right)}$ and $K_{s\left(3\right)}$ will generate runoff while unsaturated.

Frehg is run separately with the DW and SWE algorithms for each of the $K_s$ values in Table 2. The results are compared with those simulated by the CATHY and ParFlow models [45], which use the kinematic wave (KW) approximation of the SWE (see Section 3.2). The comparitive study by Maxwell et al. [46] reveals negligible difference between KW and DW models for this test problem. For low slopes and low rainfall rates inertial terms can be expected to be small, so the results from the SWE should be close as well.

Figure 3 shows the surface runoff, at $(x,z)=(400,0)$ m of the sloping plane, which corresponds to the downstream outlet in Sulis et al. [45]. Figure 3a corresponds to a saturated hydraulic conductivity of $1.16\times {10}^{-4}$ m/s, which is higher than the rainfall rate. Thus, the unsaturated subsurface domain will be filled first before runoff is generated. In Figure 3b,c, the hydraulic conductivity is less than the rainfall rate, meaning that infiltration and runoff generation occur simultaneously. Since a deeper initial water table is used in Figure 3b,c, less runoff is generated despite their small conductivity. In all scenarios, Frehg has good agreements with CATHY and ParFlow in terms of runoff. Compared with solving the SWE, applying the DW approximation leads to slightly higher surface runoff, but the difference between SWE and DW does not exceed the difference between CATHY and ParFlow, indicating both SWE and DW are reasonable methods for this test problem. For Figure 3a,b, the DW method produces the highest flow rate on the rising limb, which is likely due to the asynchronous coupling scheme [46].

The depth and velocity at the $(x,z)=(400,0)$ m are shown in Figure 4. Similar data are not reported in Sulis et al. [45], Maxwell et al. [46], but are recommended as providing insight by more recent studies [18]. For all test scenarios, the DW approximation overestimates velocity while underestimating depth. A similar finding is reported in Cea and Blade [16]. Since flow depth affects inundation area, which affects surface-subsurface exchange, these results suggest the DW approximation could influence the surface-subsurface exchange in the intertidal zone that experiences wetting/drying. This idea is investigated in Section 4.3 and Section 4.4.

Figure 5 shows the numerical dimensionless numbers (Equations (18)–(20)) as functions of time. It can be seen that $\widehat{Rf}$ is close to 1 for all tested conductivity values. Furthermore, $\widehat{Fr}$ and $\widehat{Rt}$ are orders of magnitudes smaller than $\widehat{Rf}$. Because these non-dimensional numbers represent the relative importance of terms in the momentum equations (see Section 2), these results indicate that the pressure gradient and bottom friction are well-balanced and the inertia terms are relatively unimportant for this rainfall-runoff case. Thus, the use of the DW approximation is justified and the similarity of the SWE and DW results is expected for this test case.

4.3. Inundation of an Intertidal Zone

Of particular interest in modeling a tidally-driven marshland is the behavior of the surface-subsurface exchange with successive wetting and drying of a sloping landscape. As an idealized test case, we examine SWE and DW behaviors in the test domain shown in Figure 6, which has a uniform sloping face with a bed slope of 0.001. The subsurface domain extends 3 m below the lowermost point of the sloping surface plane. This domain is 2D and driven by sinusoidal oscillation of surface elevation, which represents idealized tidal forcing with period $T_{\mathrm{tide}}$ and amplitude $a_{\mathrm{tide}}$.

Model parameters for a baseline simulation are presented in

Table 3. Common parameters for modeling inundation of an intertidal zone. Parameter values for $\varphi$, ${\theta }_r$, $\alpha$, and $\widehat{n}$ are those listed in Table 1.

Parameter	Value	Units
$\mathsf{\Delta }x$	5	m
$\mathsf{\Delta }z$	0.05	m
$S_0$	0.001	–
$\mathsf{\Delta }t$	2.0	s
$K_s$	${10}^{-6}$	m s${}^{-1}$
$q_{e\left(\mathrm{surface}\mathrm{water}\right)}$	${10}^{-7}$	m s${}^{-1}$
$T_{\mathrm{spin}-\mathrm{up}}$	30	day
$T_{\mathrm{test}}$	3	day
$a_{\mathrm{tide}}$	0.1	m

. According to Equation (1) and (2), the DW approximation favors large time scale (T) and bottom friction (n). As previously discussed, due to strong variation of spatial scales in tidal-driven marshes, $Rs$ and $Rt$ have limited practical use as criteria to assess the validity of the DW approximation. However, the basic scaling relations in these equations should still hold, i.e., the DW model should converge to the SWE model as T and n increase. Thus, the following three test cases are examined (

Table 4. Parameters varied for three test cases of modeling inundation of an intertidal zone.

Parameter	(1)	(2)	(3)	Units
$n_{\left(i\right)}$	0.03	0.03	0.003	m${}^{-1/3}$ s
$T_{\mathrm{tide}\left(i\right)}$	24	12	24	hour

): (1) baseline with 24 h tidal period ($T_{\mathrm{tide}\left(1\right)}$) and Manning’s n of 0.03, (2) reduced tidal period, and (3) reduced Manning’s n. The initial water table elevation is same as the initial tidal elevation, which is 0.1 m above the zero-elevation plane in Figure 6. Note that the initial water table is set to a constant elevation, in contrast to the constant depth below surface used in Section 4.2. A model spin-up period ($T_{\mathrm{spin}-\mathrm{up}}$) is applied, during which the horizontal conductivity is increased to $10\times$ the baseline value of Table 3 to rapidly create reasonable head and moisture fields in the subsurface domain. Data are collected and analyzed from the SWE and DW simulations over the $T_{\mathrm{test}}$ testing period at the 1/2 and 3/4 locations shown in Figure 6. The first corresponds to the tide’s neutral level (the center of the intertidal zone), and the second is the upper edge of the continuously inundated zone.

Figure 7 shows evolution of the surface water over two days of simulation. In frames (a), (b) and (c), the normalized inundation area is defined as the inundation area divided by the maximum inundation area of the baseline SWE model. Figure 7a indicates the DW inundation area is similar to the SWE model, but has a slight phase lag. In Figure 7b the phase lag is increased when $T_{\mathrm{tide}}$ is halved, and in Figure 7c the phase lag is decreased (and almost eliminated) when Manning’s n is reduced by an order of magnitude. In these figures we can also see that increased phase lag is associated with a slight damping of the inundation area. However, as illustrated in Figure 7d–f, the phase lag in different tests does not significantly affect the depth evolution at the 1/2 and 3/4 sampling locations.

Surface profiles on the rising and falling tides are shown in Figure 8 to illustrate how the phase lag affects the water distribution. These effects can be interpreted as the DW introducing an artificial dissipation that scales directly with the bed roughness and inversely with the time scale.

For insight into the behaviors illustrated above, the discrete dimensionless numbers $\widehat{Fr}$, $\widehat{Rt}$ and $\widehat{Rf}$ of Equations (18)–(20) are plotted in Figure 9 for the SWE simulations. All three cases clearly meet the DW requirement with $\widehat{Fr}\ll 1$, which is similar to that proposed by [13]. This result indicates advective inertia (the nonlinear term in SWE, neglected in DW) is unimportant. However, if we focus on cases with larger n, i.e., Figure 9a–d, we find conditions where the DW requirements of $\widehat{Rf}\sim 1$ and $\widehat{Rt}\ll 1$ are violated. In particular, we see $\widehat{Rt}\sim 1$ during several episodes when $\widehat{Rf}\ll 1$, indicating that the unsteady term dominates the friction term, hence the DW approximation is violated. Thus, the discrepancies between DW and SWE in Figure 7a,b are largely explained by neglect of the unsteady term that becomes episodically important.

The results with the lower value of Manning’s n, Figure 9e,f are intriguing. These indicate the approximations of the DW method are never satisfied, i.e., $\widehat{Rf}\ll 1$ and $\widehat{Rt}\sim 1$, which indicates that the SWE momentum reduces to a balance between pressure gradient and unsteadiness throughout the simulation. Furthermore, yet, the phase lag and inundation area errors for this case in Figure 7c are substantially smaller than those in Figure 7a,b where the DW approximation is only violated episodically. The explanation for this behavior lies in the numerical discretization approach to using the DW approximation in an unsteady solver. The DW approximation, Equation (6), is formally a steady-state relationship between $\partial \eta /\partial x$ and friction (driven by u) that is used to close the unsteady equation for mass conservation, Equation (3). Thus, although the unsteady term is not included in momentum, unsteadiness will still exist in the coupled equation set. Furthermore, it can be shown (see Appendix A) that a discrete balance of unsteady momentum and the pressure gradient is equivalent to the DW approximation for small values of Manning’s n when u is slowly varying relative to the time step. Note that the collapse of these equations under limited discrete conditions should not be taken as a reason to use the DW approximation for flows where $\widehat{Rf}\ll 1$ and $\widehat{Rt}\sim 1$. Such conditions (with $\widehat{Fr}\ll 1$) should arguably be modeled with a linearized unsteady momentum equation (e.g., a wave equation).

Figure 10 shows the surface-subsurface exchange flux in the intertidal zone and the fully inundated zone. In the intertidal zone, the exchange flux switches between infiltration (negative) and exfiltration (positive) in accordance to the tidal cycle. During falling tide, evaporation and topography-driven exfiltration leave the top layer of the subsurface unsaturated. During rising tide, a strong infiltration flux is observed that fills the unsaturated subsurface voids. In the fully inundated zone, the exchange flux is two orders of magnitude smaller and is monotonically downward during the entire simulation period. These results highlight the significant role of the intertidal zone in enhancing surface-subsurface exchange. Neglect of tidal wetting/drying and the variably-saturated subsurface flow is likely to significantly underestimate the exchange fluxes.

The value of the full SWE, as compared to the DW, can be seen in Figure 10a,c,e. The surface water phase lag illustrated in Figure 7 reappears here as a change in the timing of infiltration and exfiltration. The SWE model in Figure 10 predicts higher maximum infiltration fluxes and longer exfiltration period because of its fast response to the tidal boundary condition. These results illustrate that model simplifications (i.e., using DW rather than SWE) affect both timing and magnitude of surface-subsurface exchange, especially in the intertidal zone.

4.4. Surface-Subsurface Exchange in the Trinity River Delta

The third test problem uses terrain of a real coastal wetland—the Trinity River Delta near Houston (Texas, USA). High-resolution lidar data was used to characterize the topography at the Trinity Delta and build the Trinity Delta Hydrodynamic Model (TDHM) reported in Li et al. [48]. The TDHM used the FrehdC code [40], which was a precursor of the Frehg code presented herein. Numerical studies with TDHM highlighted the important roles of boundary conditions, subgrid-scale topography and groundwater flow on surface hydrodynamics and salinity transport [48]. However, TDHM only models fully saturated groundwater flow, which provides limited insights regarding the surface-subsurface exchange process at the intertidal zone. The Frehg model enables simulation of variably-saturated subsurface flow, which allows for a preliminary investigation of surface-subsurface exchange on a realistic wetland topography. Unfortunately, obtaining adequate field data for validation of both models remains an unresolved issue.

Figure 11 shows the bathymetry used in the present study at 150 m resolution, which is upscaled from available 1 m resolution lidar data. The coarse resolution was selected to minimize the computational costs while allowing simulation over complex topography. Admittedly, this coarse resolution would be insufficient for validating the model or conducting detailed analyses of flow, but it should suffice for the purposes of evaluating differences between the DW approximation and the SWE over a sufficiently complex system. Further details on the topographic data and the Trinity River Delta are available in Li et al. [49]. All of the elevations used herein are given relative to the NAVD88 elevation datum.

Model parameters for the Trinity River Delta simulations are presented in

Table 5. Parameters for modeling the Trinity River Delta. Parameter values for $\varphi$, ${\theta }_r$, $\alpha$, and $\widehat{n}$ are as listed in Table 1. Parameter values for n, $K_s$ and $q_{e\left(\mathrm{surface}\mathrm{water}\right)}$ are as listed in Table 3.

Parameter	Value	Units
$\mathsf{\Delta }x$	150	m
$\mathsf{\Delta }y$	150	m
$\mathsf{\Delta }z$	0.25	m
$\mathsf{\Delta }t$	10.0	s
$T_{\mathrm{spin}-\mathrm{up}}$	30	day
$T_{\mathrm{test}}$	20	day
${\eta }_{river}$	0.3	m
$u_w$	2.0	m/s

. In the surface domain, the river boundary condition on the north side of the domain is provided by a constant water level, ${\eta }_{river}$. The southern boundary is open water in Galveston Bay, to which we apply a time-varying tidal elevation as an open boundary condition. The tidal data is extracted from the NOAA record for tide & current station 8770613 (Morgans Point, Barbours Cut, Texas) for year 2018, which is available at 6 min intervals. The subsurface domain extends down to $z=-3.0$ m with a uniform vertical grid resolution. The bottom and side boundaries are all impermeable. The soil parameters are same as those in Table 3. As in the simpler study in Section 4.3, above, during $T_{\mathrm{spin}-\mathrm{up}}$ the model is run with increased horizontal conductivity to rapidly build up reasonable head and moisture profiles in the subsurface domain. Analyses (below) is based on $T_{\mathrm{test}}$ period that covers a 14-day spring-neap cycle with a 3-day buffer on either side.

A total of 4 test scenarios (

Table 6. List of the four test cases for the Trinity Delta test problem.

Name	SWE	SWE + Evap	SWE + Evap + Wind	DW + Evap
Surface flow solver	SWE	SWE	SWE	DW
Evaporation	No	Yes	Yes	Yes
Wind	No	No	Yes	No

) are established to examine sensitivity to model simplifications and environmental processes. When evaporation is included, a constant evaporation flux of $1\times {10}^{-7}$ m/s is applied to wet regions. Equation (17) is used for dry regions. When wind stress is included, a constant wind speed of $u_w=$2 m/s pointing towards the north direction is enforced. Note that all tests in Table 6 are performed with the same $\mathsf{\Delta }t$ as listed in Table 5 for consistency, but a much larger value of $\mathsf{\Delta }t$ can be used when surface flow is modeled with SWE.

Figure 12 shows the numerical dimensionless numbers at locations L1 and L2 defined in Figure 11. L1 is located in a permanently-wet lagoon distant from the open boundary with surface fluxes in both x and y directions. In contrast, L2 is located in the intertidal zone close to the open boundary and has surface fluxes only in the x direction. It can be seen that at L1, $\widehat{Rt}$ is around 1 in both directions, indicating unsteadiness is important at this location, which is consistent with results for the simpler study of inundation in the intertidal zone, in Section 4.3. Similarly, $\widehat{Fr}$ is typically less than 1. We see flow direction becoming important with $\widehat{Rf}$, which is less than 1 in the x direction, but close to 1 in y direction. At L2, $\widehat{Rf}$ is close to 1, while $\widehat{Rt}$ and $\widehat{Fr}$ are mostly much less than 1. However, a few higher values of $\widehat{Rt}$ are also observed. Detailed examination (not shown) reveals that these high $\widehat{Rt}$ values at L2 occur during high tide.

Overall, the above observations imply that due to unsteadiness, the DW approximation is generally invalid in the Trinity Delta unless the depth is small. This result seems contradictory to the simple scaling argument of Equation (2) that states increasing importance of unsteadiness at small depth. The reason might be the thin-layer drag model used in Frehg [10], which increases bottom friction at small depth and stabilizes the flow field (this drag model is deactivated in Section 4.2 for a consistent comparison to the benchmark results). It should be noted that the coarse grids used in this test problem prevents investigation on the effects of small-scale topographical features. At finer grid resolution with more complex topographical features, advective inertia could become important, which brings further question to the validity of the simplified DW equations [50].

Figure 13 shows the modeled depth, exchange flux and water content at L1 and L2. Figure 13a indicates that, although L1 is relatively far from the tidal boundary, the depth (hence water surface elevation) still oscillates at tidal frequency. The reason is that we use a relatively coarse grid resolution ($\mathsf{\Delta }x=150$ m), which neglects small-scale water-blocking structures and enhances surface connectivity [10]. At L1, the 3 scenarios with SWE produce negligible differences, but the DW model overestimates amplitude of depth oscillations.

Figure 13b shows that during spring tide, the intertidal zone experiences wetting/drying transitions. During neap tide, the intertidal zone experiences continuous dry period. The effect of evaporation in removing thin layers of water during falling tide is evident. Without evaporation, L2 never becomes completely dry during spring tide. Compared to SWE+evap, additional inundation events are found on day 1 and day 3 when wind effect is modeled. Additional inundation events are found on day 1, 3, 12, 17 and 20 when SWE is replaced by DW.

These differences in modeled inundation patterns have influence on surface-subsurface exchange. Figure 13c–f display the exchange flux and the water content of the top-most subsurface cell at L1 and L2. As long as evaporation is modeled, the peak infiltration flux at L2 is two orders of magnitudes larger than that at L1, which highlights the critical role of the intertidal zone in promoting surface-subsurface exchange. By comparing Figure 13b,d,f, the occurrence of these extreme infiltration fluxes is caused by evaporation and wetting/drying: during low tide or neap tide, evaporation leads to reduced water content in the subsurface domain (hence large negative pressure head). When the intertidal zone is flooded again, since increased head gradient is formed and plenty of void space becomes available, large infiltration flux is observed. Factors affecting the wetting/drying status of the intertidal zone, such as wind and use of the DW approximation, alters timing and magnitude of the exchange flux because they have influence on the evolution of the inundation area, which affects subsurface water content (Figure 13f). These findings and mechanisms are similar to those revealed in Figure 10, but they further illustrate the role of spring-neap cycle in depleting the subsurface domain. In permanently-wet regions such as L1, the subsurface domain is unaffected by evaporation and is always fully saturated (Figure 13e), so switching between SWE and DW has minimal influence on the exchange flux (Figure 13c). Thus, the key point for simulating surface-subsurface exchange is to accurately capture the wetting/drying status in the intertidal zone. This relies on the use of the SWE (rather than the DW), the inclusion of the variably-saturated subsurface flow and important environmental processes such as wind and evaporation. Arguably, accurately predicting the evolution of inundation area also depends on topography and grid resolution, which will be investigated in future studies.

To quantitatively illustrate the importance of the intertidal zone, we can compare the area of the intertidal zone (defined as the area that experiences wetting/drying transition at least once during the simulation period) and the exchange flux through this zone to the total area and total exchange flux. Our analysis indicates that for SWE+evap, the intertidal zone is only 3.34% of the total wet area in the Trinity Delta but provides 11.2% of the surface-subsurface exchange flux; hence the intertidal zone is of outsized importance in the exchange.

5. Discussion and Conclusions

In the present study, a coupled surface-subsurface flow model (Frehg) is developed to understand the effects of various model simplifications that are often adopted in shallow marsh simulations. The surface flow module in Frehg has been previously tested in existing shallow marsh studies. The subsurface flow module is implemented based on the predictor-corrector-allocation scheme by Li et al. [12]. The coupled model is (i) validated using a rainfall-runoff test problem in the literature, (ii) applied to simulate tidal-driven flow on a synthetic sloping plane, and (iii) demonstrated on a coarse-grid simulation of real marsh topography in Texas, USA.

Our study in Section 4.2 supports the prior literature that the DW approximation is appropriate and falls within the expected non-dimensional scales for rainfall runoff on a slope. Compared to the SWE, the DW approximation tends to underestimate flow depth and overestimates velocity but the errors are typically within acceptable limits.

However, for intertidal flows on a slope in Section 4.3, we found the DW approximation is not supported by non-dimensional analyses of simulation results. Formally, the unsteady term that is neglected in DW becomes important in the intertidal zone. The effect of neglecting the unsteady term in the setup for a time-marching solution (as presented in Section 3.2) is that the scheme develops increased dissipation, which shows up as phase lags in tidal motions. Interestingly, by decreasing Manning’s n, which is in the denominator for u in Equation (7), the effect of the dissipation might be tuned – allowing the “right answer for the wrong reason.” Such an outcome might account for successful use of the DW approximation in conditions where the strict scaling analyses might not hold.

The study of the Trinity River Delta in Section 4.4 indicates that unsteady tidal motions across the wetting/drying intertidal zone have an outsized impact on surface-subsurface exchange, with roughly 3% of the area contributing 11% of the exchange flux. Thus, models that poorly represent the intertidal zone will have their error amplified in exchange flux predictions. As our work shows that the DW approximation is not well-founded for predicting the extent, timing and duration of inundation, it follows that the SWE should be preferred in coupled surface-subsurface simulations of tidal marshes.

The present work is limited in that we have not tested effects caused by the neglect of the advective inertia in the DW approximation. Within our test cases, advective inertia is relatively unimportant with $\widehat{Fr}<1$ throughout. However, one can imagine that small-scale topography (neglected in our coarse-grid model) could enhance non-linearity of the surface flow and increase the importance of the advective term, even while $\widehat{Fr}$ remains small. This topic deserves further investigation in the future.

Although the fully-inundated area of a marsh is typically the largest area, our work shows that the smaller intertidal area plays an outsized role due to evaporation. The exposure of the bare surface during low tide leads to decreased water content (i.e., increased gradient of pressure head) in the subsurface domain, which promotes infiltration upon flooding. It follows that the succession from neap to spring tides that change the inundation patterns and extents will play a vital role. Thus, the key to model surface-subsurface exchange in shallow wetlands is to accurately capture the range and the inundation of the intertidal zone. To achieve this, the model should be equipped with solvers for the full dynamic surface flow equations (i.e., the SWE) and variably-saturated subsurface flow equations. The critical environmental processes of evaporation and rainfall cannot be neglected. Simplified simulations that adopt DW approximations or neglect processes are unlikely to provide high-fidelity results, although they may be tuned to correctly represent observations.

On a more speculative note, the perceived importance of the intertidal zone in simulation of surface-subsurface exchanges may be closely tied to the model grid scale, which we have not investigated herein. Topography of a marsh has fractal features, with increasing complexity and extent of the wetting/drying front becoming visible at different scales. As the grid scale is refined we expect to see greater variability across the wetting/drying front, which will likely further increase its importance in the overall exchange budget. However, it remains to be seen as to what scale is necessary to accurately define the temporal dynamics of the intertidal areas and understand their impact on the exchange fluxes.

The present study only reports preliminary development and application of the Frehg model with a goal of understanding whether or not a DW approximation could be substituted for the SWE for a simpler model. To that end, we have determined that the SWE is necessary. The main limitations of the present study are rooted in the lack of field data to provide more robust validation and insight. In the future, the model can be improved by including high-resolution topography, density-driven salinity transport, and time-varying meteorological inputs.