Salt Transport Under Tide and Evaporation in a Subtropical Wetland: Field Monitoring and Numerical Simulation

Tidal wetland ecosystems are sensitive to porewater salinity dynamics. However, it is unclear how salts move and distribute in these wetlands, particularly how the salts accumulated by evaporation get removed from the wetland soil, so that the salinity levels may stabilize to accommodate vegetation. We conducted a combined field and modeling study to identify the porewater flow and salinity patterns in a subtropical wetland subjected to tidal inundation and evaporation. Measured and simulated salinity contours indicated the formation of hypersaline porewater plumes in the upper intertidal zone and the fresher porewater zones with salinity close to that of seawater near the creek and in the supratidal zone. Simulations indicated the discharge of the hypersaline upper intertidal porewater to the creek with a discharge pathway developed under the fresher near‐creek porewater zone, and this was further confirmed by the field‐observed significant salinity gradient under the creek. Our model suggested that both water and salt discharge from the wetland soil occurred predominantly at the creek bank and creek bed. The porewater discharge is more intensive through the creek bank than the creek bed, while the salt discharge across both the creek bank and creek bed was comparable due to the much higher salinity level under the creek bed. Salt discharge driven by density gradients and tidal‐induced porewater circulation provides a mechanism for removing salts accumulated in the upper intertidal zone due to evaporation and could prevent salt flat formation and marsh plants dieback.

surface. Therefore, porewater flow and salt transport patterns in a tidal wetland could be significantly different from those in other nearshore environments.
Many studies, including numerical simulations and field observations, have been conducted to characterize tide-influenced porewater flow in estuarine creek-wetland systems. Porewater transport near the creek is predominantly governed by tide amplitude (Wilson et al., 2015), soil permeability and stratigraphy (Gardner, 2007;Xin et al., 2012), and macropores Xin et al., 2009). Frequent tidal fluctuation and soil submersion lead to an aerated active porewater exchange zone close to the creek, creating idealized conditions for pioneer tidal wetland marsh plants (Nuttle, 1988;Wilson & Gardner, 2006;Xia & Li, 2012;Xin et al., 2013). Toward the landward intertidal and supratidal zone, the groundwater level fluctuation is attenuated due to the tidal signal damping over the distance to the creek (Wilson et al., 2011). Porewater flow in these interior areas is more sensitive to evaporation and rainfall as they are less frequently inundated by tides (Cao et al., 2012;Montalto et al., 2006;Wilson & Gardner, 2006;Xin et al., 2011).
Compared with porewater flow, solute transport in the tidal wetland subsurface system, which includes salt distribution and movement in the unsaturated zone and the saturated aquifer, has been less studied and characterized . Several numerical simulations (Table 1) and some field investigations have attempted to trace the tidal wetland salt distribution with a focus particularly in the intertidal zone. Werner and Locking ton (2006) demonstrated that tidal oscillation in an estuarine creek reduces the extent of the horizontal seawater intrusion compared with a steady creek water level. The tide-induced salt exchange between the wetland aquifer and seawater was found to occur only through the section of the creek bank within the tidal range (Wilson & Gardner, 2006). Seawater salinity in the creek could vary temporally either due to tidal flushing or rainfall. A rainfall-induced drop of seawater salinity may result in a short-lived freshening of the creek hyporheic zone (Lenkopane et al., 2009). By adopting a tidally varied seawater salinity in their numerical study, Xiao, Li, Xia, et al. (2019) found that density variation of the creek water leads to the formation of salt fingers in the aquifer under the creek and dilution of the saltwater wedge near the creek.
Despite the previous numerical characterizations on intertidal wetland salinity, only tidal force was considered as the governing force to the porewater flow and salt transport in their systems. Salt transport and distribution could be affected by evaporation, which was however ignored in their models (Xin et al., 2022). Surface evaporation results in salt accumulation, particularly near the soil surface (Wang et al., 2007;Zhang et al., 2014). Transpiration by vegetation, such as the mangroves in tidal wetland, could also increase the porewater salinity as salts are usually excluded during vegetation water uptake (Fass et al., 2007;Passioura et al., 1992). The spatial variation of evaporation rate due to changes in micro-topography, saturation level, and soil properties may cause nonuniform surface salt distribution as evidenced by the reported salinity patches on the tidal wetland surface from previous field investigations (Albuquerque et al., 2014;Galliari et al., 2021;Miklesh & Meile, 2018;Poh, 2013;Ridd & Stieglitz, 2002).
The density contrast between the overlying higher salinity porewater and the less saline porewater could trigger unstable fingering flow, which can alter the flow and salt pattern in the groundwater and the exchange between the creek water and wetland porewater (e.g., studie 5 and 6 in Table 1, Shen et al., 2015;Shen et al., 2016). The evaporation-induced salt accumulation near the wetland soil surface, on the other hand, could form salt crusts (see the efflorescence and subflorescence in Zhang et al., 2014) and thus may change soil properties and evaporation rate. However, most of the previous tidal wetland models (Table 1) only considered evapoconcentration in a simplified form if at all. While a few simulations incorporated tidal fluctuation and nonphysically based salt evapoconcentration (i.e., no salt crust precipitation and associated evaporation rate change) in their model, they either prescribed a flat wetland surface fully inundated at high tide (study 7 in Table 1, Xin et al., 2017) or the studies (studies 4 and 8 in Table 1, Geng & Boufadel, 2015 were based on beach environment subject to strong salinity reduction by inland fresh groundwater. In a tidal wetland with shallow groundwater and limited inland freshwater input, evapocencentration and unstable salt fingering are expected to be significant (Simmons et al., 2002), and the salt transport and distribution in this environment could deviate from those in the past investigations. Therefore, a more comprehensive characterization of both intertidal and supratidal porewater salinity distribution is needed. More importantly, previous salt distribution and transport patterns were investigated through numerical simulations and are yet to be verified by field data, such as spatial and temporal porewater salinity, porewater pressure head, and local weather.
This study aims to investigate in the field and by numerical simulation the porewater flow and salt transport in a coastal wetland under the combined influence of tides and evaporation. We chose a subtropical tidal wetland for investigation because the low-lying site is not only subjected to regular seawater inundation, but also strongly influenced by salt accumulation as the vegetation-free salt flat is evident on site. We hypothesized that salt accumulation on the surface as a result of the interplay between tide and evaporation may influence the deeper wetland aquifer salinity and the salt exchange between creek water and porewater. This hypothesis was tested using a combination of spatiotemporal groundwater hydraulic head change and porewater salinity monitoring and a 2D variable-density, variable-saturation porewater flow, and solute transport model with physically based salt evapoconcentration. We also further tested the variation of the modeled salinity by running two additional simulation cases to gain insights of the salt distribution dynamics under a future higher evaporation rate and impacted by a rainfall event.

Study Site
The field site is a tidal wetland adjacent to the Hays Inlet Conservation Park (27°16'01"S, 153°02'59.5"E) near Moreton Bay, Southeast Queensland, Australia. A tidal creek extends into the wetland from the estuary and terminates on the north side of the wetland (Figure 2). The tidal creek only receives surface water from the estuary by tidal inundation and rainfall as no other surface water flows into the creek were identified. Site topography is overall low-lying with the lowest elevation at the creek bed around −0.8 m AHD (Australian Height Datum). The highest elevation is away from the creek in the supratidal area and is about 1.5 m AHD (Figure 3a). The creek bank slope is relatively steep (∼0.175), while the slope of the wetland surface is mild (∼0.002). The local climate is subtropical with an average potential evaporation rate of around 4.8 mm/day in the summer and around 3.5 mm/ day in the winter (BOM, 2021).
Plant zonation (Figure 2) is exhibited in this field site. The plant species identified in the near-creek intertidal zone (∼20 m near the creek) are Ceriop australis (Yellow mangroves) and Sporobolus virginicus (Saltcouch). The rest of the intertidal zone accommodates the succulent high-salinity resistance saltmarsh plants Sarcocornia quinqueflora and Tecticornia halocnemoides (Johns, 2006) with bare soil developed in the lower elevation intertidal zone. The supratidal zone is dominated by S. virginicus, and Casuarina cunninghamiana (Coastal river oak tree) forest occupies the inner land.  Guimond and Tamborski (2021) and Robinson et al. (2007 Local sediment characterization had been reported by Poh (2013). Soil samples were collected along the elevation gradient (from the creek to the interior) at different depths during borehole drilling. Compared with the mud-sand-layered system reported in many other tidal wetland sites (e.g., Guimond, Yu et al., 2020;Xiao et al., 2017), the sand component in the soil of this site is higher with a loamy sand layer to a depth of around 1 m below the mean sea level and a coarser sandy layer underneath ( Figure 3). The shallow 0-20 cm root zone contains more mud, and several clay patches were found near the tidal creek.

Field Methods: Groundwater and Surface Water Monitoring
Groundwater piezometers were installed along the elevation gradient ( Figure 3a) from the creek to the supratidal zone to capture groundwater table fluctuations due to tidal influence. All piezometers were constructed of 50-mm PVC pipe and screened over the lowermost 50 cm. Two sets of piezometers with the same distance to the creek Figure 2. Map of the field site located in Southeast Queensland (the red star is the site location) with the major vegetation zones, the location of groundwater piezometers (P1-P9), and the surface water monitoring point S1.

of 21
while different screen depths, namely P1, P2 (1.5 and 0.5 m, respectively below the surface) in the near-creek intertidal zone and P5, P6 (1.5 and 3 m, respectively below the surface) near the high tide line, were used to determine the vertical groundwater flow condition. Calibrated Aqua Troll 200 data loggers (In-situ Inc. USA) were deployed in the piezometers to monitor the groundwater pressure head and electrical conductivity (EC) with a 15-min recording interval. Elevation along the piezometer transect was surveyed with a rotating laser level (Leica Rugby 640G, ±2.2 mm at 30 m) by comparing the elevation difference with an elevation survey mark (SCR #93850) about 200 m away. Elevation and groundwater hydraulic head were both reported relative to AHD. Saturated hydraulic conductivity of the shallow loamy sand sediments was estimated using the auger hole method (Montalto et al., 2006) at three boreholes during the piezometer installation. Local meteorology data were obtained from the Redcliffe weather station around 2 km away (ID: 040958) (http://www.bom.gov.au/climate/dwo/IDCJDW4099.latest.shtml). Tidal creek water (surface water) level and EC were recorded at S1 by another Aqua Troll data logger installed close (∼10 cm) to the stream bed and were regularly inspected to ensure no sediment burial had occurred. Groundwater and surface water monitoring started in July 2019 with the intensive monitoring carried out during June and July 2020 when most of the piezometers were equipped with loggers.

Field Methods: Soil Water Sampling and Analysis
MacroRhizon suction porewater samplers (Rhizosphere Research Products B.V., The Netherlands) were deployed next to the groundwater piezometers at the depths of 10 cm, 25 cm, 45 cm, and 65 cm below ground level ( Figure 3a) to collect soil porewater using 50-mL syringes (Seeberg-Elverfeldt et al., 2005). Soil water samples were mostly collected in January and July 2020 at low tide when the wetland surface was not inundated by creek water. Samples after rain events were also obtained in January 2020 to identify the changes caused by rainfall. Groundwater sampled from the piezometers and the soil porewater samples were immediately analyzed for their pH, EC, temperature, and dissolved oxygen (DO) with a handhold multiparameter probe (Hanna, HI98194). The chloride concentration of water samples was further determined by discrete analyzer. It should be noted that while the porewater salinity (unit: ppt, parts per thousand) represents the total dissolved salt in the water, which is commonly converted from the EC value using the following equation: Salinity (ppt) = 0.4665 × EC (mS∕cm) 1.0878 (Williams, 1986), we separately identified the seawater-originated salinity (Salinity sw ) from the chloride concentration using the theory of constancy of seawater composition (Salinity sw (ppt) = 1.80655 × 10 −3 × Cl − (mg∕ ) ) (Hughes et al., 2012) to compare the field porewater salinity to the simulated groundwater-seawater salinity exchange.

Numerical Simulation
Numerical simulations were conducted using SUTRASET (Shen et al., 2018), a modified version of the finite element model SUTRA (Voss & Provost, 2010), to simulate the 2D variably saturated porewater flow and solute transport in soil with variable fluid density. SUTRASET is capable of implementing the occurrence of seepage, tide fluctuation, evaporation, and salt precipitation (America et al., 2020;Shen et al., 2018). A summary of the governing equations of the SUTRASET model is presented in Text S1 and Table S1 in the Supporting Information S1. The right blue boundary was prescribed with hydrostatic head with zero pressure located at 0.8 m and a salt concentration of 0.035 kg/kg. The black surface boundary is mixed with specific pressure (i.e., tide overtopping), seepage, and specific flux (i.e., evaporation) boundaries, depending on the surface inundation and saturation condition. High tide and low tide signals are prescribed at 1.30 and 0.88 m, respectively. The gray-dashed line indicates sediment stratification in both figures.

Conceptual Model
A conceptual model representing the 2D creek-normal transect from P1 to P9 was adopted in the numerical model as shown in Figure 3b. The domain consists of one loamy sand layer and an underlying sand layer. The bottom of the domain is impermeable, representing the weathered sandstone clay sediment (Poh, 2013). The domain surface includes creek bed (−20 m ≤ x < −5 m), creek bank (−5 m ≤ x < 0 m), and wetland ground surface (0 m ≤ x < 320 m). The elevation of the domain surface was designed based on the in-situ surface point elevation measurements and the elevation profile from the 1-m LiDAR DEM (Digital Elevation Model) data (Department of Resources, 2014). The creek water level fluctuates with the tide and can submerge part of the wetland surface. Evaporation occurs on the surface once it is exposed to the atmosphere and is unsaturated. Groundwater can flow through the landward boundary of the domain, and the creek center is considered as a natural flow divide (Wilson & Gardner, 2006). Vertical infiltration from the rainfall event is considered as the opposite process to evaporation, which feeds freshwater flux to the model through the exposed wetland surface while it is unsaturated. The 2D transect model, considering the local groundwater flow governed by upland pressure head, seaward tidal fluctuation, and evaporation, represents an idealized flow situation of a transient and more complex 3D flow field.

Boundary and Initial Conditions
The vertical creek center flow divide and the domain bottom were both treated as no-flow boundaries. The inland boundary was treated differently at the calibration and simulation stages. During the calibration stage, the 3 months hydraulic head monitored in P9 was prescribed at the inland boundary to validate the conceptual model and the choice of soil properties. During the long-term simulation with evaporation, constant hydrostatic pressure with the water table located at 0.7 m below the soil surface (z = 0.8 m), representing the averaged hydraulic head recorded in P9, was assigned to the boundary. The surface boundary condition varied depending on the soil saturation and submersion condition: (a) when the wetland surface is exposed to the atmosphere and the surface soil is unsaturated, a Neumann boundary condition was assigned to simulate the evaporative flux; (b) when the soil surface is below the tidally fluctuating creek water level, a pressure head being equal to the hydrostatic head with respect to the creek water level was assigned to the submerged soil surface; (c) when the surface soil is saturated while the surface elevation is above the surface water level, which occurs during the falling tide, seepage was assigned to the surface as described by Wilson and Gardner (2006). The salt concentration of the creek water was taken as 35 ppt as the average seawater salinity. The salt concentration of the groundwater flow through the upland boundary was set as 35 ppt according to the long-term salinity monitoring in a groundwater well near the upland forest (not shown in Figure 2, salinity monitoring results are shown in Figure S1 in Supporting Information S1). Rainfall-induced vertical flux has a salinity of 0 ppt with the infiltration rate, depending on the designed rainfall intensity and local infiltration capacity (Yu et al., 2017). Equations for infiltration rate determination are in Text S2 in Supporting Information S1. During the implementation of the event-based rainfall boundary, the inland boundary was assumed to be a no-flow boundary. Creek water level fluctuating signal was generated by fitting the spring-neap tide records at the monitoring station S1 by tide harmonic analysis (Text S3 and Table S2 in Supporting Information S1) as done by the previous researchers (e.g., Robinson et al., 2007;Xia & Li, 2012;Xiao et al., 2017). Table 2 were assigned using a combination of field measurement data, previously published values, and model calibration. Sensitivity analyses were conducted on the soil hydraulic conductivity , porosity , soil characteristic parameters and , and potential evaporation. Preliminary sensitivity analysis results showed that the groundwater head variations were most sensitive to parameters , , and . The implementation of evaporation rate at 4 mm/day did not cause an evident influence on the groundwater hydraulic heads ( Figure S2 in Supporting Information S1). Therefore, model calibration was conducted by adjusting the saturated hydraulic conductivity and soil characteristic parameters and to achieve the best fit to the observed groundwater hydraulic heads in P2 located in the near-creek intertidal zone, P4 in the lower intertidal zone, and P5 in the upper intertidal zone. Details of the calibration results will be introduced in the following numerical results section.

Numerical Implementation
The model domain consists of 47470 nodes and 46900 elements with the mesh size refined near the creek and the intertidal wetland surface. The vertical mesh size Δz ranges from 0.04 to 0.063 m in the unsaturated zone to better describe the water and salt gradients and the horizontal mesh size Δx = 0.5 m. It should be noted that this selection of mesh resolution is a trade-off between the computation times and the simulation accuracy of such a 2D coastal transect with a length of more than 300 m. Our mesh resolution, particularly the vertical resolution in the unsaturated intertidal zone, is similar to that used by Geng and Boufadel (2015) and America et al. (2020). Finer mesh (i.e., with 196131 nodes and 194600 elements) and a reduced time step were tested, and the simulation results were found to be independent of these numerical parameters.

Simulation Cases
The base case (Case 1) is a long-term simulation under a maximum evaporation rate of 4 mm/day, which represents the local annual-average potential evaporation value (BOM, 2021). The model was run for 22.75 years with a time step of 5 min until the simulated salinity reached a quasi-steady state. Case 2 is designed to simulate the salinity condition under the higher potential evaporation rate (maximum evaporation rate = 8 mm/day), representing the future global warming climate or the climate in the dry coastal regions. The adjustment of maximum evaporation rate was achieved by adjusting the evaporation process parameters in Table 2. An additional shortterm simulation (Case 3) applies a rainfall boundary in replacement of the evaporation boundary, simulating a 12-hr rainfall event with an intensity of 7.91 mm/hr. Case 3 aims to explore the impact of a frequent (50% Annual Exceedance Probability, see Text S4 in Supporting Information S1 for details) storm event on the base case (Case 1) steady salinity condition. Numerical model results were primarily analyzed around base Case 1 by examining the porewater and salt flow patterns of the last two spring-neap tide cycles (15 days) of the long-term simulation. The quasi-steady-state salinity level of Case 2 and the post-rainfall salinity distribution of Case 3 were briefly discussed.  (Carsel & Parrish, 1988). d Laboratory or field measurement. e Similar value as used by Shen et al. (2018). f Value used by Geng and Boufadel (2015) and Wilson and Gardner (2006). g Evaporation-related parameter values used to generate a potential evaporation rate at 4 mm/day with the consideration to represent the semiarid to subtropical environment. N.A.: Calibration not applicable.

Hydraulic Head and Groundwater EC
Groundwater hydraulic head and EC records from 28 June 2020 to 19 August 2020 were retrieved from the transducers at different locations. The local potential evapotranspiration (ET) rate and rainfall data (BOM, 2020) during field monitoring are shown in Figure 4b, where the potential ET was calculated by the Penman-Monteith equation (FAO56-PM) (Allen et al., 1998).
Creek water level exhibited a semidiurnal fluctuation pattern influenced by the spring-neap tide (Figure 4a) with the highest level at 1.34 m and the low level at −0.86 m. Groundwater hydraulic head in the near-creek intertidal zone displayed rapid and significant variations to the tide fluctuation (Figure 4c). During the spring rising tide, the wetland platform in the near-creek zone was inundated and seawater from the creek infiltrated in the soil, causing the rise of groundwater level in P1 and P2. However, the inundation only lasted for a short time (∼2.5 hr/ day) during the spring tide. During the neap low tide period, the creek water level was below the near-creek groundwater level, leading to the discharge of groundwater into the creek and the drop of groundwater level. The creek water level was all the time below the wetland platform during the neap tide period so the near-creek zone was exposed to air. However, groundwater head in the deeper aquifer (P1) did not rise to the ground surface during the inundation, and the variations were much smaller compared with that in shallow groundwater (P2). This was due to the less conductive loam and clay patch surrounding the screen of P1.
In contrast to that in the near-creek intertidal zone, the hydraulic heads in the higher wetland, including the lower and upper intertidal zone (Figure 4e) and the supratidal zone (Figure 4g), exhibited more pronounced time lag and damping of the tidal signal. The hydraulic gradient in P4, P5, and P6 remained positive for a few tidal cycles during the spring tide (Figure 4e). Despite the intermittent rainfall that might help maintain the water table, the stable water level during the spring tide suggests that the drainage in the higher wetland is very slow, particularly at the lower intertidal zone (P4), which is mainly bare soil flats.
The comparison of hydraulic heads in P5 and P6 highlights the vertical head gradient at the same location in the upper intertidal zone. During the recharge periods, as indicated by the increase of hydraulic heads in P5 and P6 (i.e., 6-14 July and 19-22 July 2020; Figure 4e), the upward hydraulic gradient increased, suggesting downward vertical flow due to infiltration. Following the recharge events (i.e., 29 June -05 July, and 14-18 July 2020; Figure 4e), the upward hydraulic gradient between P5 and P6 reduced, suggesting a weakening downward flow.
Surface inundation did not occur in the supratidal zone despite the occasional rainfall during the monitoring period with the water table in P9 varying from 0.98 to 1.18 m (Figure 4). It is worth mentioning that during the transition from spring to neap tide (i.e., 07-16 July), the groundwater level in the upper intertidal zone (P5, P6) bulged across the creek-normal transect with levels higher than adjacent groundwater and creek water. Such an elevated hydraulic head in the intertidal zone may form a groundwater wave that transports landward beyond the boundary of the creek-normal domain.
EC of surface water in the creek showed a slight tidal variation (Figure 4a) around 52 mS/cm (34 ppt), a typical seawater salinity value. The EC of the shallow groundwater near the creek (P2 in Figure 4d) was close to the surface water due to active tidal flushing. It increased slightly during the inundating spring tide period and decreased back to seawater salinity during the neap tide period. The deeper groundwater (P1 in Figure 4d) near the creek was more saline than the shallow groundwater and the surface water, indicating the presence of a vertical salinity gradient. The highest groundwater EC was observed near the high tide line (Figure 4f), which fit the previous report (America et al., 2020;Shen et al., 2018) on the formation of a hypersaline zone near the high tide level due to evaporation. However, the deep groundwater also showed a stable hypersaline condition (P6 in Figure 4f), which indicated the downward migration of the salt. Groundwater EC in the supratidal zone further decreased but was still higher than the surface water EC despite the fact that seawater rarely floods this area. The hypersaline groundwater in the upper intertidal zone and the supratidal zone indicated the accumulation of salt in these areas, and it most likely resulted from the evapoconcentration. The rainfall events only resulted in a slight decrease of the EC of the surface water and groundwater in P5. Overall, there was no significant deep groundwater salinity change observed in piezometers in response to the rainfall.

Porewater Salinity
The chloride-based porewater salinity data collected from the field at 11 January and 16 January 2020 were analyzed and presented as contour plots in Figure 5, since the porewater samples were not collected during the groundwater head observation period from June to July 2020. Rainfall of 50.6 and 38 mm occurred on 12 January and 16 January 2020, respectively. Therefore, the comparison of the porewater salinities between 11 January and 16 January highlighted the influence of rainfall.
Field porewater salinity measured on 11 January clearly showed the formation of a hypersaline plume under the high tide line (x ≈ 110 m, Figure 5a) with the highest porewater salinity being 153.5 ppt more than 4 times that of seawater salinity. Porewater salinity decreased to 73 ppt deeper in the aquifer under the hypersaline plume. The hypersaline plume tends to expand toward the creek. Porewater salinities close to that of seawater (35-65 ppt) were measured near the creek and in the supratidal zone. Near the creek where tidal flushing occurred frequently, a lower porewater salinity "wedge" was formed. Notably, the porewater under the creek bed (x = −5 m, z = −1.8 m and −2.5 m, Figure 5a) showed much higher salinity (47.5 and 66.1 ppt, respectively) than the creek water (35 ppt). This most likely resulted from the discharge of the hypersaline porewater at the creek bed and the creek bank. More details on the formation of high-salinity porewater under the creek will be discussed later in the modeling results.
The average salinity in the hypersaline plume was significantly reduced after the rainfall, and the highest surface salinity decreased from 153.5 to 110 ppt (x ≈ 110 m, Figure 5b). However, the deep groundwater salinity in the intertidal zone did not change significantly after rainfall, consistent with the groundwater EC recorded in July 2020 (Figure 4f). Therefore, the hypersaline plume still existed in the aquifer under the high tide line after rainfall but with reduced salinity. In the lower and near-creek intertidal zone, porewater salinity after the rainfall was higher than that before the rainfall, and the "wedge" formed before rainfall advanced to the creek bank (Figure 5b), indicating seaward movement of salt. The porewater salinity at the supratidal zone was overall freshened by infiltration with surface porewater salinity in some ponding areas reduced to around 15 ppt. Compared to the intertidal zone, deep groundwater salinity in the supratidal zone, particularly near the inland boundary, was evidently reduced by the rainfall. Rainfall induced infiltration saturated the soil in the wetland at the locations that would otherwise never be inundated by tide, forming a stronger hydraulic gradient between groundwater and creek water, particularly at the falling tide.

Numerical Model Calibration
The saturated hydraulic conductivity of the loamy sand layer measured in the field ranged from 5.77 × 10 −6 m/s to 1.48 × 10 −5 m/s. The calibrated saturated hydraulic conductivity values of the loamy sand layer and the sandy layer were 5.00 × 10 −5 m/s and 8.00 × 10 −5 m/s, respectively. The longitudinal dispersivity was 0.1 m and the transverse dispersivity was 0.01 m, which are the values commonly used by other coastal groundwater models in a similar field scale (e.g., Geng & Boufadel, 2015;Shen et al., 2015;Wilson & Gardner, 2006). The diffusion coefficient was 5.00 × 10 −8 m 2 /s in the same order of magnitude as that used by Shen et al. (2018). Sensitivity analysis (Text S5 and Figure S3 in Supporting Information S1) of the diffusion and dispersion coefficients indicated that increasing diffusivity and dispersion both speed up salt migration. The salt transport process is more sensitive to dispersivity, while the diffusion process is weak in the dynamic coastal environment influenced by tide (Geng & Boufadel, 2015). The model was calibrated using the field-observed hydraulic head records from July to September 2019 (starting from 29 July 2019 0:00 a.m.), which covers 1,500 hr in total, as presented in Figure 6. The model well captured the behavior of groundwater head in the cross-creek transect, where the hydraulic head is less influenced by the tide with the groundwater fluctuation amplitude decreasing as the location moves more landward. The discrepancies between the simulation and observation were evident in the landward side (P4 and P5, Figure 6), where the simulated head could be ∼20 cm higher than the observed head. These discrepancies could be resulted from (a) relatively coarse topography characterization of the creek-normal transect by the 2D model since the field site contains many micro-topography features, (b) the driving tidal signal ignores some extreme events, such as king tide, (c) the simplification of each soil layer as to be homogeneous and isotropic, and (d) the lack of evaporation.

Evaporation Rates and Soil Saturation
The soil surface evaporation rate, porewater flow, and soil saturation condition varied at different stages during a spring-neap tide cycle, according to the simulation results. The spring high tide inundated the intertidal wetland surface and resulted in strong downward and landward seawater recharge into the aquifer (Figure 7a). Surface evaporation was disabled under tide overtopping. When the tide retreated to the low tide line, the soil saturation in the lower intertidal zone (see 40 m < x < 80 m, Figure 7b) remained at a high level, causing the mounded water table, which drives groundwater flow toward both the inland and the creek boundaries. Strong groundwater seepage into the creek developed near the low tide mark on the creek bank and through the creek bed. The discharge to the creek results in an aerated zone near the creek (0 m < x < 20 m, Figure 7c). Surface soil saturation (>0.8) was highest in the lower intertidal zone and smaller in the near-creek and upper intertidal zone. However, the actual evaporation rate, estimated by Equation 1 that considers the impacts of surface saturation, salinity, and salt precipitation, was similar across the intertidal zone (∼3.8 mm/day, except the nearcreek aerated zone) as the surface moisture level is sufficient to maintain a high evaporation rate (saturation is overall greater than 0.4). During the neap tide period, evaporation occurred throughout the intertidal surface soil as the disruption from tidal inundation was absent. Surface soil saturation decreased significantly with the lowest saturation decreasing to 0.2, leading to the decline of evaporation rate in the lower and upper intertidal zones as shown in Figures 7c and 7d. The actual evaporation rate in the supratidal zone (x > 120 m) stayed close to zero all the time, and a permanent unsaturated zone with a depth up to 70 cm was formed. The deep groundwater table under the supratidal soil surface and the weak capillary forces of the loamy sand soil could not provide sufficient hydraulic connections between the groundwater table and soil surface to maintain evaporation (Shen et al., 2018).

Simulated Porewater Salinity Distribution
Tidally averaged porewater salinity contour at the quasi-steady state of Case 1 is shown in Figure 8. Salinity distribution development during the simulation can be found in Figure S4 in Supporting Information S1. Strong porewater salinization occurred due to evaporation (comparing the salinity level at 1 year and the end of the simulation, Figure S4 in Supporting Information S1). Distinct salinity zones were developed along the model domain under the influence of tide and evaporation. A hypersaline porewater plume was formed in the upper intertidal aquifer, and the surface salinity reached the highest value of 152 ppt, which is close to the highest field-measured surface salinity of 153.5 ppt. The downward flow lines in the higher intertidal zone indicated the migration of salt into the deeper aquifer, resulting in the elevated groundwater salinity to around 120 ppt. Hypersaline groundwater was then transported back to the tidal creek by the tide-and density-induced groundwater circulation. Meanwhile, the hypersaline plume was partially transported toward the inland boundary due to the hydraulic gradient along the transect by hydrodynamic dispersion and the density-driven flow.
A small pocket of less saline porewater formed between the adjacent hypersaline zones around the high tide (100 m < x < 118 m, Figure 8) line near the surface. This resulted from the short inundation time of the spring high tide (3.4 hr per 15 days) and the weak capillary forces, which neither supplied enough salt nor maintained In the near-creek intertidal zone, a lower porewater salinity "wedge" (∼45 ppt) was formed and extended 5 m below the surface due to the rapid surface water and groundwater exchange (the 45-day flow line, Figure 8). The salinity transition zone from the hypersaline plume to the near-creek "wedge" was located in the lower intertidal zone, where salt accumulation by evaporation is weaker than salt removal by tide flushing. Under the  "wedge," a salt discharging pathway developed with an increased porewater salinity, transporting the accumulated salt from the hypersaline plume in the higher intertidal zone to the creek through the creek bed (see the 409-day and 1624-day flow lines, Figure 8). Similar to salt transport to the inland boundary, this hypersaline porewater discharging pathway was formed due to the creek-normal hydraulic gradient and density contrast in the soil profile. The simulated hypersaline plume, and particularly the high-salinity groundwater underneath the creek, agrees well with the field measured salinity distribution patterns shown in Figure 5a.

Saltwater and Solute Exchange Across the Creek Boundary
The temporal variations of saltwater flux and salt exchange across creek bank and creek bed were examined as shown in Figure 9. Overall, both the saltwater and salt flux corresponded with the tide signal as recharge occurred during high tide and discharge during low tide. The saltwater recharge (positive saltwater flux) occurred predominantly at the creek bank with its peak flux dictated by the tidal amplitude (see the red-dashed line, Figure 9b). Only a limited recharge occurred through the creek bed even during the spring high tide period. On the contrary, saltwater discharge occurred simultaneously through the creek bank and creek bed with a similar maximum rate of ∼2,500 kg/day. The cumulative saltwater and salt fluxes were determined by integrating their corresponding rate over time, respectively. The cumulative saltwater fluxes suggest discharge dominated at both the creek bank and creek bed, while the cumulative discharge through the creek bank almost doubled that through the creek bed (Figure 9c). Similar to saltwater recharge, salt inflow to the soil mainly occurred along the creek bank as it is driven by the saltwater recharge across the creek boundary (Figures 8b and 8d). The salt outflow from the creek bank was comparable with that from the creek bed, while the flux peak from the creek bed was slightly higher than that from the creek bank, a phenomenon that was opposite to the saltwater discharge. Despite that the cumulative saltwater discharge across the creek bank doubled that on the creek bed (Figure 9c), the cumulative salt outflow through the creek bed was almost identical to (even slightly overtook) that along the creek bank (Figure 9e). The significant salt outflow on the creek bed was due to the high-salinity porewater in the creek hyporheic zone (beneath the creek bed), resulting from the density-driven discharge of the hypersaline plume. Figure 10a represents the quasi-steady salinity distribution obtained by Case 2 under a higher potential evaporation rate at 8 mm/day. The overall salinity spatial distribution pattern by Case 2 was similar to that by the base case (Case 1) with the highest salinity plume developed in the upper intertidal zone and a lower salinity "wedge" near the creek. The higher evaporation rate resulted in elevated salinity, reaching almost 160 ppt in the deep upper intertidal aquifer and more than 170 ppt on the surface near the high tide mark. The near-creek "wedge" zone (0 m < x < 20 m) salinity also increased around 33% (∼45 ppt in Case 1-∼60 ppt in Case 2) despite the unchanged tidal regime. The deep groundwater was strongly salinized compared to that from Case 1 (maximum groundwater salinity in Case 1 was around 120 ppt), further confirming the downward transport of the accumulated salt into the aquifer. This thus led to the elevated salinity of the groundwater under the tidal creek, promoting the salt outflow to the creek. Inland salinization was also intensified under higher evaporation. The existence of an upland groundwater head, however, prevented the complete salinization of the inland soil.

Hypersalinity Under Higher Evaporation Rate and Rainfall
Case 3 simulated the salinity change after a 12-hr rain event, following the quasi-steady state achieved in the base case (Case 1). Evident porewater salinity change was observed in the unsaturated zone particularly in the supratidal area, where the pre-rainfall high salinity formed by dispersion and diffusion (Figure 8) got diluted by rainfall recharge with salinity reduced from ∼85 to ∼10 ppt (Figure 10b). Similar salinity change also occurred in the upper intertidal hypersaline plume near the surface with the highest salinity reduced from 152 to 100 ppt. Overall, rainfall resulted in strong freshening of the shallow porewater near the surface, the ecologically important root zones of the wetland. The hypersaline groundwater was however not impacted by this short-term rain event. These findings agreed with the field porewater salinity distribution after rain events shown in Figure 5b.

Discussion
Seawater-groundwater interaction and salt exchange in a tidal wetland are more influenced by evaporation than coastal beach setting, which tends to result in salt accumulation. The field groundwater EC, measured porewater salinity, and numerical modeling all identified the existence of hypersaline plumes at this wetland. Due to the density gradient between the hypersaline porewater and seawater, and the tidally driven porewater circulation, the heavier hypersaline plumes were discharged toward the creek via the creek bank and creek bed. Such a transport pathway provides a mechanism of salt removal from the intertidal zone. Without an efficient salt removal mechanism, evaporation is likely to cause the formation of permanent salt crusts, the over-salinization of the underlying root zone, and the consequent reduction of wetland vegetation coverage (Hughes et al., 2012). Unvegetated areas are a phenomenon that is present but not widespread across the wetland surface at our site. The long-term simulation of Case 1 and Case 2 in this study both showed that the salt distribution reaches a steady-state condition with a maximum porewater salinity slightly lower than the solubility limit (∼265 ppt).
Both field measurements and numerical simulation indicated the insignificant formation of the salt crust. However, the salinity of deep groundwater estimated by the model was higher than that in the field (i.e., ∼110 ppt in the base case model vs. ∼70 ppt in-situ measurement). Such discrepancy is more likely linked to the freshwater recharge, for example, from rainfall, which was not considered in the long-term numerical simulation. Simulation Case 3 and field measurement after two rainfall events confirmed that the shallow porewater was diluted evidently in the intertidal zone and particularly in the supratidal zone (Figures 6b and Figure 10b). Though the numerical simulation only attempted to simulate the influence of a short-term rain event, regular seasonal rainfall that occurs throughout the year during the formation of the hypersaline condition should be capable of diluting or removing (via runoff) salt in shallow soils, thus leading to less intense salinization of the deep groundwater that receives the salt from above by tidal recharge and density-driven flow. On the other hand, the surface water Figure 10. (a) Tidally averaged salinity distribution simulated by Case 2 under a potential evaporation rate of 8 mm/day for 22.75 years (as compared to the base case with a potential evaporation rate of 4 mm/day) and (b) salinity distribution simulated by Case 3 after a 12-hr rain event based on the steady salinity distribution from the base Case 1. The solid black line represents the groundwater and tidal water level, and the dashed lines represent the high and low tide marks. salinity of estuarine tidal creeks can be greatly impacted by rainfall. Rainfall-induced creek water salinity drop has been reported by other studies (Lenkopane et al., 2009;Sadhuram et al., 2005;Xiao, Li, Xia, et al., 2019) and in this study site (see Figure 4a). While the numerical models in this study regarded the creek water salinity as being constant, simulations conducted by Lenkopane et al. (2009) and Xiao, Li, Xia, et al. (2019) both suggested that the intertidal porewater, particularly the porewater under tidal creek water, could be diluted when the creek water salinity is reduced. The freshened creek water, either varied by the tide or influenced by rainfall, could be an important factor in reducing wetland salinity. Nevertheless, the field measurements identified high-salinity (determined by chloride concentration) porewater under the creek bed support the modeled discharging of the hypersaline porewater.
The formation of the hypersaline porewater zones in the coastal wetland is primarily driven by evaporation. Case 2 and Case 3 simulations with evaporation implemented through the whole simulation period both resulted in strong porewater hypersalinization. High evaporation rate such as 4-8 mm/day is common in the semidry to dry coastal regions and during the summer dry seasons in the subtropical coastal areas (BOM, 2021). Hypersaline porewater condition similar to this study has been more often documented in the low-latitude regions under the similar climate, such as the Argentine Atlantic coast (Galliari et al., 2021) in South America, the subtropical coastal marshes in the Gulf of Mexico (Borchert et al., 2018;Forbes & Dunton, 2006), the South Atlantic coastal areas in South Africa (Bornman & Adams, 2010), and the subtropical to tropical regions in Northern Australia (Duke et al., 2019;Ridd & Sam, 1996;Ridd & Stieglitz, 2002). Those areas typically have high evaporation rates and very limited or seasonal rainfall, which mostly only occur in the wet season months (∼3 months per year) (Albuquerque et al., 2014;Ridd & Sam, 1996), thus allowing the formation of hypersaline porewater zones in the wetland by evaporation during the dry period as shown by the simulations herein. However, most of the previous studies, particularly the numerical modeling on tidal wetland salt transport dynamics, were based on the tidal fresh-brackish wetland with the fresh groundwater sustained at uplands. This is a common condition for the coastal wetlands in middle-to high-latitude regions, such as the wetlands along the Atlantic and Alaska coast, North America (e.g., , Guimond, Yu et al., 2020Kelsey et al., 2016;Montalto et al., 2006), and the regions with reliably regular rainfall through the year (e.g., Cao et al., 2012;Xiao, Li, Shananan, et al., 2019), as the salt accumulation by evapoconcentration is either weak or can be removed by incoming freshwater. Based on the distinctive and significant hypersalinization process formed by evaporation in this study, it is important for future studies to include the evaporation effects should they conduct investigation in regions under a similar climate as this site.
Enhanced evaporation is also expected in the future under a warmer climate (Konapala et al., 2020;Wang et al., 2018). Case 3 simulation, with evaporation rate increased to 8 mm/day, suggested that salinization could be further intensified not only in the upper intertidal aquifer but also toward the inland and near the creek. Under the future warmer climate, enhanced evaporation might be another important factor causing coastal area groundwater salinization on top of the sea level rise influences, particularly in the low-lying coastal wetland areas with the limited inland groundwater input. In the near-creek zone, the 33% increase in porewater salinity due to the rise of potential evaporation rate from 4 mm/day to 8 mm/day suggests that salt removal by current tidal flushing may not be capable of offsetting the salt accumulation by higher evaporation. The salinity increase in those areas might result in plants' zonation pattern change as the local plants have to adapt to the environmental stress from a higher salinity level (Moffett & Gorelick, 2016).
Due to the formation of the hypersaline porewater discharge zone under the creek, the salinity at the bottom of the creek may be much higher than the overlying seawater, which could influence the local ecosystem. Peterson et al. (2016) presented a similar downward salinity stratified system in Long bay, South Carolina. They identified the occurrence of hypersaline and anoxic SGD from the carbonate aquifer under the seafloor in their field site. As the result of the density contrast between these hypersaline SGD water bodies and the ambient ocean water, a salty water mass layer on the seafloor had been formed and further inhibited oceanic water mixing. This phenomenon accounts for the coastal hypoxia events (severe DO depletion in coastal water) in their study site. Likewise, the mixing of discharging hypersaline water and the overlying seawater may also affect the biota in the hyporheic zone of the creek, an active and important ecotone (Boulton et al., 1998). Future study is encouraged to explore the correlation between the creek water quality and stratification and the accumulation of salts in the adjacent coastal wetland system.

Conclusions
Previous studies of porewater flow and salt transport in the tidal wetland mostly ignored the evaporation process, which could significantly impact the porewater salinity distribution and alter local porewater flow. This study monitored the groundwater head and porewater salinity variations along a representative cross-creek transect in a subtropical tidal wetland. Numerical simulation of porewater flow and solute transport, considering variable saturation, variable density, tidal fluctuation, and evaporation-induced salt accumulation was also conducted.
Field measurements and long-term simulation confirmed the formation of hypersaline porewater plumes under the high tide line. The hydraulic connections between the groundwater table and wetland surface, which support evaporation, are key to the formation of the hypersaline plumes. The plumes further transported into the deeper aquifer and caused aquifer salinization in the upper intertidal zone. Due to tide fluctuation-induced porewater circulation and density gradients, the hypersaline porewater was discharged toward the tidal creek via a salt exporting path formed under the fresher saltwater "wedge" (with salinity close to that of seawater) in the nearcreek intertidal zone. This provides a mechanism of salt removal from the upper intertidal zone moderating salt flat formation. Quantitative analysis on salt flux across the creek boundary revealed that the creek overall is receiving a large amount of salt that has accumulated in the wetland soil by evaporation. The creek bed discharged a comparable amount of salt as the creek bank despite that most porewater exchange on the creek boundary occurred along the latter. The formation and discharge of the hypersaline porewater caused salinity stratification in the creek hyporheic zone and the near-creek intertidal aquifer.
The evaporation-induced hypersaline condition was intensified under a higher evaporation rate, leading to stronger landward groundwater salinization and higher salinity near the tidal creek. The salt outflow to the creek was thus increased as the result of a higher salinity level in the creek hyporheic zone. Despite the increased hypersaline plume salinity under a higher evaporation rate, no salt crust developed in the simulation due to the existence of the saline water discharge pathway. Field measurements on post-rainfall porewater salinity and the simulation results of a short rain event confirmed the diluting influence of rainfall, forming a thin layer of freshened porewater in the root zone while the deep groundwater still stayed in hypersaline condition. However, further future modeling and field studies are needed to investigate the effects of regular seasonal rainfall on modifying the porewater flow and solute transport dynamics in the tidal wetland system. These findings highlight the necessity of considering evaporation in the low-lying coastal wetland environment, particularly in the subtropical, and semiarid to arid areas. Porewater salinity in the creek hyporheic zone and the near-creek area could show distinct variation, and thus future field measurements should be designed to cover these areas at a higher resolution.

Data Availability Statement
Field data sets and SUTRASET model used to support this paper are available at http://www.hydroshare.org/ resource/fd4c5288e8864148bb31c70a4fe6cf0b. support from Don Gynther. The authors are grateful to the three reviewers and the editors for their constructive comments that have led to significant improvement of this paper. Open access publishing facilitated by The University of Queensland, as part of the Wiley -The University of Queensland agreement via the Council of Australian University Librarians.