## Abstract

This paper assesses the extent of saltwater intrusion (SWI) in coastal aquifers in two sub-districts of Jeneponto Regency, Indonesia, for the first time using an analytical method in the form of a sharp-interface approach and numerical dispersive simulation of a finite difference solution developed by SEAWAT, based on field observations and secondary field data. In general, the SWI extents obtained from the analytical solutions for four aquifers in the region are around 30% greater than numerical solutions. The results show that the maximum SWI extent occurs in the Binamu2 aquifer in Jeneponto Regency, where the analytical solution and numerical simulation produced SWI extents of 850.4 m and 510 m, respectively. The higher results in the sharp interface model compared to the numerical model are most likely due to the pushing seaward effect caused by the presence of circulating flow inside the saltwater wedge below the freshwater-saltwater interface, which is considered in the numerical solution but not in the analytical solution, as explained in the previous SWI study. The results also revealed that the SWI conditions in Jeneponto Regency are still within reasonable limits, but groundwater exploitation control is critical in the region to maintain this situation.

## HIGHLIGHT

The higher results in the sharp interface model compared to the numerical model are most likely due to the pushing seaward effect caused by the presence of circulating flow inside the saltwater wedge below the freshwater-saltwater interface, which is considered in the numerical solution but not in the analytical solution, as explained in the previous SWI study.

### Graphical Abstract

## INTRODUCTION

Groundwater has shown the potential to be exploited in large quantities as a water resource alternative, with the caveat that the amount extracted cannot exceed the amount of water entering the aquifer (Yusuf *et al.* 2021). As a result, groundwater use must follow the yield capacity of aquifers, and exploitation must be governed by a set of policies at the national and regional levels.

The coastal area is topographically lowland with a coastal plain morphology (Bhattachan *et al.* 2018). Composing rock is typically alluvial deposits made up of clay, sand, and gravel that result from the transport and erosion of rocks in upstream river areas (Zhang *et al.* 2019). In general, the rocks in the plain are less compact, giving rise to the concept of reliable groundwater potential. Although unconfined aquifers can be a good source of groundwater, particularly in coastal areas, suitable coastal aquifers are generally in the form of confined aquifers. The significant issues in the coastal area are the diversity of aquifer systems and the position and distribution of seawater, both naturally and artificially caused by groundwater extraction for domestic, agricultural, and industrial needs. Changes in the hydrogeology of coastal areas can cause seawater to move towards land, potentially contaminating groundwater in aquifers, a condition known as a saltwater intrusion (SWI) (Badaruddin & Mehdizadeh 2021).

SWI has historically been caused by excessive pumping or extraction of groundwater on land, resulting in significant losses in groundwater availability in coastal aquifers around the world (FAO 1997; Badaruddin *et al.* 2015). However, the effects of climate change (e.g., sea-level rise and a decrease in groundwater recharge quantities) can also cause SWI (Post 2005). As a result, the vulnerability of coastal aquifers to climate change, increased groundwater extraction volume, and sea-level rise must all be factored into groundwater management.

The SWI is a complex process involving variable-density flow, solution transport, and hydrochemical processes (Werner *et al.* 2012), making groundwater assessment difficult and costly. As a result, large-scale assessments of coastal aquifer susceptibility to SWI generally rely on qualitative methods such as GALDIT (Lobo-Ferreira *et al.* 2007) and CVI (Ozyurt 2007) only take into account a limited number of factors thought to influence SWI. Furthermore, these methods generally lack theoretical foundations and focus more on selecting only one element related to the SWI. Werner *et al.* (2012) developed an alternative prediction of large-scale SWI. This method is based on the steady-state conditions of the Strack equation (1976), which assumes that the interface of saltwater and freshwater in an aquifer is a sharp line, so it involves the physical mechanics of the SWI, albeit under ideal conditions. Whereas numerical modeling, which assumes that the contact between saltwater and freshwater in an aquifer takes the form of a dispersive mixing zone, is required to obtain conditions close to actual conditions when predicting SWI.

Over the last decade, there has been rapid population growth worldwide, including in Indonesia, particularly in Jeneponto Regency in South Sulawesi Province, which has resulted in massive groundwater exploitation (Hatta *et al.* 2020). This phenomenon has harmed groundwater quantity and quality, including decreased groundwater levels, increased fluctuations, and decreased groundwater quality, as well as SWI in several regions (Asmal *et al.* 2020). Asmal *et al.* (2020) conducted groundwater salinity investigation in several groundwater wells in Bangkala District, Jeneponto Regency, and found that mostly the water in the wells that near from the coastline were brackish. In addition, Wahyuni *et al.* (2019) also investigated the SWI in the same district in Jeneponto Regency and discovered that saline water was identified around the coastal area i.e., 100 m to 1,000 m from the coastline. As a result, a concerted effort on the government, the public, and the private sectors are required to mitigate these negative consequences.

This study used Werner *et al.* (2012) analytical methods and a numerical dispersive model to predict SWI in two sub-districts of the Jeneponto Regency. For the first time, analytical methods and 2-dimensional numerical modeling were used to determine the extent of saltwater intrusion in this area. Jeneponto Regency is one of the districts in South Sulawesi, Indonesia, whose coastal region has a high risk of experiencing an SWI due to the region's long-term extraction of groundwater for domestic and irrigation purposes (Syamsuddin *et al.* 2009). Given the conditions mentioned above, it is critical to conduct preliminary assessments of SWI conditions in this area and understand the causes of SWI, which will help local decision-makers develop groundwater resource protection policies.

## METHODOLOGY

### Conceptual model

In this study, SWI length was determined only in four cross-sections at the locations shown in Figure 1 and only under steady-state conditions. The research was conducted in two districts in Jeneponto Regency, namely Binamu and Arungkeke. Due to the limitations of available hydrological and hydrogeological data (e.g., aquifer thickness, soil stratigraphy, recharge, and specific storage), secondary data from previous studies were still considered (i.e., Syamsuddin *et al.* 2009). Due to the lack of detailed boring log data that can provide a comprehensive description of the stratigraphic conditions of the soil in the study area, it is assumed that the aquifer type in the study site is an unconfined aquifer based on Syamsuddin *et al.* (2009). Table 1 contains the hydrogeological data used in this study's SWI analysis.

Parameter . | Aquifers . | |||
---|---|---|---|---|

Binamu1 . | Binamu2 . | Arungkeke1 . | Arungkeke2 . | |

K (m/d) | 5.80 | 5.80 | 5.80 | 5.80 |

h (m) _{s} | 68 | 68 | 64 | 64 |

GWL h (m) _{f} | 3.6 | 2.0 | 9.5 | 8.7 |

x (m) _{f} | 5,000 | 4,800 | 3,300 | 3,200 |

n (-) | 0.46 | 0.46 | 0.46 | 0.46 |

S (-) _{y} | 0.32 | 0.32 | 0.32 | 0.32 |

α (m) _{L} | 4 | 4 | 4 | 4 |

α (m) _{T} | 0.4 | 0.4 | 0.4 | 0.4 |

D (m_{m}^{2}/d) | 8.6 × 10^{−5} | 8.6 × 10^{−5} | 8.6 × 10^{−5} | 8.6 × 10^{−5} |

δ (-) | 0.025 | 0.025 | 0.025 | 0.025 |

W (mm/y) _{net} | 56.70 | 56.70 | 56.70 | 56.70 |

Parameter . | Aquifers . | |||
---|---|---|---|---|

Binamu1 . | Binamu2 . | Arungkeke1 . | Arungkeke2 . | |

K (m/d) | 5.80 | 5.80 | 5.80 | 5.80 |

h (m) _{s} | 68 | 68 | 64 | 64 |

GWL h (m) _{f} | 3.6 | 2.0 | 9.5 | 8.7 |

x (m) _{f} | 5,000 | 4,800 | 3,300 | 3,200 |

n (-) | 0.46 | 0.46 | 0.46 | 0.46 |

S (-) _{y} | 0.32 | 0.32 | 0.32 | 0.32 |

α (m) _{L} | 4 | 4 | 4 | 4 |

α (m) _{T} | 0.4 | 0.4 | 0.4 | 0.4 |

D (m_{m}^{2}/d) | 8.6 × 10^{−5} | 8.6 × 10^{−5} | 8.6 × 10^{−5} | 8.6 × 10^{−5} |

δ (-) | 0.025 | 0.025 | 0.025 | 0.025 |

W (mm/y) _{net} | 56.70 | 56.70 | 56.70 | 56.70 |

### Analytical solution

The analytical solution developed by Strack (1976) for locating the contact between seawater and freshwater in fixed flow conditions forms the basis of the method described by Werner *et al.* (2012) and will be summarized here. The conceptual model for unconfined and confined aquifers is depicted in Figure 2.

The water quantities in the aquifer domain as shown in Figure 2 consists of net groundwater recharge (*W _{nett}* [L/T]), which takes into account infiltration, evapotranspiration and pumping distributions, fresh water drained into the sea (

*q*

_{0}[L

^{2}/T]), and the lateral flow in aquifers from inland (

*q*

_{b}[L

^{2}/T]). Hydraulic pressure head (

*h*[L]) (groundwater table for unconfined aquifer conditions) relates to the depth of the location of the interface line of saltwater and freshwater (

_{f}*z*[L]) explained in the Ghyben-Herzberg relation;

*z*=

*h*

_{f}/

*δ*. Here,

*δ*[-] is the dimensionless ratio of density

*δ*= (

*ρ*−

_{s}*ρ*)/

_{f}*ρ*, and

_{f}*ρ*(=1,025 kg/m

_{s}^{3}) and

*ρ*(=1,000 kg/m

_{f}^{3}) is the density of saltwater and freshwater, respectively, resulting

*δ*= 0.025. The thickness of fresh water is

*h*[L], the bottom of the aquifer is

*z*

_{0}[L] below the average sea level and the hydraulic conductivity is

*K*[L/T].

*x*[L]) is obtained as shown by Equation (3):

_{t}*x*[L]) is obtained as shown by Equation (6):

_{t}### Numerical model

For predicting the SWI under numerical-dispersive solutions, numerical modeling is usually required to validate analytical calculation results and obtain a more realistic picture. In this case, two-dimensional modeling is used with SEAWAT version 4, which is used for flow with density variations and solution transport. This program uses a finite difference method with only modeling saturated flow. The numerical method and equations used in SEAWAT are described in Guo & Langevin (2002) and Langevin *et al.* (2008).

As shown in Figure 3, the conceptualization of aquifers from the four research locations used in the numerical model is configured as separate aquifers for Binamu1, Binamu2, Arungkeke1, and Arungkeke2 aquifers and is modeled in two dimensions with cross-sections perpendicular to the coast. In this study, all representations in the numerical model, including the boundary conditions, are head-controlled (Werner & Simmons 2009), and the energy level is assumed to be constant despite the effect of pumping groundwater.

For simplification, the domain model is uniformly discretized in carrying out steady-state simulations for the four aquifers, where for the Binamu1 aquifer uses 500 vertical columns with a width of 10 m and 76 horizontal layers with a thickness of 1 m, the Binamu2 aquifer uses 480 vertical columns with a width of 10 m and 76 1 m thick horizontal layer. Arungkeke1 aquifer uses 330 vertical columns with a width of 10 m and 78 horizontal layers with a thickness of 1 m, while the Arungkeke2 aquifer uses 320 vertical columns and 78 horizontal layers with 1 m thick. This discretization is consistent with the *Peclet* number, smaller than 4, as Voss & Souza (1987) recommended to reduce numerical oscillation. The condition of a particular energy head is assumed to be at the groundwater boundary condition on land, and the constant concentration condition is assumed to be at the coastal boundary conditions, with a concentration of saltwater of 35 kg/m^{3}. Preconditioned conjugate-gradient 2 (PCG2) and general conjugate gradient (GCG) are used as a solution for solution flow and transport equations. The differential scheme is used for advection solutions with Courant numbers of 0.75. Courant numbers smaller or equal to 1 are usually needed to limit numerical dispersal occurrence to achieve more accurate results (Zheng & Bennett 2002).

## RESULTS AND DISCUSSIONS

### Groundwater level measurement in the research area

Table 2 presents data on the water level in four observation points in Binamu and Arungkeke sub-districts in Jeneponto Regency. Based on water level observation data, it is noticed that the groundwater level at the observation site is quite varied, which may be caused by heterogeneous conditions of the porous medium. In this study, soil heterogeneity of aquifer was ignored because of financial constraints in acquiring these data; therefore, simplification in aquifer heterogeneity was adopted in this study.

Aquifer code . | Sub-districts . | Coordinates . | GWL from MSL (m) . | |
---|---|---|---|---|

Latitude . | Longitude . | |||

Binamu1/SDJP273 | Binamu | 5 °39′25″ | 119 °43′52″ | 3.6 |

Binamu2/SDJP54 | Binamu | 5 °39′33.5″ | 119 °43′46″ | 2.0 |

Arungkeke1/AK1/TP1 | Arungkeke | 5 °39′28.7″ | 119 °47′42″ | 9.5 |

Arungkeke2/AK2/TP2 | Arungkeke | 5 °39′2.4″ | 119 °48′10″ | 8.7 |

Aquifer code . | Sub-districts . | Coordinates . | GWL from MSL (m) . | |
---|---|---|---|---|

Latitude . | Longitude . | |||

Binamu1/SDJP273 | Binamu | 5 °39′25″ | 119 °43′52″ | 3.6 |

Binamu2/SDJP54 | Binamu | 5 °39′33.5″ | 119 °43′46″ | 2.0 |

Arungkeke1/AK1/TP1 | Arungkeke | 5 °39′28.7″ | 119 °47′42″ | 9.5 |

Arungkeke2/AK2/TP2 | Arungkeke | 5 °39′2.4″ | 119 °48′10″ | 8.7 |

#### SWI results from analytical method

The hydrogeological parameters listed in Table 3 together with Equations (1), (2), (3), and (6) are used in estimating the amount of groundwater flow discharge into the sea and also the theoretical length of SWI in steady-state conditions at each the aquifers studied are shown in Table 3. These results represent the theoretical conditions of SWI for an extensive duration (steady-state) based on recent hydrological and hydrogeological parameters. Relative to the SWI geophysical investigation conducted by Wahyuni *et al.* (2019) in the same regency but in different district, the results found in this study were in a good agreement with their results where saltwater was identified only in the area around 100 m to 1,000 m from the coastline.

Aquifer code . | Sub-districts . | Coordinates . | SWI extents from coastline (x) (m)
. _{T} | |
---|---|---|---|---|

Latitude . | Longitude . | |||

Binamu1/SDJP273 | Binamu | 5 °39′25″ | 119 °43′52″ | 607.0 |

Binamu2/SDJP54 | Binamu | 5 °39′33.5″ | 119 °43′46″ | 850.4 |

Arungkeke1/AK1/TP1 | Arungkeke | 5 °39′28.7″ | 119 °47′42″ | 234.8 |

Arungkeke2/AK2/TP2 | Arungkeke | 5 °39′2.4″ | 119 °48′10″ | 250.6 |

Aquifer code . | Sub-districts . | Coordinates . | SWI extents from coastline (x) (m)
. _{T} | |
---|---|---|---|---|

Latitude . | Longitude . | |||

Binamu1/SDJP273 | Binamu | 5 °39′25″ | 119 °43′52″ | 607.0 |

Binamu2/SDJP54 | Binamu | 5 °39′33.5″ | 119 °43′46″ | 850.4 |

Arungkeke1/AK1/TP1 | Arungkeke | 5 °39′28.7″ | 119 °47′42″ | 234.8 |

Arungkeke2/AK2/TP2 | Arungkeke | 5 °39′2.4″ | 119 °48′10″ | 250.6 |

#### SWI results from numerical model

For the parameters listed in Table 2, the steady-state conditions for Binamu1, Binamu2, Arungkeke1, and Arungkeke2 aquifers are achieved using an extensive duration of 300 years, until the end of the SWI is in a stable position (no change in position during the duration observation). For the four aquifers observed, the value of *xt* obtained was smaller than the results obtained from analytical solutions. A more significant SWI prediction by analytical solutions (thin line method) is expected to occur (for example, from the case of Pool & Carrera 2011). The flow circulation in the interface zone between freshwater and saltwater due to mixing between them both will suppress the intrusion boundaries towards the sea, which is considered in numerical methods but is ignored in analytical solutions.

As explained earlier, the parameter values used in the SEAWAT numerical simulation are shown in Table 2. For the purposes of determining the distance of intrusion and comparison with analytical solutions, the length of SWI (*x _{t}*) is defined as the intersection of 50% isochlor concentration line and the bottom of the aquifer. However, the 5 and 95% isochlor concentrations lines are still shown as additional information. As seen in Figure 4, the value of

*x*obtained for Binamu1, Binamu2, Arungkeke1 and Arungkeke2 aquifers is 390, 510, 160, and 190 m, respectively. In contrast to Arungkeke1 and Arungkeke 2, the length of the SWI is almost the same for Binamu1 and Binamu2, although the observation locations are only about hundreds of meters apart, but the

_{t}*x*value obtained from numerical solutions and analytical solutions is notably different. This is correlated with the results of observations of water levels that are notably different between Binamu1 and Binamu2 aquifers (i.e., 1 m) as obtained from field observations, which is most likely due to the exclusion of soil heterogeneity effects on SWI.

_{t}As shown in Figure 4, the higher the difference between GWL and MSL, the smaller the magnitude of SWI produced in the SEAWAT numerical model. The same pattern of SWI was also pictured in the sharp interface modeling results. However, a significant mismatch was noted for the results of these two methods in the exact location. For example, the magnitude of SWI in the Binamu 1 aquifer from the sharp-interface model was 607 m, while 390 m of SWI length was produced in the SEAWAT numerical model. Similar conditions were also seen in the other locations where the results of the sharp interface model were more significant than the results of the SEAWAT numerical model. This phenomenon is probably due to the circulation flow inside the saltwater wedge as the system's response to the difference of salt and freshwater density (Badaruddin *et al.* 2015). This phenomenon is considered in the dispersive SEAWAT numerical model and excluded in the sharp interface model. Badaruddin *et al.* (2015) concluded that under passive SWI conditions, the circulation flow in the saltwater wedge tends to resist the inland movement of SWI, and this condition is ended when the type of SWI is turned to be active. Compared to the previous geophysical research conducted by Wahyuni *et al.* (2019) and groundwater salinity investigation by Asmal *et al.* (2020), the results of this research is in a reasonable agreement with their results where saltwater was identified only in the area near the coast (i.e., 100 m to 1,000 m from the coast line). Since there is no enough field investigation of SWI (i.e., groundwater salinity profiling) has been conducted in the respected area, it is difficult to measure the performance of the present research using MAE (Mean Absolute Error) and RMSE (Root Mean Square Error).

## CONCLUSIONS

Some conclusions can be drawn from the results. One is that varying water levels in adjacent locations; for example, GWL in Binamu1 and Binamu2 differs to 1 m despite its proximity for only about hundreds of meters. It is most likely due to the heterogeneity of the soil layer at the study site, where the exact data is still unknown. The study's findings also show that the maximum SWI occurred in Binamu1 aquifers, with analytical solutions and numerical simulations producing SWI lengths of 850.4 m and 510 m, respectively.

In general, from the study results, the prediction results of SWI are considerably different on both analytical and numerical solutions where the value of SWI from the analytical method is better, by around 20 percent, than the value of SWI from numerical models. It is due to differences in the review of the two methods, where the presence of a mixing zone between saltwater and freshwater calculated in numerical methods which are ignored in analytical solutions, giving a pushing seaward effect (due to the influence of flow circulation in the mixing zone) on SWI boundaries.

Although there are differences between analytical and numerical solutions, it was evident that the results of analytical solutions are still functional in the initial assessment of the SWI due to its rapid results. Unlike numerical solutions, it takes longer to prepare the resources needed to perform modeling simulations to predict SWI (e.g., software and domain model setup). It can also be concluded that the availability of secondary data in the study area can be considered to be very minimal, especially data on hydrogeological conditions that have not been explored in detail and thoroughly. Further investigation must obtain more accurate and comprehensive results in predicting SWI. Since there is a lack of SWI research has been conducted in the region, therefore, the results obtained in this research are significant as a preliminary SWI investigation. However, involving different boundary conditions such as sea level rise and gradual groundwater decline effects on SWI are also important to be explored and therefore, may become the future subjects on this type of research.

## CONFLICTS OF INTEREST STATEMENT

No conflicts of interest.

## DATA AVAILABILITY STATEMENT

All relevant data are included in the paper or its Supplementary Information.

## REFERENCES

*Vulnerability of Coastal Areas to Sea Level Rise: A Case Study on Goksu Delta*

*Simulation of Ground Water Level Fluctuations in the Coastal Area of Jeneponto*