Vulnerability assessment considering impact of future groundwater exploitation on coastal groundwater resources in northeastern Jeju Island, South Korea

Jeju Island is the largest island in South Korea. Recently, extensive groundwater abstraction has been reported from the shallow aquifer in the northeast region of the island. This study simulated the freshwater resources of the aquifer to estimate the sustainability of groundwater use on Jeju Island in terms of its vulnerability to seawater intrusion. Three-dimensional finite-difference numerical groundwater models were simulated using the MODFLOW-family code SEAWAT. Precise and recent groundwater level and multi-depth salinity data obtained from the study site were used for model calibration; the simulated results showed good agreement with the observed data. SEAWAT was used to delineate the current seawater-freshwater interface to quantitatively estimate the coastal fresh groundwater resources. Future stress scenarios were also simulated in response to increased pumping and various changes in the recharge. The results showed that current groundwater use in the coastal aquifer did not induce seawater intrusion in the coastal aquifer, but seawater intrusion will occur if the dry season continues for the next ten years. The vulnerability assessment based on the predicted groundwater levels and ion concentrations using numerical simulations suggests future vulnerability in the aquifer; therefore, continuous assessment and visualization of the aquifer sustainability is vital. Future projections by the integrated SEAWAT simulation and GALDIT assessment showed that an increase in groundwater pumping may escalate the vulnerability status of coastal groundwater resources from moderate to high in some areas of the study site, by inducing lateral seawater intrusion in deeper areas of the unconfined aquifer.


Introduction
Coastal areas are often characterized by urbanization, population growth, and environmental degradation. These factors can lead to saltwater intrusion (SWI). Water pollution, ground subsidence, depletion of water resources, and a decline in groundwater levels due to excessive pumping intensify SWI into coastal groundwater, such that they are considered core hazards in water resource management. Furthermore, increased groundwater abstraction wells, land-use changes, droughts, and sea-level rise due to climate change are emerging issues in coastal groundwater management (Chang et al. 2011;Feseker 2007;Masterson and Garabedian 2007;Oude Essink et al. 2010;Ranjan et al. 2006;Rozell and Wong 2010;Sherif and Singh 1999;Houben and Post 2017).
Numerical models are commonly used to estimate the spatial and temporal heterogeneity of coastal groundwater resources and predict aquifer environments (Chang et al. 2011;Feseker 2007;Masterson and Garabedian 2007;Oude Essink et al. 2010;. Mathematically computed models based on density-driven flow determine SWI based on the inland movement of the fresh/saltwater interface (saltwater wedge) in response to different natural and/or anthropogenic activities (Feseker 2007;Younes and Fahs 2014;Shao et al. 2018). Currently, field investigations combined with numerical modeling are one of the best techniques to assess aquifer conditions and establish aquifer management plans based on scientifically reliable predictions (Feseker 2007;Masterson and Garabedian 2007;Oude Essink et al. 2010;Sulzbacher et al. 2012;Yu et al. 2010;Oki 2005).
Coastal groundwater resources in island regions are more vulnerable to SWI than mainland coastal aquifers because they exhibit geological and topographic characteristics that complicate the use of surface water resources Rozell and Wong 2010;Sulzbacher et al. 2012;Oki 2005). In other words, groundwater resources are more vulnerable to SWI when a greater land area is in contact with the ocean. A three-dimensional (3-D) fresh/saltwater model of Oahu, a well-studied volcanic archipelago, was developed by applying a dispersion model, with further advancements made regarding changes in the freshwater/saltwater interface caused by geological effects as well as groundwater pumping (Oki 2005;Gingerich and Voss 2005;Oki et al. 1998).  conducted a rigorous assessment of the island's groundwater recharge and asserted the need for the protection and management of groundwater resources based on a water budget analysis. Moreover, the impacts of groundwater pumping and recharge were modeled for Manukan Island in eastern Malaysia , and it was determined that groundwater pumping should be reduced to 75% of the current volume. This would ensure sustainable groundwater use in the region with the implementation of artificial groundwater recharge. Chang et al. (2016) simulated a simplified numerical model of a barrier island in the United States to confirm whether the water from the groundwater pumping site conformed to drinking water standards for chloride, and estimated the water resources for future scenarios in which coastal water resources and anthropogenic activities had decreased due to climate change.
Coastal aquifer vulnerability assessments using precise scientific methods have been required to support decisionmaking regarding land use controls and water resource management. For example, Werner et al. (2012) used a quantitative technique to modify analytical solutions to assess the SWI vulnerability in Australian coastal aquifers. Although general groundwater vulnerability assessment techniques, such as the DRASTIC and SINTACS models, have been applied to coastal aquifers (Polemio 2016), new techniques more suited to the characteristics of coastal aquifers have also been developed including the GALDIT model (Lobo-Ferreira and Chachadi 2005). Previous studies have compared the GALDIT method with index-based methods, including DRASTIC, SINTACS, and AVI (Saidi et al. 2014;Kardan Moghaddam et al. 2016;Luoma et al. 2016;Kura et al. 2015), indicating that GALDIT can be used to assess the extent of aquifer vulnerability to SWI in coastal and island regions (Saidi et al. 2014;Kardan Moghaddam et al. 2016;Luoma et al. 2016;Kura et al. 2015;Recinos et al. 2015;Trabelsi et al. 2016;Pedreira et al. 2015).
Jeju Island is the largest island in Korea that shows the highest precipitation compared to other inland regions, with a range of 1142-1967mm/year (Jeju Special Self-governing Province, 2013. Heavy rainfall in the region leads to abundant areal recharge into the island's subsurface. Few studies have assessed the vulnerability of groundwater resources on Jeju Island. For example, El-Kadi et al. (2014) conducted a numerical study of the island to assess water resource sustainability and predicted dried springs in the island under drought scenarios. Furthermore, Chang et al. (2019) conducted a GALDIT-type vulnerability assessment for Jeju Island, revealing that groundwater vulnerability is increasing for certain parts of the island.
The objective of this study is to comprehensively understand changes in the coastal groundwater resources in northeastern Jeju Island, South Korea, due to groundwater exploitation. A new approach is used, involving three-dimensional numerical simulations and a GIS-based vulnerability assessment with two dynamic parameters, i.e., groundwater level and salinity concentration. We examine the drivers of SWI by developing a variable-density dependent model of a representative area on eastern Jeju Island. Detailed steady-state and transient SWI models are developed in regions with deep groundwater abstraction wells to describe groundwater dynamics, as well as the shape and movement of the saltwater front. Planned increases in pumping rates are incorporated to project future SWI under various climate scenarios. Subsequently, a groundwater vulnerability assessment is conducted to predict changes in the SWI vulnerability status of coastal aquifers. To the best of our knowledge, this is the first study to integrate the results of a detailed SWI numerical approach with a GIS-based SWI vulnerability assessment to explore the impact of anthropogenic activity on a coastal groundwater system. Figure 1a shows the position of Jeju Island, which is located to the south of mainland Korea. The island has an area of approximately 1850 km 2 and a maximum elevation of 1950 m at the crest of Hallasan Mountain. The terrain is characterized by coastal plains with a steeply sloping mountainous interior and diverse spatiotemporal distribution of rainfall. The Jeju Island aquifer supplies most of the island's freshwater (Koh et al. 2017). Only a few perennial streams exist on the island because of its well-drained geological structure composed of porous volcanic rocks with high permeability. Hallasan Mountain is the largest unpolluted groundwater recharge area on the island, and its low-permeability (Seogwipo) layer plays an important hydrogeological role (Yoon et al. 2005). Figure 1b shows the aquifer system on Jeju Island. The unconfined aquifer corresponds to a basaltic formation, and the strata of the Seogwipo Formation are generally classified as aquicludes or low-permeability layers (Koh et al. 2017). Koh (1997) analyzed core logging data from Jeju Island and concluded that groundwater levels were mainly determined by topography because the groundwater tends to flow seaward from the mountainous area located at the center of the island and that the distribution of the subsurface Seogwipo Formation is a secondary factor that influences groundwater levels because it acts as a barrier to groundwater flow.

Site description
Groundwater recharged into the upper and middle areas of Hallasan Mountain passes through the underground permeable layer to the coastal area, where the groundwater level is near the surface. Groundwater levels on Jeju Island range between approximately -26 and 282 m (Won et al. 2005). The water levels were found to be lower in the western and southern basins by manually measuring 300 groundwater wells from August 20 to September 27, 2002(Won et al. 2005. The island is predominantly composed of volcanic pyroclastic layers, clinker layers, and primary and secondary effective voids in the basalt (Hahn et al. 1997).
Groundwater development is active in the Sungsan and Pyoseon watersheds in the eastern and northern areas of Jeju Island, respectively (Won et al. 2005). Youn et al. (2003) and Kim et al. (2005) confirmed the location of the saltwater/freshwater interface in the Handong-ri watershed in the eastern area of the island using three observation wells. They reported a sharp increase in the electrical conductivity at 37, 31, and 64 m below sea level at the Handong 1, 2, and 3 observation wells, respectively. Youn et al. (2003) used electrical conductivity and temperature values from the doublesloped section to confirm the existence of the reservoir at 44.1, 37.2, and 74.5 m below sea level at the Handong 1, 2, and 3 observation holes, respectively. Furthermore, they found that the interface was shallower than the theoretical  Koh et al. (2017)) estimates calculated using the Ghyben-Herzberg principle (Ghyben 1888; Herzberg 1901). Kim et al. (2005) calculated the hydraulic characteristics of the coastal aquifer using a tidal response technique and concluded that tidal effects on eastern Jeju Island extend up to 3-5 km off the coast. The integrated surface water-groundwater hydrologic analysis SWAT-MODFLOW model was tentatively applied to the Pyoseon watershed to estimate groundwater recharge (Kim et al. 2008(Kim et al. , 2009. The study site is an area of 13.9 km 2 located in the northeastern area of Jeju Island. Figure 2a shows the locations of the study area on the island. In this area, brackish groundwater pumping wells are concentrated near the coast, due to which the use of brackish groundwater for heating and fishery has gradually expanded since 2000. Black circles represent the Seawater Intrusion Monitoring Network (SIMN).
There are approximately 150 wells within SIMN in Jeju, some of which have been operating since 2001. Figure 2b shows the locations of the drilling sites in the study area, and Fig. 2c shows the lithostratigraphy of the drilling core obtained at the study site. The sequences of volcanic and sedimentary rocks observed in the geological logs analyzed in this study were incorporated into the numerical model layers, in which the dominant rock types (i.e., basaltic, hyaloclastite, lower basaltic, and sedimentary) were identified. Basaltic rock is the dominant geological formation, extending from the surface to approximately 40 m below sea level at the study site, followed by hyaloclastite and thin discontinuous sedimentary formations. The hyaloclastite and lower basalt extend from 50 to 95 m below sea level, and the sedimentary formation extends from 90 to 132 m below sea level. Koh et al. (2019) analyzed the MW 2-1 drilling core and simplified its lithostratigraphy. They defined the basalt, hyaloclastite, and the following lower basalt as a hydrologically connected volcanic rock. From the well logs (i.e., the gamma-ray and conductivity logs) for the 2-101 m interval of the drilling core, they concluded that the sedimentary formation (unit F in Fig. 2b) is the Seogwipo Formation in this study. The location of the Seogwipo Formation determined herein is consistent with the locations reported in previous studies that determined the depth of the Seogwipo Formation (Chang et al. 2019;Hahn et al. 1997).
In the study area, monitoring wells are not evenly distributed. They are concentrated along the shoreline to detect the seawater intrusion near the industrial park. Two types of monitoring data were used in this study site. One SIMN monitoring well (SIMN JD-HD1), which was installed before this study commenced within the study site, provided 20 years of historical groundwater level data to calibrate areal-recharge input data based on the seasonality and longterm trend. The second type of monitoring data used in this study was obtained from new monitoring wells (MW_1, MW_2-1, MW_2-2, MW_3, MW-4-1, MW_5, MW_6-1, MW_6-2) installed since 2017. These monitoring wells provided groundwater, EC, and salinity data. The spatial distribution of groundwater levels and multi-depth salinity data at a certain time (June 2018) were used for calibrating the numerical model in this study. Figure 3 shows the temporal trend of groundwater level and EC data collected from the MW_1, MW_2-1, MW_2-2, MW_3, MW-4-1, MW_5, MW_6-1, and MW_6-2 wells. The dots in the figure represent the measured data using real-time sensors that collected the data every 10 min. The left axis represents the scale of the groundwater level and the right axis represents the scale of EC. Because the monitoring wells are close to the ocean, the groundwater level ranged between 0 and 2 m above sea level. The groundwater levels increased slightly during July and August; however, the EC trend did not show any Fig. 2 a Location of the SIMN wells on Jeju island, b location of drilling sites in the study area, and c lithostratigraphy of the drilling cores significant change. Note that quality data should consider the monitoring depth because it varies vertically. The realtime sensor is located at different depths for each monitoring well, due to which it does not enable the comparison of EC values directly. Instead, the multi-depth measurement for salinity was also performed, which will be presented in the calibration section.

Numerical analysis via SEAWAT
A detailed numerical study using the MODFLOW-family computer code SEAWAT (Guo and Langevin 2002) was conducted to evaluate coastal aquifer systems on Jeju Island. The model was discretized into uniform grids of 50 m × 50 m, consisting of 105 rows, 91 columns, and 7 layers. Figure 4 shows a conceptual diagram of the numerical model indicating the boundary conditions. The coastal area surrounded by seawater was set as a constant boundary, with a concentration of 35 kg/m 3 and a groundwater level of zero. The inland boundary was set to the no-flow condition to let the head at the boundary automatically adjusted in response to the amount of areal recharge. During the simulation period, seasonal variations in the recharge rates were applied to the model using monthly data. The wells, bounded by orange circles, were applied to future pumping increase scenarios at sites that currently pump brackish water in the mixing zone.
Aquifer tests around the study site indicate a hydraulic conductivity of 85-470 m/day. The hydraulic conductivity was adjusted to 320 m/day during calibration. The geologic units analyzed from Fig. 2c were assumed to be hydrostatically connected as one unconfined aquifer, previously defined as the Basaltic Formation in Fig. 1b. The conceptualized numerical model subdivides one uniform hydrogeological unit into seven numerical layers to obtain a better resolution of the vertical seawater-freshwater interface position. The sedimentary formation represented as unit F in Fig. 2c is considered the impermeable layer, previously defined as the Seogwipo Formation in Fig. 1b. Therefore, the bottom of the numerical model is considered the top of the sedimentary formation. The vertical hydraulic conductivity was set to 2% of the horizontal hydraulic conductivity values.
Layer 1 was assigned a specific yield (S y ) of 0.1. The specific storage (S s ) value assigned to layers 2 to 7 was 0.00001 m -1 . The porosity was set to 0.3. To minimize the effect of variance, the dispersity coefficient was set to zero. The densities of freshwater and seawater were 1000 and 1025 kg/m 3 , respectively. The concentration of saltwater was set to 35 kg/m 3 with a density slope of 0.714. The fully implicit finite difference solver package was used in all our simulations. Table 1 summarizes the model parameter values used in this study.
After adjusting the average values of historical data from 2001 to 2019, the recharge value was applied to the predevelopment model without pumping to reach a steady state of the head and salinity for the model. The spatial distribution of the simulated head and salinity from the steady-state predevelopment model were used as the initial condition for all transient SEAWAT modeling runs.
The transient model was run for 2000-2019 according to the model input data, with the pumping and recharge from 2001 to 2019. The monthly recharge rate of the Gujwa watershed (location of the study site) ranges between 35 and 65% of the precipitation when calculated based on the SWAT model (KICT 2015;Arnold et al. 1996). January presents the lowest monthly recharge rate (calculated as 35% of precipitation), whereas August presents the highest monthly recharge rate (calculated as 65% of precipitation). To use the calibrated recharge rate in the SEAWAT simulation, the average monthly watershed recharge rate was multiplied by the precipitation data from 2010 to 2019 to calculate the average amount of monthly recharge, which was input into the numerical model (Fig. 4c). Through trial and error, we found that 95% of the recharge rates reported in a previous study (KICT 2015) fit well with the observed data. The calibrated recharge values are similar to those reported by Jeju Special Self-governing Province (2013), in that the eastern sector of Jeju Island showed the highest recharge rate among all sectors, with an annual average recharge rate of 47.8% based on the latest 20-year average data from 1992 to 2011. This high recharge rate is due to the geologic character of the island that features complex and highly permeable volcanic structures with a significant amount of rainfall affected by oceanic climate. Figure 5 shows a multistage SEAWAT simulation process. The steady state of the SEAWAT model was simulated by assuming the predevelopment period. Using the simulation results of the predevelopment period as an initial condition, changes in the SWI due to groundwater development from 2029.10

SWI vulnerability assessment via GALDIT
An analysis of the SWI vulnerability was conducted for the aquifers on the northeastern coast of Jeju Island. The GALDIT technique (Chachadi 2005) is an overlay and ranking vulnerability assessment tool using GIS that visualizes the SWI vulnerability of coastal aquafers. This technique uses the aquifer characteristics of hydraulic conductivity, groundwater level, distance from the shore, current level of SWI, and aquifer depth as the main assessment factors. In this study, the vulnerability assessment was performed by collecting groundwater data and hydrogeological parameters over several years to visualize the current extent of SWI. The GALDIT index was calculated for each observation point. In this study, the calculated points were converted to the Thiessen polygon in the map using GIS. Chang et al. (2019) modified the original GALDIT index parameters to describe the characteristics of SWI on Jeju Island and reflect decreasing groundwater levels due to anthropogenic activities. Table 2 provides the weight, variable range, and importance rating of the individual parameters. The GALDIT index was determined by multiplying the weights by the rating parameters assigned to each data point as follows: (1) In this study, the salinities or saline concentrations simulated by SEAWAT should be converted to EC data and used for the parameter I in the GALDIT assessment. For the conversion, we compared the monitoring data for electrical conductivity and salinity from the monitoring well sensors used in this study with more than 100,000 data points, which were measured (monitored) once every 10 min. The EC and salinity data have a relatively linear relationship. Based on this comparison, 1000 µS/m corresponds to 0.498 kg/m 3 (or practical salinity unit, PSU), 2000 µS/m corresponds to 1.028 kg/m 3 , and 3000 µS/m corresponds to 1.575 kg/m 3 . Similarly, a previous study by Chang et al. (2019) compared the EC and chloride ion concentration, as well as EC and mole fraction to convert different indicators for representing the groundwater quality of coastal aquifers. Table 3 shows the index value range according to the vulnerability status based on the GALDIT assessment.

SEAWAT simulation
The SEAWAT model comprises three sequential processes: steady-state simulation, historic transient-state simulation, and future-scenario simulation. The first step determines predevelopment by stabilizing the salt wedge with no pumping well installations. Next, the transient model was run for 2000-2019. Figure 6a shows a vertical view of the salinity distribution in the A-A' cross-section in the model according to the transient simulation results for 2019. In the figures, red indicates saltwater, blue indicates freshwater, and yellow to turquoise indicate brackish water, i.e., mixing of saltwater and fresh groundwater. The saltwater wedge migrates over 3000 m inland.
As groundwater exhibits seasonal variations, the dry and wet seasons were compared in this study. Figure 6b shows the results of the transient simulation for 2019 for the dry season and wet seasons. White lines represent the contours of the groundwater levels. Pumping areas near the coast exhibit a decrease in the groundwater below sea level, but this groundwater level change is relatively small, despite heavy pumping of brackish groundwater. In the dry season, groundwater levels in inland areas decrease by approximately 1.0 m more than during the wet season. The colored saline distribution showed that the SWI did not fully advance in the area near the brackish pumping field.
Calibration processes in SWI studies generally follow two simultaneous steps using the groundwater level and groundwater quality (electrical conductivity, salinity, or chloride ion concentration) data. In this study, the long-term transient simulation results were compared to groundwater level data observed at SIMN JD-H1. As the sensing depth for the electrical conductivity data at SIMN JD-H1 is not known, the EC data were not compared to the simulated results.
We first explored the sensitivity of the long-term groundwater data observed at SIMN JD-H1 to various hydraulic conductivity (K) values. Figure 7 compares the simulation results with various simulation tests. Figure 7a shows the daily precipitation data at the study site for the period from 2001 to 2019, which indicates the seasonal variability characterized by high rainfall in the summer and low rainfall in the fall and winter. Figure 7bshows the total pumping of groundwater development occurring in layers 4 and 5. Pumping wells range between -40 and -70 m, which correspond to the distance between layers 4 and 5. Figure 7c shows the monthly recharge used as the input for the numerical simulations during 2001 and 2019. The subsequent calibration process used the peak altitude and lowest groundwater level to calibrate the major hydrogeological parameters, including the hydraulic conductivity (K) and specific yield (S y ). Based on these long-term observations, the hydraulic conductivity and recharge rate were calibrated. In Fig. 7d, the K value was varied from 200 to 450 m/day. In all sensitivity simulations, the other parameters were fixed at the base-case levels listed in Table 1. The data show that, when the K value was small, the simulated groundwater level was relatively higher than when the K value was larger. The maximum intrusion length was also small for smaller values of S y . When K was large (at K = 450 m/day), the seasonal change in the groundwater level was relatively small. The results show that when K = 320 m/day, the simulated data matched relatively well with the observed data, even from 2001 to 2019. The simulation satisfactorily described the magnitude of the seasonal variance between the highest peak and lowest groundwater level, by showing that the peak height mainly varied with respect to the recharge values. The difference between the simulated and observed level was approximately 0.2-0.3 m in the highest seasonal peaks. Figure 7e presents the temporal variations in the groundwater level according to the S y values. The impact of S y on the seasonal groundwater level was examined by varying S y within the range from 0.05 to 0.3. The seasonal groundwater level was relatively insensitive to S y , and the Table 3 Indication of groundwater vulnerability to SWI in terms of the GALDIT index (modified from Lobo-Ferreira and Chachadi (2005)) Index value range 2.5-5 5-7.5 7.5-10 GALDIT Vulnerability Low vulnerability Moderate vulnerability High vulnerability simulation best fit among several sensitivity tests when S y = 0.1. Additionally, sensitivity tests for specific storage values were conducted in the magnitude range from 0.00001 to 0.0000001 m -1 . The results show that the change in the groundwater level was also insensitive to S y values. Additional calibration was conducted using the transient simulation results and groundwater monitoring data from seven monitoring wells in June 2018. Figure 8 compares the results obtained from the transient numerical model to the average monthly groundwater monitoring data from the observed groundwater levels of monitoring wells. The calibration results indicate that the simulated groundwater levels were similar to the observations, except for MW 1 (Fig. 8a). It is regarded that the observation for MW 1 may possibly contain a measurement error because the observed groundwater level was higher than 1 m at a location near the ocean. Results were obtained for sensitivity tests for the groundwater monitoring well at various K and S y values, with the base case showing R 2 values of 0.96. In addition to comparing the water levels, the simulated concentration levels were compared with data from the vertical salinity profiles of monitoring wells (Fig. 8b). Figure 8b shows the vertical salinity profiles collected from selected monitoring wells. Certain salinity discrepancies occur in the saltwater concentration range of 10-90%. However, most of these discrepancies are minor and can be attributed to sharp changes in the salinity levels inside the narrow seawater-freshwater mixing zone in the model.
The simulation was extended to explore the future impacts of increased groundwater pumping after 2019. Future projections used identical boundary conditions and model parameters as the historical simulations. Seasonally repeated recharge data were averaged from the historical  Figure 9 shows these recharge patterns and pumping patterns. SEAWAT simulations were allowed to evolve for ten years under four scenarios. The future scenario simulated the effects of increased pumping and climate from November 2019 to October 2029. The increased pumping rate applied to the water development scenarios was calculated according to actual plans for groundwater expansion in the region. The scenario was configured to simulate excess groundwater demand due to the expected growth of industrial water use at the study site. Assuming this relationship is applicable to the water demand in scenario 1, the pumping rates increased from 1000 to 10,000 m 3 /day for one well from November 2019, followed by an increase to 20,000 m 3 /day from November 2020. The well "industrial park" that was expected to see an increase in pumping is located relatively inland, away from the other brackish groundwater pumping wells. Assuming this relationship is also applicable to the groundwater demand in scenarios 2 and 3, these future simulations were compared with the baseline scenario, which simulated groundwater changes under current groundwater pumping conditions over ten years. Scenario 2 explored the impacts of a dry climate, in addition to increased pumping. A dry climate recharge scenario was obtained using the recharge pattern for 2017, where the total amount of recharge for the next ten years of the simulation period was only 44% that of the base-case scenario, which is the lowest value in the historical recharge data from 2000 to 2019. Scenario 3 explored the impacts of a wet climate, in addition to increased pumping. A wet climate recharge scenario was obtained using the recharge pattern for 2012, where the total recharge was 40% higher than that in the base-case scenario, which is the highest value in the historical recharge data from 2000 to 2019. Figure 10 shows a cross-sectional view of the saline concentration profile across the south-north A-A' cross-section, exhibiting the detailed SWI pattern simulated under the four scenarios. The continuous lines in the figure represent the 50% concentration profile in the saltwater-freshwater mixing zone, defined as a saltwater wedge. To determine the position of the saltwater wedge, a sectional cut in the model in column 50 was selected. The x-axis represents the distance from the ocean boundary in the south-north direction of the A-A' cross-section. The pumping well is located 1750 m from the left (sea) boundary. This well increased the pumping rate in future scenario 1. The spatial extent of the intruded saltwater wedge is generally regarded as one of the representative indicators of SWI because its behavior, such as intrusion or recession, would depend on hydrological stress, geological aquifer parameters, or anthropogenic activities including recharge rate, hydraulic conductivity, pumping, and so on. The results show that the wedge is located at a depth of 10 to 65 m below sea level in the curved shape. This study analyzes how the wedge moves vertically and observes how the vertically upward movement or slope change of SWI is expected to also affect the horizontal intrusion in this numerical model. The graph showed that the simulation results for the base-case scenario are similar to the results for the historic simulation, indicating that the aquifer exhibits no additional SWI during the base case scenario, and the current pumping condition under an average climate does not have a substantial impact on SWI movement over ten years. In the shallow region of the aquifer, the results show that the interface moves upward under scenario 1, indicating the impact of increased pumping under an average climate. Compared with the results for scenario 1, we can clearly see the effect that a dry climate has on the SWI in scenario 2. The entire wedge interface moved upward under scenario 2, indicating the impact of increased pumping under a dry climate. In the deep inland region of the aquifer, the results for scenario 3 show that the interface moved downward with increased pumping for a wet climate. The results showed that the slope of the wedge for scenario 3 is less steep than that for the other cases and that SWI would occur more into the inland. Based on the projection of the slope, it was approximately calculated that the dry climate would induce more than 2 km of horizontal SWI compared to the base case at 100 m below sea level, which is the bottom in the model.

GALDIT SWI vulnerability assessment
Six GALDIT parameters were indexed according to the categories described in Table 2. The GALDIT assessment allocates an importance rating of 2.5, 5, 7.5, or 10 according to the range of each rating. Figure 11 shows the spatial distributions of the importance ratings of the GALDIT parameters, representing (a) groundwater occurrence, (b) hydraulic conductivity, (c) distance from the shore, (d) saturated aquifer thickness, and (e) impact of the existing seawater intrusion status. As a higher weight indicates a larger impact on the SWI, the height of the groundwater level (L) and distance from the shore (D) have a weight of 4, such that they have the largest impact among the six GALDIT parameters. Based on Table 2, the testbed of this study was an unconfined aquifer, and all sites were assigned an importance rating of 7.5 for the groundwater occurrence parameter (G). The horizontal hydraulic conductivity (320 m/day) corresponded to an importance rating of 10 for the aquifer hydraulic conductivity parameter (A). The depth distribution considers that the majority of the aquifer is deeper than 10 m, corresponding to a saturated aquifer thickness (T) importance rating of 10. As the study site is relatively small and does not show marked variations in geological formations, a spatial map of T ratings is homogeneous. The salinity levels are the lowest inland and higher near the shore. Based on the relationship between the EC and salinity data previously observed in Fig. 6, the average monthly salinity was determined to range from 8950 to 53,000 uS/cm for all sites in the model, corresponding to an importance rating of 10, which represents the impact of the existing SWI (I) status. Consequently, most areas in the study site are assigned an importance rating of 10, which shows the impact of the existing SWI. As the SWI previously occurred at the bottom of the aquifer and did not result in any substantial changes in the concentration in low areas of the aquifer in future scenarios, the importance rating of 10 remains constant when simulated via the SEAWAT based on future scenarios. Table 4 lists the predicted groundwater level from SEA-WAT, depicting what the values were converted to when they were indexed according to the GALDIT parameter ranges. The values in the parenthesis in the table represent the change in the groundwater level for the scenarios compared with the base-case scenario output. In the base case, the spatial distribution of the monthly groundwater level (L) above sea level ranges from -0.12 to 1.00 m, corresponding to importance ratings of 5-10. Under Scenario 1, there was an approximately -0.01 to -0.15 m decrease in the groundwater level, indicating that increased pumping development-induced changes in the groundwater level. Scenario 2 was the least favorable scenario for recharge prediction, where the groundwater levels at all sites further decreased due to climate factors, as compared with Scenario 1. However, in Scenario 3, wet Bars represent the recharge rates (mm/month), and circles represent the pumping rate in the industrial park (m 3 /month) climate factor predictions compensated for the impact of increased groundwater development, such that the groundwater level returned to levels similar to those in the base case. Figure 12 shows the spatial distribution of the GALDIT assessment results under different future scenarios calculated using Eq. (1) and Table 2. These future GALDIT values correspond to the temporal variance of the groundwater level and salinity based on the pumping scenarios with different recharge patterns. As the SEAWAT simulation indicates that the base-case scenario did not exhibit severe SWI compared with the 2019 conditions, as previously observed in Fig. 10, the SWI vulnerability status for 2019 was assumed to be similar to that of the base-case scenario. Figure 12a shows the future SWI vulnerability based on SEAWAT predictions under the base-case scenario. According to Eq. (1), the computed GALDIT values range from 6.5 to 9.8, indicating moderate to high vulnerability. Wells located inland generally exhibit the lowest index values (6.5). Figure 12b shows the future SWI vulnerability based on groundwater levels predicted by SEA-WAT under future scenario 1, which has an average climate with increased pumping. As the SEAWAT simulation under scenario 1 indicates aquifer conditions with a decrease in the groundwater level compared with the base case, the SWI vulnerability status worsened west of the model center due to increased pumping. Thus, additional groundwater abstraction increases the vulnerability of the coastal aquifer due to the decline in the groundwater levels. Figure 12c indicates that an equal amount of increased groundwater development occurred under a dry climate. Figure 12d shows the GALDIT index map of scenario 3, which has the same amount of increased groundwater development under a wet climate, with respect to the aquifer vulnerability. The vulnerability showed a similar status as future scenario 1. The vulnerability map showed no additional deterioration in terms of SWI under future dry and wet climate scenarios compared with the future average climate, indicating that the recharge condition did not have a substantial impact on the SWI of the aquifer under a future groundwater development plan with an increase in the pumping from 1000 to 20,000 m 3 /day. The reason why the severe change is not exhibited in the future dry climate is because a high concentration of salinity already exists in the area, so only the groundwater level has an impact on the vulnerability status change. A further SWI of approximately 2 km inland is expected, according to the SEAWAT simulation results. Even though we did not extend the map, we could expect an escalation in the vulnerability status in the dry climate scenario at the south of the model boundary.

Conclusions
The sustainable development of groundwater resources is challenging under conditions of rapid economic growth, particularly in island areas that are susceptible to SWI. Although coastal groundwater resources have been extensively exploited in the study area, previous studies have not yet analyzed the effect of groundwater development due to pumping wells installed in the saltwater-freshwater mixing zone. Furthermore, plans for this area involve increasing pumping by 20,000 m 3 /day, which presents the possibility of severe SWI in the future in the eastern area of Jeju Island. This study conceptualized the study site through field investigations and simplified hydrogeological units in the numerical model. Numerical models were used to simulate SWI in the coastal aquifers to describe current aquifer conditions and predict future conditions based on different water-use scenarios. In this study, we assessed aquifer changes by Fig. 11 Spatial mapping of the importance ratings for the GALDIT parameters for a groundwater occurrence, b hydraulic conductivity, c distance from the shore, d saturated aquifer thickness, and e impact of existing seawater intrusion status predicting the groundwater level and salinity distribution. Parameters selected using a reliable method were integrated into a GIS-based SWI vulnerability assessment tool. Moreover, effective methods were devised to transform numerically simulated data into a vulnerability assessment. In this manner, numerical modeling was successfully integrated with a GIS-based SWI vulnerability assessment, providing important insights into the management of the unconfined aquifer system in Jeju Island. The results from the integrated SEA-WAT simulation and GALDIT assessment indicated higher vulnerability due to future pumping in some areas of the study site; however, most of the areas would be relatively sustainable for the various future climatic conditions and pumping plans. The approaches used in this study can serve as a baseline to integrate more advanced technologies in future field studies. This study has limitations in the calibration process, because the number of monitoring wells is limited, and calibration was not conducted for the predevelopment period. Further investigations are required to develop a more detailed calibration process to describe the characteristics of SWI. Ongoing research will focus on more sophisticated scenario-based estimations to describe land use and land cover changes, as well as climate change scenarios. Considering  . 12 GALDIT map considering groundwater development and climate effects for the a base-case scenario: average climate with current pumping conditions, b future scenario 1: average climate with increased pumping, c future scenario 2: dry climate with increased pumping, and d future scenario 3: wet climate with increased pumping the regional impacts of climate change, one of the possible further studies is to develop an integrated method of estimating how the GALDIT parameter, 'the height of groundwater level above sea level (L)' can reflect the projection of sea-level rise due to future climate change. In addition, the impact that groundwater development in the saltwater-freshwater mixing zone has on coastal aquifers should be further assessed using multiple methods. Geochemical studies or analyses based on temperature indicators, which were not addressed in this study, are also an interesting topic for future research.

Compliance with ethical standards
Conflict of interest The authors declare no conflicts of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.