A 2-D coastal wetland model was developed to study the salinity change and solute transport of aquifers with variable saturation and density using the COMSOL finite element software. The simulation of flow with variable saturation and density was based on Richard’s equation (Equation 1). Solute transport was based on the advection–dispersion–reaction equation (Equation 2) in porous media (Shuai et al., 2017):

where ρ is fluid density (kg m^-3), C_i is the concentration of species i (mM), P is the pressure (Pa), and Z is the elevation head (m). The value of k_s is the saturated hydraulic permeability (Equation 3), while k_r stands for the relative permeability (m s^-1) (Equation 4), u is the Darcy velocity (m s^-1), θ is the water content (-) (Equation 6), g is the acceleration of gravity (9.81m s^-2), Q_m in this equation represents a stress source term (kg m^-3 s^-1), w_m is the specific moisture capacity (m^-1) (Equation 5), D is the hydrodynamic dispersion coefficient (m² s^-1), R_i is the reaction rate for species i (mM s^-1), S_e is the effective saturation (-) (Equations 7, 8), and S is water storage (P_a ^-1). The constitutive relationship between relative permeability and pore pressure obeys the van Genuchten model (van Genuchten, 1980):

where θ_s and θ_r is the saturated water content and relative water content (-), respectively. K is the hydraulic conductivity (m s^-1), where n_i , which is related to m in Equation 8, and a is the fitting parameter which describes the shape of both the moisture (m^-1) and relative permeability functions (-) obtained by van Genuchten (van Genuchten, 1980).

The model is a 2-D beach aquifer cross-section ( Figure 1 ), where the left and bottom boundary (BAE) and the upper boundary (CD) are both set as no flow and zero solute flux boundaries (Xin et al., 2010; Shen et al., 2018). The right boundary (DE) is set as the inland boundary, while the ocean boundary (BC) is the permeable layer that is used to represent the interface between the beach surface and seawater or atmosphere and assigned a semi-pervious layer or seepage face boundary. The process is realized by the following formula (Equations 9–11):

where R_b is the conductance term (s^-1), defined as the ratio of the saturated hydraulic conductivity (k_s ) with a coupling length scale (L, m). It was set at a high value allowing water to readily move in and out of the interface. H_b is the external head representing sea level (m), H is the total head (m), and C_sea is seawater salinity (ppt).

In order to deal with the boundary conditions reasonably, two conditions have been considered. When the aquifer is saturated and the external head (H_b ) of the sea level is higher than the elevation head, the head difference is determined by the difference between the overlying seawater and the beach interface. However, if the total head (H) is higher than the sea level (H_b ), and the pressure head is equal to or higher than the atmospheric pressure. The boundary is a seepage surface. If the value of R_b is 0 when the aquifer is unsaturated, the boundary is a no-flow boundary. The salinity is determined by the groundwater velocity at the node. If the velocity is inward (to the aquifer), the salinity is estimated by seawater salinity (C_sea = 35 ppt). If the groundwater flows from an aquifer, the zero concentration gradient was specified at the nodes. Previous studies have demonstrated that this modeling approach is useful in describing the water level movement and residence times within the transition zone (Xin et al., 2010; Shen et al., 2019).

A time-varying head (H(t)) Equation 12 was implemented at selected cells in this zone by:

where H(t) is tidal head (m), determined by time; H_msl is the mean sea level (10 m); A is tidal amplitude (0.5 m); ω = 2 π / T is the angular frequency; and T is the tidal cycle (semi-diurnal tide, 12 h). These parameters were used in many studies (Anwar et al., 2014; Shuai et al., 2017).

The model domain represented a homogeneous and isotropic coastal aquifer with a thickness of 12 m and a sloping beach boundary (slope of 0.1; Figure 1 ), which was comparable to a typical sandy coastal aquifer system (Carsel and Parrish, 1988) with a hydraulic conductivity of 15 m d^-1, longitudinal dispersion coefficient (a _L) of 0.2 m, transverse dispersion coefficient (a _T) of 0.02 m (Anwar et al., 2014), and porosity of 0.45. Moreover, seawater salt concentration and density were set to 35 ppt (mass fraction) and 1,025 kg m^-3, respectively. The freshwater salt concentration of 0 ppt and a density of 1,000 kg m^-3 (see Table 1 ) were both adopted from Robinson et al. (2014) and Shen et al. (2019).

The preferential flow had a width of 0.16 m and was assumed to be a homogeneous medium with high permeability and porosity (the permeability was 100 times higher than that of the matrix and the porosity was set as 1). This technology enabled simulations of the preferential flow rapid responses to tidal-water-level fluctuations without affecting the simulation results (Xin et al., 2009). Similar techniques and parameters have been used in groundwater models to simulate pore water flow and solute transport induced by macro-porosity (Xin et al., 2009; Xiao et al., 2019). The preferential flow was set upstream of the intertidal zone, and this was adopted from other studies where it was reported that preferential flow is mainly concentrated in the upper supratidal zone and intertidal zone (Fanjul et al., 2008; Ying, 2021).

For all simulations, the mesh density increased in the preferential flow zones and at the beach. The values for the Courant number and numerical Péclet number did not exceed 1 and 4, respectively, which satisfied the stability criterion and avoided numerical oscillation. During the simulation, a constant temperature was maintained, resulting in constant dynamic viscosity.

The simulation first conducted a no-preferential flow model under the tidal action. Subsequently, to investigate the influence of depth and quantity of preferential flow, we set different depths (2.3, 3.3, and 4.3 m) and quantities (from 1 to 3) of preferential flow conditions under the tide action (five groups of experiments were conducted). The specific parameters are listed in Table 2 .

To study the transformation of marine nutrients, a reaction network of four reactive species was applied to the reaction model. The model considered nitrification, denitrification, aerobic respiration, and DOC degradation.

The reaction model simulated the first-order reaction that required mixing for the inflow of DO, NO₃ ^-, NH₄ ⁺, and DOC. The study omitted the production of NH₄ ⁺ from DOC degradation and nitrification and only focused on the transformation of NO₃ ^- from marine-derived sources. The reactions and kinetic rate expressions are shown in Table 3 . Kim et al. (2017) validated this reaction network by comparing the simulation results with field data from the Cape Shores beach adjacent to Delaware Bay.

Based on the studies of Bardini et al. (2012) and Spiteri et al. (2008), the parameters of the reactions are shown in Table 1 . This study determined the influence of preferential flow under tide influence on the chemical transformations in a nearshore aquifer. The sensitivity analyses on the specific kinetic parameter values adopted were not performed, although these parameter values vary in different coastal settings. Similar to salt transport, the nutrient concentration was set along with the aquifer–ocean interface for the reactive solutes. The boundary conditions of the nutrients are presented in Figure 1 , and the specific nutrient values are shown in Table 1 . The TOE (the saltwater wedge toe) was used to evaluate the length of the saltwater wedge, and the NO₃ ^- removal efficiency (R_N ) Equation 13 was quantified by:

where Ω is the domain area, t is the tidal cycle time (12 h), f_b is the boundary flux, R NO 3 − is the rate of NO₃ ^- reaction, and l is the length of the boundary layer.

We performed four steps: (i) groundwater flow and salinity distributions for all cases at the static water level were run to steady state, (ii) the hydraulic gradient and salinity distributions obtained at the static water level were used as the initial conditions under tidal action, (iii) the salinity distributions obtained in step (ii) were used as the initial conditions for the chemical solute transport and were run to steady state; (iv) and the results of the previous step were run as the initial conditions for the reaction network until they stabilized. The simulations were run for 700 days to allow the reactive solutes to reach equilibrium.

To verify the reliability of the model, we developed a numerical model consistent with the lab experiment of Xie et al. (2023), with an aquifer length of 4 m and width of 0.7 m and porosity of 0.45; the parameters a and n_i were set to 5.9 m^-1 and 2.68, respectively, and the longitudinal dispersion coefficient (a _L) was 0.005 m and transverse dispersion coefficient (a _T) was 0.0005 m [for the specific parameters, please refer to Xie et al. (2023) for details]. We consider three case conditions: non-preferential flow, preferential flow at x = 0.8 m, z = 0.3 m, and preferential flow at x = 1.8 m, z = 0.3 m, and the preferential flow scales and parameter settings are also consistent with those of Xie et al. (2023).

From the experimental results in Xie et al. (2023) and numerical simulation results ( Figure 2 ), the simulated results matched the experimental in all cases well ( Figure 2 ). We can find that tidal forcing formed a USP and SW in the intertidal zone where saltwater infiltrated near the high tide mark and was discharged near the low tide mark. A fresh groundwater discharge zone (FDZ) separated the circulation cell from the lower interface located at the base of the beach. As expected, under the impact of preferential flow, when the preferential flow location is at x = 0.8 m, it has some effect on the salinity distribution and reduces the saltwater wedge deed extension. When the preferential flow location was at x =1.8 m, it has little effect on the salinity distribution, which is consistent with the salinity distribution pattern of the no-preferential flow condition.

The experimental and numerical results consistently demonstrated that the preferential flow could affect the salinity distribution in terrestrial environments, and the 2-D numerical model could be used to predict the water–salt exchange processes in aquifers.

The salinity distribution under two different conditions is presented in Figure 3 . Preferential flow modified the salinity distribution under tidal influence. The average salinity (average value of the whole aquifer) in the aquifer decreased from 7.21 to 7.03 ppt, while the TOE decreased from 24.2 to 23.09 m. Under the tidal condition ( Figure 3 ), the area of the mixing zone (MZA) (the area between 10% and 90% salinity contour) was 53.87 m², the area of USP (USPA) was 10.37 m², and the area of SW (SWA) was 156.0 m². However, under the preferential flow, the mixing zone and USP were increased, and SW was decreased, i.e., the MZA is 64.20 m², the USPA is 12.98 m², and the SWA is 148.79 m². This indicated that the preferential flow enhanced freshwater–saltwater mixing and USP ( Figure 3E ); the preferential flow weakened the length of SW in the aquifer and reduced the salinity.

To determine the effect of preferential flow on groundwater flow, saltwater circulation, and saltwater–freshwater exchange under tidal conditions, the travel time for specified particles in the aquifer was calculated [readers may refer to Xin et al. (2010) for more details]. Three particles were released both near the inland boundary and the beach surface (refer to Figure 3 , Table 4 for the specific release location).

Under tidal action ( Figure 3 ), a USP developed in the intertidal zone, with particles released from inland along the moving path of freshwater bypassing the USP and SW discharge to the sea. The deeper particles took longer than those in the shallow, e.g., the one starting from the shallow (x = 75 m, z = 9 m) took 257.1 d, while that from the deep inland (x = 75 m, z = 3 m) took 267.8 d ( Figure 3 ), increasing by 4.2%. The difference in travel time among the three particles at the seaside varied dramatically as the particles released to the deeper beach moved further landward and took longer paths, e.g., the particle starting from x = 13 m, z = 8.8 m took 176.7 d, while that from x = 3  m, z = 7.42 m took 1,026.9 d ( Figure 3 , Table 4 ).

Preferential flow leads to pore water flow and an increase in USP. Therefore, the particles released from the inland shallow layer experienced USP resistance across the preferential flow, with the travel time increased slightly, e.g., it took 257.1 d for particles released from x = 75 m and z = 9 m and 258.2 d for preferential flow ( Figure 3 ). The travel time of particles released from the beach was consistent with the travel time under non-preferential flow, decreasing with an increase in USP, e.g., the particle released from x = 3 m and z = 7.42 m needed 1,026.9 d with the non-preferential flow but only 908.4 d in the presence of preferential flow.

To further discuss the influence of preferential flow on groundwater flow and circulation, Figure 4 presents the velocity distribution of rising tide, high tide, falling tide, and low tide without and with preferential flow. Only the velocity distribution in the area of preferential flow is shown for comparison. As shown in Figure 4 , the velocity distribution differed with the seawater level. When the tide rose, salt water flowed into the aquifer through the intersection of the tidal level and the aquifer. When the tide fell, the tidal water retreated, and the pore water flowed to the ocean ( Figure 4 ). Under the preferential flow condition at the rising tide, a relatively large pore water velocity appeared and rapidly flowed into the aquifer, thereby providing a preferential gateway for the pore water. When the flood completely flowed over the marsh platform (t = 3 h), the saltwater level reached the maximum (10.5 m), and a sharp rise in the pore water velocity near the preferential flow model was observed [ Figure 4 —(2B)]. The direction of the pore water flow under preferential flow was vertically downward, and the magnitude increased considerably. This indicated that the over-topping water rapidly enters preferential flow and diffuses into the surrounding marsh soil. Furthermore, preferential flow enhances the mixing and dilution of saltwater–freshwater, thus reducing the salinity of groundwater and increasing the quality of drinking water.

The distributions of DO, NO₃ ^-, NH₄ ⁺, and DOC under tidal [ Figure 5 —(1A), (2A), (3A), and (4A)] and preferential flow conditions [ Figure 5 —(1B), (2B), (3B), and (4B)] are presented in this section. There were differences in the distribution of nutrients. DO was mainly concentrated in the shallow layers of the aquifer and NO₃ ^- was in the shallow and medium layers of the aquifer (Note: we used shallow, middle, and deep zones to differentiate the characteristics of the distribution area of DO, NO₃ ^-, NH₄ ⁺, and DOC. In this study, these three intervals are a relative region varying with the preferential flow, so we are only distinguishing between the same model) [ Figure 5 —(1A), (1B), (2A), and (2B)], while NH₄ ⁺ [ Figure 5 —(3A) and (3B)] and DOC [ Figure 5 —(4A) and (4B)] were distributed in the entire aquifer with a relatively high concentration in the shallow layer and in the middle and shallow layers, respectively, e.g., the average concentration of DOC in the region z > 4 m is 0.058 mM larger than that at z< 4 m (0.025 mM), and the average concentration of NH₄ ⁺ in the region z > 4 m is 0.012 mM larger than that at z< 4 m (0.011 mM). The nutrient concentration increased along preferential flow in the reaction process [ Figure 5 —(1B), (2B), (3B), and (4B)].

The reaction rate distributions of nitrification [ Figure 6 —(1A) and (1B)], denitrification [ Figure 6 —(2A) and (2B)], respiration [ Figure 6 —(3A) and (3B)], and DOC degradation [ Figure 6 —(4A) and (4B)] are described in this section.

Nitrification and respiration occurred primarily in the surface layers of the aquifer, while denitrification was dominant under the influence of anaerobic bacteria in the middle and deep layers where oxygen is depleted ( Figure 6 ). DOC degradation occurred in the USP and SW of the whole aquifer under tidal action and preferential flow. However, the reaction rate in the middle and upper layers was relatively high [ Figure 6 —(4A) and (4B)], e.g., the average rate of DOC degradation in the region z > 4 m is 1.49 × 10^-5 mM d^-1 larger than that at z< 4 m (6.54 × 10^-6 mM d^-1). Compared to the non-preferential flow condition, preferential flow resulted in the extension of nitrification, respiration, and DOC degradation along with the preferential flow in the reaction process. However, denitrification consistently took place in the deeper USP.

We quantitatively analyzed the NO₃ ^- removal efficiency (R_N ). Under non-preferential flow condition, the R_N is 7.9%. When we insert the preferential flow, the NO₃ ^- removal efficiency (R_N ) was modified, e.g., under the preferential flow action, the R_N is 8.97%. These results suggest that preferential flow increases the NO₃ ^- removal efficiency, thus decreasing the salinity distribution.

A USP promotes the retreat of SW under tidal action, which renders the SW closer to the ocean and limits saltwater intrusion. Preferential flow enhanced ( Figure 3 ) this effect, driving SW further seaward, and this impact gradually increases with the depth of preferential flows ( Figures 7 , 9 ). These results are consistent with those of Gao et al. (2023), but they consider the under spring-neap tide action. The high hydraulic conductivity of the preferential flow increased the salt transport, resulting in a larger USP than the non-preferential flow; this could control SWI to some extent. Several previous studies have been performed to investigate saltwater–freshwater mixing dynamics—for example, Xie et al. (2023) investigated the effect of fracture characteristics on salinity distribution and groundwater flow through experiments and numerical simulations. They reported that the vertical fractures had a limited impact on most SWI properties. This seems to be different from our results. This is mainly due to the effect of fracture on SWI properties depending on their relative position to the saltwater. In our study, we focused on the longitudinal rift located inside the USP and connected to it. In addition, they noted that fractures can increase the mixing zone, which is consistent with our results. Preferential flow acted as a drain for the shallow layer during the falling tide and for a recharge well during the rising tide. This was the main reason for the increase in USPs and the decrease in the salinity of the aquifer. Additionally, Xiao et al. (2019) found that the salinity in crab burrows in the intertidal zone at high tide was slightly higher than in the soil matrix, e.g., the salinity of crab burrows reached 34 PSU and larger than the nearby soil matrix (33.4 PSU). This suggested that the presence of crab burrows can greatly enhance salt transport in salt marshes.

The nutrients from the sea were initially input in the numerical model with the salinity simultaneously. Nitrification primarily occurred in the surface layer ( Figure 6 , Supplementary Figures S2 , S4 ) since nitrification, respiration consumption of deeper DO, and denitrification were dominant in the middle and deep layers of the aquifer ( Figure 6 , Supplementary Figures S2 , S4 ). Shuai et al. (2017) had the same report, and they assumed that DO, NO₃ ^-, NH₄ ⁺, and DOC primarily originate from rivers. They noted that nitrification occurs in the shallow layers, while denitrification occurs in the deeper anaerobic layers. In the presence of preferential flow, nitrification mainly occurred in the surface layer of the aquifer and extended along with the preferential flow ( Figure 6 , Supplementary Figures S2 , S4 ). Nitrification and respiration consumption of deeper DO ( Figure 5 , Supplementary Figures S1 , 3 ) and the preferential flow increase in the area of sediment air/water interface increased the exchange between surface water and groundwater and the concentration of DO in the aquifer. This resulted in the extension of nitrification and respiration along the preferential flow ( Figure 6 , Supplementary Figures S2 , S4 ) and decreased the carbon stock of the sediment. Preferential flow increased the concentration of nutrients at the USP by dramatically enhancing nutrient transport and reducing SW as the USP pushed the SW seaward. However, total nitrification increased under preferential flow conditions and gradually increased with depth and amount of preferential flow; however, it decreases after a certain quantity.

In the presence of tidal action, the distribution of nutrients in the aquifer formed a USP and a SW. The travel and residence times were longer in the deeper layers of the aquifer for nutrients, nitrification, and respiration consumption of deeper DO ( Figure 5 , Supplementary Figures S1 , S3 ), creating anaerobic conditions for denitrification. Subsequently, the NO₃ ^- produced was used up during denitrification, resulting in the reduction of NO₃ ^- in the deep layer ( Figure 5 —(2B), NO₃ ^- distribution). Denitrification occurred in the deeper USP as nitrification in the upper USP not only provided the reactant for denitrification but also offered an anaerobic environment. The NO₃ ^- removal efficiency increased compared to that in non-preferential flow, and this effect increases progressively with the depth and number of preferential flows; however, it decreases after a certain quantity. It can be observed that the existence of macro-porosity (such as crab burrows and invertebrate nests) facilitates the removal of NO₃ ^- from sea sources. Some previous studies demonstrated that denitrification also increased with the size of the mixing zone (Spiteri et al., 2008; Heiss et al., 2017; Gao et al., 2023). Heiss et al. (2017) reported that saltwater–freshwater mixing could promote denitrification in the intertidal zone.

In addition, Gao et al. (2023) investigated the effect of macropores on the salinity distribution and denitrification of the aquifer under spring-neap tide action. They pointed out that macropores can lead to NO₃ ^- removal efficiency increase. They investigated the reaction between DOC from the sea source and NO₃ ^- from the land source. In this study, we considered the relationship between DOC from the sea source and NO₃ ^-, which is the main reason for the difference, and it is also noteworthy that they also suggested that macropores can increase the mixing of saltwater and freshwater. This is consistent with our results. Thus, our study provides new insights into the relationship between macropores and salinity distribution and denitrification in the aquifers.

This study elucidated the effect of preferential flow on salinity change and solute transformation under tidal conditions. It determined the importance of the depth and quantity of preferential flow, which guide the mitigation of saltwater intrusion and contaminant removal. However, the location and density (Fanjul et al., 2008; Ying, 2021) distribution of preferential flow should be further studied. The preferential flows are based on the actual and generalized, and the real aquifer preferential flow distribution pattern and specifications are extremely complex. Future research can focus on the real aquifer preferential flows on the saltwater intrusion characteristics as well as on the impact of solute transport transformation. Waves drive seawater recirculation (Anwar et al., 2014), and seasonal groundwater level drives groundwater fluctuations (Michael et al., 2005; Liu et al., 2016). This leads to the more complex nearshore mixing dynamics of freshwater and seawater, which, in turn, complicates NO₃ ^- transformation. In addition, other factors may also impact the biogeochemical processes in aquifers, including precipitation, evaporation, and long fluctuating tides (e.g., spring-neap tides), and sediment heterogeneity (Heiss and Hichael, 2014; Geng and Boufadel, 2017; Kreyns et al., 2020; Gao et al., 2023; Zheng et al., 2023), which need to be considered in the future work.