Effects of Multiconstituent Tides on a Subterranean Estuary With Fixed-Head Inland Boundary

While tides of multiple constituents are common in coastal areas, their effects on submarine groundwater discharge (SGD) and salinity distributions in unconfined coastal aquifers are rarely examined, with the exception of a recent study that explored such effects on unconfined aquifers with fixed inland freshwater input. For a large proportion of the global coastline, the inland areas of coastal aquifers are topography-limited and controlled by constant heads. Based on numerical simulations, this article examines the variation of SGD and salinity distributions in coastal unconfined aquifers with fixed-head inland boundaries at different distances from the shoreline (i.e., 50, 100, 150, and 200 m). The results showed that the fluctuation intensity of freshwater input was enhanced as the inland aquifer extent decreased, e.g., the range of tide-induced fluctuations in freshwater input increased by around 5 times as the inland aquifer extent decreased from 200 to 50 m. The frequency spectra of the fluctuations of SGD and salinity distributions showed that the coastal aquifer of a shorter inland aquifer extent smoothed out fewer high-frequency tidal constituents but enhanced interaction among different tidal constituents. The interaction among tidal constituents generated new low-frequency signals in the freshwater input and salinity distributions. Regressions based on functional data analysis demonstrated that the inland freshwater input and salinity distributions at any given moment were related to the antecedent (previous) tidal conditions weighted using the probability density function of the Gamma distribution. The influence of the antecedent tidal conditions depended on the inland aquifer extent.


INTRODUCTION
Ocean tides are identified as an important driving force that controls the hydrodynamic and biogeochemical processes in coastal aquifers (Werner et al., 2013;Robinson et al., 2018). In subterranean estuaries (i.e., mixing zones between terrestrially derived fresh groundwater and recirculating seawater in coastal aquifers), tides significantly affect submarine groundwater discharge (SGD) and salinity distributions (Werner et al., 2013;Robinson et al., 2018). Tidal forcing, in combination with the density contrast between terrestrial fresh groundwater and nearshore seawater, can lead to the formation of two saltwater recirculation zones in unconfined aquifers: a lower saltwater wedge (SW) and an upper saline plume (USP) (Boufadel, 2000;Robinson et al., 2006).
Recently, Yu et al. (2019b) numerically examined the effect of multiconstituent tides on both SGD and salinity distributions in a subterranean estuary. They found that the response of water exchange (e.g., SGD) across the aquifer-sea interface to tides is almost instantaneous. As such, the present (at any given moment) water exchange is mainly influenced by the present (i.e., at the same time) tidal conditions. However, the response of salinity distributions to tides is delayed; thus, salinity distributions are dependent on the antecedent (previous) tidal conditions. Yu et al. (2019b) found that coastal aquifers behave as a lowpass filter smoothing-out tidal signals, particularly highfrequency tidal constituents (e.g., semidiurnal and diurnal tidal signals). Furthermore, the combination of high-frequency tidal constituents generates new low-frequency signals when the tidal signals propagate in the aquifer, as determined by Li et al. (2000) using analytical solutions for groundwater tidal propagation. These generated low-frequency signals, combined with the primary tidal fluctuations created directly by the ocean tide, jointly affect the salinity distributions (e.g., SW and USP) in unconfined aquifers (Yu et al., 2019b). With the antecedent tidal conditions weighted in the form of convolution using the probability density function of the Gamma distribution, Yu et al. (2019b) demonstrated that the prolonged and cumulative effects of antecedent tidal conditions on the present salinity distributions in unconfined aquifers could be quantified. Yu et al. (2019b) studied multiconstituent tides within a coastal unconfined aquifer with a fixed-flux inland boundary and inland extent of 150 m. Both fixed-flux (Michael et al., 2005;Robinson et al., 2007c;Kuan et al., 2012;Evans and Wilson, 2016) and fixed-head (Li et al., 1997;Ataie-Ashtiani et al., 1999;Heiss and Michael, 2014;Geng and Boufadel, 2015) inland boundaries have been widely adopted in previous studies of SGD and salinity distributions. The former assumes that freshwater discharge to the sea is constant, whereas the latter applies to topographycontrolled aquifers (Werner and Simmons, 2009;Michael et al., 2013). In low-elevation coastal areas, the watertable is often shallow, leading to the overlying unsaturated zone being thin. Under this condition, the rainfall-induced infiltration would be inhibited, i.e., once the soil is saturated, surface water runoff would occur and the watertable tends to be fixed to the local surface elevation. This likely forms a fixed-head inland boundary (Michael et al., 2013;Sawyer et al., 2016). Michael et al. (2013) suggested that approximately 70% of the world's coastlines are topography-controlled. Coastal aquifers with fixed-head landward boundaries are recognized to be more vulnerable to seawater intrusion, as salinity distributions are more responsive to changes to ocean boundary conditions (Werner and Simmons, 2009;Michael et al., 2013). Despite the commonality of fixedhead inland boundaries, the question remains open as to how they regulate the SGD and salinity distributions in subterranean estuaries subjected to multiconstituent tides.
In this study, we focused on the effects of multiconstituent tides on SGD and salinity distributions in unconfined aquifers with fixed-head inland boundaries at different landward positions, under the assumption that the yearly averaged inland freshwater input is the same in aquifers of alternative landward extent. The following key questions were explored based on numerical simulations: 1) How does a change in the inland aquifer extent affect the SGD and salinity distribution responses to multiconstituent tides? 2) How do multiconstituent tides affect inland freshwater input in inland fixed-head aquifers? 3) How does the inland aquifer extent influence the role of antecedent tidal conditions on salinity distributions?

Conceptual Model and Simulations
Similar to Yu et al. (2019b), a 2D cross-shore section of a coastal unconfined aquifer was considered in this study, as shown in Figure 1. Boundary DEF was set as the time-variant head boundary (aquifer-sea interface) affected by multiconstituent tides. To investigate the response of the aquifer to multiconstituent tides, we used six years (2011-2016) of real tide data measured from Brisbane Bar tidal station (website: http://data.qld.gov.au/dataset/brisbane-bar-tide-gauge-archivedinterval-recordings). The measurement interval was 600 s. Four dominant tidal constituents were detected in this dataset based on Boundaries BC and CD were set to no-flow boundaries. Rainfall and evaporation were excluded and thus AF was also set as a no-flow boundary. Different from Yu et al. (2019b), the landward boundary AB was set to a fixed-head boundary. The fresh groundwater recharged from the inland boundary AB and discharged into the ocean through the aquifer-sea interface DEF. It is worth noting that to make the simulation results comparable, inland watertable elevations (H f ) were adjusted to get a similar yearly averaged freshwater input (around 2.1 m 3 /m/d). We focused on four cases with the fixed-head inland boundary located at different locations as follows: (1) Fixed-head and L I 50 m (L I is the inland aquifer extent, which measures the horizotal distance from the inland boundary to the shoreline), with the seaward aquifer extent L S 50 m, the seaward aquifer thickness H S 27 m, and the inland aquifer thickness H I 33 m (note that the same L S , H S , and H I were set in all the cases). In addition to these four cases with fixed-head inland boundaries, one case with fixed-flux inland boundary (L I 150 m, similar to that used in Yu et al. (2019b)) was also considered for comparing the behavior of aquifers with different inland boundaries. We tested the effect of inland aquifer extent (L I ) on salinity distributions in the aquifers with fixed-flux inland boundaries. The results showed that the inland aquifer extent did not remarkably affect the variation of the USP and SW in the aquifers. Therefore, only the results from one case with fixed-flux inland boundary (fixed-flux and L I 150 m) were included in this article. The aquifer properties were similar to those used in Yu et al. (2019b). That is, the model setup was based on a homogeneous sandy beach with hydraulic conductivity K s 10 m/d and porosity ϕ 0.45. The beach slope was 0.1. The longitudinal and transverse dispersivities were set to 0.5 and 0.05 m, respectively (Hunt et al., 2011). The residual water saturation S Wres was set to 0.1 and water retention parameters α and n were set to 14.5 m −1 and 2.68 (typical values for sand (Carsel and Parrish, 1988)) for the van Genuchten (1980) equations.
The numerical simulations were conducted using SUTRA (Voss and Provost, 2010), which simulates variably saturated and density-dependent pore-water flow and associated solute transport. In the inland subdomain, a mesh with Δx 1.5 m and Δz 0.5 m was used, while in the intertidal and seabed subdomain, the mesh was refined (Δx ≈ 0.33 m and Δz ≈ 0.1 m). For all the simulations, the model was initially run to the steady states (pressure and salinity unchanged) with a static sea level (0 m). The measured tidal conditions (2011)(2012)(2013)(2014)(2015)(2016) were then introduced as the boundary conditions in the model with a time step of 300 s. Here, we focused on the analysis of the last two years (2015-2016). As suggested by Yu et al. (2019b) and will be discussed later, after 4 years (2011)(2012)(2013)(2014), the behavior of the system was independent of the initial conditions subjected to the static sea level.

SGD and Salinity Distribution Fluctuations Induced by the Multiconstituent Tides
As mentioned earlier, for the five considered cases (four with fixed-head inland boundaries and one with fixed-flux inland boundary), the averaged inland freshwater input was similar. The two-year (2015-2016) averaged salinity distributions in the aquifers were similar for the five cases, e.g., the 50% isohaline of the seawater (17.5 ppt) largely overlapped ( Figure 2). This suggests that the yearly averaged inland freshwater input FIGURE 1 | Conceptual diagram of a subterranean estuary of an unconfined aquifer. Flow processes considered in this article are as follows: 1) density-driven seawater circulation, 2) tide-induced seawater circulation, and 3) terrestrial groundwater discharge. The colors represent the salinity (yellow for freshwater and red for saltwater). Boundary ABCDEF represents the model domain. The x-z coordinate origin was set at the mean shoreline. L I and L S are the inland and seaward aquifer extents, respectively. H I and H S are the inland and seaward aquifer thicknesses, respectively. H f is inland watertable elevation at the inland boundary.
Frontiers in Environmental Science | www.frontiersin.org January 2021 | Volume 8 | Article 599041 rather than the inland aquifer extent significantly affected the salinity distributions at a long-time scale. For the four cases with fixed-head inland boundaries, the calculated yearly averaged amounts of salt stored in (per unit width) saltwater wedge (SM SW ) and upper saline plume (SM USP ) were also close to those for the case with fixed-flux inland boundary (8,425 and 872 kg/m, respectively) ( Table 1). The averaged tide-induced (around 1.4 m 3 /m/d) and density-driven (around 1 m 3 /m/d) seawater circulation fluxes were also similar for all the cases ( Table 1). Since the yearly averaged inland freshwater input was similar, the inland aquifer extent did not remarkably affect the averaged SGD and salinity distributions in the aquifers subjected to the multiconstituent tides. We further investigated the variation of fluxes and salinity distributions (both water flux and amount of salt stored) in response to the multiconstituent tides by presenting the variation in the last two years (2015)(2016). To quantify the fluctuations for the different cases, the averaged value (μ) and the standard deviation (σ) over each principal semidiurnal lunar tidal cycle (12.42 h) in these two years were calculated.
The calculated standard deviations for the tide-induced and density-driven influx were similar among the four cases with fixed-head inland boundaries and differed slightly from that for the case with fixed-flux inland boundary ( Table 1). The SGD, as the total efflux across the aquifer-sea interface, varied intensively in response to the multiconstituent tides ( Figure 3). The averaged values of the SGD over each principal semidiurnal lunar tidal cycle for the five cases largely overlapped. Consistently, the calculated standard deviations were quite similar and around 0.62 m 3 /m/d. This suggests that, for the simulated cases with similar averaged freshwater input, SGD was predominantly affected by the present tidal conditions. The response of SGD to the tidal fluctuations was mostly instantaneous; thus, the effect of the inland aquifer extent was not obvious.
In contrast, the fluctuations of freshwater input differed case by case ( Figure 4). As the inland aquifer extent decreased,   the averaged values of the freshwater input over each principal semidiurnal lunar tidal cycle started to fluctuate around that for the case with fixed-flux inland boundary (2.1 m 3 /m/d). As L I reduced from 200 to 150, 100, and 50 m, the standard deviation increased from 0.08 to 0.12, 0.21, and 0.47 m 3 /m/ d, respectively. It is worth noting that the equal reduction in L I led to unequal increase in the standard deviation of the freshwater input, leading to an accelerating trend. The inland freshwater input was controlled by the lateral hydraulic gradient. As setup in the simulations, a similar yearly averaged inland freshwater input was produced ( Table 1). It is expected that the instantaneous hydraulic gradient was continuously modified by the inland watertable and seaward tidal level. The multiple tidal constituents propagated into the aquifer and thus affected the hydraulic gradient in a complex fashion. In particular, the low-frequency signals would propagate further inland and obviously affect the freshwater input (details in Analysis on Amplitude Spectrums).
For both the cases with fixed-head and fixed-flux inland boundaries, SM SW and SM USP fluctuated remarkably ( Figures  5, 6). Overall, the cases with fixed-head inland boundaries were related to intensive fluctuations in the averaged value over each principal semidiurnal lunar tidal cycle ( Table 1). As L I reduced from 200 to 150, 100, and 50 m, the standard deviation for SM SW increased from 162 to 167, 180, and 218 kg/m, respectively, whereas the standard deviation for SM USP increased from 90 to 95, 107, and 154 kg/m, respectively. Clearly, the fluctuation intensity for SM SW was more obvious than that for SM USP , suggesting that SM SW tended to be more affected by the lowfrequency tidal signals (details in Analysis on Amplitude Spectrums).

Analysis of Amplitude Spectra
The variations of the SGD and salinity distributions could be reflected by amplitude spectra. Similar to Yu et al. (2019b), the instantaneous results were selected with a 2 h time interval (X). We also used the following standardized and nondimensional variables (X*). . (1) As expected, four dominant frequencies (M 2 , S 2 , K 1 , and N 2 ) appeared in the tidal series ( Figure 7A) and also in the SGD variations ( Figure 7B). As pointed out by Yu et al. (2019b), in the SGD, two new signals with frequencies of 7.83 × 10 -7 Hz (spring-neap tidal cycle, 14.8-day period) and 4.24 × 10 -7 Hz (27.3-day period) were generated by the combinations of M 2 and S 2 tidal constituents (f (S 2 )−f (M 2 )) and M 2 and N 2 tidal constituents (f (M 2 )−f (N 2 )), respectively. These two new signals also appeared in the freshwater input, SM SW and SM USP (Figures 7C-E).
Different from the tidal and SGD series, the low-frequency signals started to play important roles in the freshwater input, SM SW and SM USP (Figures 7C-E). The contributions (defined as the ratio of a specific nondimensional amplitude to the sum of all the nondimensional amplitudes) of these signals changed in different ways as L I reduced. For example, for the freshwater input, the contribution of the signal of the 365-day period decreased from 7.21% to 6.10%, 4.51%, and 1.83% as L I reduced from 200 to 150, 100, and 50 m, respectively. The contribution of the signal of the 365-day period kept unchanged for SM SW , whereas for SM USP , it increased from 2.44% to 3.43%, 5.35%, and 9.33% as L I reduced from 200 to 150, 100, and 50 m, respectively.  As mentioned earlier, the low-frequency signals likely propagated further landward and thus interacted more with the inland boundary. The boundary reflections would damp these signals. The extent of damping depended on the inland aquifer extent and the tidal frequencies. With respect to the variation of the freshwater input, the low-frequency signal of the 365-day period was not apparent with L I 50 m, where the highfrequency signals (M 2 , S 2 , K 1 , and N 2 ) were relatively strong. Thus, the fluctuation intensity of the freshwater input was significantly enhanced in the case with L I 50 m.
Similar to Yu et al. (2019b) considering the fixed-flux inland boundary, the combination of M 2 , S 2 , and N 2 tidal constituents (3f (M 2 ) -2f (N 2 )f (S 2 )) generated a new signal with a frequency of 6.50 × 10 -8 Hz (178.1-day period) in SM USP ( Figure 7E). However, for the cases with fixed-head inland boundaries, this signal of the 178.1-day period competed with that of the 365-day period, both affecting SM USP ( Figure 7E). The contribution of the former decreased, whereas that of the latter increased as L I decreased. This was significantly different from the case with fixed-flux inland boundary, in which no signal of the 365-day period was detected for SM USP .

Quantification of SGD and Salinity Distributions Affected by Multiconstituent Tides
From the spectrum analysis, it can be seen that the variations of the SGD and salinity distributions were more complex for the cases with fixed-head inland boundaries, in comparison with the case with fixed-flux inland boundary. The question of whether these variations can be quantified remained interesting, and the regression modeling approaches in Yu et al. (2019b) were tested.
Similar to the results of Yu et al. (2019b), the response of SGD to the multiconstituent tides was instantaneous. As such, the present SGD could be quantified by a linear regression model without considering the effect of the antecedent tidal conditions.
where Y is the output of the regression model. X σ and X M are the standard deviation of tidal level and the averaged tidal level over each semidiurnal lunar tidal cycle (12.42 h), respectively, which denote the extent of tidal fluctuations over the tidal cycle and the averaged sea level. σ and M indicate the standard deviation of   S2). With increased inland aquifer extent, the weight of the averaged tidal level (a M ) in the regression model decreased and the performance of the regression model slightly declined. As expected, the response of the freshwater input, SM SW and SM USP , to the tides was delayed, which failed the regression model in Eq. 2 excluding the antecedent tidal conditions. Similar to Yu et al. (2019b), we defined the cumulative effect of antecedent multiconstituent tides (X σ and X M ) in the regression models based on functional data analysis (Ramsay and Silverman, 2005). The cumulative effect was weighted in the form of the following convolution: where DTL i is the parametric regressor, i.e., the weighted antecedent tidal conditions combined in a cumulative fashion. The subscript i indicates the independent value. t is the present time and t − jΔt is the past time associated with the quantified antecedent tidal condition. The increment Δt of each time step was set to 12.42 h. X i,t−jΔt is the variable at that past time. The subscript j indicates the past time period considered in the convolution, with the minimum and maximum values of n and m, respectively (set to 1 and 1,411 (two years), respectively). ζ i, j is a time-dependent weighting factor described by the following probability density function of the Gamma distribution: where α i and β i are the shape and scale factors, respectively. The ratio α i /β i is the mean of the Gamma distribution, which indicates a tail of the distribution (a large α i /β i indicates a long tail). In addition, it reflects the overall contribution of the antecedent tidal conditions. The freshwater inputs, SM SW and SM USP , were then quantified using a linear combination of DTL σ and DTL M . Firstly, we tested the freshwater input using the regression model as follows: For the cases with fixed-head inland boundaries, this regression model well captured the variation of the freshwater input for the cases with L I 200, 150, and 100 m ( Table 3 and Supplementary Figures S3, S4; adjusted R 2 0.99, 0.99, 0.96, respectively). However, the regression result (adjusted R 2 0.65) for the case with L I 50 m was not as good as the other cases with longer inland aquifer extents. The freshwater input was affected by the lateral hydraulic gradient, which was a combination of the inland watertable and tidal level. Given a short inland aquifer extent, the lateral hydraulic gradient was more responsive to the tidal fluctuation as it took less time for tidal signals to propagate into the inland boundary. The signal reflection would thus moderate the inland freshwater input. In this way, the contribution of the antecedent tidal conditions would decrease, whereas that of the present tidal conditions would increase. This was clearly shown by α i /β i values for both DTL σ and DTL M ; for example, with L I reduced from 200 to 150, 100, and 50 m, α i /β i for DTL σ decreased from 7.02 to 4.48, 2.38, and 1.44 days, respectively (Table 3 and Figure 8). As the regression model considered the antecedent tidal conditions only, it became poor for the cases with a short inland aquifer extent (e.g., that with L I 50 m).
Overall, α i /β i values for both DTL σ and DTL M were less than 10 days, suggesting that the recent tidal conditions were affected more on the freshwater input (Figure 8). Both the variation of α i and α i /β i were monotonic. In particular, for all the cases with fixed-head inland boundaries, α i values were larger than 1, suggesting that the maximum effect of the antecedent tidal conditions was delayed. With L I reduced from 200 to 150, 100, and 50 m, α σ for DTL σ decreased from 2.97 to 2.73, 2.24, and 1.52, respectively, whereas α M for DTL M decreased from 8.43 to 4.87, 2.44, and 2.25 days, respectively. These trends were consistent with the tidal propagation; that is, it took less time for tidal signals to reach the inland boundary and moderate the freshwater input with the reduced L I .
Secondly, we examined SM SW . The regression model for the freshwater input (Eq. 5) well captured the variation of SM SW for the cases with fixed-head inland boundaries (Table 4 and Supplementary Figures S5, S6) as well as that for the case with fixed-flux inland boundary (all the R 2 values were larger than 0.95). As L I increased, both α i and α i /β i values increased monotonously (Figure 9), suggesting that the effect of the antecedent tidal conditions on SM SW was further delayed and prolonged as the inland aquifer extent increased.   a L I is the inland aquifer extent; a σ , a M , and a 0 are the coefficients of the regression model; DTL σ and DTL M are the parametric regressors relevant to X σ and X M ; α i and β i are, respectively, the shape and scale factors of the Gamma distribution function; the ratio αi /β i is the mean of the Gamma distribution function.
Frontiers in Environmental Science | www.frontiersin.org January 2021 | Volume 8 | Article 599041 9 Finally, we checked SM USP using the regression model given in Yu et al. (2019b) for the cases with fixed-flux inland boundaries.
where ω is the angular frequency related to the low-frequency of 6.50 × 10 -8 Hz (period of 178.1 d). a A and δ are the amplitude and the phase shift of this signal, respectively. This regression model failed to capture the variation of SM USP , particularly under the conditions of L I 50 and 100 m ( Table 5 and Supplementary Figures S7, S8). The adjusted R 2 values for L I 50 and 100 m were, respectively, 0.23 and 0.47. This further demonstrated that the hydrodynamic process in the aquifers with fixed-head inland boundaries was more complicated than that with fixed-flux inland boundaries. As discussed earlier, in the former, the inland freshwater input was variable and affected by both the standard deviations and the averaged values of the antecedent tidal conditions. Therefore, we introduced the following new regression model with DTL M included to account for the varying freshwater input.
Then, the variation of SM USP was well captured in both trend and magnitude for the cases with fixed-head inland boundaries ( Table 5 and Supplementary Figures S9, S10). The fitted results were improved in the new regression model (adjusted R 2 of 0.23, 0.47, 0.63, and 0.70 for the old model (Eq. 6) in comparison with adjusted R 2 of 0.92, 0.88, 0.87, and 0.87 for the new model (Eq. 7). It is interesting that α i /β i values for both DTL σ and DTL M monotonically increased as L I increased ( Figure 10). However, α i values for DTL σ and DTL M showed opposite trends. The former decreased but the latter increased. This suggests that as  a SM SW is the amount of salt stored in the saltwater wedge (per unit width aquifer); L I is the inland aquifer extent; a σ , a M , and a 0 are the coefficients of the regression model; DTL σ and DTL M are the parametric regressors relevant to Xσ and XM; αi and β i are, respectively, the shape and scale factors of the Gamma distribution function; the ratio αi/β i is the mean of the Gamma distribution function. the inland aquifer extent increased, the effect of the mean tidal level on the USP was further delayed, whereas the effect of standard deviation was further enhanced. The USP developed in the seaside aquifer and was altered by both the fluctuation intensity of tides and the inland freshwater input. The opposite trends of the enhanced tidal effect and the reduced effect of the inland boundary on the USP are consistent with the previous studies (Werner et al., 2013;Robinson et al., 2018).

DISCUSSIONS
There are generally two types of inland boundary conditions for coastal unconfined aquifers: topography-controlled and constant freshwater discharge, to which the fixed-head and fixed-flux inland boundaries could be applied, respectively (Werner and Simmons, 2009). Given the similar freshwater input, the yearly averaged SGD and salinity distributions remained similar in the  aquifers subjected to the multiconstituent tides. However, with the fixed-head inland boundary, the multiconstituent tides led to more remarkable fluctuations of the USP and SW, in comparison with those for the case with fixed-flux inland boundary. This reflects intensive mixing between freshwater and seawater, which was shown to intensify as the inland aquifer extent decreased. As chemical compositions of seaward seawater and land-sourced fresh groundwater are apparently different, these enhanced mixing between freshwater and seawater would favor biogeochemical reactions in coastal aquifers (Robinson et al., 2018).
With the fixed-head inland boundaries, the lateral hydraulic gradient varied in response to the multiconstituent tides. This led to a varying inland freshwater input, which was more remarkable as the inland aquifer extent decreased. However, the SGD, as the total water efflux, kept a similar fluctuation trend for all the cases and was not significantly altered by the inland aquifer extent. Based on the principle of mass conservation, it is expected that the varying inland freshwater input would moderate the storage of both the freshwater and land-sourced chemicals stored in the aquifer.
The spectrum analysis of the SGD and SM SW showed that the importance of different signals was similar for the aquifers with fixed-head or fixed-flux inland boundaries. The spectrum of freshwater input demonstrated the significant effect of lowfrequency signal (365-day period) in the fixed-head inland aquifers and its influence increased as the inland aquifer extent increased. This suggests that the signal of the year-long period in the tidal series was enhanced in the inland freshwater input. The influence of this signal also increased with the increased inland aquifer extent for the size of USP in the aquifers with fixed-head inland boundaries. Consistently, the 365-day period signal did not appear in the spectrum of SM USP for the aquifer with fixed-flux inland boundary. These results were overall consistent with previous studies that coastal aquifers with fixed-head inland boundaries are more responsive to changes to ocean boundary conditions as groundwater discharge and salt distributions can be easily changed (Werner and Simmons, 2009;Michael et al., 2013). Our results further highlighted the importance of low-frequency tidal fluctuations. Not only the SW but also the USP was affected by low-frequency signals. In particular, the USP was affected more by lowfrequency signals as the inland aquifer extent increased. At a longer-time scale, the location of the fixed-head inland boundary in a natural aquifer depends on the local topography and mean sea level. A relative sea-level rise, as a combination of the two, is expected to reduce the inland aquifer extent. In this way, the highfrequency tidal signals would play a more important role in affecting salinity distributions (both the SW and USP).
SGD is extensively recognized for delivering considerable land-sourced chemicals to the sea and thus affects marine ecosystems such as inducing coastal eutrophication. Salinity distributions in aquifers affect terrestrial ecosystems and fresh groundwater use in coastal areas. Based on our results, the response of the SGD to multiconstituent tides is mostly instantaneous. This suggests that an intensified oceanic forcing condition may instantaneously lead to large input of chemicals to the sea. The variations of the SW and USP are insensitive to the present tidal conditions but related to the antecedent tidal conditions. Oceanic forcing conditions would affect the salinity distributions in aquifers in a prolonged fashion. This suggests that the management of coastal groundwater resources such as control of seawater intrusion could carefully consider antecedent oceanic forcing conditions and long-lasting impacts of human activities (e.g., groundwater extraction and recharge). The regression models presented in this study provide a way for predicting salinity distributions in coastal aquifers under different scenarios and would help guide decision-making.

CONCLUSIONS
This article has examined the SGD and salinity distributions in coastal aquifers with fixed-head inland boundaries and subjected to the influence of multiconstituent tides. The numerical simulations with different inland aquifer extents represented the typical topography-controlled unconfined coastal aquifers. The following conclusions can be drawn: • Given a similar yearly averaged freshwater input, the SGD responded instantaneously to the multiconstituent tides. The effect of the inland aquifer extent on SGD was minor. However, the USP and SW fluctuations were enhanced as the inland aquifer extent decreased. • Inland freshwater input was affected by the combination of the inland aquifer extent and multiconstituent tides. Filtering effect of the aquifer on the high-frequency tidal signals was weakened as the inland aquifer extent decreased and that led to enhanced fluctuations of the freshwater input. The varying freshwater input in turn affected remarkably the salinity distributions in the aquifer. • The regression modeling approaches developed in Yu et al. (2019b) were shown to be able to quantify the variation of the SGD and SW but did not work well for the USP. The cumulative influence of the antecedent lateral hydraulic gradients (as a combination of inland watertable and tidal level) was further included in the new regression model for the USP.
In natural systems, the inland condition is more complex than that of either fixed-head or fixed-flux inland boundary as the inland freshwater input is commonly irregular in response to local rainfall, evapotranspiration, and anthropic groundwater use. While our regression models well captured the SGD and salinity distributions in both aquifers with fixed-head and fixed-flux inland boundaries, they remain to be tested for real aquifer systems. It is worth noting that most natural aquifers are also affected by high-frequency waves and episodic rainfall. All these forcing factors would affect SGD process and salinity distributions in a prolonged fashion (Xin et al., 2014;Yu et al., 2017). Long-term investigations are necessary for applying our regression models in real coastal aquifers. With available longterm datasets, possible quantitative analysis on the interaction among these forcing factors would help to uncover complex SGD process and salinity distributions. Nevertheless, the present study has provided guidance for future research, including field investigations that aim to establish the principle link of SGD and salinity distributions to hydrological forcing factors.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

XY:
methodology, formal analysis, data curation, conceptualization, and writing the original draft. PX: conceptualization, methodology, supervision, and writing, review and editing. CS: methodology and supervision. LL: methodology and writing, review and editing.