Identification of groundwater mean transit times of precipitation and riverbank infiltration by two‐component lumped parameter models

Groundwater transit time is an essential hydrologic metric for groundwater resources management. However, especially in tropical environments, studies on the transit time distribution (TTD) of groundwater infiltration and its corresponding mean transit time (mTT) have been extremely limited due to data sparsity. In this study, we primarily use stable isotopes to examine the TTDs and their mTTs of both vertical and horizontal infiltration at a riverbank infiltration area in the Vietnamese Mekong Delta (VMD), representative of the tropical climate in Asian monsoon regions.

unique data set of stable isotope records for a tropical region. We quantified the contribution that the two sources contributed to the local shallow groundwater by a novel concept of two-component lumped parameter models (LPMs) that are solved using δ 18 O records.
The study illustrates that two-component LPMs, in conjunction with hydrological and isotopic measurements, are able to identify subsurface flow conditions and water mixing at riverbank infiltration systems. However, the predictive skill and the reliability of the models decrease for locations farther from the river, where recharge by precipitation dominates, and a low-permeable aquitard layer above the highly permeable aquifer is present. This specific setting impairs the identifiability of model parameters. For river infiltration, short mTTs (<40 weeks) were determined for sites closer to the river (<200 m), whereas for the precipitation infiltration, the mTTs were longer (>80 weeks) and independent of the distance to the river.
The results not only enhance the understanding of the groundwater recharge dynamics in the VMD but also suggest that the highly complex mechanisms of surface-groundwater interaction can be conceptualized by exploiting twocomponent LPMs in general. The model concept could thus be a powerful tool for better understanding both the hydrological functioning of mixing processes and the movement of different water components in riverbank infiltration systems. time (mTT) describes the average time that water particles spend travelling through a system (K. J. McGuire & McDonnell, 2006), and the transit time distribution (TTD) describes the whole spectrum of transit times of those water particles transported through the system (Maloszewski, 2000). Knowledge of the TTD and its corresponding mTT is essential for groundwater resources management, for example, when installing groundwater extraction systems for water supply (Hiscock & Grischek, 2002), calibrating numerical groundwater transport models (Bethke & Johnson, 2008), evaluating the security of drinking water supplies (e.g., Darling, Morris, Stuart, & Gooddy, 2005;Eberts, Böhlke, Kauffman, & Jurgens, 2012), or understanding the sources of contamination (e.g., Morgenstern et al., 2015).
Most mTT estimation methods are based on the lumped parameter model (LPM), pioneered by Małoszewski and Zuber (1982). The LPM does not require detailed hydrological information and thus can be used for initial investigations of little known systems (Mook & Rozanski, 2000) where data are insufficient, for example, in developing countries or ungauged basins (K. J. McGuire & McDonnell, 2006).
It is based on a lumped convolution integral (a black box model), in which an input signal is related to a specific transfer function (or a TTD) to obtain an appropriate output signal. A common assumption for the application of LPM is the time invariance of TTDs (see Maloszewski & Zuber, 1996;K. J. McGuire & McDonnell, 2006).
Groundwater modelling in the Vietnamese Mekong Delta (VMD) was pioneered by Haskoning B.V., DWRPIS (Boehmer, 2000) who set up a regional groundwater model. The author pointed out that groundwater recharge in most of the delta ranges from 0.01 to 1 mm/day and is dominated by (a) infiltration of precipitation and irrigation water, (b) downward leakage through the semipermeable layers of the Holocene aquifer, and (c) seepage from rivers, streams, and lakes. The water balance analysis suggested that recharge from rainfall and irrigation is significant (for the Plain of Reeds, see Figure 1) and slightly (for the whole VMD) smaller than that from the Mekong river branches and the canal system. Also, the hydraulic connection between shallow and deep groundwater is insignificant, except for the dune area along the coast of the eastern VMD (Boehmer, 2000).
In this study, we primarily used stable isotope (δ 18 O) time series to identify the mTTs of shallow groundwater and the optimized TTDs best describing the subsurface flow conditions when applying the two-component LPMs. The Plain of Reeds, serving as the seepage area for groundwater infiltration in the north-west region of the VMD (Boehmer, 2000), was chosen as a pilot area. For the sampling campaigns, the study site An Long in the Plain of Reeds (Figure 1) was selected, primarily considering logistic constraints. Although this area may not be representative for the entire Mekong Delta, the site is likely indicative of the general nature of the near-stream subsurface flow dynamics because of the low variation of lithology and topography (Nguyen, Ta, & Tateishi, 2000) throughout the region. The information on subsurface mixing processes, the preferred flow pathways, and the transit times of water infiltration at riverbank areas resulting from long-term isotopic records could enhance the understanding of the groundwater dynamics as well as groundwater vulnerabilities in the VMD.
The general objective of this work was to test the applicability of two-component LPMs to examine TTDs and their mTTs in a riverbank infiltration system where two distinct water components are present.
Our specific objectives were to (a) identify the dominant TTDs that best describe the subsurface flow conditions, (b) quantify the subsurface mixing processes and the time-variant mTTs of water near the riverbank, and (c) determine the uncertainty associated with parameter identification and model performance.

| STUDY AREA
The study area is located in the Plain of Reeds, in the VMD between latitudes 10 42 0 7 00 N to 10 48 0 9 00 N and longitudes 105 22 0 45 00 E to 105 33 0 54 00 E (Figure 1). The average elevation ranges from 1 to 4 m above sea level. The average annual rainfall is 1,400-2,200 mm, characterized by a distinct seasonal distribution (GSO, 2016;Renaud & Kuenzer, 2012). The annual average temperature is 27 C with monthly averages ranging from 25 C to 29 C. The annual average relative humidity ranges from 77% to 88%. The monthly evaporation F I G U R E 1 The Plain of Reeds in the Vietnamese Mekong Delta (right) and the study site (zoomed in) of An Long (left). The screening depths of Wells A, B, and C are 15, 12, and 14 m, respectively. The distance from Wells A, B, and C to the Mekong river is 140, 190, and 660 m, respectively. The pond with an area of approximately 500 m 2 and a depth of 2 m is used for fish farming. The distances from the pond to Well C and the Mekong river are around 40 and 700 m, respectively rate ranges from 67 to 80 mm and 76 to 109 mm in the rainy and dry season, respectively (GSO, 2016;Renaud & Kuenzer, 2012).
The hydrogeological units in the VMD are classified according to their geological formation: the Holocene, the Pleistocene, the Pliocene, and the Miocene aquifer systems. At the study site, the target aquifer is the Holocene sediment sequence. It is characterized by the uppermost layer of silt and clay (low-permeable aquitard), overlying a layer of fine to coarse sands (high-permeable aquifer). The average depths of aquitard and aquifer are approximately less than 11 m and 30 m below ground level, respectively (Minderhoud et al., 2017;Wagner et al., 2012). The grain sizes defined for clay, silt, fine and coarse sand, and fine gravel are 1-4 μm, 4-63 μm, 63 μm-2 mm, and 2-8 mm, respectively (Wentworth, 1922). The hydraulic conductivities of aquitard and aquifer are 0.02-0.2 m/day (for silt and clay) and 12-200 m/day (for fine to coarse sands), respectively (Boehmer, 2000). The effective porosities of clay and sand layers are 0.5 and 0.2, respectively .

| Water sampling and isotopic analysis
Precipitation and river water were sampled at An Long station.
Groundwater was sampled at three wells (A, B, C) closed to the An Long station and the Mekong river ( Figure 1). The distances from Wells A, B, and C to the Mekong river are 140, 190, and 660 m, respectively. The screening depths of these wells are 15, 12, and 14 m, respectively. These wells are used for household water supply only. Following the classification of aquifer systems by Wagner et al. (2012), the groundwater samples were collected from the Holocene aquifer and representative of the shallow groundwater in the VMD. A pond, located 700 m away from the Mekong river and used for fish farming, was included to provide information about the isotopic fractionation of local surface water by evaporation, characteristic for the floodplains during the monsoonal floods. The total number of samples and the schedule of water sampling are summarized in Table 1.
To avoid evaporation effects, we stored the collected samples in 30-ml plastic sample bottles with tight screw caps and kept all samples in the dark before the laboratory analysis. The stable isotope samples were analysed at the Alfred-Wegener-Institute (AWI) in Potsdam, Germany. The measurements were performed with a Finnigan MAT Delta-S mass spectrometer using equilibration techniques to determine the ratio of stable oxygen ( 18 O/ 16 O) and hydrogen ( 2 H/ 1 H) isotopes. Analytical results were reported as δ 2 H and δ 18 O (‰, relative to Vienna Standard Mean Ocean Water) with internal 1σ errors of less than 0.8‰ for δ 2 H and 0.1‰ for δ 18 O. The detailed measuring procedure is described in Meyer, Schönicke, Wand, Hubberten, and Friedrichsen (2000).

| Hydrological measurements
Groundwater and river water levels were recorded every 15 min between June 2015 and July 2017 by pressure sensors (HOBO U20 Fresh Water Level Data Logger). River water levels were monitored at An Long station, about 2 km upstream of the wells. Groundwater levels were observed at two additional monitoring wells screened at a depth of 15 m below ground level, in order to avoid disturbance of the level records by water extraction. The first monitoring well was installed between Wells A and B, and another one was located between Well C and the pond (Figure 1). The distance from the first and second monitoring wells to Wells A and C are 20 and 25 m, respectively. Sediment samples taken during the installation of these monitoring wells indicated that the upper aquitard layer (from 8 to 10 m below ground level) is dominated by clay and silty clay, whereas the aquifer layer below consists mainly of coarse sand. A terrestrial survey was carried out in June 2016 to reference all recorded water level measurements to the gauge at An Long, a national water level monitoring station. All water levels are reported as meters above sea level.

| Two-component LPMs
LPMs are based on the lumped convolution integral approach (Małoszewski & Zuber, 1982) to transform the tracer input signal (C in ) into the tracer output signal (C out ), considering a distribution of transit times according to a transfer function. M. K. Stewart and McDonnell (1991) introduced a more robust approximation by adding flow weights (w) to the isotopic composition of the input so that the outflow composition reflects the mass flux of tracer leaving the system where τ is the transit time, t is the time of exit from the system, and (t − τ) represents the time of entry into the system; C in and C out are the input and output tracer signature, respectively; and g(τ) is the transfer function representing the assumed TTD of the  Maloszewski, & Zuber, 1984;Maloszewski & Zuber, 1996) as the weighting term w(t) can include any appropriate factor such as rainfall rates, throughfall rates, or effective rainfall (K. J. McGuire & McDonnell, 2006).
The recharge sources of shallow groundwater in the VMD are mainly river water and precipitation (Boehmer, 2000;Ho et al., 1991;Wagner et al., 2012). The output isotopic composition (δ 18 O) thus stems from the water of two different sources with likely different transit times. Therefore, the two-component LPMs were chosen for transit time modelling. We excluded a deep groundwater component because the connection between shallow and deep groundwater is insignificant Boehmer, 2000;Ho et al., 1991). During the calibration, integrated parameters were adjusted to fit the measured isotopic records for each investigated well, following Weiler, McGlynn, McGuire, and McDonnell (2003).
The model was rewritten in a modified two-component convolution equation: where C inR and C inP are the input tracer signature of the river and precipitation infiltration, respectively; C Well is the output tracer signature of an investigated well; g R (τ)and g P (τ) are TTD functions of the river and precipitation infiltration, respectively; p and (1 − p) are the fractions of the river and precipitation infiltration in the investigated well, respectively; and the weighting terms w R and w P are defined as follows: where 1 N is number of measurements; Q i is river discharge (m 3 /s); and P i is rainfall amount at An Long station (mm); 2 river infiltration coefficient: α R = 1 and α R = 0 for periods when the river water level is higher (losing streams) and lower (gaining streams) than the groundwater level, respectively. In this sense, α R = 1 (or α R = 0) indicates mass flux (or no mass flux) of tracer infiltrated from the river to a well; 3 precipitation infiltration coefficient: α P = 1 for months with precipitation infiltration (months in the rainy season) and α P = 0 for months without precipitation infiltration (dry season). We assumed that rainfall infiltrating to the shallow groundwater in the dry season is not significant due to the thickness (~9 m) of the upper low-permeable aquitard (see Section 2) and the small rainfall amount during the dry season.
The fractions of the river (p) and precipitation (1 − p) infiltration in each well can be derived by a linear mixing equation using long-term averages: where C well is the mean δ 18 O value of the investigated well and C inR and C inP are the weighted mean δ 18 O values of river and precipitation infiltration (weighted by the weighting terms w R and w P , respectively).
Hence, the optimization of two-component LPMs can be considered as an experimental approach to test the application of the selected TTDs.
Typically, each TTD function requires one or two fitting parameters. Table 2 summarizes the equations, the fitting parameters, and the predefined parameter ranges of TTDs used. The initial parameter ranges were assumed to be bounded, uniform distributions: • The range of the mTT (τ m ) was limited to a maximum of 250 weeks, equivalent to approximately the maximum 5 years that mTTs can be determined using stable isotopes of water (Maloszewski & Zuber, 1996 • The dispersion parameter (P D ) in the AD model should not exceed 2 for a constant tracer input (Maloszewski & Zuber, 1996). To improve the convergence of Monte Carlo simulations, we limited P D to 1 (cf. Cartwright & Morgenstern, 2016), which is appropriate for kilometre-scale flow systems (Mook & Rozanski, 2000).
• The range of the shape parameter (α) in the gamma model was limited to 10, following Timbe et al. (2014) and M. K. Stewart et al. (2017).
To set up the two-component LPMs, the free combination of two TTDs (e.g., the exponential combined with the dispersion model) would yield a large number of possible set-ups. This approach would require a multitude of assumptions, increase computational cost, and model uncertainties. For computational reasons, as well as to keep the model as simple as possible, we assumed that the TTDs of precipitation and river infiltration are of the same type (e.g., exponential or linear piston-flow TTDs). The number of models to be tested is thus equal to the number of selected TTDs. Theoretically, unreasonable models, as for example, the double exponential model, which includes an exponential TTD for precipitation, are also included in order to test the theoretical limitations with observations in the model fitting. The hypothesis is that theoretically unreasonable models should be rejected during the model fitting.

| Model performance and uncertainty analysis
Monte Carlo experiments were used to find the best parameter sets. The generalized likelihood uncertainty estimation (GLUE) methodology (Beven & Binley, 1992)  and RMSE as likelihood measures (see Equation (6)). The best 5% solutions in terms of D E were selected as behavioural models, from which we constructed 90% confidence intervals of the estimated mTTs (see Mosquera et al., 2016;Timbe et al., 2014). We also examined the dotty plots to check that the selected solution provides a reasonably wide range of behavioural parameter set. In this study, D E = 1 indicates a perfect fit.
3.6 | Data preparation Because data were collected at different temporal resolutions, from fortnightly to subweekly, all data were aggregated to weekly mean T A B L E 2 The TTD functions and their parameter ranges (assumed uniform distribution) for the GLUE analysis Abbreviations: P D , dispersion parameter (dimensionless); α, shape parameter (dimensionless); β, scale parameter (dimensionless); Γ, gamma function; η, parameter indicating the contribution of each flow type (dimensionless), expressed as the total volume/volume with exponential (or linear) TTD; τ m , subsurface mTT (weeks).
values for consistency. The early fortnightly river water samples were repeated in order to obtain a weekly time series. This simplification does not affect the mTT analysis, as these data were used for model warm-up only. The actual model fitting was performed for the period of subweekly sampling of both river water and precipitation.
Considering the different length of the input time series (9 years of river water and 3 years of precipitation; cf. Table 1), the time series of precipitation was repeated back to January 2009 using monthly weighted means from the 3 years of available data. This procedure implies a stable inter-annual variation of precipitation isotopes, a reasonable assumption in the VMD (see Duy, Heidbüchel, Meyer, Merz, & Apel, 2018). The approach does not change the results of the mTT estimation while giving the models more room to find stable results (Hrachowitz, Soulsby, Tetzlaff, & Malcolm, 2011).

| Modification of the input function
The isotopic fractionation of precipitation before infiltration (e.g., due to evaporation or mixing processes) should be considered during the calibration of LPMs. Precipitation falling on the ground likely mixes with local surface water (e.g., ponds, rice paddies, and irrigated, inundated, or wetland areas) and partly evaporates before infiltrating.
We considered the vertical infiltration as a mix of local surface waters and precipitation, both affected by evaporation. This assumption is in line with the suggestion that both river and evaporated-surface-water sources recharge groundwater along the Mekong river (Lawson et al., 2013;Lawson et al., 2016;Richards et al., 2018). Therefore, the input of precipitation infiltration was corrected, considering isotopic enrichment caused by evaporation.
We modified the input functions by adding a correction factor (Δ) to the isotopic composition of precipitation. This factor accounts for the isotopic enrichment due to the evaporation and mixing processes before the infiltration. This value was assumed to be constant and was derived by accounting for the potential evaporation in the region.
However, in order to quantify the uncertainty that is introduced by a constant correction factor, a sensitivity analysis of the model results to the isotopic correction of precipitation infiltration was conducted (cf. Section 3.8). All modified input functions are referred to as precipitation infiltration. The input of river infiltration was not isotopically corrected, implying no isotopic fractionation before and during the infiltration process.

| Model set-up
In order to get deeper insights into the model behaviour, parameter identifiability, and uncertainties, three LPM set-ups were defined (Tests 1, 2, and 3). The uncertainties were attributed to errors (1) the modified input functions (correction factor Δ), (2) the mass balance analysis (value of p integrated into the twocomponent LPMs), and (3) the assumed nonstationary or steady-state conditions during the calibration.
We assumed steady-state conditions to be dominant in the groundwater system and estimated time-invariant mTTs in Tests 1 and 2. For Test 3, the two-component LPMs were applied in a moving-window approach (e.g., Heidbüchel et al., 2012;Hrachowitz et al., 2009) to examine the time-variant TTDs and their corresponding mTTs. The detailed set-ups of these tests are as follows.

| Test 1: Modified input functions
We varied the correction values (Δ var ) in the range between 0‰ and 5‰ (with increments of 0.2) to create 26 modified input functions.
Δ var = 0‰ indicates no isotopic enrichment, and Δ var = 5‰ (mean value + standard deviation) represents the likely maximum isotopic enrichment before infiltration. The upper limit is derived from the distribution of differences between isotopic content of rainfall and pond water. It represents the mean value + 1 standard deviation, that is, the 84% quantile. Because we could not use another tracer (e.g., Cl) to independently assess the uncertainties of the mass balance analysis, the fraction of river infiltration was considered a fitting parameter (p cal ) and calibrated accordingly. Depending on the assumed TTD function, three or five parameters were fitted during the calibration with each of the 26 modified input functions. In this test, 9-year records of δ 18 O were used (the first 6 years for a warm-up and the last 3 years for analysis).

| Test 2: Mass balance
We fixed a correction value (Δ fix = 1.81‰), defined by the isotopic difference between the arithmetic mean value of the pond water and the weighted mean value of precipitation. The aim was to match the mean values of the modified input function and the pond water, implying that the precipitation ponding on the ground is mixed with preexisting local surface water and partly evaporates before infiltrating to the groundwater. In this test, the fraction of river infiltration was predefined by Equation (5) and used as a fixed parameter (p fix ).
The motivation for this test was that the contribution of river (or precipitation) infiltration should be identical independent of the selected TTDs. This was tested with this approach. Consequently, two or four fitting parameters (depending on the assumed TTD) were fitted with the modified input function created by Δ fix . Similar to Test 1, we used 9-year records of δ 18 O.

| Test 3: Nonstationarity
Nonstationarity is implicitly considered by the weighing according to discharge (Equation (3)) and the alpha coefficient determining the seasonal variation between losing and gaining stream conditions.
where n is the number of observations, k is the number of fitting parameters, and SSE is the sum of squared errors.

| Surface-groundwater interaction
Considerable dynamics of surface-groundwater interaction were observed at the study site, as depicted by the similar water level variations up to 2 m annually in the river and groundwater (Figure 2).
Seasonal changes in groundwater levels observed between Wells A and B mostly lay between those found in Well C and the river. The groundwater level observed at wells closer to the river (Wells A and B) exhibited a higher seasonal variation than the site farther from the river (Well C).
Gaining or losing stream conditions were defined as river water penetrating into the groundwater system or groundwater seeping out into the river, respectively (Fetter, 2001). Losing stream periods (higher monthly river water level than groundwater level) were detected mainly during the flood season from July to November, whereas gaining stream periods were observed primarily during the end of the dry season from April to June (Figure 2). From December to February, the differences between river water and groundwater levels were insignificant. These months were considered as the transition period between losing and gaining stream conditions.   Figure S1. The fitting accuracy generally increased with an increasing correction factor (Δ var ) up to 1.4‰, remained stable around the peak for Δ var between 1.4‰ and 2.6‰, and decreased slightly after that. The best performances were identified for the LPF model (Figure 4d), followed by the AD model ( Figure 4e). Notably, the optimum correction value (Δ var = 1.8‰), providing the best LPF model performance, was very close to the isotopic difference between the F I G U R E 2 The daily groundwater and river water levels at An Long. Grey background indicates losing stream periods when the monthly river water level is higher than the monthly groundwater level arithmetic mean value of the pond water and the weighted mean value of precipitation (Δ fix = 1.81‰). This supports the validity of the assumed fixed correction factor.

| Stationary transit time modelling
To evaluate the sensitivity of the parameter identifiability to the isotopic input correction, we focused on the LPF models providing the best-fit accuracies. The behavioural solutions (90% F I G U R E 4 Sensitivity of the efficiency of two-component lumped parameter models to isotopic correction values (Δ var ) measured by using the Euclidean distance (D E ) between (1 -Kling-Gupta efficiency [KGE]) and root mean square error values. The green reference line indicates the isotopic difference (Δ fix ) between the arithmetic mean value of the pond water and the weighted mean value of precipitation. Δ var-best corresponds to the best-possible model performance and reasonable parameter identifiability (see Figure 5). Negative likelihood measures (D E < 0) are not shown F I G U R E 3 The isotopic data at An Long in a dual isotope plot. The regression line (RL) for groundwater is derived using all groundwater samples from all the wells. GW, groundwater; LMWL, local meteoric water line; VSMOW, Vienna Standard Mean Ocean Water confidence bound of the GLUE analysis) corresponding to a threshold of 5% of the best predictions by LPF models are shown in Considering each investigated well, the estimations of river mTTs were identical for the identified range of acceptable isotopic corrections (Figure 5a-c, blue boxes in the grey-shaded areas).
Shorter mTTs of river infiltration were determined consistently for sites closer to the river. The optimized river mTTs ranged approximately from 15 to 20 weeks for Well A, from 35 to 40 weeks for Well B, and from 25 to 240 weeks for Well C. The parameter identifiability of river mTTs was better for the sites close to the river (e.g., Wells A and B).
F I G U R E 5 Sensitivity of the parameter identifiability of the linear piston flow model (i.e., the best-performing model) to the correction values (Δ var ) for all tested wells. The box plots indicate the 90% confidence intervals of fitting parameters given by the generalized likelihood uncertainty estimation analysis. Each box shows the interquartile range, the whiskers show the minimum and maximum values associated to 1.5 times the interquartile range. p (calibrated parameter) and p fix (predefined by Equation (5)) are the fractions of recharge from river (or precipitation) infiltration; mTT is the tracer's mean transit time; η is a parameter indicating the ratio of total volume/volume with linear transit time distribution. Δ var-best indicates the range for reasonable parameter identifiability. The explanation of Δ fix is similar to Figure 4 The behavioural solutions of precipitation mTTs were identical for Wells B and C, independent of the isotopic correction and the distance of wells to the river. For Well A, the behavioural mTTs were also very similar, with a tendency to a larger range for higher acceptable correction factors. Compared with riverbank infiltration, the precipitation infiltration exhibited longer mTTs. Optimized precipitation mTTs were between 75 and 110 weeks with low uncertainties for all tested sites (Figure 5a-c). However, the LPF model provided poor constraints of parameter η, indicating different ratios of linear or piston flows at different investigated sites (Figure 5d-f). The uncertainty bounds of river infiltration fraction (p cal ) were quite narrow and increased with higher Δ var (Figure 5i-k). Analogously, the fraction of precipitation infiltration (1 − p cal ) decreased with higher Δ var .
Results in Test 2 were compared with the best-fit results in Test 1 (e.g., using Δ var = 1.8‰), considering the model efficiency (Figure 6), the fractions of water components contributing to the shallow groundwater ( Figure 5 i-k), and the behavioural solutions of optimized mTTs (Figure 7). The best-fit results (LPF models) of the tests are reported in Table 3. The model performances (corresponding to the best-matching likelihood measures) of these tests were comparable ( Figure 6a,b,d,e) and relatively good (>0.7) in terms of the KGE statistic for all investigated wells. Although the model parsimoniousness is slightly better within the set-ups in Test 2, illustrated by lower BIC values (see Figure 6c,f and Table 3), the dotty plots indicate that the parameter identifiability of two of these tests is comparable (see Figures S3 and S4). Better parameter identifiability of river mTTs is observed for sites close to the river (e.g., Well A). Conversely, parameter identifiability of precipitation mTTs is much better for sites farther from the river (e.g., Well C).

| Identification of best-suited TTD
Out of all 18 models (three tested sites and six models per site), the threshold of model acceptance (KGE > 0.5) was fulfilled in 13 cases (Figure 6a,d). The five poor models with KGE < 0.5 were the exponential (E), the linear (L), and the gamma (G) models at two sites located farther from the river (Wells B and C). The EPF, the LPF, and the AD models provided satisfactory performances for all sites. Unsurprisingly, the more complex models (EPF, LPF, AD, and G) performed better than the simpler models (E and L),  (Figure 7),  Figure S4). Figure 8 shows the best-fit modelled δ 18 O of LPF models and the uncertainty interval (90 % confidence bound of the GLUE analysis). Considering each investigated well, the best-fit mTTs of river infiltration were identical (stable in time) with acceptable parameter uncertainties ( Figure 9). The best-fit mTTs of river infiltration were approximately between 16 and 24 weeks for Well A, between 36 and 48 weeks for Well B, and between 17 and 55 weeks for Well C. The uncertainties of river mTT were better constrained for sites close to the river. Regarding precipitation infiltration, the best-fit mTTs were relatively similar (mainly around between 85 and 125 weeks) for all investigated wells. However, the uncertainties of precipitation mTTs were poorly constrained for the Well A close to the river compared with Well C farther from the river. In general, the uncertainty of river mTTs increased with distance to the river, whereas the uncertainties for precipitation infiltration mTTs decreased with distance to the river.

| Time-variant transit time modelling
Also, stationary and time-variant mTTs estimated by LPF models were comparable, in terms of both best-fit and behavioural solutions (Table 3).

| Mechanisms and sources of groundwater recharge
Similar seasonal fluctuations between groundwater and river levels over the monitoring period ( Figure 2) suggested a good hydraulic connection between surface water and groundwater along the Mekong river. The semi-annual reversal of gradients between the river and the groundwater indicates groundwater recharge, that is, that the river loses water to bank infiltration and recharges the Holocene aquifer during flooding. In contrast, groundwater is released from the aquifer to compensate for the small amount of river water at the end of the dry season. Considering the difference of the hydraulic conductivity F I G U R E 7 Estimated mean transit times (mTTs) of river (top) and precipitation (bottom) infiltration in Tests 1 (left) and 2 (right). The error bars indicate the 90% confidence intervals of mTT given by the generalized likelihood uncertainty estimation analysis and the elevation of the aquitard and aquifer layers, the shallow groundwater is mainly in horizontal hydraulic contact with the river via bank infiltration at the highly permeable aquifer (characterized by the medium and coarse sand layer), instead of the low-permeable aquitard (characterized by the silt and clay layers).
Note. The simulated results and model efficiencies are presented for the best-matching likelihood measures (D E ). The model set-ups (including modified input functions, calibrated parameters, and assumed conditions) of these tests are described in Section 3.8.
a Mean values of the 29 best-fit LPF models are reported.
isotopes observed in the groundwater with increasing distance from the river (Figure 3) might be explained by the mixing of riverbank infiltration with a more significant contribution from an evaporated recharge source to the shallow groundwater. The isotopic signatures indicate the different importance of the two sources for the groundwater recharge at the investigated sites, as one would also expect F I G U R E 9 Model efficiencies (Kling-Gupta efficiency [KGE]) and the mean transit times (mTTs) of river and precipitation infiltration corresponding to the best-matching linear-piston flow models. The shaded areas represent the best behavioural solutions (90% confidence) of mTT predictions by generalized likelihood uncertainty estimation analysis. Results of time-variant mTTs are shown for every 2 weeks F I G U R E 8 The observed and modelled δ 18 O plotted with the behavioural solutions (90% confidence bound of generalized likelihood uncertainty estimation analysis) corresponding to a threshold of 5% of the best prediction by the linear-piston flow models. Error bars indicate the analytical reproducibility of the δ 18 O measurements. VSMOW, Vienna Standard Mean Ocean Water from a hydraulic point of view due to the different distances to the river. With increasing distance, the exchange between groundwater and river is usually dampened due to infiltration length and associated longer transit times.
The groundwater regression line deviated significantly from the LMWL exhibiting a less steep slope, whereas it compared well with the regression line of the pond water ( Figure 3). This indicates that the local surface water, being affected by evaporative fractionation processes, is likely a second source recharging the shallow groundwater. The slopes derived from the regression lines of groundwater and pond water in our study are comparable with those in Lawson et al. (2013), Lawson et al. (2016), and Richards et al. (2018), who also suggested an evaporated source of surface recharge (e.g., from the wetland and ponds) to the groundwater in Cambodia. Such comparable results support the assumption that the precipitation is mixed with preexisting local surface water and evaporates before infiltrating down the unsaturated zone towards the shallow groundwater in the VMD. Considering the offset of isotopic samples of the river from the local precipitation (Figure 3), it is unlikely that the river was recharged by the local precipitation but is preferably sourced from upstream of the VMD (i.e., its run-off stems almost exclusively from the Mekong basin).
Generally, the analysis of water level fluctuations and isotopic signatures suggests that the shallow groundwater is recharged by two distinct water components via different flow pathways, justifying the use of two-component LPMs, beyond simple fitting considerations.

| Sensitivity of modelling results to isotopic correction
The sensitivity analysis (Test 1) illustrates that using modified input functions within the calibration of two-component LPMs can provide better fitting accuracies (Figure 4) without altering the estimated mTTs ( Figure 5). Better performance was also reported when LPMs were calibrated with input functions modified by canopy interception (e.g., Stockinger et al., 2014), evapotranspiration (e.g., Stumpp, Stichler, & Maloszewski, 2009), or the correction of the isotopic mass balance (e.g., Viville, Ladouche, & Bariac, 2006). In this study, the correction was necessary because of (a) isotope enrichment in ponding water and (b) theoretical preconditions for the application of the mixing models, specifically the assumption that the water at the wells is the product of mixing of both river water and precipitation. has to be noted that the quantification of contributions from the two sources is sensitive to the isotopic correction of precipitation, depicted by the increase of p cal with higher Δ var .
A comparison between Tests 1 and 2 illustrates that the predefined mixing process (e.g., using Equation (5)) and fixed correction for the δ 18 O precipitation (e.g., Δ fix = 1.81‰ to produce the same mean value of precipitation infiltration as the ponding water) provide reasonable modelling results. In this context, it is noteworthy that the best model results are obtained with Δ var = 1.8‰, which is almost identical to Δ fix . This similarity supports the validity of the definition of Δ fix by the mean difference between precipitation and ponding water isotopic content. The fixing of the correction factor not only fulfils real-world and theoretical constraints but also improves the parsimoniousness of the applied model. A similar correction of δ 18 O of precipitation (1.4‰) was applied by Calderon and Uhlenbrook (2016) for a hydrograph separation in Nicaragua, a climatic environment somewhat similar to this study.

| Dominant subsurface flow conditions
The relatively high performance (KGE > 0.7 and RMSE < 0.4‰,  with the related linear models (e.g., LPF model; see Figure 6). This suggests that the subsurface transport of water is better characterized by a linear distribution rather than by an exponential distribution. This confirms that the exponential model is inadequate to represent recharge to groundwater collected at larger depths below the ground surface (Zuber et al., 2011). Secondly, mTT modelling results suggest a considerable fraction of piston flow in the TTDs. For example, considering the results of Test 2 from the LPF model for Well A (see Table 3), the best-fit value of η = 1.82 for river water implies that 55% of the volume of the river water infiltration passes through the aquifer as linear flow, whereas 45% can be characterized by piston flow behaviour. Accordingly, the value of η = 1.90 for precipitation implies a 53% of volume portion of linear flow and a 47% volume of piston flow in the TTD of precipitation infiltration.
The statistical findings indicate that the subsurface flow condition at the study site is likely best described with a linear distribution accounting for the infiltration along the river followed by the hydraulic replacement of groundwater caused by pressure gradients that adds the piston flow component to the model. The explanation is consistent with the hydrogeological setting at the study site characterized as a partially confined aquifer that does not create a phreatic system.
This situation can be represented by a linear-piston model or a dispersion model according to the five hydrogeological settings described in Małoszewski and Zuber (1982). The accordance of the statistical finding with these theoretical considerations serves as a corroboration of the ranking of models (see Section 4.4) by the model fitting and the GLUE analysis.
Compared with the gamma model, the EPF and AD models provide better performances ( Figure 6). The results agree well with the dominance of these TTDs in riverbank infiltration studies (e.g., Kármán et al., 2014;Maloszewski et al., 1992;Stichler et al., 1986;Stichler et al., 2008) and/or in groundwater studies (e.g., Cartwright & Morgenstern, 2016;Stewart et al., 2017), whereas the gamma model has been frequently used in catchment studies (cf. K. J. McGuire & McDonnell, 2006;Hrachowitz et al., 2010). Although the LPF distribution has been introduced early (see Maloszewski & Zuber, 1996;Małoszewski & Zuber, 1982), it has rarely been tested, to our best knowledge, within the lumped parameter approach. Our study agrees well with Timbe et al. (2014), who suggested that the LPF model could be a reliable method to determine water transit times in southern Ecuador.

| Two-component LPM reveals recharge mechanism
Comparisons of the results (e.g., estimated mTT, parameter identifiability, and model efficiency) of three tests (  Table 3)  to the river (Wells A and B), but cannot be constrained for sites farther from the river (Well C). Notably, using stable isotopes alone cannot provide reliable longer mTTs, despite the used LPMs (Seeger & Weiler, 2014;Kirchner, 2016a), indicating the contribution of older water (>5 years) to the groundwater system of the VMD. The mTTs of precipitation infiltration were independent of the distance to the river, depicted by the relative similarity of mTTs (82-95 weeks) for all investigated wells. The fact that the estimated mTTs of precipitation infiltration are longer than the ones of river infiltration is attributable to the hydraulic conductivity of the soil. The horizontal infiltration from the river takes place mainly via the highly permeable aquifer, resulting in short mTTs (<40 weeks) for the investigated wells located close to the river (e.g., <200 m). Meanwhile, the vertical infiltration of precipitation (after ponding on the surface) takes place primarily via a lowpermeable overlying aquitard, resulting in considerably longer mTTs (>80 weeks) for all investigated wells.
Overall, the results follow the general understanding of groundwater hydraulics and are reasonable from physical and hydrological points of view, corroborating the applicability of two-component LPMs to identify groundwater mTTs at riverbank infiltration systems.
However, in the given lithological setting, the predictive skill and particularly the reliability of the models decrease for locations farther from the river, where recharge by precipitation dominates and a lowpermeable aquitard layer above the aquifer is present. This specific setting impairs the identifiability of model parameters in this case. In other settings, for example, without an overlying aquitard, better model performance and parameter identifiability can be expected even for larger distances to the river.
F I G U R E 1 0 Conceptual model of subsurface flow conditions at the study site. The mean transit times (mTTs; arrows) and recharge contributions (pie charts) of river water and precipitation shown here result from the linear-piston flow model. The characteristics of the aquitard and aquifer layers (e.g., the thickness, type of soil, porosity, and vertical hydraulic conductivity) are referenced from Boehmer, 2000;Minderhoud et al., 2017)

| Limitations and wider implications
Although the identified results are hydrologically plausible, corroborating the validity of the model concept, we acknowledge the limitations of this study. Here, we point out possible reasons related to the inevitable errors of two-component LPMs and shortcomings of data that could add uncertainties to the mTT analysis.
First, the lumped convolution modelling approaches rely on steady-state conditions and assumed nonstationary TTDs (Tests 1 and 2). Such assumptions are probably less problematic for groundwater than they are for surface water systems (Christian Birkel et al., 2016), yet they are rarely met in any real hydrologic setting (Rinaldo et al., 2011 Fourth, the study would have benefited from higher sampling frequencies (e.g., daily) to provide better insights into short-term system responses (see Hrachowitz et al., 2010;Christian Birkel et al., 2016).
Higher sampling resolution could also improve model conceptualization and calibration (e.g., C. Birkel et al., 2010) and reduce potentially misleading insights (Hrachowitz et al., 2011) and uncertainties of mTT modelling (Timbe et al., 2015). However, given the practical constraints and costs of isotope sampling, higher sampling frequencies are difficult to realize in general.
Fifth, the reconstruction of the precipitation record (see Section 3.6) could be a potential source of error. However, the approach of looping precipitation isotopic signature has been common practice when input time series are too short to constrain mTT estimates adequately (e.g., Timbe et al., 2014;Christian Birkel et al., 2016;Mosquera et al., 2016;Muñoz-Villers et al., 2016). For more reliable palaeoclimate reconstructions of precipitation isotopes in Asian monsoon regions, model-based statistical approaches (e.g., the combination of global climate models with statistical analyses) could be applied (see Duy et al., 2018, and references therein).
Sixth, choosing the same types of TTDs to combine in the twocomponent LPMs cannot provide an entire picture of all possible combination of selected TTDs. Although mixing different kinds of TTDs (e.g., the combination of AD model for precipitation infiltration and LPF model for river infiltration) could improve the model performance, we expect that this approach invalidates the estimated mTTs, because most of the TTD types are relatively flexible and tend to accommodate themselves to the data. However, this approach should be considered in further study.
Finally, although the isotopic correction due to the evaporation process is in line with the hydrological setting (discussed in Section 3.7) and essential to fulfilling the theoretical constraints (discussed in Section 5.2), adding a constant value to the isotopic signature of rainfall to estimate the isotopic signature of vertical recharge unavoidably introduces uncertainties. This approach assumes a stable isotopic fractionation process for the whole study period, which is probably quite unrealistic. The sensitivity analysis revealed that although the mTTs are relatively insensitive to the correction factor, the identifiability of the contribution of the different sources to recharge is impaired by the correction factor. Moreover, if the water does pond on the surface before it infiltrates, the isotopic signature may be attenuated prior to recharge. This will result in mTTs being overestimated, because the presented approach does not account for the time required for evaporative changes in isotopic composition, but not for the time required for this. With the available data, it is impossible to analyse the isotopic enrichment of local surface water (e.g., Skrzypek et al., 2015) or independently assess actual contributions of infiltrated water components in this study. This could, for example, be achieved by using another tracer (e.g., Cl).

| CONCLUSIONS
This study investigated groundwater transit times and subsurface flow conditions at the riverbank infiltration areas in the VMD.
Precipitation, river, groundwater, and local surface water were sampled on a subweekly to weekly basis for different periods between 2009 and 2017 and analysed for stable isotopes. The applicability of two-component LPMs (allowing different TTDs for different recharge components) in conjunction with hydrological and isotopic measurements to identify subsurface flow conditions and the contribution to groundwater mixing was tested. The proposed method proved to be able to identify the TTDs and their corresponding mTTs of both river and precipitation infiltration to shallow groundwater using δ 18 O records.
LPMs based on the LPF distribution were able to capture isotopic variations in shallow groundwater in response to the modified input function. Although the exact contribution of the water components infiltrating to the groundwater system remains uncertain, the dynamics of the surface-groundwater interaction could be identified. River water infiltrates horizontally mainly via the highly permeable aquifer, resulting in short mTTs (<40 weeks) for locations close to the river (<200 m). The vertical infiltration from precipitation takes place primarily via a low-permeable overlying aquitard, resulting in considerably longer mTTs (>80 weeks). The outcomes are hydrologically plausible, corroborating the validity of the applied approach. Our findings enhance the understanding of the shallow groundwater recharge dynamics and may serve as a baseline for future groundwater studies using environmental isotopes in the VMD. Groundwater resources management needs to consider the different recharge mechanisms and mTTs (mainly controlled by the distance to the river), resulting in different management options for different areas in the delta.
Our study suggests that the highly complex mechanism of surface-groundwater interaction and subsurface mixing processes at riverbank infiltration systems can be conceptualized by exploiting two-component LPMs. Regardless of the restrictions associated with certain errors of LPMs and the use of stable isotopes, the model concept can be transferred to other locations. Therefore the proposed model concept with the associated model selection procedure could provide a comprehensive hydrological tool for the analysis and understanding of groundwater recharge by different sources.

ACKNOWLEDGMENTS
We acknowledge the German Academic Exchange Service (DAAD) for the scholarship that enabled the first author to pursue a PhD degree at the University of Potsdam, Germany. We thank the staff at An Long DT/14-19/C11).