Effects of Unsaturated Flow on Salt Distributions in Tidally Influenced Coastal Unconfined Aquifers

Unsaturated flow influences both the seawater extent under steady‐state conditions and the propagation of tides in coastal aquifers. However, its effects on salt distributions in tidally influenced coastal aquifers are little investigated. The present study used numerical simulations and data from laboratory experiments to analyze the effects of unsaturated flow on density‐dependent solute transport in coastal unconfined aquifers. The effects of the inland boundary condition (i.e., constant‐head or constant‐flux) were tested. Compared to a stable sea level, the results show that unsaturated flow has a more pronounced influence on salt distributions in coastal unconfined aquifers when tides are considered, regardless of the type of inland boundary condition. Neglect of unsaturated flow effects leads to expansion of the upper saline plume (USP), shrinkage of the saltwater wedge (seaward movement of saltwater wedge), and overestimation of water and salt exchange across the aquifer‐ocean interface. This is caused by a lower head in the nearshore area during high‐tide periods with the unsaturated zone effects removed. Thus, without the unsaturated zone, stronger head gradients within the nearshore aquifer occur at high tide, leading to stronger tidally driven seawater infiltration and hence a larger USP. Counterintuitively, ignoring unsaturated flow effects leads to greater average inland head over a tidal period, which shifts the saltwater wedge seaward. It is concluded that unsaturated zone effects should not be neglected for modeling tide‐affected seawater intrusion, especially if quantification of near‐shore conditions is important.

Generally, the amplitude of tidally driven watertable fluctuations decreases, and the phase lag increases with distance from the shoreline (Ferris, 1952). Unsaturated flow enhances the inland propagation of tidally driven watertable fluctuations in coastal unconfined aquifers (e.g., Barry et al., 1996;Jeng et al., 2005;Kong et al., 2013Kong et al., , 2015Li et al., 2000). This occurs because neglecting unsaturated flow effects leads to a larger specific yield that weakens the tidal propagation (Luo et al., 2023). Thus, it is expected that the effects of tides on salt distributions are greater when the unsaturated zone is included. Although Geng et al. (2014) used numerical modeling to show that unsaturated flow shortens the transit time of land-driven pollutants in coastal unconfined aquifers subject to regular waves, their analysis neglected the density contrast between seawater and fresh groundwater that is a strong control on coastal groundwater flow processes.
While the unsaturated zone is often incorporated into studies of salt distributions in coastal aquifers (e.g., Ataie-Ashtiani et al., 1999;Kuan et al., 2012;Shen et al., 2020;Xin et al., 2010;Yu et al., 2017), until recently the effects of the unsaturated zone on estimates of salt distributions in coastal aquifers were unexplored. Luo et al. (2022) derived an analytical solution to examine unsaturated flow effects on salt distributions in coastal unconfined aquifers. They indicated that neglecting unsaturated flow leads to a more seaward saltwater wedge in the absence of tidal fluctuations. This occurs because the unsaturated zone transmits part of the freshwater flux towards the sea that otherwise flows entirely below the watertable in models without an unsaturated zone, leading to a lower watertable and hence a more seaward saltwater wedge in models that include the unsaturated zone. It is likely that the inclusion of tidal fluctuations complicates the effects of the unsaturated zone on salt distributions relative to the non-tidal findings of Luo et al. (2022).
This study aims to evaluate the effects of unsaturated flow on salt distributions in tidally influenced coastal unconfined aquifers, building on the analysis of Luo et al. (2022) concerning unsaturated zone effects in non-tidal coastal aquifers. The approach adopted here first compares numerical models of variably saturated, variable-density tidal propagation with existing laboratory experiments. Subsequently, numerical modeling is used to address two questions: (a) To what extent does unsaturated flow affect salt distributions in tidally influenced coastal unconfined aquifers? (b) What are the key factors controlling unsaturated flow effects on salt distributions under tidal conditions? Figure 1 shows the conceptual model of an idealized coastal unconfined aquifer. As mentioned previously, the USP develops above the saltwater wedge because of tidal fluctuations, causing fresh groundwater discharge The arrows indicate flow directions. One of two options is considered as the inland boundary condition: Constant head (h l ) or constant flux (q f ). The distance from the mean sea level to the aquifer base is denoted as d. Green and white arrows refer to freshwater flow and saltwater flow, respectively. Note that for the sake of conciseness, the unsaturated zone indicated here includes the capillary fringe occurring above the watertable.

Conceptual Model
10.1029/2022WR032931 3 of 18 (FGD) to occur between the wedge and the USP. The origin of the x-z axes is placed at the junction (O) of the inland boundary OA and the impermeable base OF, with the x axis pointing seaward and the z axis directed upward. AB is an atmospheric boundary condition. Previous studies showed that inland boundary conditions play an important role in controlling salt distributions in coastal unconfined aquifers (e.g., Luo et al., 2022;Werner & Simmons, 2009). Therefore, when conducting numerical simulations, either a constant hydraulic head h l [L] or a constant inland flux q f [L 2 T −1 ] was specified at OA, referred to respectively as head-controlled and flux-controlled conditions hereafter following Werner et al. (2012). BEF is the sea boundary subject to tidal fluctuations. The mean sea level is located at height d [L] above the impermeable base. During the tidal ebb period, a seepage face appears between the contact point of the watertable at the shoreline and sea level. Consistent with previous studies (e.g., Ataie-Ashtiani et al., 1999;Kuan et al., 2019;Robinson, Li, & Barry, 2007;Xin et al., 2010), the following assumptions are made: (a) the aquifer is isotropic and homogeneous (b) any recharge (e.g., rainfall), or discharge (e.g., pumping and evapotranspiration) through AB is neglected, and (c) hysteresis effects within the unsaturated zone are ignored. Shen et al. (2020) conducted laboratory experiments of freshwater-seawater interactions in a sand flume with dimensions of 7.7 m (length) × 0.16 m (width) × 1.0 m (height) (Figure S1a in Supporting Information S1). In their experiments, a constant hydraulic head of 0.74 m was maintained at the inland boundary, whereas either a steady sea level or a sinusoidal tide (with amplitude A of 0.075 m and period T of 240 s) was considered at the sea boundary. The mean sea level was 0.7 m above the impermeable base. The sand used in the experiments had an effective porosity (n e ) of 0.4 and a saturated hydraulic conductivity (K s ) of 3.0 × 10 −3 m s −1 .

Laboratory Experiments
Experiments displaying freshwater-seawater interactions under tidal fluctuations were also conducted by Kuan et al. (2019). They used a sand flume with dimensions 3.4 m (length) × 0.02 m (width) × 0.72 m (height) (Figure S1b in Supporting Information S1). A constant flux of 2 × 10 −5 m 2 s −1 was imposed at the inland boundary. In addition, experiments to examine salt distributions under a constant sea level and a sinusoidal tide (A = 0.043 m and T = 62 s) were undertaken. The mean sea level was 0.55 m above the impermeable base. The sand used in the experiments had n e of 0.46 and K s of 5.26 × 10 −3 m s −1 .
In the experiments of Shen et al. (2020) and Kuan et al. (2019), the salinity (C) at the inland and sea boundaries remained constant at 0 and 35 ppt, respectively. The salt distribution was identified with red food dye and was captured by digital photography. In the experiments of Shen et al. (2020), the thickness of the unsaturated zone (i.e., the depth to the watertable from the upper surface) varied due to tides from 0.225 to 0.375 m, whereas it varied from 0.127 to 0.213 m in the experiments of Kuan et al. (2019).

Numerical Simulations
The current study applied SUTRA (Voss & Provost, 2010), a well-known variable-saturation and variable-density flow and transport model, to simulate the experiments of Shen et al. (2020) and Kuan et al. (2019). The numerical setups were in accordance with the experiments. With reference to the schematic shown in Figure 1, AB and OF were treated as no-flux boundaries, while OA was set to a constant-head boundary in simulating the experiments of Shen et al. (2020), whereas it was a constant-flux boundary for modeling the experiments of Kuan et al. (2019). Following the approach of Xin et al. (2010), the shoreline boundary (BEF), which was subjected to tidal fluctuations, was represented in SUTRA using a time-varying boundary condition that included seepage face development. Specifically, submerged nodes were assigned hydrostatic pressures reflecting the sea level. For other nodes on BEF, the soil saturation was taken from the previous time step. The local pressure of saturated boundary nodes above sea level was set to atmospheric pressure. Otherwise, boundary nodes were assumed to be no flow in the direction perpendicular to the shore.
The widely used van Genuchten-Mualem model (Mualem, 1976b;van Genuchten, 1980) was implemented in SUTRA to describe the soil water retention curve and relative hydraulic conductivity function. The van Genuchten parameters α [L −1 ] and n of the sand adopted in the experiments of Shen et al. (2020) were respectively 7.5 m −1 and 6. These values were obtained by fitting the van Genuchten-Mualem model to the soil moisture-capillary head measurements of Xin et al. (2018) (Figure S2 in Supporting Information S1). Consistent with Kuan et al. (2019), the van Genuchten parameters used in simulations were 5.9 m −1 and 2.68, which were obtained from Carsel and Parrish (1988). Other numerical parameters are listed in Table 1.
In addition to laboratory-scale models, unsaturated zone effects at the field scale were examined. The dimensions of field-scale models were 154 m (length) × 20 m (height) with a shoreline slope of 1:8 (Figure S1c in Supporting Information S1). The adopted slope and aquifer depth are similar to several published case studies of coastal settings (e.g., Heiss & Michael, 2014;Mao et al., 2006;Raubenheimer et al., 1999;Vos et al., 2020;Xiao et al., 2019;Zamrsky et al., 2018;Zhang et al., 2017). The aquifer was assumed to consist of fine sand (Rehovot sand; Mualem, 1976a) with n e = 0.4 and K s = 1.33 × 10 −4 m s −1 (Röper et al., 2013;Vandenbohede & Lebbe, 2006. The van Genuchten parameters α and n were 4.5 m −1 and 3.26, respectively, adopted from the study of Luo et al. (2019) (Figure S2 in Supporting Information S1). A sinusoidal tide with amplitude of 1 m and period of 12 hr was applied at the sea boundary. As with the laboratory-scale models, two types of inland boundary conditions were considered in the field-scale models. To compare unsaturated flow effects on salt distributions for head-controlled and flux-controlled aquifers, the magnitude of q f for the flux-controlled aquifer was determined as an average based on the head-controlled aquifer. That is, taking h l = 14.8 m for the constanthead boundary condition, the averaged inland flux over a tidal period was calculated to be 1.04 × 10 −5 m 2 s −1 at quasi steady state with unsaturated flow considered. Other scale-dependent parameters are given in Table 1.
The numerical mesh for each simulation is described in Table 1. The grid resolution satisfies the Péclet number stability criterion (4α L > ΔL, where α L [L] is the longitudinal dispersivity and ΔL [L] is the local distance between element sides along the flow direction; Voss & Souza, 1987). In all numerical simulations, the salinity was set to 0 and 35 ppt at the inland and sea boundaries, respectively. The corresponding densities were 1,000 and 1,025 kg m −3 for fresh groundwater and seawater, respectively. For laboratory-scale cases, the time step was 4 s for the experiments of Shen et al. (2020) and 1 s for those of Kuan et al. (2019). For field-scale cases, numerical simulations were first run by assuming a steady sea level with a time step of 120 s. Initially, the pressure distribution was hydrostatic, assuming the watertable was at mean sea level and the aquifer contained only freshwater. After reaching steady state, the pressure, velocity and salinity were used as the initial condition for tidal simulations, which was then run with a smaller time step of 30 s (Table 1) Shen et al. (2020). e Experiments conducted by Kuan et al. (2019) where the van Genuchten parameters are compiled from Carsel and Parrish (1988). f Both tidal and non-tidal cases are considered in the experiments. SUTRA cannot directly switch off unsaturated flow effects. To avoid uncertainties from adopting different numerical codes for variably saturated flow (i.e., SUTRA) and saturated-flow only (e.g., SEAWAT; Guo & Langevin, 2002), unsaturated flow effects were reduced in SUTRA by adjusting the van Genuchten parameter α. According to previous studies, α is approximately equal to the reciprocal of the capillary fringe height (van Genuchten & Nielsen, 1985), which is a key characteristic of unsaturated flow effects, particularly in relation to tidal propagation (Geng et al., 2017. That is, a larger α creates a smaller capillary fringe and weaker unsaturated flow effects. To reduce unsaturated flow effects, α was assigned a value of 55 m −1 (a capillary fringe height of around 0.018 m), while other parameters remained unchanged. For the three soils considered, the soil water content and corresponding hydraulic conductivity decrease to the residual soil water content and near-zero, respectively, over short distances when α = 55 m −1 ( Figure S2 in Supporting Information S1).
The ILU-preconditioned orthomin solver (Voss & Provost, 2010) was used in all numerical simulations, with convergence tolerances of 10 −12 kg m −1 s −2 and 10 −12 for solver iterations during pressure and transport solutions, respectively. All tidal simulations were run to quasi-steady-state conditions, whereby the position of the 17.5 ppt contour (50% of seawater), averaged over a tidal cycle, remained the same (within 0.5%) in consecutive tidal cycles. Moreover, testing of alternative mesh sizes and time step lengths was conducted to ensure convergent numerical results. The predicted 17.5 ppt contour location averaged over a tidal cycle was almost identical (i.e., location difference within 0.5%) when the mesh resolution was doubled and the time step lengths were halved relative to those listed in Table 1. For example, numerical tests for the experiments of Shen et al. (2020) and Kuan et al. (2019) are presented in Figure S3 in Supporting Information S1.

Comparison With Existing Laboratory Experiments
Figure 2 compares model results from the current study to the experimental results of Shen et al. (2020) for the situation of a constant-head inland boundary. As expected, there was no USP and only the saltwater wedge was apparent with a constant sea level (Figure 2a). Introducing tides caused a USP to form above the saltwater wedge, regardless of whether unsaturated flow was included or not ( Figure 2b). For the head-controlled aquifer with a constant sea level, the difference between the predicted results with and without unsaturated flow effects was almost indiscernible and both simulations slightly overestimated the measured saltwater wedge (Figure 2a). However, when a tidal fluctuation was added to the sea boundary, the simulated saltwater distributions with and without unsaturated flow effects showed more substantial deviations from the laboratory results ( Figure 2b). Specifically, neglecting unsaturated flow effects resulted in a more extensive USP and a smaller saltwater wedge. In combination, these led to FGD occurring at a lower shoreline elevation. This suggests that in comparison with the non-tidal case, tidal fluctuations intensify the role played by unsaturated flow in affecting salt distributions and the outflow face of FGD. Compared to the experimental results, both simulations with and without unsaturated zone effects overestimated the extent of the USP and underestimated the size of the saltwater wedge, although the numerical model with unsaturated zone effects was a closer match to both the USP and the saltwater wedge.
A USP was also apparent under tidal conditions for aquifers involving a constant-flux inland boundary whether or not unsaturated flow effects were included in the numerical model ( Figure 3). In contrast to the head-controlled results (Figure 2a), there was a noticeable difference between the non-tidal numerical results with and without unsaturated flow effects under flux-controlled conditions (Figure 3a). This is because unsaturated flow has a more significant impact on the hydraulic head of flux-controlled aquifers than that of head-controlled aquifers (Luo et al., 2022). The numerical result without unsaturated flow effects matched well with the non-tidal measurements, while the inclusion of unsaturated flow effects led to a slight overestimation of the saltwater wedge extent.
Consistent with the head-controlled case, the inclusion of tidal fluctuations in the flux-controlled case increased the deviation between the predicted results with and without unsaturated flow effects (Figure 3b). The numerical model including unsaturated flow effects slightly underestimated the saltwater wedge extent, while a close prediction of the USP was obtained. Ignoring unsaturated flow effects caused a substantial deviation from the tidal measurements, including a more extensive USP and more seaward saltwater wedge. Thus, unsaturated flow effects on salt distributions appear to be intensified when tidal fluctuations are considered regardless of the inland 10.1029/2022WR032931 6 of 18 boundary condition. The numerical model performed much better in predicting both the saltwater wedge and USP when taking unsaturated flow into account, consistent with the head-controlled results.
To decrease the uncertainty from assuming a large α to reduce unsaturated flow, the experiments of Shen et al. (2020) and Kuan et al. (2019) were further simulated by SEAWAT with the parameters listed in Table 1. As shown in Figures S4 and S5 in Supporting Information S1, the results of SEAWAT were consistent with those of SUTRA with α = 55 m −1 for both experiments, both with and without tides. This confirms that setting α to a large value is a valid approach for ignoring unsaturated flow.
The above results indicate that accounting for unsaturated flow effects significantly improved the numerical model performance for both head-controlled and flux-controlled experiments when tidal fluctuations were included. Nevertheless, there were some deviations between the numerical and experimental results, even when the unsaturated zone was represented. At least in part, this is likely due to unavoidable experimental uncertainties. For example, the van Genuchten parameters used in the experiments of Kuan et al. (2019) were compiled from Carsel and Parrish (1988), instead of values specific to the sand used in those experiments. There may also be unsaturated zone hysteresis effects in the experiments that were ignored in the numerical model (e.g., Werner & Lockington, 2003). Moreover, Zheng et al. (2022) pointed out that the original Richards' equation may not predict watertable dynamics accurately when the unsaturated zone is heavily truncated during watertable fluctuations due to the capillary fringe encountering the land surface. However, this is the case in the sloping beach areas for both experiments of Shen et al. (2020) and Kuan et al. (2019) with tides. Thus, truncation of the unsaturated zone in the sloping beach area may have introduced errors in predicting the watertable dynamics and salt distributions with SUTRA based on the original Richards' equation. Overall, despite differences between the numerical and experimental results, the numerical model with unsaturated flow effects allows a better replication of the laboratory results.

Sensitivity Analysis of Unsaturated Flow Effects on Salt Distributions: Laboratory Scale
The effects of unsaturated flow on salt distributions in both head-controlled and flux-controlled aquifers were explored through sensitivity analysis based on simulations using different values of the van Genuchten parameters α and n. As mentioned earlier, these parameters were not measured in the experiments by Kuan et al. (2019). Therefore, numerical models were established based on the experimental setup of Shen et al. (2020). For the head-controlled aquifer, α values of 5.5, 7.5 (measured), 11,15,22,35,45, and 55 m −1 and n values of 3, 6 (measured) and 10 were evaluated, while other parameters remained unchanged relative to those of the respective laboratory experiments. Note that when conducting the sensitivity analysis for α, n is kept at the measured value, and vice versa. In addition to varying α and n, the constant head assigned to the inland boundary of the experiment was replaced by an equivalent constant flux to assess sensitivities under flux-controlled conditions. As done above, the flux assigned to the inland boundary was determined from the corresponding head-controlled simulation. For Case L1 with tides, the averaged inland flux over a tidal period was calculated to be 9.95 × 10 −7 m 2 s −1 at quasi steady state with unsaturated flow considered. This flux was assigned to q f for the constant-flux boundary condition. Note that phase-averaged results (averaged over a tidal cycle) are presented for all numerical simulations with tidal fluctuations hereafter. Figure 4 and Figure S6 in Supporting Information S1 show the sensitivity of the salt distribution and saturation distribution to α for the tidally influenced head-controlled aquifer, respectively. Salt distributions differed depending on the value of α, highlighting the importance of the capillary fringe thickness, and unsaturated flow more generally. Figure 5a represents the corresponding USP area (calculated as the area between the 17.5 ppt contour and BE in Figure 1) and the saltwater wedge toe location (the distance from EF in Figure 1 to the most landward 17.5 ppt contour at the base of the aquifer). An increase of α from 5.5 to 55 m −1 led to larger USP areas from 0.11 to 0.35 m 2 , while the saltwater wedge toe retreated from 1.45 to 1.00 m (a reduction of 31%) (Figure 5a). The USP expansion and saltwater wedge retreat caused the FGD outflow face to shift further below the low-tide mark about 0.1 m. The USP and saltwater wedge changed substantially when α increased from 5.5 to 22 m −1 , while further increases in α led to only mild changes to the salinity distribution, reflecting nonlinear interactions between salt distributions and unsaturated flow (Figure 5a). The opposite responses of the USP and the saltwater wedge to reductions in unsaturated zone effects is due to the more extensive USP pushing the saltwater wedge seaward to maintain the outward flow of FGD. Nguyen et al. (2020) observed the same interactions between the USP and the saltwater wedge in response to freshwater-saltwater temperature contrasts. In their study, movement of the saltwater wedge caused a reduction in the USP whereas the opposite appears to have occurred in the current analysis of unsaturated zone effects. The corresponding numerical results with a constant sea level but with varying α are presented in Figure S7 in Supporting Information S1. A comparison of Figure 4 to Figure S7 in Supporting Information S1 indicates that unsaturated flow has a greater effect on salt distributions when tidal fluctuations are considered in head-controlled aquifers.
Consistent with the head-controlled aquifer, changing α significantly affected salt distributions in flux-controlled aquifers (Figure 6), as well as the saturation distribution ( Figure S8 in Supporting Information S1). An increase of α from 5.5 to 55 m −1 induced an expansion of the USP with the area increasing from 0.14 to 0.40 m 2 , along with a reduction in the size of the saltwater wedge with its toe location decreasing from 1.11 to 0.69 m (Figure 5b). The corresponding numerical results with a steady sea level showed that removing unsaturated flow effects led to a seaward retreat of the saltwater wedge ( Figure S9 in Supporting Information S1). This is consistent with the study of Luo et al. (2022), which showed that the unsaturated zone transmits part of the freshwater flux toward the sea, leading to a lower watertable and hence a more seaward saltwater wedge when the unsaturated zone is included. Again, unsaturated flow noticeably affected salt distributions in flux-controlled aquifers subject to tidal fluctuations ( Figure 6 and Figure S9 in Supporting Information S1). The salt distribution changed only slightly, relative to those resulting from changes to α, when varying n from 3 to 10 for both the head-controlled and flux-controlled aquifers ( Figure S10 in Supporting Information S1). This indicates that α is the key parameter influencing the effects of unsaturated flow on coastal flow patterns. In addition, the sensitivity analysis of α and n confirms that unsaturated flow effects are largely removed and hence can be treated as the case without unsaturated flow when setting α = 55 m −1 for the sand adopted in Case L1.
To further explain why neglecting unsaturated flow induces an expansion of the USP and seaward retreat of the saltwater wedge, hydraulic heads over one tidal period at x = 6.6 m and z = 0.7 m (indicated by the white circle in Figure 2b) were extracted, as shown in Figure 7. Note that only simulated results with α = 7.5 and 55 m −1 are presented, corresponding to the earlier numerical cases with and without unsaturated flow effects. As can be seen, unsaturated flow has a major influence on the predicted hydraulic head time series of the observed point over one tidal period for both head-controlled and flux-controlled aquifers. During the flood period (0-60 s), seawater infiltrated into the head-controlled aquifer because of the head difference between the observed point and the tidal forcing (Figure 7a). By comparison, ignoring unsaturated flow effects led to lower hydraulic heads and hence greater landward hydraulic head gradients during high tide. This is consistent with previous findings that removing unsaturated zone effects enhances tidal attenuation (e.g., Barry et al., 1996;Kong et al., 2013Kong et al., , 2015. Consequently, a larger amount of seawater infiltrated into the aquifer and the USP expanded. The flow direction gradually reversed during the ebb period, causing groundwater discharge from the aquifer. In contrast to the flood period, the predicted hydraulic heads were generally higher in the absence of unsaturated flow effects during the ebb period, especially for the flux-controlled case (Figure 7b). It follows that the inclusion of unsaturated flow caused more intensive drainage of the aquifer (and therefore a larger head decline), despite the larger head difference between the observed point and the sea without unsaturated flow effects. This is consistent with findings of Kong et al. (2016), who concluded that unsaturated flow accelerates drainage of hillslope unconfined aquifers. Counterintuitively, ignoring unsaturated flow effects leads to greater average shoreline head over a tidal period, which shifts the saltwater wedge seaward.   Both the flux-controlled and head-controlled cases show the above general trends, although differences between cases with and without unsaturated flow are accentuated in the flux-controlled case (Figure 7b).

Upscaling From Laboratory to Field Scale
The laboratory-scale simulation was upscaled to evaluate the effects of unsaturated flow in tidally influenced field-scale aquifers ( Figure S1c in Supporting Information S1, Table 1). Figure S2 in Supporting Information S1 shows the soil water retention curve adopted for field-scale numerical simulations. Compared with the experiments of Shen et al. (2020) and Kuan et al. (2019), the soil in the field-scale simulations has stronger capillarity effects due to a smaller α.
There are clear differences between salt distributions predicted with and without unsaturated flow effects for the field-scale, head-controlled aquifer ( Figure 8). As α changes from 4.5 to 55 m −1 , the USP area increases from 264 to 373 m 2 (an increase of 41.3%), whereas the saltwater wedge toe retreats 4.1%, from 70.3 to 67.4 m. As already noted, this is because the head in the nearshore area is lower during the high-tide period without unsaturated flow. The effect of the unsaturated zone on the salt distribution is weaker in the field-scale case compared to those of the laboratory-scale models (Figures 2b and 8c). This arises because the ratio of the capillary fringe height to the aquifer thickness is larger in the laboratory-scale simulation. Thus, the laboratory-scale cases have greater relative changes in the saturated thickness when the unsaturated zone is excluded.
For the flux-controlled, field-scale aquifer (Figure 9), a similar trend in salt distributions to that of the head-controlled aquifer was obtained. As α increased from 4.5 to 55 m −1 , the USP area increased 23.8% (from 260 to 322 m 2 ), while the saltwater wedge toe retreated 3.0% (from 70.3 to 68.2 m). These results highlight that unsaturated flow is also important when predicting salt distributions in tidally influenced aquifers at the field scale.
Following previous studies (e.g., Robinson, Li, & Barry, 2007;Xin et al., 2010Xin et al., , 2011, flow paths and associated transit times within the aquifer were further calculated based on the mean velocity field (white lines with a number shown in Figures 8 and 9). Regardless of whether head-controlled or flux-controlled conditions were assigned to the inland boundary, there were noticeable differences between the results calculated with and without considering unsaturated flow effects. In both cases, unsaturated flow effects caused USP shrinkage and saltwater wedge extension, and hence transit times of particles entering the aquifer through the sea boundary increased within the saltwater wedge and decreased within the USP. However, transit times of particles within the fresh groundwater regime differed between head-controlled and flux-controlled cases when unsaturated flow was included. In the former, transit times were substantially reduced as a result of USP shrinkage and the associated shorter flow paths (Figure 8d), whereas the transit time slightly increased for the flux-controlled case. This is because the transit time is the balance of the trajectory length and flow velocity. In the flux-controlled case, part of the horizontal flux occurs within the unsaturated zone leading to reduced velocities and hence the longer transit time despite the shorter flow path length (Figure 9d).
In terms of the water and salt exchange across the aquifer-ocean interface, there was still an obvious difference between the predictions with and without unsaturated flow for both the head-controlled and flux-controlled aquifers ( Figure 10). Neglecting unsaturated flow effects leads to a wider variation of both water and salt fluxes and overestimation of the peak of water and salt fluxes. For the head-controlled aquifer, the peak water influx increased by 25.3% from 1.98 × 10 −3 kg m −1 s −1 to 2.48 × 10 −3 kg m −1 s −1 , whereas the corresponding salt flux increased by 24.4% from 6.94 × 10 −5 kg m −1 s −1 to 8.63 × 10 −5 kg m −1 s −1 (Figures 10a and 10c). For the flux-controlled aquifer, the peak water influx increased by 22.2% from 1.98 × 10 −3 kg m −1 s −1 to 2.42 × 10 −3 kg m −1 s −1 , whereas the corresponding salt flux increased by 22.8% from 6.92 × 10 −5 kg m −1 s −1 to 8.50 × 10 −5 kg m −1 s −1 (Figures 10b and 10d). Thus, neglecting unsaturated flow effects leads to expansion of the USP because of increased salt infiltration.

Sensitivity Analysis: Field Scale
Numerical simulations were further conducted based on Cases F1 and F2 (Table 1) by varying K s (0.83 × 10 −4 , 1.33 × 10 −4 , 1.83 × 10 −4 , 2.33 × 10 −4 m/s), h l (14.8, 15.0, 15.2, 15.4 m) and q f (1.04 × 10 −5 , 1.25 × 10 −5 , 1.50 × 10 −5 , 1.75 × 10 −5 m 2 s −1 ), which vary significantly in practice. For both the head-controlled and flux-controlled aquifers, significant differences between the simulated results with and without unsaturated flow effects were obtained for the range of K s values adopted ( Figure 11). This indicates that unsaturated flow is important in predicting salt distributions in a wide variety of aquifer systems. For the head-controlled aquifer, the USP area decreased and the saltwater wedge toe location increased with an increase of K s , regardless of whether unsaturated flow effects were considered (Figure 11a). The opposite occurred in the flux-controlled aquifer cases,  in which the USP increased in extent and the saltwater wedge toe location shifted seaward with increasing K s . This is a consequence of complex relationships between the USP size, freshwater and saltwater fluxes, and the aquifer properties. For example, increasing the saturated hydraulic conductivity enhances seawater infiltration during the rising tide for both the head-controlled and flux-controlled aquifers. However, the fresh groundwater discharge increases only in the head-controlled aquifer, while fresh groundwater discharge is independent of K s in the flux-controlled aquifer. By comparison, unsaturated flow effects on salt distributions (the difference between solid and dashed lines for the same color) increase at small saturated hydraulic conductivity values for the head-controlled aquifer, and at large saturated hydraulic conductivity values for the flux-controlled aquifer.
When varying h l for the head-controlled aquifer or q f for the flux-controlled aquifer, the USP area and the saltwater wedge toe location predicted with unsaturated flow greatly differed from those without unsaturated flow ( Figure 12). For the head-controlled aquifer, the USP area decreased and the saltwater wedge toe location increased with increasing h l , regardless of whether unsaturated flow effects were included (Figure 12a). Increasing h l leads to a greater fresh groundwater discharge, which inhibits the USP development. By comparison, the deviation between the results with and without unsaturated flow effects decreased with increasing h l , suggesting weakened unsaturated flow effects for a larger inland head. This is because of a decrease in the ratio of the capillary fringe height to  the aquifer depth. In addition, increasing q f slightly decreased the deviation between the results with and without unsaturated flow effects for the flux-controlled aquifer, which is consistent with increasing h l for the head-controlled aquifer (Figure 12b).
Increasing the inland head or inland flux weakens unsaturated flow effects on salt distributions, which implies that an increase in the saturated aquifer thickness weakens unsaturated flow effects on salt distributions. However, this does not mean unsaturated flow effects are unimportant for aquifers with a large thickness. According to Robinson, Li, and Prommer (2007), USPs are expected at sites with a large tidal range and beach slope but intermediate aquifer thickness, inland hydraulic gradient and aquifer dispersivity. This means that the USP can appear in aquifers with a large thickness by varying other factors. As long as there is an USP, unsaturated flow effects are expected to play an important role in affecting salt distributions.
It is clear from both laboratory-scale and field-scale simulations that unsaturated flow affects salt distributions, as well as aquifer-ocean material exchange in coastal unconfined aquifers. Previous studies noted that the USP is a hotspot for biogeochemical reactions because of intensive fresh groundwater and seawater mixing (e.g., Anwar et al., 2014;Heiss et al., 2017;Kim et al., 2020). These two water bodies generally contain different components (e.g., dissolved oxygen, dissolved organic carbon and nutrients) and thus their mixing would alter biogeochemical reactions that are either beneficial or deleterious to the marine ecosystem. Since including unsaturated flow effects leads to a shrinkage of the USP, it is expected that considering unsaturated flow in coupled groundwater flow and reactive transport models with variable density will impact biochemical reactions.
In addition, the FGD is important for land-sourced pollutant loading to the ocean (Luijendijk et al., 2020;Moore, 2010;Robison et al., 2018;Slomp & Van Cappellen, 2004). For instance, nitrogen fertilizer that cannot be absorbed by crops is leached to the unconfined aquifer in the form of nitrate (Min et al., 2022). Consideration of unsaturated flow leads to changes in the transit time of inland particles and material exchange across the aquifer-ocean interface, which could further affect land-sourced nitrate loading to the ocean and hence coastal water quality. The effects of unsaturated flow on reactive pollutant transport and transformation require further investigation.

Conclusions
Combining numerical simulations with laboratory experiments, this study examined the effects of unsaturated flow on salt distributions in tidally influenced coastal unconfined aquifers with constant-head or constant-flux inland boundary conditions. The results lead to the following conclusions: • Compared with the aquifer subject to a constant sea level, unsaturated flow affects salt distributions more significantly in tidally influenced coastal unconfined aquifers, regardless of the inland boundary condition. As a result, an incorrect representation of unsaturated flow will cause errors in estimates of tide-affected seawater intrusion, especially if near-shore conditions are a focus of the modeling effort. • Neglecting unsaturated flow effects leads to an expansion of the USP and shrinkage of the saltwater wedge. This is because of the lower head in the nearshore area during high-tide periods, which facilitates seawater infiltration into the aquifer and hence formation a larger USP. Counterintuitively, ignoring unsaturated flow effects leads to greater average shoreline head over a tidal period, which shifts the saltwater wedge seaward. • With unsaturated flow, USP shrinkage and saltwater wedge extension causes the transit time of particles from the sea boundary to increase within the saltwater wedge and decrease within the USP. Meanwhile, water and salt exchange across the aquifer-ocean interface decrease. • The effect of the unsaturated zone on salt distributions increases with the ratio of the capillary fringe height to the aquifer thickness. • Unsaturated flow effects on salt distributions increase with decreasing/increasing saturated hydraulic conductivity for the head-controlled/flux-controlled aquifer. • Unsaturated flow effects on salt distributions reduce with increases in both the inland head and the inland flux.
The aquifer considered here was assumed to be isotropic and homogeneous, whereas aquifers are usually anisotropic and heterogeneous. This affects the extent of seawater within the aquifer (Geng & Michael, 2021). In addition, the sea boundary signal is usually a combination of tides and waves, instead of a single-frequency tide (Xin et al., 2010). Although uncertainties due to these factors could provide further insights, this study shows that unsaturated flow significantly affects salt distributions in coastal unconfined aquifers.

Data Availability Statement
The experimental data used to produce Figures 2 and 3 are compiled from Shen et al. (2020) and Kuan et al. (2019), respectively. Simulation input files can be found at https://zenodo.org/record/7480816#.Y6el5nbP1Q9. Shumei Zhu and Huiqiang Wu for helping with SEAWAT simulations. The authors greatly appreciate constructive comments from the Editor, associate Editor and three reviewers (L. Stoeckl and two anonymous reviewers), which led to significant improvement of the paper. Open access funding provided by École Polytechnique Fédérale de Lausanne (EPFL).