VD models are based on the spatial variability of groundwater density, which ranges from saline water density to freshwater density. The density changes because of the mixing of freshwater and seawater, which results in the existence of a transition zone across the entire coastline. The width and exact position of the zone depends on the aquifer parameters and the pumping regime. In the current paper, thermal and viscosity effects are neglected and the density changes are attributed only to concentration effect. The mathematical formulation of the VD model has been thoroughly analyzed in several works (e.g., Frind, 1982; Kolditz et al., 1998; Younes et al., 2009; Pool and Carrera, 2011; Kourgialas et al., 2015; Fahs et al., 2022). In this study thesis the notation proposed by (Guo and Langevin, 2002; Langevin et al., 2008) is followed in order to be consistent with the SEAWAT numerical code, which is utilized for the VD simulations.

In general, the flow and solute-transport equations are used to describe mathematically the VD model. The two equations form a coupled differential equation system, which could be expressed as follows

where ρ is fluid density, q is the Darcy velocity vector, ρ_s is the density of water entering from a source of leaving through a sink, q_s is the volumetric flow rate per unit volume of porous medium representing sources and sinks, S_f is the specific storage, h_f is the freshwater head, n is the porosity, C is the solute concentration, D is the hydrodynamic dispersion tensor, v is the fluid velocity vector, and C_s is the solute concentration of water entering or leaving through sources and sinks respectively. Since solute reaction is not considered, fluid density is only a function of the solute concentration C, according to the following equation

in which ρ_o is the freshwater density, ϵ is the density difference ratio (Equation 3), C_o is reference concentration, and C_s is the maximum concentration. In this paper, the following values are used for the parameters of Equation (3): ρ_o = 1.000kg/m³, C_o = 0kg/m³, and C_s = 35kg/m³.

The density difference ratio is expressed as

where ρ_s stands for the maximum seawater density. In this study, we consider ρ_s = 1025kg/m³.

The Darcy flux term q of Equation (1) for constant viscosity and freshwater properties could be expressed as

where q_x, q_y, and q_z are the components of the specific discharge in the principal directions, K_fx, K_fy, and K_fz are the components of the freshwater hydraulic conductivity in the same directions and ρ_f is the freshwater density.

Equations (1)–(5) are the mathematical representation of VD approach of seawater intrusion. The well-established SEAWAT code is used to solve numerically the above equation set. SEAWAT is a modular finite-difference computer code created by USGS, which couples MODFLOW and MT3DMS, in order to solve iteratively the fluid flow and solute transport equations (Guo and Langevin, 2002).

The SI model employed in the present work is based on the single flow potential formulation proposed by Strack (1976). The Strack's approach considers three basic assumptions: (1) the Dupuit-Forchheimer approximation, (2) Ghyben-Herzberg approximation, and (3) Hydrostatic conditions in the saltwater zone. For steady state conditions, the incorporation of the three assumptions result in a simplified one-fluid, 2D SWI model, which, for the general case of irregular-shaped aquifers with spatially variable parameters (e.g., hydraulic conductivity and recharge), could be solved numerically using a groundwater flow code (Mantoglou et al., 2004). The differential equations used to calculate the flow potential in an unconfined coastal aquifer has the following form

where K = K(x, y) is the hydraulic conductivity of the aquifer, which is usually considered as a function of the spatial coordinates (x, y), N is the aquifer recharge, and Q(x, y) represents the distributed well pumping rates. The flow potential ϕ_τ at the location of the toe in an unconfined aquifer is calculated from the following equations (Mantoglou, 2003).

where d denotes the thickness of the aquifer from the horizontal bottom to the sea level. From the flow potential value field we are able to calculate the exact location of the sharp interface toe, which is of considerable interest in coastal aquifer management problems.

The SI model proposed by Strack (1976) and adopted in several studies (e.g., Werner et al., 2012; Ketabchi and Ataie-Ashtiani, 2015) have been widely used to simulate SWI in coastal aquifers or pumping optimization problems in coastal aquifers (e.g., Mantoglou, 2003; Mantoglou et al., 2004). However, it has been observed in several previous studies (e.g., Dausman et al., 2010; Koussis et al., 2015; Kopsiaftis et al., 2019a) that for a specific set of pumping rates the SI model tends to overestimate the extent of the seawater wedge compared to the variable density approach, which is considered as a reference result for the estimation of seawater intrusion. Llopis-Albert and Pulido-Velazquez (2014) further examined the efficiency of the Ghyben-Herzberg based SI model in producing comparable results with the VD model, within a range in the values of critical hydrogeological parameters and defined several limitations in the validity of the SI model.

In order to eliminate the discrepancy between the two models, several correction methods have been recently proposed (e.g., Pool and Carrera, 2011; Koussis et al., 2012, 2015; Lu et al., 2016; Werner, 2017; Koussis and Mazi, 2018; Christelis et al., 2019), which aim to incorporate the effect of the dispersion mechanism and improve the accuracy of the SI approach. Pool and Carrera (2011) were the first to suggest a novel formula to extend Strack's equations to the case of the mixing zone. In particular, they proposed an multiplication factor for the saltwater-freshwater density ratio, which is expressed as an empirical equation of specific physical parameters of the aquifer. The multiplication factor has the following form

where ϵ^* is the modified buoyancy factor or saltwater-freshwater density ratio and a_T is transverse dispersivity. This specific correction factor has been widely used in a number of studies (e.g., Lu et al., 2012; Christelis and Mantoglou, 2013; Lu and Werner, 2013; Lu and Luo, 2014; Koussis et al., 2015). Lu and Werner (2013) proposed a modification of the (Pool and Carrera, 2011) correction factor. Specifically, they performed a significant number of variable density simulations and suggested that the 1/6 exponent of Equation (8) should be replaced with 1/4. In a more recent study, Christelis and Mantoglou (2016a) observed that in highly pumped aquifers, the value of the density ratio that allows the SI to resemble the salinity spatial distribution depends not only on the aquifer physical parameters included in Equation (8), but also on the well pumping rates. In order to incorporate this effect, they suggested a dynamic adjustment of the density ratio during the pumping optimization process. In this case, the correction is based on the modification of the density ratio in an intermediate optimization step, which minimizes the discrepancy between the VD and SI models. Kopsiaftis et al. (2019a) compared the SWI extent based on five models: (i) the VD model, (ii) the Strack SI model, (iii) the Pool and Carrera (2011) correction, (iv) the Lu and Werner (2013), and (v) several combinations of the above SI models, in the form of ensemble SI model prediction. In the same study, the authors used all the examined models in a pumping optimization model and concluded that the ensemble solutions provide promising results. The ensemble SI prediction could be expressed as follows

where Ŷ_ens is the ensemble prediction of the toe, P is the number of the individual SI models, w_i is the weight corresponding to Y_i toe prediction of the i-th SI model. In the present work, three SI models are examined for the ensemble prediction of the toe: (1) the original Strack SI model Strack (1976), (2) the modified Strack model proposed by Pool and Carrera (2011), and (3) the modified Strack model proposed by Lu and Werner (2013). The contribution of the three SI models is considered equal, thus the weights are w₁ = w₂ = w₃ = 1/3 and the Ŷ_ens is calculated as the average value of the individual SI model predictions.

In this work an implicit correction of the Strack (1976) SI model is investigated using the position of a specific isohaline, which represents the extent of SWI based on the variable density approach. In the multi-fidelity (MF) related literature, the proposed correction is denoted as additive correction (e.g., Alexandrov and Lewis, 2001; Zhou et al., 2016; Giselle Fernandez-Godino et al., 2019) and, in general, it could be described as an additive shift of the low-fidelity response y_LF(x) to yield the high-fidelity response y_HF(x). The correction has the following general mathematical form Giselle Fernandez-Godino et al. (2019).

where ŷ_MF(x) is the MF surrogate model, δ^(x) is the surrogate constructed to model the discrepancy between the y_HF(x) and y_LF(x) at data point where we have model computations for both the HF and LF models. For the groundwater application examined in the present work, the input vector x corresponds to the vector Q of the pumping rates.

In this study, RF method was utilized as the discrepancy function δ^(x). RF algorithm is an ensemble learning technique for solving both classification and regression problems. RFs can be used to predict continuous target variable based on a set of input features. Each decision tree in the forest is built using a random subset of the available features, and the final output is obtained by averaging the predictions of all the trees. This results in a more stable and accurate model than using a single decision tree, as it reduces the risk of overfitting and the impact of noisy or irrelevant features. Ho (1998) was the first to embed a random subspace method in decision forest. Breiman (2001) further extended the algorithm and introduced the term “Random Forests”. RFs have been widely used in several components of the water resources field under different context. For example, Papacharalampous et al. (2018, 2019) aimed at forecasting several hydrological processes and parameters (e.g., temperature and precipitation time-sires) and they compared the outcomes with stochastic methods. A popular implementation of RF refers to the predictive modeling of several groundwater pollution sources, which is crucial in areas with significant agricultural or urban growth, especially in arid and semi-arid climates (e.g., Wheeler et al., 2015; Canion et al., 2019; Knoll et al., 2019; Messier et al., 2019; Ouedraogo et al., 2019; Lahjouj et al., 2020; Pham et al., 2021; He et al., 2022). Other noteworthy applications of RF are the prediction of groundwater levels (e.g., Wang et al., 2018; Saha et al., 2022), assessment of drought impact on groundwater potential (Masroor et al., 2021), classification of surface-groundwater interaction (Yang et al., 2019), identification of potential sites for groundwater artificial recharge (Norouzi and Shahmohammadi-Kalalagh, 2019; Naghibi et al., 2020), investigation of the spatial and temporal variations of soil salinity (Fathizad et al., 2020) and modeling of the shallow water table at high spatial resolution to calculate the risk for groundwater-induced flood events (Koch et al., 2019).

It should be noted that in the present work, in order to capture the complex relation between the well pumping rates and the proposed additive correction factor, the 0.1 kg /m³ isohaline and the corresponding SI toe are divided in several segments and a RF model is trained for the edge points of each segment. Thus, the proposed correction demonstrated good adaptive capabilities, which are further analyzed in Section 4.3.

The SWI models described in Sections 2.1–2.3, as well as the optimization frameworks of Section 2.1 are applied on a hypothetical unconfined aquifer of a simplified rectangular shape. The horizontal dimensions of the examined aquifer are L×W = 7.000m×3.000m and the base of the aquifer is considered horizontal and is located at a depth of -25 m below sea level. Along the coastline, a hydrostatic boundary condition is set with a constant salinity concentration of 35 kg /m³. Also, a specified flux is applied along the entire inland boundary, while the remaining boundaries—lateral and aquifer bottom—are considered impermeable. The aquifer is considered to be homogeneous and anisotropic with regard to the hydraulic conductivity and is pumped through ten fully penetrating wells. Finally, a uniform recharge replenishes the aquifer along its upper boundary. It should be noted that for both the recharge and the well pumping rates we used constant average values throughout the entire 50-year simulation period. Figure 1 is a conceptual representation of the coastal aquifer model, while Table 1 includes all the model's parameters. Each VD simulation requires approximately 75 s in a i7-4770 CPU of 3.7 GHz, 16 GB RAM to complete. The duration of the corresponding SI simulation is in the order of 0.3 s. The significant difference in computational time between the two SWI models highlights the benefit of the SI approach and its usefulness in applications which require a large number of simulations, such as a simulation-optimization method in groundwater management problems.

As already mentioned in Section 4.1, there is a clear advantage to utilizing the SI model rather than the VD model for groundwater management applications, due to its computational efficiency. In the present study, we further examine the correction of the SI model discussed in Section 2.3 by searching for any potential spatial dependencies of the density ratio, which would lead not to a unique correction factor for the entire flow field—as in previous studies—but to a locally adapted modification of the SI toe. To this end, we create a large dataset of 5, 000 pumping scenarios using the Latin Hypercube Sampling (LHS) statistical method. The pumping scenarios are applied on both the VD and SI models and subsequently the locations of the 0.1 kg /m³ isohaline and the SI toe are extracted from the concentration results and the flow potential field, respectively. Figures 2A, B depict the advancement of SWI for the VD model and the four variations of the SI model which include: (i) the original Strack model, (ii) the Pool and Carrera (2011) correction, (iii) the Lu and Werner (2013) correction, and (iv) the ensemble SI model described in Section 2.3. The results indicates that for low pumping rates the model proposed by Strack (Strack, 1976) outperforms the other SI models. However, in high pumping rate scenarios—which is of practical interest for groundwater management problems—the Strack model overestimates the extend of SWI, while the Lu and Luo (2014) models seems to capture more accurately the benchmark VD solution. The ensemble model toe closer resembles the shape of the 0.1 kg /m³ isohaline, though shifted several hundred meters inland. Overall, the results demonstrate a considerable dependency of the SI models' efficiency on the total pumping rates, confirming the findings of previous studies (e.g., Christelis and Mantoglou, 2016a).

In order to define possible spatial dependencies of the correction factor both the 0.1 kg /m³ isohaline and the SI toe based on Strack's model were divided in an equal number of segments and the correction of the latter was examined at the endpoints of each segment. In the hypothetical aquifer of the present study we selected a dense segmentation of the lines to achieve a relatively comprehensive and complete analysis. In particular, following the finite different discretization of the numerical models, we divided both lines into 59 segments. For each one of the segments' endpoints and for each pumping scenario we calculated an optimal density ratio that minimizes the distance between the two points. This problem is equivalent to calculating the potential flow value at each point of the 0.1 kg /m³ isohaline through an interpolation method, and then solve Equation (7) for ϵ. Figure 3 presents the optimal values for the 5, 000 pumping scenarios, for three segments of the a the critical isohaline, representing both the northern and the southern part of the aquifer, as well as the middle area. The figure leads to the following two key conclusions

The findings and observations of Section 4.2 suggest that the correction factor proposed in earlier studies for the Strack SI model (e.g., Pool and Carrera, 2011; Lu and Werner, 2013), enhance the determination of the toe position considerably, aligning with the benchmark VD approach. However, it has certain limitations and its overall precision is contingent on the pumping regime. Furthermore, applying a uniform correction factor across the entire flow field does not have the capability to adapt to the specific conditions of distinct areas within the aquifer, which—for the case of homogeneous aquifer, like the hypothetical aquifer examined in the present study—are mostly related to the relative position of the wells and the applied pumping rates.

Following the line segmentation approach of Section 4.2, a total number of 60 points is selected from both the 0.1 kg /m³ and the SI toe based on Strack model, as depicted in Figure 4. The arrows in the figure denote the intended correction for indicative points of the line. In order to investigate a spatially variable correction factor, we utilize a ML method, and particularly the RF algorithm described in Section 2.3. The training sample consists of 4, 000 variable sets, representing the 80% of the total sample. Each variable set has 11 parameters, namely the pumping rates of the 10 extraction wells and the initial position of the examined point relative to the coastline. The output set consists of a single variable, that is the distance of the examined SI toe point from its relative point lying on the 0.1 kg /m³ isohaline. Figure 5 schematically presents the overall concept of the proposed method in contrast to a typical approach in coastal aquifer managements problems. Regarding the RF based correction of the SI, the following points should be mentioned

In our experiments, we used a k-fold cross-validation technique to mitigate possible overfitting and achieve a more robust evaluation of the RF models' generalization ability. The train dataset is divided into k = 5 subsets (folds) of equal size. The RF hyperparameters were optimized using Bayesian optimization over these 5k-fold cross-validation sets. Table 2 summarizes the optimizable variables and the corresponding value ranges. The number of splits falls within the range of 1 to 3995.

For the evaluation of the RF models performance we use as metrics the mean absolute error (MAE), the root mean squared error (RMSE), and the R-squared (R²), which are defined as follows:

where N is the number of variable sets in the train or test dataset, y_i represents the true values of the modeled parameter, ŷ_i stands for the estimated values and ȳ_i is the mean value of the parameter. These metrics are applied in both training and test sets.

In order to ensure the accuracy and robustness of the final RF models, the entire training and validation process was repeated five times for all the examined points of the SI toe. The performance metrics were computed for each one of the five output values datasets. The repetition of the process allowed the approximate quantification of the uncertainty of the proposed method. The uncertainty is provided as the standard deviation of the metrics (e.g., MAE) over the five repetitions. In particular, the calculations for the MAE are based on the following formula:

Table 3 contains the values of the statistical measurements used to evaluate the performance of the proposed method, namely the mean absolute error (MAE), the root mean square error (RMSE) and the R². For the MAE we also provide the corresponding standard deviation. It should be noted that the values listed in the table are the averages over the five repetitions of the calculation process for half of the 60 points selected from the SI toe. Given the magnitude of the examined parameter and the size of the finite difference cell, the MAE and RMSE values are considered satisfactory. At the bottom of the table, we also provide the average values of all the metrics over the entire set of points. Regarding the high values of the R², the use of the cross-validation approach and the repetition of the process minimize the possibility of overfitting and ensure the model's generalization ability. The proposed methodology significantly improves the location of the SI toe, with respect to the 0.1 kg /m³ isohaline.

In order to test the usefulness of the proposed LF model in pumping optimization problems, we apply the optimization framework described in Section 3 on the unconfined hypothetical aquifer. For consistency, we solved the optimization problem not only for the HF model, but also for the original Strack model. The parameters of the EAS optimization algorithm are listed in Table 4 and their values are set based on the work of Efstratiadis and Koutsoyiannis (2002), Tsoukalas et al. (2016) and Christelis and Mantoglou (2019). The optimization is terminated if at least one of the following criteria is met: (i) the convergence criterion ε reaches its preset value, and (ii) the number of objective function evaluations exceeds the maximum number n_max.

Table 5 contains the optimal pumping rates for the examined HF and LF models. As expected, the Strack's SI model provides a relatively low estimation of the total pumping rates, which does not exceed 75% of the VD optimal results. Using the proposed correction of the SI significantly improved the results, allowing the extraction of a freshwater amount that corresponds to 95% of the VD estimation. Furthermore, it is worth mentioning that the evaluation of the optimal pumping rates with the VD optimization framework confirms that it is a feasible solution, since none of the constraint functions was violated. Figure 6 shows the 0.1 kg /m³ isohaline corresponding to the optimal pumping rates for the two examined SI models. The RF-based correction of the Strack SI model provides a solution which approaches the benchmark VD solution, without violating the constraint set.