Multiple Seepage‐Faces in Tidal Flat With Very Gentle Slopes

Large‐scale seepage‐faces occur on small‐slope tidal flats. All previous studies assume that the seepage‐face has only one single exit point. Here we show via numerical simulations of tidally‐influenced groundwater flow that, in a two‐dimensional vertical, homogeneous transect of a tidal flat with gentle beach slope of 1‰, multiple seepage‐faces may occur with at most four unsaturated beach surface segments which separate four seepage‐faces. Salinity‐variation induced density‐dependent flow leads to this complex phenomenon. While the seepage‐faces are the groundwater discharging zones on the beach surface, the unsaturated zones are the recharging zones. The whole aquifer beneath the tidal flat is almost occupied by seawater and forms a wall blocking the horizontal seaward discharge of inland fresh groundwater so that inland freshwater discharges mainly occur near the high tide mark. This is in great contrast with the traditional results that inland freshwater discharge occurs mainly near low tide mark.

with time are discussed. The reasons for leading to these multiple seepage-faces are analyzed. To this end, two different salinity values are set for seawater to compare the numerical results. One is 35 g L −1 , the standard salinity for seawater, and the other is zero salinity (i.e., fresh water). Numerical simulations results show that the behaviors of seepage-faces of these two cases have distinct differences and can help us to understand the formation and development of multiple seepage-faces in case of existence of density-dependent groundwater flow. Seepage-faces are kind of important boundary for coastal groundwater discharge. In order to accurately model the groundwater flow and solute transport in intertidal aquifers, the complex nonlinear boundary conditions both for the groundwater flow and solute transport on the seepage-faces should be strictly and accurately implemented. This paper made such an attempt.

Model Setup
Consider a two-dimensional vertical cross-section model domain of a coastal aquifer perpendicular to the coastline shown in Figure 1. The domain is horizontally 5 km long with the middle part lying in the intertidal zone. The aquifer has a horizontal impermeable bottom, and its depth is 35.5 m at the landward boundary and 29 m at the seaward boundary. The tidal flat has a horizontal extension of 3.5 km with a slope of 1‰. To model the effects of a slope break, the tidal flat also includes a 300 m long steeper left part with a slope of 1%. Such a landward slope increment is common in reality (e.g., Hou et al., 2016;Ma et al., 2015;Robinson, Gibbes, et al., 2007;Xia & Li, 2012). The mean sea level is set to be the elevation datum, that is, z = 0.0 m.
The governing equations of the two-dimensional density-dependent, variably-saturated groundwater flow and solute transport are given in Supporting Information S1. They are generalization of the original model given by MARUN (Boufadel et al., 1999) in the sense that the tidal loading effects are considered.
Boundary conditions are summarized in Table S1 of Supporting Information S1. On the beach surface or seabed submerged by the tidal seawater (Γ Sea ), the tidal moving boundary conditions Xia et al., 2010) are applied, namely, specified water pressure determined by the tidal seawater depth above the boundary point, and flow-direction dependent salinity (constant seawater concentration where flow is entering the domain, zero dispersive flux normal to the boundary where flow is leaving the domain). The tidal signal reflecting diurnal tide is defined by the following equation: where h tide [L] is the tidal level, A is the amplitude [L], and T p is the tidal period equaling 24 hr to simulate diurnal tide.
The ground or tidal flat surface higher than sea level (EF + FG + GA, see Figure 1) comprise either the unsaturated zone(s) Γ U or saturated zone(s) Γ S (i.e., seepage-faces). Here both Γ S and Γ U may comprise several segments separated by each other, as schematically shown in Figure 1. On Γ S , the boundary conditions for groundwater flow are given by (e.g., Xiao et al., 2017;Xin et al., 2017;Wilson & Morris, 2012), namely where ψ is the equivalent freshwater pressure head [L], q n is the outward Darcy velocity [LT −1 ] normal to the boundary. Equation 2 describes the fact that the saturated seepage-face is exposed to the air so the pressure head there is zero. The inequality (Equation 3) describes the fact that the groundwater flow velocity normal to the seepage-face boundary should be outward. On the exposed unsaturated surface Γ U no-flow boundary condition is used for groundwater flow, which implies neither rainfall infiltration no evaporation, namely For salinity the following zero-dispersion boundary condition is used on both Γ S and Γ U : where c is the groundwater salinity [ML −3 ], the term D [L 2 T −1 ] represents the hydrodynamic dispersion tensor (see Equation A10 in Supporting Information S1), ∇ = ( ∕ ∕ ) , ⃖ ⃑ is the outward unit vector normal to the seepage-face boundary. On Γ U , using Equation 4 the boundary condition defined by Equation 5 can be further simplified into The challenges to deal with boundary conditions defined by Equations 2-6 on Γ S and Γ U are to find the unknown interface point(s) (I 1 , I 2 , I 3 ,…, etc. in Figure 1b) of the seepage-face(s) and unsaturated zone(s). Iterative numerical method to find these points will be given later.
For the initial salinity condition a sharp salt-freshwater interface was recommended so that landward of the high tide mark the salinity was zero and seaward, standard seawater salinity of 35 g L −1 . Different initial salinity conditions may also work but need longer computation time to obtain periodic numerical solution with respect to salinity. The initial condition for the pressure head was a vertically hydrostatic pressure distribution (Li et al., 2008).
Typical model cases were designed to analyze the developments of seepage-faces. Model parameter values for the base case and freshwater case are shown in Table S2 of Supporting Information S1. These parameter values were chosen based on the literature investigations of coastal aquifer systems beneath tidal flats (Batu, 1998;Domenico & Schwartz, 1990;Freeze & Cherry, 1979;George et al., 2015;Hou et al., 2016;Kuang et al., 2020;Morris & Johnson, 1967), the details are given in Supporting Information S1.

Numerical Implementation
Numerical simulations were conducted with the 2-D density-dependent, variably saturated groundwater flow model MARUN (Boufadel et al., 1999). MARUN code has been successfully applied in many previous studies, such as simulations of temporal variation of one single seepage-face scales in a dimensionless beach (Li et al., 2008), flow and its ecological effects in a brackish marsh (Xiao et al., 2017). The boundary conditions defined by Equations 2-5 incorporated in the original MARUN code assumes that the seepage-face boundary Γ S and the unsaturated tidal flat Γ U has only one single interface point I 1 , which may be unrealistic in certain circumstances. Here the MARUN code was improved to assume that Γ S and Γ U may have several interface points (I 1 , I 2 , …, I m ) as shown in Figure 1. The following iterative algorithm was used in the numerical simulation to determine the locations of these interface points: 1. Assume there is no seepage-face at the initial time, and assume no-flow and zero-dispersion boundary conditions on all nodes higher than the tidal sea level to obtain the numerical solution; the set of the interface point of Γ S and Γ U of the numerical solution in this case is empty, that is, P (1) I = {} . If the set of the interface point of the last time step is not empty, use the set of the last time step as the initial set of this time step. 2. To update the boundary conditions of for groundwater flow (or equivalently ψ) on all nodes higher than the tidal sea level according to the newly-updated numerical solution. For a given node, a) If ψ > 0, then let ψ = 0 on this node by implement Dirichlet boundary condition, b) If ψ < 0, then no-flow boundary condition is implemented. c) If ψ = 0, then the normal velocity q n is further considered, which may be either nonnegative or negative. If q n ≥ 0, then regard this node as seepage-face and let ψ = 0 on this node by implement Dirichlet boundary condition, otherwise, if q n < 0, then regard this node as unsaturated and no-flow boundary condition is implemented. For all the cases (a-c) above the boundary condition for salinity is zero-dispersion defined by Equation 5. 3. Use the new boundary conditions in (2) to update the numerical solution. Identify and record the set of interface points of Γ S and Γ U of the numerical solution where m is the number of interface points and k is the iteration number of numerical solution (i.e., number of update). 4. Convergence judgment: for three interface sets P ( −2) I , P ( −1) I and P ( ) I in three successive iterations, if they meet each of the following conditions, then the iterations with respect to seepage-face boundary conditions are regarded as converged: a) The sets P ( −2) I , P ( −1) I and P ( ) I have the same number of interface points, b) For each interface point ( ) in P ( ) I , either ( ) = ( −1) , or ( ) = ( −2) holds (i = 1,⋯,m). Here ( ) = ( −2) is used to describe convergence since the iterative accuracy of the location of the interface points depends on the distances between adjacent nodes on the tidal flat. If the locations of an interface point in three successive iterations changes between two adjacent nodes, it is regarded that the iteration is converged. 5. If the iterations converge, go to other calculation tasks such as next time step (if no solute transport is considered) or solving the solute transport equation (if solute transport is considered), otherwise, go to step 2 for next iteration.
The domain in Figure 1 was discretized into a mesh with 48,421 nodes and 94,400 triangular elements. The mesh was relatively fine near the beach surface in the intertidal zone (triangles with a length of 3.5 m and height of 0.05 m), and became gradually coarse toward the bottom boundary. The height of the triangular elements at the bottom was the largest (1.575 m). The numerical models were running for 100 years of model time so that both the salinity and pressure head becomes periodical with time.

Seepage-Faces
The numerical simulation results of the base case related to seepage-faces are shown in Figure 2. For the sake of comparison, the results of the freshwater case (the salinities of seawater and initial groundwater are set to zero) are also shown. Figure 2a shows how the elevations of the upper and lower ends of seepage-faces change with time during one tidal cycle. During the whole tidal cycle, 6 different seepage-faces can be observed for the base case but only one seepage-face for the freshwater case, and the tidal level is always act as a lower end for the seepage-faces. For example, the highest seepage-face lasts for about 2.5 hr (from t = 0-2.5 hr). At low tide (t = 12 hr), four different seepage-faces occur for the base case. The pressure head and normal velocity on the beach surface at low tide are shown in Figure 2b, and the spatial distributions of salinity and flow field at low tide are shown in Figures 2c (full scale), 2d (local scale near high tide mark), and 2e (local scale near an unsaturated tidal flat zone) for the base case, and in Figure 2f for the freshwater case. In Figure 2b, the values of pressure heads for both the base case and freshwater case equal zero whenever the Darcy velocities normal to the beach surface are positive, indicating the correct and accurate numerical applications of the seepage-face boundary conditions defined by Equations 2 and 3.
The following observations can be made by comparisons of the results of freshwater case and base case in Figures 2a and 2b. (a) Multi-seepage-faces occur only for the base case while there is only one seepage-face for the freshwater case. (b) The normal velocity on the seepage-face presents unstable oscillating around a mean value for the base case, but becomes smooth for the freshwater case. The oscillating normal velocity is caused by unstable downward density-gradient which will be described below.
The reasons for the occurrence of the multi-seepage-faces can be analyzed by comparing the base case and freshwater case. Since there is only one single seepage-face for the freshwater case, obviously the multi-seepage-faces are caused by the density difference between seawater and freshwater. More precisely, such density effects can be classified into two aspects. On the one hand, for the base case the whole aquifer beneath the tidal flat is occupied by seawater and there forms a "salty wall" with high density under near the high tide mark (Figures 2c and 2d). This wall essentially inhibits the further horizontal seaward moving of inland fresh groundwater (Figures 2c and 2d; Figure S1 in Supporting Information S1) and most freshwater discharges out of beach surface from the landward side of the wall. Due to the entrainment caused by dispersion in the salt water-freshwater mixing zone, some salt groundwater originated from seawater also discharges together with the freshwater, thereby forming a clockwise seawater-groundwater circulation under near the high tide mark. As a result, significant downward seawater recharge occurs landward of the high tide mark, which eliminates the seepage-face there (Figures 2a  and 2c).
On the other hand, Landward of the above-mentioned clockwise seawater-groundwater circulation, the groundwater flow direction is generally horizonal and seaward. Thus, there forms a groundwater divide near x = 545 m. As the location of divide changes with time, groundwater with lower salinity landward of the divide will move seaward, therefore the salinity on the divide is always less than that of the seawater and decreasing with depth (Figures 2c-2e). This kind of salinity distribution extends seaward, forming a slight downward salinity gradient (or equivalently downward density gradient, see Figures 2c and 2e) in the shallow tidal flat aquifer. This downward salinity gradient may induce unstable oscillations of the normal velocity on the beach surface and sometimes even downward groundwater flow, leading to forming the unsaturated zone on the beach surface with a "salinity contour ridge" underlying (e.g., x = 2,000 m, Figures 2c and 2e). As a result, the downward groundwater flow below the unsaturated zone will squeeze the groundwater in the deep part of the aquifer and enhance the outflow velocity on adjacent seepage-faces resulted in the formation of groundwater discharge hotspot (e.g., x = 1,500 m and x = 2,300 m, Figures 2b and 2c).
Although there are many previous studies that conducted either field motoring or numerical simulations of the seepage-faces on various tidal flats (Li et al., 2008;Turner, 1993Turner, , 1995Vousdoukas, 2014), the beach slopes considered in these studies are not less than 1.75% and no muti-seepage-faces are found. To the authors' knowledge, this is the first study that multi-seepage-faces occur on a homogeneous tidal flat with very gentle-slope of 1‰. These facts imply that gentle-slope of 1‰ of the tidal flat is an important factor for the occurrence of multi-seepage-faces.

Groundwater Flow and Salinity
Tidally-averaged groundwater flow field and salinity distribution of the base case and freshwater case are shown in Figure 3. The following observations can be made.
1. For the base case, the whole aquifer below the intertidal zone is almost occupied by seawater and forms a vertical salty wall near the high tide mark where salinity increases dramatically from 10% seawater salinity to 99% seawater salinity and forms a clockwise circulation under near the high tide mark (Figure 3). The reasons for the formation of this circulation are explained above. 2. Seaward of the above-mentioned clockwise seawater-groundwater circulation, there is an anti-clockwise seawater-groundwater circulation, where seawater enters the aquifer from the intertidal zone and discharges from the subtidal zone. The groundwater divide of the two reverse circulations is located near x = 550 m.
In this anti-clockwise seawater-groundwater circulation, the magnitude of the groundwater Darcy velocity decreases seaward approximately exponentially. 3. For the freshwater case, there is no clockwise seawater-groundwater circulation near the high tide mark (Figures 3b and 3d). There is only one anti-clockwise seawater-groundwater circulation in the whole aquifer.
When the beach slope is relatively steep (>1.75%), previous studies on beach tidal groundwater flow (e.g., Boufadel, 2000;Fang et al., 2021Fang et al., , 2022Greskowiak, 2014;Heiss, 2020;Heiss & Michael, 2014;Heiss et al., 2020;Li et al., 2008;Robinson, Li, & Barry, 2007;Shen et al., 2019;Xiao et al., 2019) identified that in the aquifer beneath the beach surface, if the seaward direction is chosen to be rightward, there are usually an upper saline plume related to the intertidal zone with anticlockwise seawater-groundwater circulation, a lower saltwater wedge with clockwise seawater-groundwater circulation below the subtidal zone, and a fresh groundwater discharge tube between them. The salinity distribution, groundwater flow fields and seawater-groundwater circulation patterns for the base case here are completely different from those in these previous studies where the beach slopes considered are not less than 1.75%, which is much greater than 1‰. Both the groundwater discharge tube and the lower saltwater wedge with clockwise seawater-groundwater circulation below the subtidal zone disappears. Almost all the intertidal aquifer and the whole subtidal aquifer are occupied by a salt plume very close to the seawater salinity with an anticlockwise seawater-groundwater circulation (Figure 3a). A new anticlockwise seawater-groundwater circulation forms under near the high tide mark. Thus, we believe that the beach slope has significant impacts on the formation of multi-seepage-faces, and consequently on the tidal groundwater flow and solute transport beneath the beach surface.

Conclusions
Numerical simulations of tidally-influenced groundwater flow and solute transport are conducted in a two-dimensional vertical, homogeneous transect of a tidal flat with gentle beach slope of 1‰. It is found that multiple seepage-faces may form on the tidal flat. The reasons for the occurrence of these multi-seepage-faces are explained by (a) comparisons of the numerical simulation results for the base case (standard seawater with salinity of 35 g L −1 ) and freshwater case (a hypothetical case setting zero salinity for the seawater) and (b) analyses of the groundwater flow field and salinity distributions. The following conclusions can be drawn: 1. Multiple seepage-faces occur only for the base case, and there is only one seepage-face for the freshwater case. This indicates that the multi-seepage-faces are induced by salinity-variation induced density-dependent groundwater flow. 2. A clockwise seawater-groundwater circulation forms under near the high tide mark. Seawater infiltrates the aquifer near seaward of the high tide mark, flows downward and slightly landward, mixes with the freshwater from inland, moves upward and slightly seaward, and finally discharges near landward of the high tide mark. 3. Seaward of the above-mentioned clockwise seawater-groundwater circulation, there forms another anti-clockwise seawater-groundwater circulation. Seawater infiltrates the intertidal zone and discharges from near the low tide mark and the subtidal zone. 4. In the shallow aquifer of the intertidal zone, the salinity decreases slightly with depth unevenly, and forms heterogeneous density gradient, which may induce unstable oscillations of the normal velocity on the beach surface and sometimes even downward groundwater flow, forming the unsaturated zone on the beach surface. The downward groundwater flow below the unsaturated zone will squeeze the groundwater in the deep part of the aquifer and enhance the outflow velocity on adjacent seepage-faces. 10.1029/2023GL104189 9 of 10 5. As a result of (2), (3) and (4), multi-seepage-faces occurs on the tidal flat. 6. The salinity distribution, groundwater flow fields and seawater-groundwater circulation patterns in an aquifer beneath a tidal flat with very gentle slope of 1‰ are completely different from those in these previous studies where the beach slopes considered are not less than 1.75%, which is much greater than 1‰. Small beach slopes and density-difference between seawater and freshwater have determinant impacts on the formation of multi-seepage-faces, and consequently on the tidal groundwater flow and solute transport beneath the beach surface.

Data Availability Statement
The input files and output results of the numerical models are available online (https://doi.org/10.5281/ zenodo.7858053).