Particle tracking as a vulnerability assessment tool for drinking water production

1 Introduction

Drinking water contamination is a major matter of concern for worldwide public health, with worldwide economic and social effects (Daud et al., 2017; Sharma and Bhattacharya, 2017; Gwimbi et al., 2019; Turner et al., 2021). This emphasizes the need for decision-support modeling tools that could provide 1) reliable predictions of the risk associated with water contamination in prospective scenarios and 2) optimized production settings to mitigate the impact of reported contamination. To this end, models may be of great interest provided a series of conditions are satisfied (Doherty, 2015; Hermans, 2017; Doherty and Moore, 2020). First, that models should properly account for the processes governing the outputs of interest. Second, that an appropriate and reliable observational dataset constrains the parameters of importance. Third, that the model is “practical”, which implies reasonable computation time and effective tools to automatize time consuming operations (Bakker et al., 2016; White et al., 2020a). Such conditions may not be easy to fulfill, in particular when dealing with the complex process of contaminant transport in heterogeneous aquifers.

In such a context, a thorough analysis of the simplicity-complexity trade-off becomes a critical step for which guidelines are now available (Hrachowitz et al., 2014; Guthke, 2017; Schwartz et al., 2017; Hugman and Doherty, 2022). In practice, appropriate choices have to be made on 1) the parameterization of hydraulic and transport properties and 2) the simulated physio-chemical processes driving the propagation of contaminant. For the parameterization of hydraulic properties, the options range from homogeneous equivalent hydraulic conductivity zones to heterogeneous fields with discrete features (de Marsily et al., 2005; Carniato et al., 2015; Pool et al., 2015). For the selection of processes governing transport, the options range from simplified models that rely on assumptions such as steady-state and dominant processes (e.g., advection), to complex, advection-dispersion reactive transport models to provide a detailed description of the physics of the problem (Anderson et al., 2015). Several studies highlighted that the best balance in terms of robustness, efficiency and reliability may be achieved with relatively simple “surrogate” models based on simplified representations of the physics (Razavi et al., 2012; Asher et al., 2015; Burrows and Doherty, 2015). Surrogate transport models present fast run times allowing for thousands of model executions of transport calculations, which are necessary for history matching of highly parameterized models. The interest of including a diverse observational data types has been highlighted by Hunt et al. (2006) and concentration in particular has been recently discussed by Schilling et al. (2019) and Knowling et al. (2020). Furthermore, fast run time allows the use of more robust but demanding algorithms for uncertainty quantification (Rajabi et al., 2018).

During the last decades, parameter estimation and uncertainty quantification algorithms have been made available with an ever growing variety of approaches (Doherty, 2016; White et al., 2020b). Their use is now facilitated by Python interfaces (White et al., 2016). More recently, studies provided repeatable script-based workflows which greatly facilitate the replication of the presented approaches (White et al., 2020a; Fienen et al., 2022). Though facilitated, the interfacing of parameter estimation and uncertainty quantification algorithms with complex models remains difficult for transport models, which are characterized by long computation times and numerical instabilities.

Focusing on the widely reported risk of contaminant transfer from contaminated rivers to groundwater production units, we present a framework based on particle tracking as a fast and effective surrogate model for contaminant transport. The approach, initially presented by (Cousquer et al., 2018), is made available in a newly developed Python module, TrackTools which will facilitate its replication. The module provides particle seeding capabilities and post-processing options that are valuable for analyzing drinking water vulnerability of production wells or drains to pre-defined sources of contamination. The script facilitates exploratory parametric analysis, which can be useful to investigate the system response to different configurations. The interfacing with the PEST + + suite (White et al., 2020b) is detailed on a real-world case study which paves the way for history matching, hypothesis testing, and optimization of decision variables, which are essential to support decision related to the definition of wellhead protection area and the optimization of production settings.

The theoretical background and numerical tools related to the developed Python module are described in Section 2. A simple synthetic model is presented in Section 3 and a parametric study is conducted to investigate the driving factors of mixing ratios computed with this approach. In Section 4, the interest of the method is illustrated on a drinking water production site with a fully script-based approach from model setup to pumping optimization through parameter estimation.

2 Methodology

A common practice in vulnerability analysis is to investigate the origin of water withdrawn at a groundwater sink (well or drain) originating from one, or a series of potential or effective contaminant sources. This can be conducted from the analysis of the flow contributions of each source to the discharge rate of the sink. The methodology described hereafter is an extension of the approach described by Cousquer et al. (2018) initially designed for a single river reach. The method may now be applicable to multiple weak or strong source boundary conditions (e.g., fixed head, river, general head boundary condition).

Consider a groundwater flow model with a groundwater sink (well or drain) subject to a potential are effective contamination originating from one of the model boundary conditions. A set of N starting points is seeded around the model cells where the sink condition is applied. Given the flow field, the origin of flow to the sink can be described by backward particle tracking.

The sink discharge rate Q can then be decomposed as follows: $Q={\sum }_{i=1}^N{q}_i$, where q_i is the contribution of the i-th particle to the discharge rate at the sink. This can be rewritten considering particle velocities, v_i: $Q={\sum }_{i=1}^N{v}_i\cdot S_i$ where S_i is the area of the surface crossed by q_i. When particles are evenly distributed around the model cell, it can be assumed that $S_i=\frac{S}{N}$ for all i.

Assume that from the N particles seeded around the sink, n_j originate from the j − th boundary condition, subject to contamination. The contribution of this boundary condition to the sink contamination can be written as follows:

where β_k is the mixing ratio in the source model cell where the k-th particle ends. β_k = 1 for “strong” sources and 0 ≤ β_k ≤ 1 for “weak sources” where a mixing occurs (Pollock, 2016). In the latter case, β_k can be quantified from the cell water budget at the endpoint cell (Cousquer et al., 2018). This allows the method to account for mixing in the source aquifer cell, where the boundary condition (typically a river) is applied.

The approach is now integrated in the TrackTools Python module. Our approach relies on the USGS MODFLOW6 (Langevin et al., 2017) and MODPATH7 (Pollock, 2016) programs to simulate fluid flow and transport processes, respectively. We rely on the FloPy module (Bakker et al., 2016) for the processing of model input and output files.

Pre-processing capabilities of the TrackTools module are provided in the ParticleGenerator class, which can handle:

• Particle seeding around groundwater sinks, which can be described by a Python geometry or ESRI^TM shapefile,

• Adding, removing and merging particles into groups for FloPy and MODPATH7.

After the simulation of flow with MODFLOW and particle tracking with MODPATH, the TrackingAnalyzer class of TrackTools can be used to process particle pathline data and derive the values of mixing ratios at each groundwater sink. The MODFLOW cell-by-cell water budget file is also required to derive the source cell mixing ratio β_k. Results can be provided in the form of a data frame and plotting options are proposed for the description of the origin of groundwater withdrawn in the series of sinks where particles were seeded.

3 Synthetic case

3.1 Model description

In order to illustrate the presented Python module, a synthetic 2D case is considered with a production well in an unconfined aquifer in interaction with a stream. The domain (3,750 × 5,000 m) is discretized with a 2D regular mesh (25 × 25 m cells) with a local refinement around the stream and the well (6.5 × 6.5 m cells). A Cauchy-type boundary condition is applied to the upper limit of the model domain with a stage of 40 m and conductance of 10^–3 m² s⁻¹, while a Dirichlet-type condition is prescribed to the lower boundary with a stage of 20 m. They are simulated with the General Head Boundary (GHB) and Fixed Head (FH) MODFLOW packages, respectively. The left and right domain boundaries are considered as impermeable. The stream is simulated with a head-dependent flux (Cauchy-type) boundary condition with the dedicated MODFLOW river package (RIV), with a head ranging from 35 m to 25 m from the upper to the lower boundary condition. The stream conductance is set to 10^–3 m² s⁻¹. Steady-state flow conditions are considered and the transmissivity is assumed to be independent of water level fluctuations, which corresponds to the Boussinesq assumption. The well is considered as fully penetrating and is located at position (x, y) = (1212, 1363).

We consider both homogeneous and heterogeneous hydraulic conductivity fields, the heterogeneous cases being generated with an isotropic exponential variogram for log 10(k) defined by a sill of one and a nugget of 0.1 for two different ranges of 200 and 500 m as presented in Figure 1. These fields are generated using the geostatistical python package GSTools (Müller et al., 2021). During the following simulations, the mean hydraulic conductivity value will evolve in heterogeneity patterns (Figures 1B,C).

FIGURE 1. (A) Synthetic case model structure and boundary conditions with a General Head Boundary Condition (GHB) of 40 m to the upper boundary and a Fixed Head (FH) of 20 m to the lower boundary. No flow is imposed at the left and right domain boundaries. The stream is simulated with a head-dependent flux (Cauchy-type) boundary condition and the well position is represented with a red dot. The hydraulic conductivity field patterns correspond to an exponential viariogram with a sill of 1, a nugget of 0.1 and a range of (B) 200 m and (C) 500 m.

3.2 Parametric study

The synthetic case previously described is used to illustrate the behavior of the mixing ratio α depending on key parameters such as the hydraulic conductivity K and pumping discharge Q. β from Eq. (1) is quantified from the river cell water budget and can vary from 0 to 1. Figure 2 shows the distribution of α for parameters K and Q ranging from 10^–5 to 10^–3 m s⁻¹ and from 0 to 200 m³ h⁻¹, respectively, considering a homogeneous distribution of K. Boundary conditions remain the same throughout this parametric study. We also present in Figure 3 three examples of the spatial distribution of the hydraulic head H with K set to (A) 10^–5, (B) 10^–4 and (C) 10^–3 m s⁻¹ and the same configuration but with a pumping discharge of 90 m³ h⁻¹ for D, E and F. These values have been chosen as representative of the three main tendencies that are observed for α: 1) the yellow region in which example D is located corresponding to the highest values of α, 2) the dark blue region in which example E is located corresponding to the lowest values of α and, 3) the green region in which example F is located corresponding to intermediate values of α. This shows that increasing K for a given value of Q results in decreasing (from D to E) and then increasing (from E to F) α, except for very small values of Q for which α directly increases (from E to F). These different behaviors are explained by the different sources of water that contribute to the pumped water, which are characterized by the particle paths represented on the hydraulic head distributions in blue. For small values of K (example D), most of the particles comes from the portion of the river that is located right next to the well, implying that the river provides a substantial proportion of water to the pumping well. Increasing K (example E) results in reducing the impact of the pumping on the natural hydraulic head distribution, implying that the top boundary condition contributes to the pumping. When we keep increasing K (example F), we reach configurations that do not depend on the pumping rate (green triangular region on the right side of Figure 2). For these cases, the pumping does not modify the natural hydraulic head distribution and the particle paths are fully driven by the boundary conditions, resulting in particles coming from the top of the river and thus high values of α. For small values of Q and K (black region of the bottom of the figure), the impact of the pumping on the natural hydraulic head distribution is still negligible but the values of K imply that there is no contribution of the river to the pumping (α = 0).

FIGURE 2. Distribution of the mixing ratio α depending on the hydraulic conductivity K and pumping discharge Q.

FIGURE 3. The spatial distribution of the hydraulic head H(x, y) for hydraulic conductivity K set to 10^–5 m s⁻¹ (A,D), 10^–4 m s⁻¹ (B,E) and 10^–3 m s⁻¹ (C,F). Cases (A,B,C) are without pumping, while a pumping rate of 90 m³ h⁻¹ is set for cases (D,E,F). The red dot represent the pumping well locations, the vertical bold black line represents the stream and the blue lines in cases (D,E,F) are the path of the particles from the backward particle tracking.

Figure 4 shows the distribution of α when considering the heterogeneous hydraulic conductivity fields provided in Figure 1. The three main tendencies described above for the homogeneous case are also observed for the heterogeneous cases, showing that these behaviors are mostly driven by the boundary, river and pumping conditions. However, small differences are noticeable, in particular for small values of Q and large values of $\bar{K}$, which correspond to configurations with a small impact of the pumping on the natural hydraulic head distribution. For example, the extent of the dark blue region observed at the bottom of Figure 3 is reduced and increased in Figures 4A,B, respectively, showing that the river does not contribute to the pumping for different values of Q and $\bar{K}$. We also observe the presence of yellow areas for large values of $\bar{K}$ in the heterogeneous cases, corresponding to high values of α and thus a strong contribution of the river to the pumping. These different behaviors are due to different flow paths induced by the heterogeneities of the hydraulic conductivity fields. Note that, as mentioned before, the impact of these heterogeneities is reduced when increasing the pumping rate.

FIGURE 4. Distribution of the mixing ratio α depending on the average hydraulic conductivity $\bar{K}$ and pumping discharge Q for the hydraulic conductivity patterns presented in (A) Figures 1B and (B) Figure 1C.

3.3 Derivative analysis

This synthetic case is also used to set some rules in order to obtain robust derivatives which is essential for parameter estimation and model analysis based on linearization methods. The sufficient number of particles corresponding to this synthetic case is evaluated with the derivative calculation of $\frac{\partial \alpha }{\partial K}$ which constitutes the Jacobian matrix used for model linearization. The use of gradient based methods for parameter estimation (e.g., the Gauss-Marquardt-Levenberg method on which PEST and PEST++ are based (Doherty, 2015; White et al., 2018)) requires that model outputs are differentiable or at least continuous with regard to model parameters (Doherty, 2010). Doherty and Hunt (2010) mentioned that this assumption often does not hold because of poor model performance and the lack of differentiability can be due to, among others reasons, a low number of particles in the MODPATH particle tracking model. In order to avoid this issue, we explore the differentiability of model output regarding particles number. Let the mixing ratio α be the output of interest, the hydraulic conductivity K is the parameter and ΔK is the parameter increment. The sensitivity of α to K is approximated as follows:

The estimate of $\frac{\partial \alpha }{\partial K}$ may be biased when 1) ΔK is too small, due to numerical noise, 2) when ΔK is too large due to model non-linearity, or 3) when the number of particles is not high enough to correctly simulate the mixing ratio. These three points are evaluated by calculating the value of the derivative for several numbers of particles and different values of K_i considering two values of ΔK set to 0.01 × K_i and 0.1 × K_i (Figure 5).

FIGURE 5. Exploration of derivatives quality behaviour regarding the number of particles with the hydraulic conductivity K_i ranging from 10^–5 to 10^–3 m s⁻¹ and the increment ΔK set to (A) 0.01 × K_i and (B) 0.1 × K_i.

For both ΔK = 0.01 × K_i and ΔK = 0.1 × K_i, Figure 5 shows changes in the derivative sign that are visible when the color representing the derivative goes from blue (negative values) to red (positive values) and vice versa. For a given line, these changes depend on the number of particles N_p used in the model, showing wrong estimations of the derivative when N_p is too small (due to wrong estimations of the mixing ratio). This results in instabilities of the derivative along K for small number of particles (N_p < 200−500, columns that are located on the left side of Figure 5). For parameter estimation and inversion purposes, these instabilities could be interpreted as local minima and lead to wrong estimates and optimized values. When N_p is high enough (columns on the right side), the derivative values are stable in the sense that they do not change with the value of N_p. In this case, changes in $\frac{\partial \alpha }{\partial K}$ along K are due to the non-linearity of the model, which needs to be taken into account for parameter estimation and optimization.

Furthermore, comparing Figures 5A,B shows that the choice of the increment of K influences the stability of the derivative. The number of particles that is required to observe stable values of $\frac{\partial \alpha }{\partial K}$ is lower when increasing ΔK.

Although they are expected, these results clearly show that the derivative calculation requires to carefully choose the parameters N_p and ΔK with 1) a sufficient number of particles to correctly simulate the mixing ratios and its derivative, and 2) an increment that is neither too small to avoid that the derivative calculation is tainted by numerical problems, nor too large to take into account the non-linearity of the model. In all cases, increasing the number of particles is a systematic method to improve the definition of the derivative whatever the value of the discretization step of the derivative.

4 Case study

4.1 Model description

The method described above is embedded in a workflow designed for model-based decision making and applied to a complex real world case study. We describe the script-based modeling workflow and interfacing with the PEST + +, which provides an illustration of the interest of the method and facilitates the replication of the approach to other cases. All scripts are available on a GitHub^® repository.

The well field is located in South-West France and is critical to the water supply for the city of Bordeaux (Cousquer, 2017; Valois et al., 2018; Delbart et al., 2021) as it accounts for about 20% of the needs. The aquifer lies in Oligocene limestone overlain by sandy and alluvial deposits. Limestones are locally subject to karstification, leading to strongly heterogeneous hydraulic properties. Two wells (W1, W2) and two horizontal drains (D1, D2) are used for groundwater extraction. Due to industrial activities, the main stream crossing the well field from west to east is prone to contamination. The vicinity of pumping wells and drains to the stream is favorable to stream-aquifer exchanges. The end-purpose of this modeling exercise is to maximize groundwater production while meeting drinking water quality standards.

For this model purpose, the geographical data describing the domain and geometry of boundary conditions was first processed with QGIS (QGIS Development Team, 2022). The model setup was fully script-based with the FloPy Python library (Bakker et al., 2016) for the pre- and post-processing of MODFLOW6 and MODPATH7 input and output files. A 2D single layer was considered under the Dupuit-Forchheimer approximation to represent the aquifer of interest. The model domain was discretized by Gridgen (Lien et al., 2015) with a quadtree grid with four horizontal refinement levels so that square cell dimensions range from 200 m to 12.5 m (Figure 6). External boundaries have been defined from the regional groundwater levels as no-flow or 3^rd type head-dependent flow boundary condition. The streams and drains are simulated with head-dependent (Cauchy-type) flux boundary conditions and the wells are represented by sink terms in the aquifer cells corresponding to their location.

FIGURE 6. Model domain with boundary conditions inferred from the regional groundwater flows. Mesh refinement is conducted in the area of interest around the production and observation wells and drains (see inset). The parameterization of hydraulic conductivity has been conducted with pilot points.

In order to consider contrasting settings while avoiding the burden of transient simulations, pseudo steady-state flow conditions are considered (Haitjema, 2006; Moore and Doherty, 2021). This hypothesis is supported in the present case since the permeable aquifer responds quickly to changes Cousquer (2017). The TrackTools module is used to seed particles around the two wells and the two drains and to derive the mixing ratio of the water withdrawn. It was found that simulated mixing ratios were relatively stable with 500 particles seeded around each production unit (wells and drains). As detailed in Section 2 and references herein, the mixing ratios at source cells where particles stop (β) is inferred from the cell water budget. This allows to account for mixing of aquifer and river water in model cells where a river boundary condition is applied. Note that only the main river is prone to contamination, water originating from the tributaries is not considered as contaminated.

4.2 Parameter estimation and optimization

The model interfacing with the PEST + + suite has been performed by means of PyEMU (White et al., 2016). The hydraulic conductivity field was parameterized with a dense set of pilot points. Hydraulic conductivity values at pilot points lying at the center of the study area were all considered as independent, but they were tied in the outer portion of model domain (Figure 6). This is motivated by the lack of observations in the outer zone and it is unlikely that model predictions are sensitive to details in this remote area.

A series of 10 field surveys, spread over 4 years (2014–2018) have been considered for history matching, considering pseudo steady-state conditions for contrasting stream levels and operating conditions. Though it was challenging for a well field in activity, we have made our best for operating conditions (well discharge rates and drain levels) to remain relatively stable on the period preceding each of the surveys. The observation data set is composed of hydraulic heads, drain discharge rates, and mixing ratios. The latter were derived by end-member analysis from HCO³⁻ and Ca²⁺ concentrations (Delbart et al., 2021).

In order to reduce run times, all model files for each of these surveys were setup in separate folders (Figure 7). Doing so, file operations during parameter estimation were limited to writing parameter files and reading output files. These operations were all conducted with the dedicated methods of the PyEMU Python package.

FIGURE 7. Directory tree structure used for parameter estimation and optimization workflow.

The parameter estimation was conducted with the widely used Gauss-Levenberg-Marquardt Algorithm (GLMA) as implemented in PEST + + suite (White et al., 2020b). Parameter increments (DERINC) values were adjusted by trial and error. Best results were obtained with relative parameter increments of 15% for hydraulic conductivities and 10% for all the other parameters. Both zero order (preferred value) and first order Tikhonov regularizations were employed (Doherty, 2015).

After parameter estimation, the resulting model was used to optimize the total abstraction rate. The optimization problem consists in maximizing the total abstraction rate Q while verifying the drinking water quality standards after mixing. It can be expressed as follows:

where the decision variables Q_i and h_i correspond to well discharge rate and drain level, respectively, and the constraint on water quality is expressed as a critical stream-aquifer mixing ratio α_crit.

For this illustrative exercise, the optimization is solved with a sequential version of linear programming as implemented in PESTPP-OPT (White et al., 2018) considering a 50% risk (maximum likelihood) configuration. This algorithm is fast but sensitive to model non-linearity.

4.3 Results

The history matching conducted with the GLMA converged in approximately 10 iterations (Figure 8A). The fit between simulated values with their observed counterparts can be considered as satisfying for heads (Figure 8B) and reasonable for discharge rates and mixing ratios (Figures 8C,D). Measurements of drain discharge rates can be considered as reliable and the misfit may rather be explained by the static, linear relation that is considered with the MODFLOW drain package. In contrast, measurement of mixing ratios derived from concentrations by end-member analysis are prone to uncertainties (Delbart et al., 2021). The structural error of this simplified model is therefore strong and surely contributes to the misfit, but uncertainties in the observation data and historic operating variables are also important. In this context, it is unclear whether more detailed process-modeling could improve the predictive capacity of the model.

FIGURE 8. Results of the parameter estimation with the GLMA. (A) Evolution of the measurement objective function throughout GLMA iterations. Simulated versus measured hydraulic heads (B), drain discharges (C), and mixing ratios (D). Numbers in (B,C,D) refer to the surveys; colors in (B) refer to the observation wells.

The estimated hydraulic conductivity field is highly heterogeneous (Figure 9A), as expected for partly karstified limestones. Attempts to smooth these contrasts with stronger regularization constraints lead to an important increase of the misfit with historical data. Karst conduits are expected to be narrow and of limited extension but can be at the origin of the high values of estimated hydraulic conductivity (Cousquer, 2017). On the opposite, low values can be explained by fine deposits in the riverbed (Delbart et al., 2021).

FIGURE 9. (A) Hydraulic conductivity field after parameter estimation. (B) Simulated groundwater table level and particle tracking trajectory for the initial values of parameters Q and α. Particle trajectory from the contaminated stream are colored in red and those from uncontaminated boundaries in blue. Particle tracking was conducted from production wells W1, W2 and drains D1, D2 shown in black. Pumping from wells Wx1 and Wx2 was considered in the flow model but they do not pertain to the studied drinking water facility.

We followed a deterministic approach, which can be considered as a first exploratory step for more robust methods. The “best” calibrated parameter set was therefore used to conduct the optimization of well and drain discharge given a quality constraint expressed as a global mixing ratio α_crit=25%. This corresponds to a 4-fold reduction of contaminant concentration with respect to river contamination levels. The optimization algorithm was run from a situation close to the current operating configuration over three iterations of the sequential linear programming algorithm (Figure 10). Compared to the initial settings (iteration 0), the third and last iteration lead to an increase of production rate by + 25% for a similar contamination level. The particle pathlines for this configuration are presented in Figure 9B.

FIGURE 10. Optimization results with (A) total discharge rate and mixing ratio for each optimization iteration of PESTPP-OPT and (B) discharge rate (Q), mixing ratio (MR), and drain stage (H) for each optimization iteration. it0 represents the initial parameter values, it1, it2 and it3 the values for the first, second and third iterations, respectively.

5 Discussion and conclusion

In the context of drinking water contamination issues, we provided a decision-making tool to investigate the vulnerability of groundwater production units. The TrackTools Python module offers pre- and post-processing functions of particle tracking data simulated with MODFLOW and MODPATH. Mixing ratios are inferred from flow contributions of each contaminant source to the discharge rate of a groundwater sink considering particle velocities and mixing in source cells. This can be appropriate where a production well or drain is vulnerable to a boundary condition such as a river. However, the method is not relevant when the source of the contamination is diffuse, such as agricultural contaminants driven by groundwater recharge. In the latter case, forward particle tracking is more appropriate [see e.g., Fienen et al. (2022)].

The synthetic case presented in this work illustrates that contrasting behaviors can be observed even for simple configurations. It highlights the interest of the method to explore the sensitivity of mixing ratios to operating variables (pumping rates, drain levels) and model parameters. The impact of the number of particles seeded around each sink, and the size of the increment used for model linearization was also evaluated through the analysis of the derivatives of α. Derivative instabilities are mitigated when increasing the number of particles or increasing the parameter increment for estimating model derivatives (DERINC). However, adding particles tends to increase the computation time, so that a reasonable particle number should be considered.

Relatively fast and easy to implement, this approach will facilitate the use of contaminant transport models for history matching, uncertainty quantification or optimization. It also presents the advantage to be didactic and can be used to represent contaminant transfer in a visual and potentially interactive manner, which is important for end-users.

The method was implemented on a real world case study and is provided with a series a script for the interfacing with the PEST + + suite. The parameter estimation with the GLMA and optimization by linear programming illustrated the interest of the approach. Difficulties in reproducing observations can be explained by errors in the historic dataset (observations and operating variables) and model structural error. In spite of our efforts to obtain a robust linear version of the model, the methods based on the Jacobian matrix have shown their limitations. This advocates for the use of more robust ensemble methods such as the Iterative Ensemble Smoother (IES (White et al., 2018)) for parameter estimation and uncertainty quantification, and evolutionary algorithms for optimization [PESTPP-MOU (White et al., 2022)]. The presented workflow may easily be extended to these methods recently available in the PEST + + Suite.

The method was illustrated on 2D horizontal models assuming steady state flow conditions. This simplifies the implementation and leads to particularly fast run times, but 3D models and transient conditions may have to be considered in other contexts. The method may be extended to these configurations with some additional processing. The implementation of the method on 3D models will be straightforward so long the wells and drains penetrate a single model layer. Otherwise, particles should be seeded in each of the intersected layers and mixing in the multi-layer well should be accounted considering the respective flow contributions provided by the dedicated Multi-Aquifer Well package (MAW) from MODFLOW6 (Langevin et al., 2017). The impact of transient flow dynamics on the mixing ratio at a specific instant will be easy to consider with backward particle tracking over historic conditions. However, it would be more challenging to quantify mixing ratios averaged over a period as it would require multiple runs of the particle tracking algorithm.

Among the simplifying assumptions of this approach, we assumed the respective flow contributions of particles evenly placed around the sink to be proportional to their velocities (Eq. 1). The validity of this approach has been investigated by Cousquer et al. (2018) with an advective-dispersive model on a synthetic, homogeneous model. Results were encouraging but the performance of the simplified model based on particle tracking may be put into question in other contexts. The use of a paired simple-complex model approach (Doherty and Christensen, 2011) can be suggested, where the complex model would be based on the resolution of the advection-dispersion equation.

Future work will also focus on improving the reliability of the simplified transport models by considering additional processes, such as dispersion. For this matter, various recent particle-based methods (Noetinger et al., 2016; Gouze et al., 2020; Roubinet et al., 2022) and machine learning techniques (Kang et al., 2021; Zhou et al., 2021) could be considered in order to consider more realistic configurations while keeping the low computational cost of the forward models.

As mentioned in the introduction, the optimal complexity level is to great extent, case- and purpose-specific, as it will depend on the purpose of the modeling exercise and the quality of the available observation dataset. In this perspective, it is likely that simple and effective methods such as the one presented in this study remain of interest in numerous situations.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://github.com/tracktools/.

Author contributions

AP: methodology, software development, case study, writing and editing. PM: software development, data processing and visualization. YC: methodology, synthetic case, writing DR: synthetic case, writing, reviewing, editing.

Funding

This work was initially conducted in the framework of the Mhyqadeau project supported by Suez Environnement (LyRE) and the French Aquitaine regional council.

Acknowledgments

The authors are thankful to the editors for handling this special issue and to the two anonymous reviewers for their corrections and valuable suggestions.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References