Impacts of Subsurface Dam Construction on Downstream Groundwater Levels and Salinity in Coastal Aquifers

Subsurface dam is a promising engineering technology for groundwater resources development. However, the possible impacts of these dams on the groundwater environment have been a major concern. Here, we used a three‐dimensional (3D), variable‐density, unsaturated‐saturated groundwater flow model to explore how a groundwater‐storage‐type subsurface dam, built in the freshwater domain of an unconfined coastal aquifer, affected groundwater levels and salinity in the downstream area. Model results suggested that, after subsurface dam construction, groundwater levels in the downstream area showed intensified fluctuations in terms of phase advances, greater amplitudes, and higher frequencies following heavy rainfall events. Numerical simulations with variable subsurface dam scenarios indicated that the groundwater level fluctuations were further intensified with a higher crest elevation or a shorter distance from the coast. Moreover, during the recharging period of the subsurface reservoir, sea water in the downstream area intruded landward from its initial location, which can at least temporarily threaten water quality near the coast. A higher dam crest elevation prolonged the duration of sea water intrusion, while a dam positioned closer to the coast induced sea water intrusion with a greater horizontal extent. General implications are discussed with respect to improving assessment methodologies and engineering designs of subsurface dams.


Introduction
Subsurface physical barriers are a promising engineering solution for groundwater resources development.The barriers can be generally classified into two types according to their purpose: groundwater storage and sea water intrusion control (Japan Green Resources Corporation 2004).The former is designed to raise groundwater levels so that the impounded groundwater can provide a stable water supply, and the latter is constructed in coastal aquifers to prevent lateral sea water intrusion when excessive groundwater resources are abstracted.There are usually three different barrier structures (Chang et al. 2021(Chang et al. , 2022;;Zheng et al. 2022): (1) a "subsurface dam" in contact with the aquifer basement forcing groundwater to discharge over the dam crest; (2) a "cutoff wall" embedded into the upper part of the aquifer allowing groundwater to discharge through the opening below the wall bottom; and (3) a "fully penetrating barrier" penetrated the entire aquifer depth.
Subsurface dam is a widely used structure for constructing subsurface physical barriers in coastal aquifers.Many dams were built nearby or within the freshwater-sea water mixing zone with their crest elevations higher than sea level, thereby achieving both groundwater storage and sea water intrusion control (Nishigaki et al. 2004;Yoshimoto et al. 2013).Some groundwater storage dams were constructed completely in the freshwater domain of coastal aquifers, examples of which include the ones in India (Senthilkumar and Elango 2011), Croatia (Roje-Bonacci and Bonacci 2013), China (Kang and Xu 2017), and the Ryukyu Islands of Japan (Japan Green Resources Corporation 2004) including Miyakojima Island which was studied in this paper as an example.
Despite many benefits of subsurface dams, a major concern is the potential negative impacts on the groundwater environment, especially in the upstream area.For example, subsurface dams cutoff natural groundwater flow, which may cause the accumulation of anthropogenic contaminants like NO 3 -N in the subsurface reservoir (Ishida et al. 2006;Yoshimoto et al. 2013;Sun et al. 2019).After building a subsurface dam in the freshwater-sea water mixing zone of coastal aquifers, residual sea water trapped in the upstream area may pose a long-term threat to groundwater quality (Nawa and Miyazaki 2009;Abdoulhalik et al. 2017;Zheng et al. 2021).Moreover, elevated groundwater levels in the upstream area may induce groundwater flood problems during heavy rainfall events (Milanović 2004;Ishida et al. 2011;Yoshimoto et al. 2013).Nevertheless, some researchers started to recognize the possible disturbances in downstream area.For example, numerical models predicted that the subsurface dam planned in a coastal aquifer of the Palar River, India, would lead to a decrease in downstream groundwater levels of 0.1 to 0.2 m, enhancing the risks of sea water intrusion during dry seasons (Senthilkumar and Elango 2011).Laboratory and numerical experiments discovered that constructing subsurface dams in the freshwater-sea water mixing zone can induce upward sea water intrusion in the downstream area (Chang et al. 2020).
Subsurface dam technology is highly dependent on regional geological conditions.Geological surveys are conducted to collect information such as hydraulic conductivity, effective porosity, stratigraphy, and bedrock depth, which are all necessary for selecting appropriate construction methods and locations (Japan Green Resources Corporation 2004).However, numerical studies of subsurface dams often adopt a two-dimensional (2D) model and assume a flat aquifer bottom, which may not always be representative of real world situations.In fact, three-dimensional (3D) models are necessary to account for the effects of subsurface dams on the radial flow fields and the overall water budget of the aquifer (Kaleris and Ziogas 2013;Wu et al. 2020).Especially, groundwaterstorage-type subsurface dams are designed to effectively raise groundwater levels, for example, by over ten or tens of meters in unconfined aquifers (Japan Green Resources Corporation 2004).How and to what degree such groundwater level rises alter groundwater flow conditions depend on the 3D aquifer geometry.
The dam height and location are two key design criteria.Several laboratory and numerical experiments have demonstrated that a subsurface dam should be designed to be higher than a critical height at a specified position to achieve effective control of sea water intrusion (Luyun et al. 2009;Chang et al. 2019;Zheng et al. 2021).However, a higher dam height would not only increase the construction cost, but also reduce fresh groundwater discharge in a head-controlled aquifer, which favors the accumulation of land-based pollutants in the subsurface reservoir (Chang et al. 2019;Sun et al. 2019;Ke et al. 2021).These studies contribute valuable insights with respect to improving the environmental performance of subsurface dams.However, how the dam height and position affect environmental changes that occur in the downstream area is so far not well understood.
Unlike many studies focusing on the upstream impacts, this work explores how subsurface dam construction affects downstream groundwater conditions.The Sunagawa Subsurface Dam, a groundwater-storage-type dam built in the freshwater domain of the coastal aquifer in Miyakojima Island, Japan, was studied as a realistic example.Numerical models were developed to solve variable-density, unsaturated-saturated groundwater flow and mass transport processes taking into account the 3D aquifer geometry, which has been often oversimplified in previous modeling work.Model results demonstrated the possible impacts on downstream groundwater level fluctuations and salinity distributions and showed how these impacts were controlled by different dam crest elevations and positions.Implications were attained for improving assessment methodologies and design criteria to achieve the environment-friendly usage of subsurface dam technology.

Site Description
Miyakojima Island is located in the southern Ryukyu Islands (24 • 40 N and 125 • 20 E, Figure 1).The island has a land area of approximately 159 km 2 with its northsouth distance of 21 km and west-east distance of 23 km.This subtropical island has a warm and humid climate.The average annual temperature was 23.5 • C, and the highest and lowest monthly average temperatures were 28.7 • C in July and 18.0 • C in January, respectively (Japan Meteorological Agency 2022).The average annual rainfall was 1970 mm with the two rainy months in May and September, and the dry seasons in August and winter (Japan Meteorological Agency 2022).
Miyakojima Island and the other Ryukyu Islands are located in a tectonically active zone where the Philippine Sea Plate is being subducted beneath the Eurasian plate.This leads to the uplift of coral reef deposits called the Ryukyu limestone.Miyakojima Island is composed of several levels of the uplifted Ryukyu limestone terraces.These terraces form a moderately flat topography of the island, where the highest elevation is 114.5 m amsl (above mean sea level) (Yamanaka 2012).
As shown in Figure 1, the geology is mainly characterized by the Pleistocene Ryukyu limestone with a thickness of 20 to 50 m, overlying the Late Miocene-Pliocene siltstone and mudstone of the Shimajiri Group (Yazaki and Oyama 1984).Some surface depressions are filled with 0 to 10 m thick Holocene Ohnokoshi clay which has a hydraulic conductivity of approximately 10 −6 m/s (Mori et al. 1997;Yamanaka 2012).The Shimajiri Group has a thickness of approximately 2000 m with a hydraulic conductivity of less than 10 −7 m/s (Mori et al. 1997), which acts as a hydrogeological bedrock.A highly permeable unconfined aquifer is found in the Ryukyu limestone with a heterogeneous hydraulic conductivity reported in the range of 10 −5 to 10 −2 m/s and a porosity of 0.1 to 0.36 (Imaizumi et al. 1988).In some parts of the aquifer, clay materials flowed into and filled the pores of the limestone, which reduces hydraulic conductivity and porosity (Mori et al. 1997;Ishida et al. 2015).
Northwest-southeast oriented faults are located in the island (Figure 1).As shown in the W-E crosssection (Figure 1), the cliffs of bedrock along the faults act as vertical low-permeability boundaries dividing the island into many groundwater basins.Along the southern coast of Miyakojima Island, the bedrock topography has several valley-shaped structures between the faults, and these bedrock valleys are suitable sites for constructing the groundwater-storage-type subsurface dams (Ishida et al. 2003).
Groundwater is the only source of water supply for Miyakojima Island.In 1988, the Japanese government proposed the construction of subsurface dams in the southeastern parts of some basins (Figure 1) in order to harvest groundwater for agricultural uses (Ishida et al. 2006).Among these dams, the Sunagawa Subsurface Dam was the first mega-subsurface dam, which has a length of 1677 m, a height of 50 m, and a crest elevation of 31 m amsl (Ishida et al. 2006(Ishida et al. , 2015)).The dam was embedded into the cliff of bedrock along the fault in the west and into the bedrock in the bottom and east.
The dam was constructed mainly by the cast-in-place deep mixing method, and the hydraulic conductivity of the dam material was lower than 5 × 10 −7 m/s (Japan Green Resources Corporation 2004).The construction of the central part of the dam was completed on November 22, 1993, the full construction was completed on March 31, 1994, and the subsurface reservoir became fully recharged on September 30, 1995 (Mori et al. 1997).Since 2001, groundwater has been abstracted from 147 tube wells to be used for agricultural purposes (Ishida et al. 2006).

Numerical Modeling
The 3D groundwater models were developed using the numerical code FEFLOW (Diersch 2014) which solves density-dependent groundwater flow and mass transport processes in unsaturated-saturated porous media.For concentration-related density-dependent groundwater flow, FEFLOW assumed a linear relationship between the solute concentration and the fluid density.Variably saturated flow was simulated using Richards' equation.Density-dependent salt transport processes were modeled by solving the advection-diffusion-dispersion equation.The corresponding governing equations are more thoroughly explained by Diersch (2014) and are briefly provided in Appendix S1 in Supporting information.
The 3D model domain covers the basin of the Sunagawa Subsurface Dam, including both the land and seafloor areas (Figure 1).No-flow boundaries were assumed along the two faults (A-B and C-D) and the bedrock ridge (A-D).After several trials, the boundary at the seafloor (B-C) was decided at a horizontal distance of approximately 1000 m from the coastline, which was confirmed to be sufficient for the model to capture all important freshwater-sea water interactions throughout the whole simulation stages.The model top was set to the land surface (Geospatial Information Authority of Japan 2022) and the seafloor (Yazaki and Oyama 1984) and the model bottom was set to the bedrock surface obtained from Mori et al. (1997) acting as a no-flow boundary.
A uniform rainfall recharge flux q f (L/T) with a constant mass concentration C = 0 kg/m 3 was assigned at the land surface.The daily input values of q f were calculated based on daily rainfall data (Japan Meteorological Agency 2022) and seasonal actual evapotranspiration data measured by Hiyane et al. (2004).The calculated q f resulted in an average groundwater recharge rate equivalent to 49% of the annual rainfalls, which is close to the ratio of 45% estimated by Fazal et al. (2005).At the seafloor and the boundary B-C (Figure 1), a constant sea water head h sw (L) of 0 m amsl and a constant salinity C = 35 kg/m 3 (NaCl) was applied.The sea water head was converted to equivalent freshwater head.Strictly speaking, these boundary conditions did not perfectly represent the real condition near the coast.Indeed, a system-dependent strategy of boundary setup was required, that is, at the seafloor, C = 35 kg/m 3 should be set for sea water inflow from the sea to the aquifer while no constraints should be applied to groundwater discharge from the aquifer to the sea, and at the beach surface, rainfall recharge should be applied except for the locations where groundwater discharges out.However, this complex setting was not adopted in this 3D model because it caused serious numerical oscillations and tremendously raised the computation time.In fact, our vertical 2D model confirmed that the system-dependent boundary strategy showed limited influences on the modeled salinity distribution (Appendix S2 in Supporting information).
The model was first run with a constant average rainfall recharge (867 mm/yr) for a period of 200 years which ensured a quasi-steady state.The obtained quasisteady state was used as the initial conditions for the subsequent simulations.Then, the transient simulations were divided into three stages: (1) January 1, 1986, to December 31, 1987, was the model warm-up period; (2) January 1, 1988, to November 21, 1993, was the period prior to subsurface dam construction; and (3) November 22, 1993, to December 31, 2000, was the period of after subsurface dam construction.The third stage began on November 22, 1993 when the construction of the central part of the dam was completed.The dam construction period (November 22, 1993to March 31, 1994) was not included in the model because this period was relatively short and the influence on the model results was sufficiently small for our modeling analysis (Appendix S3 in Supporting information).Over the entire simulation period, groundwater abstraction activities were not modeled because the amount of pumping was negligible (Fazal et al. 2005).Major groundwater abstraction activities began in 2001 (Ishida et al. 2006) which is after the simulation period.
The model was discretized into a total of 131,328 triangular-based prismatic elements with 18 layers in the vertical direction.The layer thickness increased with depth and was in a range of 0.1 to 6.7 m for the subsurface reservoir area and 0.1 to 3.1 m for the sea floor area.The average element size in the landward side is 205 m.For the sake of avoiding numerical oscillations and improving the accuracy of modeling results, a higher degree of refinement (element size l ≈ 10 m) was adopted in the coastal zone (Figure 2a) to ensure that the Peclet number Pe m ≈ ( l /α L ) < 4 (Voss and Souza 1987) where α L (L) is the longitudinal dispersivity.The automatic time-step control with the first-order predictor-corrector methods was applied for the temporal discretization (Diersch 2014).The initial and maximum time step lengths were set to be 1 × 10 −6 and 0.01 days, respectively, and the growth factor between subsequent time steps was restricted to be 1.1.
Based on the available geological information (Yazaki and Oyama 1984;Mori et al. 1997;Ishida et al. 2015), in some parts of the aquifer, clay materials flowed into and filled the pores of the limestone, which reduced hydraulic conductivity K (L/T) and effective porosity ε (−).Therefore, the model domain was divided into three geologic units: (1) clay; (2) limestone with clay; and (3) limestone without clay (Figure 2b).The materials of each unit were assumed to be homogeneous and isotropic.The values of K and ε for each geologic unit were manually calibrated (on a trial-and-error basis) against the measured groundwater levels between 1988 and 1996 (Mori et al. 1997).The calibrated values (Table 1) generally agreed with the reported ranges except for the calibrated values of ε for limestone with clay (ε = 0.2) and limestone without clay (ε = 0.4) which were larger than reported ones of the Ryukyu limestone ranging in 0.1 to 0.15 (Ishida et al. 2006; Tokashiki and Aydan 2011).This may be attributed to two factors, the uncertainty in the bedrock geometry data and the fact that some cavities have been developed in the limestone.Both factors influenced the water level changes following rainfalls and subsurface dam construction.The details of the calibration process and parameter sensitivity are presented in Appendix S4 in the Supporting information.As shown in Figure 3, the calibration model results generally matched the overall spatial and temporal distributions of groundwater levels well.Hence, the flow model was considered well calibrated.van Genuchten model (van Genuchten 1980) was adopted (Figure 2c) and the typical unsaturated parameter values of each geologic unit were obtained from The Unsaturated Soil Hydraulic Database (Nemes et al. 2015) and were listed in Table 1.For the dam materials, K was set to be 1 × 10 −7 m/s (Japan Green Resources Corporation 2004).The longitudinal dispersivity α L (L) and the transverse dispersivity α T (L) were set to 10 and 1 m, respectively (Gelhar et al. 1992).
Other parameters were based on typical literature values  (Diersch 2014): the density ratio of sea water and freshwater χ s (-) was 0.025; the dynamic viscosity μ (M/[L T]) was 0.001 kg/(m s); and the molecular diffusion coefficient D 0 (L 2 /T) was set to be 1 × 10 −9 m 2 /s.

Variable Subsurface Dam Scenarios
Table 2 lists the scenarios with different configurations of the dam.For comparison, scenario S0 assumed a natural condition without subsurface dam construction.S1 mimicked the realistic configurations of the Sunagawa Subsurface Dam.S2 to S7 represented different choices of the dam crest elevation E dam (L) and the horizontal distance between the dam and the coastline L dam (L).L dam was measured along the NW-SE dashed line in Figure 1.The dam locations with L dam of 500, 1500, and 2500 m (Figure 2) were all within the freshwater domain of the coastal aquifer.When L dam was set to be 500 m, E dam of 31 m amsl or 40 m amsl was not included in the scenarios because they were higher than the land surface elevation at the dam.For all scenarios, the same rainfall recharge was applied because the groundwater table is sufficiently deep below the ground surface for all scenarios.

Evaluation Indicators of Sea Water Intrusion
Lateral sea water intrusion induced by subsurface dam construction can be quantitatively evaluated based on the toe position of the freshwater-sea water mixing zone.Here, the toe position (L toe [L]) is defined as the horizontal distance from the coastline to the sea water toe measured along the NW-SE dashed line in Figure 1.The toe refers to the location of 50% sea water salinity at the aquifer bottom.Our preliminary simulations indicated that the lateral sea water intrusion occurred since the time of subsurface dam construction (t dam [T]) and reached its maximum level at the time when the subsurface reservoir became fully recharged (t full [T]).Hence, the extent of lateral sea water intrusion can be quantified by the horizontal movement of the toe position L toe (L) between t dam and t full .In addition, the duration of sea water intrusion T SWI (T) was defined as the period from t dam until the induced sea water intrusion was fully recovered (L toe of S * became the same with S0).Examples of L toe and T SWI are illustrated in Figure 7.

Impacts on Downstream Groundwater Level Fluctuations
Figure 3 shows the temporal changes in groundwater levels obtained from observation wells (Mori et al. 1997) and the model results for the construction of the Sunagawa Subsurface Dam (scenario S1).During the recharging period of the subsurface reservoir (from t dam to t full ), groundwater levels gradually increased and exceeded the level of E dam in the upstream area (Wells 1 to 4 in Figure 3b1 through 3b4), while dropped 2 m in the downstream area (Well 5 in Figure 3b5 and 3b6).This was caused by the fact that groundwater flow was temporarily cut off by the fill-up of the subsurface reservoir.After t full , groundwater levels in the downstream area started to recover as groundwater flowed over the dam crest (Figure 3b6).Our models indicated that this temporal reduction of groundwater levels in the downstream area caused lateral sea water intrusion, and the details of which are presented in the next section.
Figure 4 showed how subsurface dam construction influenced the groundwater level fluctuations responded to heavy rainfalls (rainfalls greater than 0.1 m/d).These influences included three aspects: the phase, amplitude, and frequency of groundwater level fluctuations.
First, after heavy rainfall events, groundwater levels in both the upstream and downstream areas showed earlier peaks in S1 than in S0 (Figure 4).Such a trend was also captured in Figure 5a which shows that the steadystate travel time was generally shorter in S1 compared with that in S0.This is related to the fact that elevated groundwater levels in S1 reduced the vertical thickness of the unsaturated zone in the upstream area (Figure 5c), thus shortening the travel time from the land surface to the water table.This also caused faster dam overflow and resulted in faster discharge from the upstream to downstream areas.Note that the travel time shown in Figure 5a was calculated by conducting steady-state simulations of S0 and S1, respectively, based on a constant recharge rate of 867 mm/yr.Therefore, the steady-state travel time indicates the difference in the flow behavior between S0 and S1 on an annual average basis but instead of the groundwater age or travel time under certain rainfall events shown in Figure 4a.
Second, the amplitude of groundwater level fluctuations in the upstream area was generally lower in S1 than in S0 (Figure 4b), while it was opposite in the downstream area (Figure 4c).This phenomenon is mainly related to the valley-shaped structure of the bedrock surface (Figure 1).In the upstream area, as groundwater levels increased in the case of S1, the water table expanded its surface area in the horizontal direction compared with that in S0 (Figure 5b).On the other hand, in the downstream area, the amplitude of groundwater level fluctuations was strongly affected by the flow rate of groundwater discharging from the upstream area as explained below.In the upstream area in the case of S0, a significant portion of the bedrock surface beneath the eastern part of the aquifer was exposed to the unsaturated zone (Figure 5b left), resulting in two distinct traveling routes for rainfall infiltration.One is a direct route from the land surface downwards to the water table while the other is a rather longer route consisting of a downward flow and then a lateral flow along the slope of the bedrock surface (Figure 5c left).This resulted in a wide distribution of the time required for rainfalls to infiltrate through the unsaturated zone in the upstream area and then travel to the downstream area (e.g., see the distribution of steady-state travel time at the SW-NE cross-section in Figure 5a left).Therefore, the discharge rate of groundwater to the downstream area in S0 was smoothed out.On the contrary, in the case of S1, some parts of the originally unsaturated bedrock surface became saturated in the upstream area (Figure 5c right), thereby narrowing down the distribution of the travel time of rainfall infiltration to a lower range (Figure 5a right).This resulted in sharper peaks of the flow rate of groundwater to the downstream area and therefore induced larger amplitudes of groundwater level fluctuations in the downstream area in S1 (Figure 4c).
Third, more frequent fluctuations of groundwater levels were observed in S1 than in S0 in the downstream area.For example, groundwater levels showed clearer increases following heavy rainfalls in S1 than in S0 (Figure 4c).This is related to the elevated groundwater levels in S1 causing a lower and more narrowed distribution of the travel time of rainfall infiltration in the upstream area (Figure 5a right), and consequently, groundwater level fluctuations in the downstream area became more sensitive to heavy rainfall events.
These impacts in the downstream area varied under different subsurface dam scenarios.Comparing S3 with S1 (Figure 4c), higher E dam caused greater amplitudes and phase advances.This is because higher subsurface dams would further reduce the thickness of the unsaturated zone in the upstream area.When comparing S6 with S1, a larger L dam generally showed less disturbance to groundwater level fluctuations.This is because a larger L dam caused groundwater buildup within a smaller portion of the entire aquifer, and hence, the influences on the downstream area became smaller.

Impacts on Downstream Groundwater Salinity
Model results showed that after t dam , the freshwatersea water mixing zone showed landward movement.Specially, L toe advanced approximately from 90 m at t dam to 250 m at t full (Figure 6).Note that the Sunagawa Subsurface Dam was built completely in the freshwater domain without any direct contact with the freshwatersea water mixing zone.The sea water intrusion was caused by the significantly lowered groundwater levels in the downstream area of the dam.After t full , the intruded sea water began to retreat as groundwater flowed over the dam crest into the downstream area, which moved the freshwater-sea water mixing zone back to seaward.This sea water intrusion process is illustrated by the temporal changes in L toe in Figure 7. Compared with S0, S1 showed an increase of L toe since t dam , and L toe reached its maximum value of approximately 250 m ( L toe = 160 m) at t full , and then began to return to a level similar to that of S0.The modeled duration of this sea water intrusion ( T SWI ) was 925 days or approximately 2.5 years.After t full , L toe showed seasonal changes with larger amplitudes in S1 than those in S0.This is because the groundwater level fluctuations in the downstream area had larger amplitudes in S1 than those in S0 (Figure 4c).
Variable L dam scenarios showed an inverse relationship between L toe and L dam (Figure 8a).This is explained as during the recharging period of the groundwater reservoir, groundwater flow from the upstream was blocked, thus the groundwater levels in the downstream area were NGWA.orgmaintained only by rainfall recharge.A shorter L dam resulted in less rainfall recharge area in the downstream side, thus induced sea water intrusion over a larger landward distance.Figure 8b shows that a higher E dam caused a longer T SWI .This is because a higher E dam represented a larger storage capacity of the subsurface reservoir, of which the recharging period became prolonged.Moreover, T SWI was longer for shorter L dam (Figure 8b).This is because shorter L dam led to more significant sea water intrusion (Figure 8a) that required a longer time for its recovery.
Although our model results show that the sea water intrusion process did not threaten the subsurface reservoir, the risks of groundwater salinization in the downstream area cannot be ignored if groundwater nearshore is being used in other coastal sites.Furthermore, during the recharging period of the subsurface reservoir, reduced groundwater discharges and landward shifts of the sea water-freshwater mixing zone would disturb salinity distributions and water mixing patterns.These disturbances could lead to the degradation of coastal ecosystems (Werner et al. 2013).For example, many tropical and subtropical islands surrounded by rich but fragile coral reefs are extremely vulnerable to changes in the water environment induced by inland anthropogenic activities (Lubarsky et al. 2018;Luijendijk et al. 2020).Therefore, further research efforts are required to evaluate the potential linkages between inland subsurface dam construction and coastal ecosystem health.
Our model results showed that the sea water intrusion induced by subsurface dam construction was a temporal phenomenon and started to recover after t full (Figures 6  and 7b).This is mainly because the studied aquifer was a flux-controlled system (i.e., groundwater discharge was controlled by rainfall recharge), which kept the total freshwater discharges to the sea before t dam and after t full very similar (assuming no significant groundwater abstractions after t full ).However, in the case of a head-controlled system (i.e., groundwater discharge is controlled by large surface water bodies), subsurface dams would lower the hydraulic gradients in the upstream area and permanently decrease groundwater discharge to the sea (Chang et al. 2019;Sun et al. 2019;Ke et al. 2021).Under such situations, subsurface dam construction would cause irreversible sea water intrusion.
Our model results showed that the construction of the Sunagawa Subsurface Dam could have induced sea water intrusion over a distance of 160 m (Figure 6), which was not significant compared with the basin scale of approximately 5 km.However, this result is sitespecific.The relatively shallow bedrock topography is inclined towards the sea (Figures 1 and 6), limiting the horizontal extent of the freshwater-sea water mixing zone migration.In addition, the abundant rainfalls in the subtropical climate have maintained groundwater levels although groundwater flow was cut off by the Sunagawa Subsurface Dam during the recharging period of the groundwater reservoir.This may suggest that, for coastal aquifers where the bedrock is far deeper and the climate is much drier than Miyakojima Island, they could be more vulnerable to sea water intrusion induced by subsurface dam construction.

Roles of Bedrock Structures
The 3D bedrock structures have been often oversimplified in numerical models of subsurface dams.In fact, the valley-shaped bedrock structure has been generally regarded as a preferred geologic condition for subsurface dam construction that increases groundwater storage more effectively (Ishida et al. 2003;Nishigaki et al. 2004).Our model results showed that, after the construction of subsurface dams, part of the originally unsaturated bedrock became immersed by the elevated groundwater levels, which altered the rainfall infiltration paths through the unsaturated zone in the upstream area and resulted in severer groundwater level fluctuations in the downstream area (Figure 5).In addition, the relatively shallow bedrock surface near the coast in our study site may have limited the landward shift of the freshwater-sea water mixing zone, as discussed in the previous section.Therefore, the bedrock geometry should be fully evaluated due to its impacts on both fresh groundwater flow behavior and salinity distributions in aquifers with subsurface dams.

Implications for Environmental Monitoring and Subsurface Dam Design
Many studies demonstrate the importance of environmental monitoring programs for mitigating the potential degradation of groundwater quality upstream area subsurface dams (Ishida et al. 2006;Luyun et al. 2009;Abdoulhalik and Ahmed 2017;Sun et al. 2019;Zheng et al. 2021).Our model simulations indicate that subsurface dam construction affects the groundwater environment in the downstream area as well.Unfortunately, monitoring data downstream of subsurface dams are currently insufficient.Data from continuous monitoring of groundwater levels, salinity, pollutant concentrations, and ecosystem abundance in the downstream area before and long-term after subsurface dam construction would enhance our understanding of the potential impacts on the groundwater environment in the downstream area.
Designs of subsurface dams, including the dam geometry and the selection of dam location, are usually determined based on the required water storage capacity (Imaizumi et al. 1988), optimal effectiveness of sea water intrusion control (Luyun et al. 2009;Chang et al. 2019;Wu et al. 2020), and the minimization of accumulations of contaminants in the upstream area (Sun et al. 2019).Our model results revealed the impacts of the dam crest elevation and the dam location on groundwater level fluctuations and sea water intrusion in the downstream area.

Limitations and Future Research
Model calibrations in this study were conducted using only the groundwater level data from five monitoring wells, which resulted in uncertainties regarding the modeled groundwater levels and salinity in the study site.Unfortunately, no continuous monitoring data of groundwater salinity were available for the studied period, which prevented direct confirmation of the sea water intrusion processes observed in our models.To obtain a more comprehensive understanding of the impacts of subsurface dam construction, it is necessary to collect more abundant monitoring data in both the upstream and downstream areas.
Certain simplifications were made to boundary conditions and parameterization.For example, the tidal effects were not considered for the boundary condition along the seafloor, which likely affects the modeled salinity distributions.The material properties in each geologic unit were assumed homogeneous and isotropic, which could affect the modeled groundwater flow and freshwater-sea water mixing zone.Selected literature values were adopted for the unsaturated zone properties, which resulted in uncertainties of the water travel time analysis.Groundwater abstraction activities, which play key roles in resulting in sea water intrusion, have been widely studied for many coastal aquifers, including those conceptual aquifers with subsurface dams (Wu et al. 2020).Our studies focused on 1988 to 2000 which was before the major groundwater abstraction operation started in 2001.The model results showed that subsurface dam construction, before groundwater abstraction started, induced sea water intrusion.

Conclusions
While subsurface dams have been widely studied for their environmental impacts in the upstream area, this paper focused on the influence in the downstream area.Numerical models were developed for the Sunagawa Subsurface Dam built in the freshwater domain of the coastal aquifer in Miyakojima Island, Japan.The results showed that the subsurface dam caused intensified groundwater level fluctuations and sea water intrusion in the downstream area.Although these impacts do not likely have significant threats to the local water resource, observation data (groundwater levels, salinity, pollutant concentrations, and ecosystem abundance) in the downstream areas of other subsurface dam may add new insights into the sustainable application of subsurface dam technology.
The height and location are key design criteria of subsurface dams.Our models discovered that the groundwater level fluctuations in the downstream area were further intensified if the subsurface dam crest elevation was higher or if the dam was located closer to the coast.A higher dam crest elevation prolonged the duration of sea water intrusion, and a dam located closer to the coast induced a greater horizontal extent of sea water intrusion.
Numerical hydrological models of subsurface dams often assume a flat aquifer bottom, the 3D bedrock structures were found to play important roles in how subsurface dam construction affects rainfall infiltration routes and salinity distributions in this study.Therefore, it is recommended to collect reliable information on bedrock structures in order to obtain more accurate results when assessing the impacts of subsurface dams on groundwater flow and mass transport processes.
The study focused on a specific coastal aquifer without consideration of significant groundwater abstraction activities.Further observational information and numerical simulations covering a wide range of different hydrogeological settings and anthropogenic activities help understanding the full impacts of subsurface dams on groundwater conditions and are beneficial to achieve more effective usage of subsurface dam technology.

Figure 1 .
Figure 1.Geologic map of Miyakojima Island.Water tables are shown in the W-E cross-section, and the bedrock surface contour is also shown in the Sunagawa Subsurface Dam area.

Figure 2 .
Figure 2. The 3D model domain and settings and water retention curves: (a) the spatial discretization, (b) the geologic fence diagram showing the three geologic units, and (c) modeled water retention curve using the van Genuchten model.

Figure 3 .
Figure 3.Time series data for scenario S1 and scatter plot of model calibration: (a) daily rainfall, (b) modeled groundwater levels compared with the observed data, and (c) scatter plot of model calibration.In b1 through b5, the dashed line indicates the crest elevation of the Sunagawa Subsurface Dam (31 m amsl), and the gray zone indicates the recharging period of the subsurface reservoir.

Figure 4 .
Figure 4. Time series data of 1998 to 2000: (a) daily rainfall, (b) modeled groundwater levels in the upstream area (Well 3), and (c) modeled groundwater levels in the downstream area (Well 5).Vertical dashed lines mark the timing of groundwater level peaks in S0.

Figure 5 .
Figure 5. Comparisons of (a) travel time and streamlines, (b) groundwater levels, and (c) groundwater levels and rainfall infiltration routes before (left panels) and after damming (right panels).

Figure 6 .
Figure 6.NW-SE cross-sections (vertically exaggerated by a factor of 10) showing the modeled salinity distributions and groundwater levels at different times in S1.See the cross-section location in Figure 1.

Figure 8 .
Figure 8.The effects of E dam and L dam on (a) L toe and (b) T SWI for scenarios S1 to S7.