Coupled three-dimensional modelling of groundwater-surface water interactions for management of seawater intrusion in Pingtung Plain, Taiwan

Study region: A coupled framework, linking subsurface flow and surface hydrodynamics, is developed and applied to a real-world case study of Pingtung coastal aquifer in southwest of Taiwan, in East Asia. Study focus: FEFLOW is adopted to develop a 3-D variable density and transient groundwater model of the Pingtung shallow aquifer lying 250 m below mean sea level (MSL). This model is coupled with a 1-D river network model, comprised of the main river and its two tributaries, using MIKE 11 through the IFM MIKE 11 coupling interface. The model is capable of analysing the relationship between rainfall, surface water and groundwater recharge lag time. Also, the analysis of potential river inundation and maximum river discharge enable the model to choose the best location to apply artificial recharge as a management scenario to mitigate the effect of seawater intrusion (SWI). To the authors ’ knowledge, the developed coupled model is the first detailed integrated framework analysing the interaction of surface and subsurface water, with the capability to contribute to the restoration, rehabilitation, and management of the river network. New hydrological insights for the region: The rainfall ratio in the wet season to dry season is significant in this plain comparing with the rest of Taiwan. Also, southern Taiwan experiences the largest sea and river interaction, while Kaoping River playing as a pathway role for inland lead of seawater intrusion.


Introduction
In recent decades, there has been significant spatial and temporal instability in hydrological conditions due to climate change (Simonovic and Li, 2003) increasing the vulnerability of water resources (Chang et al., 2011).Groundwater plays a vital role in providing stable regional water resources (Tsur, 1990).Fresh groundwater in coastal aquifers is vulnerable to salinization by upconing and seawater intrusion (Post, 2005).An integrated framework considering surface and subsurface water as a single resource is essential for improving the reliability of water supply (Sophocleous, 2002;Montazar et al., 2010).Numerical models can be used to investigate the behaviour of surface-subsurface water systems in coastal areas for different climate change scenarios (Yang et al., 2015).Surface water and groundwater are inseparable elements in a hydrological system.From climatic and physiographic viewpoints, they interact in different ways such that, e.g., contamination of one directly affects the other (Winter, 1998).Studies in catchment-scale have illustrated the relative connection between surface and subsurface water resources.For example, there is clear relationship between the availability of water and the effect of groundwater over-abstraction on rivers or seepage of contaminants from rivers towards groundwater (Shaad, 2015).Integrated water resource management is one of the cornerstones of the global water framework (Rahaman et al., 2004) increasing the demand for modelling tools and methodologies to investigate integrated fully dynamic surface-subsurface models at the catchment scale.
Numerical modelling of aquifer-stream interactions has become an imperative tool for management of water resources.For example, the coupled models can be used to evaluate the performance of infiltration wells in cutting off the runoff of extreme rainfall Fig. 1.Pingtung plain study map.
M. Dibaj et al. and attenuating flood peaks (Sommer et al., 2009).The groundwater can interact with streams in three ways; inflow of groundwater through the streambed, outflow of groundwater along the streambed, or both losing and gaining stream in some stream reaches (Irvine et al., 2012;Storey et al., 2003).Also, abstraction of groundwater can cause local exchange of water between streams and adjacent shallow aquifers (Sophocleous, 2002;Brunner et al., 2011).Several factors control the groundwater and surface water interaction in a large scale landscape.These include the distribution of hydraulic conductivity in both alluvial plain sediments and channel, the relationship between groundwater level and adjacent stream stage, and the position and geometry of the stream channel located in the alluvial plain (Woessner, 2000).The direction of surface and subsurface water exchange is variable and is influenced by changes in hydraulic head due to the seasonal pattern and climate change (Scibek et al., 2007).The exchange flow between stream and aquifer is a function of the water head difference between aquifer and river stage.
In this study an integrated 3-D groundwater model and 1-D river model was developed for coupled transient analysis of groundwater flow, river flow and river inundation considering the interactions between surface water and groundwater.Although there are some coupled models simulating the flow exchange between surface and subsurface water, to the authors' knowledge, this is the first numerical model that also simulates the effects of rainfall occurrence and duration on groundwater-surface water interaction lag time.The developed model is also used to select appropriate locations and sources for applying recharge wells in coastal aquifer, based on resulted river discharge and inundation.The model is validated against experimental results and applied to study the interaction between the Kaoping River and the sea with the aquifer system in Southern Taiwan (Wang et al., 2009).This model is capable of computing the monthly variation of salinity amount and the exchange of mass flux in the mixing zone.This is the first numerical model for simulating the changes in salinity of the mixing zone for this specific region.

Study area
Pingtung plain, with an area of 1210 km 2 , is approximately 20 km wide from east to west and 60 km long from north to south, and is located in the southwest of Taiwan.The plain is bounded by foothills and river valleys in the north, low hills along the west, a mountain range extending along the east and the Taiwan Strait along the southern boundary (Dibaj et al., 2020).The mountainous regions surrounding the plain are consisted of black slate embedded with sandstone in western part of mountain range, shale and sandstone in north and sand stone with sandy shale in North West.The black slate with very low water-bearing texture and low porosity, causes high rate of run-off.Consequently, water is stored only in local joints and cleavages (Ting et al., 1998a).These characteristics of the plain turn most of the precipitation to surface run off and flow towards the streams and rivers (Tu et al., 2011).
The elevation of the Plain is less than 100 m above mean sea level (MSL) and it is a relatively flat area drained by three major rivers, the "Kaoping", "Donggang", and "Linbian".These rivers divide the plain into three areas of proximal fan, mid fan and distal fan (Dibaj et al., 2020;Hsu et al., 2007) (see Fig. 1).These rivers make a regular discharge base for the regional groundwater in this basin.The 171 km long Kaoping River covers a catchment area of about 3257 km 2 , and is the second longest river in Taiwan, originating from Yushan Mountain and wandering south-westward, receiving the streamflow of four main tributaries (Laonong River, Cishan River, I-Lio River, Zhuokou River) before discharging into the Taiwan Strait (Wang et al., 2009).The Laonong and Zhuokou Rivers are out of study domain area and not included in Fig. 1 (Capart et al., 2010).The estuary of Kaoping River has the largest sea and river interaction in southern Taiwan (Wang et al., 2009).The entire aquifer underlying the plain is formed by unconsolidated sediments.The direct penetration of rainfall due to the unconsolidated sediments through the foothills along the north and east of plain is the principle mechanism for the groundwater recharge.Furthermore, the river seepage in the higher areas of the plain is another source for groundwater recharge (Chang et al., 2004).The rivers in Taiwan possess huge sediment discharge compared to major rivers of the world (Liu et al., 2008;Hung et al., 2004).This is mainly because of the general characteristics of rivers in this country, with a short length, small basin area, steep gradient and brittle geology (Wang et al., 2009).The upper and middle stream of Kaoping River flows through mainly sandstone, shale and metamorphic rocks, while the downstream portion is mostly formed by recent alluvium deposits according to the several stratigraphic units exposed in its drainage basin (Wang et al., 2009;Yang, 2001).The Cishan River watershed has an area of 842 km 2 with 118 km river length (Chiang and Liu, 2013;Chen et al., 2011).The upstream of the watershed is composed of sand and shale and the mid-and downstream are composed of gravel and sand.The I-Lio River watershed has an area of about 642 km 2 and a river length of 68.5 km (Yeh and Chang, 2019), and its watershed is composed of gravel and sand with high permeability.The tortuous up-and mid-stream has good potential for subsurface flow development (Lee and Wang, 2011).According to the literature, groundwater aquifers supply almost 60 % of tap water demand in Taiwan highlighting the importance of this water resource for this country (Liang et al., 2019).Pingtung alluvial plain is one of the major agricultural and aquacultural production regions in Taiwan (Yu and Lin, 2015).About 50.5 % of the Pingtung plain land is used for agricultural activities and 5.5 % for aquaculture (Liang et al., 2020).Pingtung Plain aquifer is the second largest aquifer in Taiwan, playing a primary role in supplying the major water resources for over 50 % of population in Pingtung County as well as neighbouring areas (Yu and Lin, 2015) while surface water only supplies 20 % of water resource demand in the plain (Lin et al., 2011).This highlights the dependency of this plain on groundwater resources.A recent study showed that every year, approximately 25 million tons of groundwater is abstracted, of which 63 %, 8 % and 29 % is used for agriculture, aquaculture and other demands, respectively (Liang et al., 2020).
The permeability of the soil decreases from north east to south west.The proximal fan is mainly comprised of gravel and sand (20 % sand and 60 % gravel) with storage factor of 6.5 × 10 − 3 and a transmissivity of about 9000 m 2 /day (Hsieh et al., 2010).It is followed by finer grain size in the mid fan zone in the centre of the plain (40 % gravel and 40 % sand) with a 9.5 ×10 − 4 storage factor and Fig. 2. Developed subsurface and surface models.2300 m 2 /day transmissivity, reaching the distal fan in the south west that is comprised of silt and clay with 0.00005 storage factor and 1200 m 2 /day transmissivity (Dibaj et al., 2020).The indicated storage factor values used by previous studies were only monitored for the top aquifer and not defined for other layers.Due to the lack of data, the default value of this parameter in FEFLOW was applied with 0.0001 m − 1 as specific storage (S s ) for whole aquifer system.
The precipitation distribution in Pingtung is closely related to the monsoon, consequently increasing with a rise in topographical elevation (Fan and Hsieh, 2010;Hsu and Guo, 2005), forming a higher precipitation rate along the north and north east of the plain.The location of the I-Lio River watershed in the north east of the domain gives it the potential to have a maximum average annual precipitation varying between 3400 and 5000 mm.This rate decreases to the range of 3000− 3400 mm and a maximum of 3000 mm in the Cishan River watershed and along the Kaoping River, respectively (Lin et al., 2011).Most of the rainfall occurs due to regular thunderstorms and typhoons over the wet season (May-October), with a little rain over the dry season (November-March) (Dibaj et al., 2020).The temporal distribution of precipitation is uneven such that the precipitation in the wet season is 9 times higher than in the dry season (Hsu et al., 2007).The summer season (May to October) has the highest precipitation concentration in the Kaoping River basin, while the rest of the year is almost dry due to prevailing winds along the northwest of the plain (Chu, 2012).Accumulation of discharge from June to September forms about 78 % of the total annual discharge occurring in the Kaoping River.The average annual water discharge in this river has been reported as 6.33 × 10 6 m 3 per year based on the observations from 1999 to 2004 (Wang et al., 2009).

Finite element modelling
A three dimensional variable density finite element groundwater model of Pingtung's multilayer aquifer system was developed, including its three aquifers and two aquitards (Dibaj et al., 2020;Huang and Chiu, 2018).The groundwater was abstracted from 25 pumping wells, and monitored through 95 observation wells over the period from 2011 to 2015 (Groundwater level (GWL) stations for top aquifer are shown in Fig. 1).First the 2D domain of the plain was generated as a polygon shapefile captured by ArcGIS, and then three river polylines and 25 pumping well points were added in the FEFLOW (Finite Element subsurface FLOW) software (Dibaj et al., 2020;Trefry and Muffels, 2007a).The mesh generation was performed as an intermediate stage representing inner and outer geometries.A triangular prism mesh was used to discretise the domain.This type of mesh has the capability to discretise irregular geometries including complex polygons or polylines with narrow sections along their boundaries, variable layer thickness, and well fields or coastlines (Dibaj et al., 2020;Dhi-Wasy, 2010).It also allows a minimum angle for all finite elements as well as providing the means for local mesh refinement with maximum element size.The finite element mesh design was implemented with a special attention to minimise the distorted elements through application of local mesh refinement in the coastal area, polygon domain, pumping well points and river polylines.The local mesh refinement and mesh smoothing would allow to build a higher mesh quality than later elemental subdivisions.To eliminate the difficulties arising from highly distorted elements in the transport simulation, the Peclet criterion was calculated and updated at each time step for determining the required mesh density during simulation (Trefry and Muffels, 2007a).The available elevation (X, Y, Z) data from 95 observation wells were used for extending the 2D model into a 3D multilayer aquifer system (Dibaj et al., 2020).The nearest neighbour interpolation method was used for allocating the elevation to the adjacent node or batch of nodes adjoining the original locations of the observation wells.The aquifer with six layers and seven slices was achieved after the interpolation.Vertical mesh disorders were prevented by discretising each layer into sub-layers, leaving the model with 42 layers (Dibaj et al., 2020) (see Fig. 2b).
The model includes 39,221 elements per layer and 20,352 nodes per slice.The elevation of the model surface varies between 2.3 m and 89 m above the MSL, with the top layer having a depth of 10− 16 m.This layer overlies the top aquifer, with 10− 121 m thickness.A semi-pervious layer with 10− 81 m depth is located beneath the first aquifer as the first aquitard, overlying the middle aquifer with 10− 96 m thickness.This (second) aquifer overlies the second aquitard (with 10− 63 m depth) and the bottom aquifer (with 10− 103 m depth).The base of the bottom aquifer reaches to around 197 m below MSL.Once the multilayer aquifer system with dedicated mesh generation and model dimensions were set up, its properties were then asigned as an unconfined free surface fully saturated water table modelling approach with moveable (variable) mesh.Therefore, the mesh can be adapted to the groundwater level variation and all unsaturated or variably saturated parts can be excluded from the model avoiding the dry cell problem.In the mass transport cases, the dry cells can also freeze the solute and make it unable to transfer through the water table, yielding an incorrect model.All the slices below the free surface were set as dependent except for the bottom slice which should always be fixed.The dependent slices underlying the free surface are able to move if necessary.The layers located under confined layers (regardless of the hydraulic head being below or above the layer) were treated as fully saturated during the entire simulation.This is also in good agreement with previous studies in this region based on the available soil permeability and precipitation data (Huang and Chiu, 2018;Dibaj et al., 2020).The unconfined layer was given a default residual water depth of 0.01 m to avoid any dry cells in top slice during simulation.
The precipitation was monitored daily at 24 rainfall stations in different locations of the plain over a period of 5 years, from 1st January 2011 to 31st December 2015 based on Taiwan's Water Resource Bureau (MOEA).The average monthly rainfall distribution data assuming 60 % run off and evaporation (MOEA) were applied to the top layer of the model as a spatial inflow using the delogarithmic Kriging interpolation method (Dibaj et al., 2020).An automatic time-step control technique based on a predictor-corrector scheme was selected for time discretisation by choosing an appropriate length for time-step in line with the changing of the groundwater head variations between the time steps (Trefry and Muffels, 2007a).This feature of FEFLOW, among others, persuaded the present study to choose this software for its detailed simulation using automatic and flexible time stepping techniques.Although this model could also be developed using physical-based models such as WASH123D (Yeh et al., 2006), it may not have the high capabilities of FEFLOW.For example, WASH123D is also capable of simulating 3D subsurface flow models coupled with surface flows, but in WASH123D the time step size for both domains should be fixed and cannot change during simulation.Besides, it has limitations for local mesh refinement that may result in simulations being computationally expensive (England et al., 2006).FEFLOW applies default second order accurate, semi-implicit Forward Adams-Bashforth/backward trapezoid rule (AB/TR) for the flow and transport models.The details of average monitored rainfall data can be found in our recent study as well as the sensitivity of the inland encroachment distance of seawater due to 30 % increase and decrease in rainfall which has been analysed as well (Dibaj et al., 2020).The results showed the model was much more under influence of pumping rate.
The groundwater hydraulic head was monitored via 95 observation wells in different layers of the aquifer on a monthly basis over the same period of time.The data for the first month (January 2011) were applied as the initial head in the model using the delogarithmic inverse distance interpolation method for each layer of the model.Due to lack of data for the abstraction rates, the depths of the wells and their exact locations in the plain, the model used the pumping rates of the year 1996 (the only available data), for 25 townships of the plain as multilayer pumping wells screened at the entire depth of the aquifer system, as input transient pumping data over the simulation period.Detailed information about townships with their allocated pumping rates in Pingtung plain is available in a recent study (Dibaj et al., 2020).Based on the available data (MOEA), the hydraulic conductivity was measured only in the horizontal direction (K x = K y ) in different layers of the aquifer, ranging from 1 to 520 m/day, 5 to 90 m/day and 1 to 100 m/day in the top, middle and bottom aquifer, respectively.Considering the leading directions of the axis-parallel anisotropic hydraulic conductivity in coincident with Cartesian coordinates (K x , K y , K z along the X, Y and Z axes) (Dibaj et al., 2020) yields the following relationship, K x = K y = 10 K z , as an initial model set up for each aquifer in the model using the Kriging logarithmic interpolation method for spatial distribution in the 3D model.A constant value of 0.0001 m/day was assigned as the hydraulic conductivity of the aquitards (Dibaj et al., 2020).An automatic time stepping method with a predictor-corrector scheme and initial time step of 0.001 day was selected for the simulation (Dibaj et al., 2020).Hydraulic conductivity along X, Y and Z axes was considered for calibration due to the sensitivity of the model to this parameter and its significance for natural groundwater recharge from precipitation.
The southern border (the coastal line) was assigned first type (Dirichlet) head and mass concentration (seawater salt concentration is 35,000 mg/l) boundary conditions, and the remaining borders were considered as no fluid flux boundary, due to the low permeability soil and faults with impervious bedrocks surrounding the plain from north, east and west (Dibaj et al., 2020;Ting et al., 1998a).The only contaminant was salt and the adsorption was assumed negligible (Dibaj et al., 2020).The groundwater flows westward into Kaoping and Donggang Rivers and south-westward into the Taiwan Strait.Detail of the boundary conditions, governing equations and the parameter definitions can be found in our pervious study (Dibaj et al., 2020).
The model was simulated in steady state and then transient conditions.The initial salinity data were obtained through a steady state analysis of the model due to the lack of data.The results showed that the entire depth of the aquifer system was intruded with seawater while from the plain view, the middle and west of the plain were affected in greater areas than the east, which is in good agreement with the previous studies (Huang and Chiu, 2018;Dibaj et al., 2020).In this study, the molecular diffusivity is assumed as zero.The hydrodynamic dispersion is considered as a function of velocity and dispersivity is assumed to be larger in the flow direction than vertical and lateral directions (Ogata, 1970;Voss, 1984).FEFLOW approximates the longitudinal dispersivity based on the mean elemental size in a way that a good convergence is also guaranteed (Trefry and Muffels, 2007a).The longitudinal dispersivity varies between 6 and 150 m for most of studies (Gelhar et al., 1992;Sherif et al., 2012).In this study, the mean elemental dimension was 70 m which was set as longitudinal dispersivity.The transverse dispersivity was selected as 10 percent of this value based on the literature and FEFLOW parameter definition (Huang and Chiu, 2018;Sherif et al., 2012;Bear et al., 1999).According to the literature, this model also set the value of the horizontal transverse dispersivity the same as longitudinal dispersivity (Huang and Chiu, 2018;Dibaj et al., 2020).The model convergence and accuracy of meshing were monitored based on maximum interior angle of triangles, which was as regular as possible with the mean angle of 124.4 degree for the entire model.Also the numerical stability of advection-dispersion transport process was investigated based on Peclet number with value less than 2 along the shoreline boundary (Trefry and Muffels, 2007a).

Surface water model development
MIKE 11 is an advanced hydroinformatics engineering software with the capability of simulating 1D flow in rivers, estuaries, channels, irrigation systems, and other water bodies.MIKE 11 solves the Saint-Venant equations using a finite difference scheme (Madsen et al., 2003;Kamel, 2008) and can reflect fully dynamic, diffusive or kinematic conditions.Some assumptions are involved in the solution of these equations; firstly, assuming water as homogeneous and incompressible leads to neglecting density variations.Secondly, if the water depth is smaller than the wave length, then the flow is assumed to be parallel everywhere to the bottom level, leading to vertical acceleration being ignored and assuming that the hydrostatic pressure varies vertically.Finally, the flow is assumed to be sub-critical because of restrictions associated with applying super-critical flow models (Kamel, 2008;Thompson et al., 2004).The river data were implemented and edited separately with different editors and linked together via the simulation editor.
The Kaoping, Cishan and I-Lio Rivers were developed using point digitisation in MIKE 11 network (see Fig. 2(c)).The river centre points (Hpoints) were achieved using the centre of the right and left banks of each cross-section data using their X and Y coordinates.This network was comprised of 146 centre points.The distance between each point from the nearest point represents the chainage, which was calculated automatically starting with 0 km in the furthest upstream part and increasing towards the downstream of the river.The Kaoping River, Cishan River and I-Lio River stretch from chainage 0− 49120 m, 0− 10793 m, and 0− 11785 m, respectively (see Fig. 2(c)).From the collected data (MOEA), the water level (Hydrological stations W.L in Fig. 1), streamflow and rainfall data that were monitored from 1st January 2011 to 31st December 2015 were used as input data for surface water simulation in MIKE 11.The streamflow was monitored in five stations, namely Shanlin Bridge, Lilin Bridge, Wannta Bridge, Da-Jin Bridge, and Liouguei Bridge (see.Fig. 1).The water levels measured at 22 stations, with only two of them located at I-Lio River upstream and Kaoping River downstream, were applied as open-end river boundary conditions.The monitored water levels in Qishan Bridge, Lilin Bridge, Gaobing Bridge, and Wannta Bridge (See Fig. 1) were used for calibration and validation of the model.The daily time series of river water level and streamflow starting from 1st January 2011 and ending on 31st December 2015 were applied for Kaoping downstream and upstream, respectively.Both upstream and downstream of Cishan River were assigned as daily streamflow time series boundary conditions but as the downstream of Cishan River was discharging into the Kaoping River, and was not an open-end boundary, hence a lateral discharge was assigned as a point source at this location.The water level time series was applied for the upstream of the I-Lio River.The rivers were assumed to be freshwater and assigned with a constant 0 mg/l as the initial salinity over the simulation time.The model did not simulate the saline river intruding into groundwater, but only mass transports from groundwater to river water.The raw cross-sectional data (X, Z), based on the distance from the far left side and elevation from datum, were applied for each river centre point representing the river profile.Then the model used the raw data to calculate cross-sectional area, level, hydraulic radius, and hydraulic resistance as processed data (see Fig. 2c, d).
The steady state of the river profile was computed automatically at the dedicated start time by using the given boundary conditions via the simulation editor.The actual cross-section data, river nodes (Hpoints) gave the head values (H) due to processed data (level, cross sectional area, hydraulic radius and flow width) and the Q (flow) values were simulated through the Manning equation assuming uniform or critical flow.This study assumed uniform flow using the following equation (Thompson et al., 2004;Panda et al., 2010;Mike, 2003): where Q(h) is the discharge amount (m 3 /s) depending on the water depth h (m) calculated in the processed cross-section data.The bed slope is represented by S, while A(h) (m 2 ) is the level relevant to the area resulting from the cross-section processed data, R is the hydraulic radius (m) and the Manning roughness coefficient n (s/m 1/3 ) represents the resistance of the river bed.

Coupling of the groundwater and surface water models
The coupling of Mike 11 and FEFLOW was not an iterative process and each time step was only calculated once by either FEFLOW or MIKE 11.For defining the actual water levels in coupled boundary nodes, FEFLOW took MIKE 11's water levels at the end of the last time step.FEFLOW calculated the exchange fluxes (q) of each single boundary condition (third type -Cauchy) between the surface and groundwater separately using the transfer coefficient (ɸ h ) (1/d) as the main element for controlling this flux (Monninkhoff and Hartnack, 2009a): where q is the Darcy flux (m/d) of fluid (positive flow is towards the groundwater) and h gw and h ref represent the groundwater and river heads (m), respectively.The river was laterally coupled to the top two slices of the FEFLOW model which is illustrated in supplementary material as Data 1.The total discharge at each node, Q (m 3 /d), was calculated through multiplication of the flux (q) at every single finite element boundary node by the related exchange area, in m 2 , at the end of each time step in FEFLOW.The obtained values (budgets) were transferred to the Hpoints of Mike 11 as single point source inflow (lateral) boundary condition (Monninkhoff and Hartnack, 2009b).Mike 11 calculated as many as internal time steps as required to reach the actual time of FEFLOW.Transferring the resultant water level in the final FEFLOW time step from Hpoints of MIKE 11 to the boundary nodes of FEFLOW was the last stage of the coupling process (Monninkhoff and Hartnack, 2009a).An interface controlled the internal time steps of MIKE 11, which could be either constant or fit to the model's dynamic.On the other hand, FEFLOW controlled the groundwater time steps, in which IFM MIKE 11 automatically integrated the spatial overlay of both meshes.Every node of the FEFLOW boundary was assigned to the adjacent Hpoint of the MIKE 11.Therefore, the possibility of assigning one Hpoint to more than one boundary node in FEFLOW could be increased.Coupling of mass transport between Mike 11 and FEFLOW is a straightforward process and once a time step is completed by FEFLOW, it is passed to Mike 11 to exchange the values, taking a time step in turn.For each coupling point, FEFLOW must pass the values of flow (m 3 /s) and mass flux (g/s) of each chemical species (inflow into river) every time a Mike 11 time step is calculated.In return, Mike 11 passes the water level, concentration (mg/l), and mass flux (g/s) of each chemical species (outflow of a river) at the Hpoints.The flow boundary nodes where the flow has been coupled are automatically used for mass boundary node generation by FEFLOW.To prevent overshooting, a lumped mass matrix was selected in FEFLOW (Trefry and Muffels, 2007b).The river boundaries were assigned with the first kind (Dirichlet) mass transport boundary condition in convective form due to its high degree of stability, especially along outflow boundaries.

Calibration and validation of the model
Providing reliable information on different features of real situation requires a well stablished calibration of hydraulic model.The Manning roughness coefficient was subjected to calibration to be adjusted for achieving the least possible gap between measured and computed water levels.For this purpose, the uniform values of this parameter, ranging between 0.031 and 0.039 m − 1/3 s, were specified through the river network system.These values were selected based on the flow of the river through the alluvium deposit (Arcement and Schneider, 1989).The observed water levels in Qishan Bridge, Lilin Bridge, Gaoping Bridge, and Wannta Bridge (see.Fig. 1) were used for comparison with the corresponding simulated values in these locations during the calibration period of 1st January to 31st December 2011.The Nash-Sutcliffe efficiency coefficient, R 2 was used to assess model performance for each station (ASCE Task Committee on Definition of Criteria for Evaluation of Watershed Models of the Watershed Management Committee, I. and D. Division, 1993).The daily calibration results for the first year of the simulation period are presented in Fig. 3. Based on the available data, the Qishan Station monitoring water levels at Cishan River upstream showed a maximum of 39 m over May to June.By using Manning numbers of 0.039, 0.037 and 0.035 m − 1/3 s, the maximum water levels were over predicted as 44, 43, and 42 m, respectively.On the other hand, the computed water levels using Manning numbers 0.031 and 0.033 m − 1/3 s were underestimated by 2 and 1 m (maximum of 37 and 38 m) respectively.Lilin Station, monitoring the water levels of the Cishan River in its downstream, indicated a maximum amount of 26 m in August and September and a minimum of 23 m from January to April.The computed water levels using Manning numbers 0.031 and 0.033 m − 1/3 s were underestimated, reaching a minimum of 19 m and maximum of 21 m during the same periods.An overestimation of about 1, 4 and 7 m (24, 27 and 30 m) was yielded due to Manning coefficients 0.035, 0.037 and 0.039 m − 1/3 s, respectively.The observed water levels at Gaoping Station showed its minimum level, fluctuating between 8 and 9 m during January to May, rising to above 9 m until the end of July and reaching its peak value of over 11 m in August and September.It then followed a decreasing trend from October to December, varying between 9 and 10 m.Among the computed levels, all of their peak values determined in August and September, were slightly below 12 m.However, the computed values with Manning numbers 0.031 and 0.033 m − 1/3 s underestimated the water level during other months, except the peak times.Applying uniform Manning roughness of 0.035 m − 1/3 s provided rather better estimation compared to other values.The last stage of calibration was accomplished at Wannta Station, where the observed water levels from January to mid-July fluctuated between 4 and 5 m, rising by up to 1.5 m from August to early November, and jumping to the maximum level of about 7.5 m by the end of November.Therefore, the Manning factor of 0.035 m − 1/3 s was applied to all rivers within the Pingtung plain.The R 2 values for each simulation are provided in Fig. 4 with their overall average values.
A reliable steady state solution is a primary need for simulation of flow and solute transport problems.To achieve a stable condition, both river and groundwater head and mass boundary conditions were kept unchanged over time using a time stepping method with time step length of 0.001 day (Voss, 1984;Diersch and Kolditz, 1998) to simulate the initial transient model.The model reached a steady condition using a predictor-corrector scheme over a long time of 8395 days, which was not unexpected for such a large-scale numerical model.Due to the lack of observed mass concentration data, the resulting salinity values obtained as the output of the steady-state simulation were applied as input mass data for the transient model after checking their validity with previous studies.The results indicated Kaoping River mouth and estuary as the origin of inland encroachment of seawater following its passage into the aquifer system, and horizontal and vertical expansion in all layers of the aquifer system at different rates.This was in good agreement with previous studies showing the middle and left sides of the domain as being more vulnerable to seawater intrusion compared to the right side (Dibaj et al., 2020;Huang and Chiu, 2018;Ellis, 1999;Jang et al., 2016).Also, the top aquifer showed more intrusion than the bottom one and the second aquifer was less intruded, which was in good agreement with the permeability distribution in different layers of the aquifer system (Hsieh et al., 2010).One of the limitations of using FePEST Calibration software is its inability to calibrate the coupled models at the same time.To overcome this issue, the IFM Mike 11 interface was deactivated during calibration and the effect of the rivers and their boundary conditions were applied into the FEFLOW model in the form of Cauchy head and Dirichlet mass boundary conditions.The hydraulic conductivities of the three aquifers were selected for adjustment in calibration using FePEST software (Doherty and Hunt, 2010).The entire calibration process was set up the same as the previous study (Dibaj et al., 2020) and accomplished for the year 2011 against all observed hydraulic heads.In this process, the hydraulic conductivity along X, Y and Z axes in each aquifer was considered as adjustable parameter with the assumption K x = K y = 10 K z .Based on the observed permeability data (MOEA) in each aquifer, K x and K y were bounded between 0.1 and 600 m/day and K z could range between 0.01 and 100 m/day as lower and upper bounds.Details of the observed permeability data can be found in (Dibaj et al., 2020).Both K y and K z were tied to the change of K x during calibration using the Kriging interpolation method for 50 random pilot points.The calibrated values for each aquifer are illustrated in the supplementary material as Data 2.
The scatter plot comparing the observed and computed hydraulic heads is illustrated in Fig. 5 using root mean square (RMS) for all 95 observation wells.The hydraulic heads were calibrated with RMS of 0.987 m as the best fit, with less than 1 m deviation between the observed and computed heads (see Fig. 5a).The Mike 11 River model was run over the period of 1st January 2012 to 31st December 2015 using a uniform Manning coefficient of 0.035 m − 1/3 s based on the output of calibration.In this case, the observed and computed results were compared in the four stations of Qishan, Lilin, Gaoping and Wannta.The values of R 2 had an overall minimum value of 0.68 and maximum of 0.86 at the Lilin and Gaoping stations, respectively.Most of the observed peak values were also predicted by the developed computational model (see Fig. 4).
This model was coupled to the transient groundwater model for the same time period (1st January 2012 to 31st December 2015) employing the resulting mass concentration at the end of 2011, the time variable monthly hydraulic head of January 2012 as initial head, calibrated conductivities and monthly precipitation data over 2012− 2015.The southern boundary of the domain was assigned first type constant head (Dirichlet) and mass concentration as well as no fluid flux condition bounding the north, east and west of the plain, as with the previous study (Dibaj et al., 2020).The model was validated using a scatter plot of computed against observed hydraulic heads with RMS of 3.012 m at the end of the simulation period (see Fig. 5b).

River-precipitation-groundwater interaction
Generally the interactions between groundwater and precipitation are represented by large scale procedures depending on landscape geometry and spatial distribution of hydraulic conductivity (Gurdak et al., 2007).The interactions can generally vary with time under influence of soil saturation conditions resulting in non-linear relationships between groundwater level and rainfall time series M. Dibaj et al. (Yu and Lin, 2015).The precipitation distribution is the main parameter affecting availability of groundwater resources in Pingtung plain while the soil permeability plays an important role in the recharging pattern (Hsu et al., 2007;Hsieh et al., 2010).The proximal fan extending along east and north east of the domain compromised by sand and gravel, has high potential for recharge and infiltration by river flow and precipitation (Hsieh et al., 2010).The main focus of this study is to investigate the interaction between river network and underlying groundwater under influence of the rainfall.This interaction is mainly due to the water level differences between river and groundwater which is controlled by the hydraulic conductivity of the interaction point (Hpoint).This will help to estimate the surface water-groundwater interaction and the mass exchange and water flux in the river network and and its adjacent nodes.The rigours calibration and validation of this model in a previous study (Dibaj et al., 2020) showed that hydraulic conductivity is the main parameter for groundwater recharge and its lag time to rainfall.
The simulated groundwater (Sim GW) and river water (Sim River stage) levels interacting with each other based on the average monthly rainfall rate, are presented at six selected river Hpoints during the simulation period (see Fig. 6).The exchange flux (Q fe ) (see Fig. 7) shows groundwater discharge (gaining stream) as positive and groundwater recharge (losing stream) as negative flow.The duration and amount of rainfall are effective parameters in the interaction of the groundwater and rivers.Although the groundwater recharge generally lagged by one to two months after peak rainfall, this period varied in different areas of the plain.From the simulated results, the Hpoints located at the proximal fan, like Kaoping River upstream and I-Lio River, showed two to three months delay in the groundwater level's reaction to the occurrence of peak rainfall (see Fig. 6a, b), which was in good agreement with the literature (Yu and Lin, 2015).Considering the river bed as a reference interaction point, the groundwater head was always lower than the river level but it reached to bed level or higher during September and October.The groundwater fluctuated between its minimum and maximum levels with up to 20 and 15 m deviation at Kaoping upstream and I-Lio River, respectively (see Fig. 6a, b).The precipitation, river stage and groundwater head in the proximal fan showed highest interannual variation over the simulation time, reaching their maximum level in the summer of 2012.While 2013 and 2014 experienced typical flow conditions, 2015 showed a lower level of groundwater and river stage due to lower average precipitation.Groundwater had its lowest level from March to April and reached its maximum level during September and October after the recharge period between June and August.In 2012, the lag time of groundwater recharge was the shortest, leading to a fast (up to one month) reaction to the peak rainfall (typhoon).The interaction of the river and groundwater increased by rising groundwater level during recharge time and was at its maximum level over September to October.The flux exchange in the upstream of the Kaoping River reached its peak value by around − 0.5, +0.5 m 3 /s in June and August of 2012 respectively by discharging of the river into the underlying groundwater and vice versa (see Fig. 7a).This amount was not significant compared to the abstraction rate (at least 20 times less than abstraction) at the pumping wells adjacent to the Hpoints (Dibaj et al., 2020).It can be said that the groundwater was highly under the influence of precipitation and abstraction, leading to its fluctuation and interaction with the river water.I-Lio River showed almost the same trend as Kaoping River upstream for groundwater and river water levels.However, the resulted exchange flux was totally different; the values of the river flow were significantly higher, and the groundwater had a high potential to be recharged by the river with maximum of about 6.3 m 3 /s after peak rainfall in June 2012 (see.Fig. 7b), which was in good agreement with previous studies (Ting et al., 1998a).This, together with a lower abstraction rate and higher precipitation, led to the right side of the plain having a higher groundwater level (Ting et al., 1998b).
The groundwater recharge lag time to the peak rainfall at all Hpoints located at the mid fan of the plain was about one to two months, which was in good agreement with the literature (Yu and Lin, 2015).Compared to the proximal fan, the groundwater had a lower interannual fluctuation in the mid fan and the maximum head difference between the groundwater and river was about 4 m (in April) and 2-3 m (in May) at Cishan River upstream and downstream, respectively (see Fig. 6d, e).Increasing the rainfall led to decrease in this deviation, reaching its minimum gap over September and October.This trend was exceptional in 2012 when heavy rainfall occurred earlier than usual in June.The computed exchange flow at Cishan River upstream showed that the river discharges into the groundwater at a higher rate than the groundwater discharges into the river.This rate reached its maximum level in June 2012 with around 1.5 m 3 /s towards groundwater (see Fig. 7d).The conjunction of the Cishan River and Kaoping River occupied the largest cross-sectional area in the river network, with more than 130,000 m 2 with a maximum hydraulic radius of 46 m covering almost 3000 m storage width.This gave it the potential to be in connection with the underlying groundwater.In addition, the abstraction of   groundwater in this area with less than 50 MLD for a large area of Qishan township (Dibaj et al., 2020), helped the groundwater interact with the river at higher levels and have a lower head deviation.Similar to the river upstream, the maximum exchange flux at this Hpoint occurred in June 2012 with 0.8 m 3 /s (see.Fig. 7e).Also, this Hpoint showed the maximum river discharge (see supplementary material, Data 3), which was in good agreement with available data reported by MOEA. Fig. 6(c) represents Hpoint located at about 30 km from the Kaoping River upstream at mid fan area.The groundwater head reached its highest level over October and November 2012 with largest interannual oscillation of about 3.5 m between April and November.The groundwater level showed a sharp increase one month after the peak rainfall in June 2012; it was recharged by the river with maximum of about 1.5 m 3 /s (see.Fig. 7c) which led the interaction of groundwater with river by exceeding the river stage by up to 1 m over August 2012 and February 2013.This interaction reached its lowest level during 2014 due to increasing the gap between groundwater and river levels growing the gap between groundwater and river.By moving towards the river downstream, the interaction between the river and the underlying aquifers decreased due to increased groundwater abstraction rate fading the groundwater head, meaning that the rivers were losing stream almost all the time.In the downstream of the Kaoping River, water flowed over a confined aquifer consisting of silt and clay with lower permeability located at the distal fan.The results show that the groundwater had a faster reaction (less than one month) to rainfall events with insignificant head fluctuation, which was in good agreement with previous literature (Yu and Lin, 2015).The grain size becomes finer and finer from the Kaoping upstream in the proximal fan towards the river mid and downstream in the mid and distal fans.In upstream portions (before the conjunction of the three rivers), the river beds are located between slices 1 and 2 of the model.However, by moving towards the Kaoping midstream, the bed of the river lowers and is located between slices 2 and 3 with lower permeability, deprecating the interaction effect and exchange flow rate.The most downstream Hpoint of the Kaoping River was selected to present the resulting surface and subsurface water interaction after coupling (see Fig. 6f).Although the typhoon event in June 2012 caused both the groundwater and the river levels to increase, there was a limited rise of about 0.6 m with almost no delay in the level of the groundwater.High rainfall caused the river and groundwater to almost maintain their own levels from June to September.The river level declined during October and November and then it maintained its level by the end of February 2013 (see Fig. 6f).Accordingly, the groundwater head started to decrease from October and reached its lowest value in January 2013 moving towards the dry season.The maximum level of the river was almost the same as the minimum, about 0.5 m higher than mean sea level, in August 2013 and 2015.

Seawater intrusion (SWI)
A cross-section along the Kaoping River from its origin in the northeast to its discharge point into the sea in the southwest of the plain has been chosen to show the effect of river discharge on seawater intrusion (see Fig. 8a,c) at the end of 2011 and 2015.The crosssection view has been magnified 100 times in the vertical direction to present a clearer view.Also, a detailed view of this encroachment is presented along 4 km of the downstream portion of the Kaoping River cross-section (see Fig. 8b,d).The results showed that seawater Fig. 8. SWI along Kaoping River cross-section (a&c) and 4 km from the sea (b&d) at the end of 2011 and 2015 with 12 times magnification.
has intruded in the entire depth of the aquifer at different inland distances in different layers of the aquifer system.The Kaoping River estuary acted as an interaction zone, mixing the river water with seawater.This is one of the largest river and sea interaction zones in the world (Wang et al., 2009).The high river discharge pushed the seawater back into the sea in the surface layer (top slice of the model) of the aquifer system.However, the river bed (which is located between top aquifer and first aquitard) in contact with the top aquifer, played a pathway role for inland movement of seawater towards the top layers of the aquifer system.Mixing of the fresh water of the Kaoping River (rivers were assumed as fresh water) with the saline seawater in the estuary turned the saline water (18, 000 < TDS < 35,000 mg/l) into brackish water (10,000 < TDS < 18,000 mg/l) before discharging into the sea with a circulation pattern through the top layer of the aquifer (Dibaj et al., 2020).The saline water found its pathway through Kaoping River downstream from the top aquitard and first aquifer towards the upper layers with a cone shape encroachment at the end of 2011 (see Fig. 8a,b).The bottom aquifer along this cross-section was further intruded by the saline water highlighting the significant effect of groundwater abstraction in this aquifer which is far from the river and its effect.
However, the middle aquifer located between two confining layers showed less intrusion.This pattern was generally in good agreement with previous studies (Dibaj et al., 2020;Huang and Chiu, 2018).The semi-brackish water (1000 < TDS < 10,000 mg/l) extended further inland in the top aquifer compared to the other two aquifers highlighting the effect of fresh river discharge in diluting the salinity concentration at the top layers.After four years, at the end of 2015, the tip of the saline water cone moved towards the lower layers, extending between the first aquitard and second aquifer (see Fig. 8c, d) highlighting the significant effect of the Kaoping River bed location in the salinization of lower aquifers through the river mouth besides the effect of groundwater abstraction.In addition to this cross-section along the Kaoping River, the inland length of seawater intrusion with 10,000 < TDS < 35,000 mg/l along cross sections A-A, B-B and C-C (see Fig. 10a) at right, middle and left sides of the domain, respectively is presented at the end of 2011 and 2015 (see Table 1).Distance of inland encroachment is measured as the distance between the shoreline (contour line 35,000 mg/l) and contour line 10,000 mg/l (see Fig. 10b).This allows to compare the results of this model with those of a pervious study that did not consider the effect of the river (Dibaj et al., 2020).The simulated mass concentrations at the end of 2011 and 2015 remained almost unchanged from the predicted encroachment at the end of 2010 and 2015 (Dibaj et al., 2020) in all the layers of the aquifer along the A-A cross-section, with nearly 1790 m and 1820 m intrusion, respectively.The seawater had a forward pattern towards the land along all the layers of the aquifer through cross-section B-B, intruding the top, middle and bottom aquifers with the inland distances of 2380, 2350 and 2400 m, respectively, with around 2 % difference from the previous study at the end of 2010 (Dibaj et al., 2020).At the end of 2015, this encroachment exceeded more than 180 m, 160 m and 90 m through the top, middle and bottom aquifers, respectively.The resulted SWI at the end of 2015 showed about 100 m less encroachment through middle aquifer while the intrusion along the top aquifer decreased about 50 m and it was almost unchanged at the bottom aquifer.The first and second aquifers showed at least 100 m less encroachment along cross-section C-C at the end of 2011, highlighting the effect of fresh river water on shrinking of the salinity zone.However, this decline did not last long, with about 600 m inland progression of seawater along the first aquifer at the end of 2015 due to the expansion of brackish water (10,000 < TDS < 18,000 mg/l) (see Table 1).This transition was slowed down through the middle aquifer and kept growing through the bottom aquifer almost at the same distance as the top aquifer.At the beginning of the simulation, mixing of fresh river water with saline seawater at the river downstream reduced the salinity of water through the top layers of the aquifer, but with time, the Kaoping River estuary played a pathway role, guiding seawater inland through the river mouth due to the hydraulic head gradients caused by groundwater over abstraction.

Mass interaction in Kaoping River estuary
An experimental investigation (Wang et al., 2009) has shown the variation of salinity in the Kaoping River estuary in different monitoring stations during winter and summer, which is in good agreement with the predicted salinity changes in the developed coupled model.Our study aims to give a clear understanding of the effects of transient flow and coupling on the salinity variations in the mixing zone before discharging to the sea.The results of mass flux exchanges between the river and groundwater are shown in Fig. 9.These results highlight the consequent effect of rainfall events on the river water discharge and change of salinity in the mixing zone of the two water bodies (river and the sea) before discharging into the sea.The results of the coupled model showed only transfer of mass from the groundwater towards the river.The rivers were assumed as freshwater bodies and the exchange of mass between river and groundwater was investigated at 5 points along 4 km of Kaoping River downstream (with 1 km distance between each two points), before discharging into the sea (see Fig. 8b, d).Due to lack of water quality data and focus of present study on groundwater salinisation, the mass exchange was simulated from groundwater towards the rivers.

Table 1
Inland distance (m) of seawater intrusion along four cross-sections in 2011 and 2015.The mass concentration at the first point (P1), located at 1 km from the sea, showed minimum average monthly salinity of 1121 mg/l in June 2012 and maximum 7700 mg/l in March 2015.The mass flux was determined by multiplying the concentration (mg/l) of the river Hpoint by the exchange fluid flux (m 3 /s) between the river and groundwater.Due to the fast and high variation of river discharge (Ning et al., 2001) compared to the salinity change between wet and dry seasons, the mass flux between these two water bodies was dominated by the changing of fluid flux.For example, in point (P1), during the dry season (January-May), the mass flux from the underlying groundwater towards the river was slow and less than in the wet season, and consequently a lower amount of fluid was exchanged between these two water bodies followed by lower river discharge, yielding salinity rise in their interaction (see.Fig. 9a).Generally, the results for (P1) showed higher mass flux from the river to the sea due to this area's position adjacent to the sea and the high salinity of underlying groundwater.The discharged salt into the sea varied between a minimum of 150 kg/d and a maximum of 9000 kg/d in January 2014 and September 2013, respectively.By moving further inland, the second point (P2) showed slightly less salt transferring towards the river, which reached its maximum of 8000 kg/d at the end of July 2013 (see.Fig. 9b).Although the maximum fluid exchange occurred in June 2012 due to heavy rainfall, the peak mass exchange happened at the end of July 2013 due to three consecutive high exchange flows during July and August 2013.The salt concentration varied between 1000 and 7000 mg/l at this point during high rainfall and dry periods respectively, by exchanging salt at the river bed from the underlying aquifer during the simulation time, but this was not a considerable change, as the range was still in the semi-brackish water range.The mass flux reached the highest level twice in 2015 due to an increase in the salinity of water in May and September.At the third point (P3), the salinity concentration changed from freshwater (TDS < 1000 mg/l) in wet months to semi-brackish water (1000 < TDS < 10,000 mg/l) in dry periods.The peak of fluid exchange occurred with heavy rainfall in June 2012, yielding the peak amounts of salt to be exchanged from top aquifer towards the river bed (see.Fig. 9c).The fourth and fifth points, located 4 and 5 km from the sea, showed lower salinity concentration, varying between 800 and 3400 mg/l and 600 and 1800 mg/l between wet and dry seasons, respectively in the top layer of the aquifer system.As with the other points, the fluid flowed downward at a maximum of 0.05 and 0.06 m 3 /s in June 2012 at P4 and P5 respectively, yielding maximum salt exchange at the same time (see.Fig. 9d, e).The salinity variation between wet and dry seasons in these two points was not as high as the three previous points, causing the fluid flux to be the driving force in the salt transition.The fluid flux in these points (except in 2012) was on average at the same level, which gave a similar trend for the rest of the simulation time.From the vertical and plan views of the model, the surface layer and top aquifer were under the influence of river effects.The first aquifer, in contact with the river bed, had about 20-30 percent more salinity than the top surface layer, which was in good agreement with the previous study (Wang et al., 2009).The salinity of the underlying layers (beneath the first aquifer) did not change significantly along the river and almost maintained the trend, as shown in a previous study (Dibaj et al., 2020).

Seawater intrusion mitigation
Artificial recharge of groundwater is one of the effective methods for SWI mitigation (Javadi et al., 2015).For this purpose, high quality water (e.g., surface water, treated waste water, rainwater or desalinated water) can be used to recharge the aquifer, yielding a rise in inland piezometric heads and retaining the seaward gradient of the system.Selecting appropriate location(s) for recharge is very important for the recharge efficiency.Previous studies (Ting et al., 1998a;Tu et al., 2011) used pilot tests for estimating the groundwater recharge rate by selecting random locations as ponding areas in Pingtung plain.Four locations were selected in parallel with the groundwater flow direction running from I-Lio River at north east towards the Kaoping River at southwest to evaluate the soil permeability and texture impact on the recharge rate.The results showed the upper (first location) and lower part (second location) of I-Lio River with coarse texture have higher rate of recharge (Tu et al., 2011).This rate was decreased by moving towards the middle of the plain before reaching to Kaoping River (third location) in mid fan and the lower parts of Kaoping River (fourth location) which was in good agreement with textural change between these areas.Although the fourth location along the Kaoping River possesses coarse texture allowing high input of rainfall and irrigation water through the land surface, it showed lower contribution through the underlying aquifers.However, this location was likely located in river discharge area where the infiltrated water was possibly lost due to shallow drainage to adjacent canals or ditches (Tu et al., 2011).This location was also intruded by seawater which was not a good option for surface recharging (ponding).The present study used recharge wells for recharge purpose using the excess of river water based on the river flow and river inundation results in agricultural fields.Based on the river model results (see Fig. 11) and available land use data collected from National Land Surveying and Mapping Centre (Fig. 1), the conjunction point of Kaoping and Cishan Rivers, indicated the maximum river discharge; this is where a weir has already been constructed to regulate the river flow (MOEA).The ecological base flow at downstream of the Kaoping weir is 37 m 3 /s (MOEA).This location was selected as first artificial recharge using river as water source.Maintaining the ecological base flow in rivers is crucial for heath of river ecological system by not exceeding 37 m 3 /s for abstraction.Although artificial recharge could reduce seawater intrusion, its efficiency could be increased by choosing a recharge source adjacent to the toe of the saltwater wedge (Hussain et al., 2019).For this purpose, the Kaoping River downstream (around chainage 43,295 m), which had the potential for maximum river inundation (see Fig. 11) was selected as the second source of water for artificial recharge.The details of maximum river discharge and inundation area values in these two chainages along Kaoping River profile are illustrated in supplementary material as Data 3.The results showed that an appropriate rate for recharge would be between 2 × 10 5 and 8 × 10 5 m 3 /d (see Table 2).The rates less than 2 × 10 5 m 3 /d showed negligible effects on the transition of the saltwater wedge.This study applied injection wells (direct recharge) screening at the entire depth of aquifer system which were the best option (compared to surface infiltration) for this area considering the high rate of evaporation (Yeh and Chen, 2004), permeability of soil and underlying confined layers, especially in the distal fan.On the transition of saline water (seawards distance) with TDS > 10,000 mg/l is illustrated in Table 2 along the dedicated cross-sections.The results show that the toe of the saline water interface in the south east of the plain (along cross-section A-A) had less than 100 m seaward movement through both recharge sources (see Table 2).However, by moving towards the middle of the domain, the saline water zone (10,000 < TDS < 35, 000 mg/l) showed more seaward transition along cross section B-B, reaching a maximum of about 100 and 300 m through the aquifer layers due to artificially recharging from Cishan River downstream (first location) and Kaoping River downstream (second location), respectively.
Recharging near the toe location (around Kaoping River downstream) was more effective in moving the saline water towards the sea compared to the recharge from the upstream point (the conjunction point of Kaoping and Cishan Rivers) (see Table 2).This was also in good agreement with the pervious literature showing the lower part of Kaoping River had more potential for recharge (Tu et al., 2011).Compared to the top and bottom aquifers, seawater was less intruded in the middle aquifer.Also, the middle aquifer showed a Fig. 11.The water level along Kaoping River profile at heavy rainfall event (18-6-2012).
slower reaction to the recharge compared to the other two aquifers.From the results, it can be concluded that seawater intrusion can be mitigated by choosing Kaoping River downstream as a recharging source.Also, the results in Table 2 show that recharge rates between 2 × 10 5 and 6 × 10 5 m 3 /day lead to a sharper decline in inland intrusion of seawater compared to rates between 6 × 10 5 and 8 × 10 5 m 3 /d, which show a gradual pattern.Choosing rates higher than 8 × 10 5 m 3 /d did not cause a significant change in the results.

Conclusion
A 3-D transient density-dependent finite element model for groundwater flow and solute transport was developed.The developed model was coupled with a 1-D river network model for integrated simulation of groundwater flow, surface water flow and seawater intrusion in Pingtung plain.A meticulous calibration and validation was conducted on the model using real data.The developed coupled model is capable of analysing the flow and transient interaction between the aquifer system and rivers.Also, the effects of occurrence time, intensity and duration of precipitation on the groundwater-surface water interaction and lag time were investigated.The results showed that the groundwater in the proximal fan has the slowest reaction time to heavy rainfall whilst the lowest delay was observed in the distal fan.It can be concluded that the top layer of the aquifer system was less intruded due to fresh river water discharge pushing back the saline water seaward at the beginning of the simulation.However, with time, the Kaoping River mouth acted as the main pathway, letting seawater inland through river bed towards the top aquifer at the end of the simulation time.The middle aquifer located between two large confining layers was less intruded compared to the first and third aquifers, while seawater at the bottom aquifer showed farther inland encroachment.The inland distance of seawater was examined along 4 cross-sections in the right, middle, along Kaoping River, and left side of the plain.The east side of the plain showed shortest intrusion distance, and it was almost unchanged during the simulation period due to natural groundwater recharge and higher groundwater level.However, by moving towards the middle and west side of the plain the inland encroachment of seawater extended significantly.The results of monthly variation of salinity of water in the mixing zone of Kaoping River, and Taiwan Strait illustrated the exchange of mass flux between river and groundwater.It highlighted the effect of precipitation on exchanging fluid flux between river and the underlying aquifer during wet and dry seasons.The developed model was able to determine the appropriate locations for artificial groundwater recharge for SWI mitigation.The conjunction of Kaoping River and Cishan River, where the river flow reached its peak, was selected as the first location to apply injection well using the excess of river water.The Kaoping River downstream which was vulnerable to inundation was selected as the second location for recharge.The resulted seaward movement of seawater highlighted the effectiveness of applying recharge wells near toe location.
The developed model has some limitations.The storativity parameter for the entire model was assigned as the default value of the software due to lack of data.The river water was assumed as fresh water and the exchange mass flux was considered only from groundwater towards the river.Also, the rate of pumping wells was only available for one specific period of time with no information about their screening depth.The same rates were applied for the consecutive years, considering the screening in the entire depth of the

Table 2
Inland distance (m) of saline water (10,000 < TDS < 35,000 mg/l) along four cross sections due to different recharge rates at Cishan and Kaopingf River downstream.

Fig. 4 .
Fig. 4. Results of calibration and validation of the surface model based on the four stations.

Fig. 6 .
Fig. 6.Rivers and groundwater levels after coupling based on average rainfall.

Fig. 9 .
Fig. 9. Exchange flow and transient mass between Kaoping River and groundwater from P1 to P5.

Fig. 10 .
Fig. 10.(a) Cross-sections A-A, B-B, C-C, (b) an example of measuring the inland distance of seawater intrusion along Cross-section C-C from shore line node to the node on the isoline 10,000 (mg/l).