Caffeine vs. carbamazepine as indicators of wastewater pollution in a karst aquifer

This paper presents the analysis of caffeine and carbamazepine transport in the subsurface as a result of wastewater release in the Sorek creek over the outcrops of the carbonate, Yarkon-Taninim, aquifer in Israel. Both caffeine and carbamazepine were used as indicators of sewage contamination in the subsurface. While carbamazepine is considered conservative, caffeine is subject to sorption and degradation. The objective of the study was to quantify differences in their transport under similar conditions in the karst aquifer. Water flow and pollutant transport in a “vadose zone–aquifer” system were simulated by a quasi-3-D dual permeability numerical model. The results of this study show that each of these two pollutants can be considered effective tracers for characterization and assessment of aquifer contamination. Carbamazepine was found to be more suitable for assessing the contamination boundaries, while caffeine can be used as a contaminant tracer only briefly after contamination occurs. In instances where there are low concentrations of carbamazepine which appear as background contamination in an aquifer, caffeine might serve as a better marker for detecting new contamination events, given its temporal nature. The estimated caffeine degradation rate and the distribution coefficient of a linear sorption isotherm were 0.091 d⁻¹ and 0.1 L kg⁻¹, respectively, which imply a high attenuation capacity. The results of the simulation indicate that by the end of the year most of the carbamazepine mass (approximately 95 %) remained in the matrix of the vadose zone, while all of the caffeine was completely degraded a few months after the sewage was discharged.

1 Introduction

Sewage infiltration into the subsurface can cause groundwater pollution. Carbonate aquifers present a higher risk for groundwater quality contamination due to the presence of preferential flow paths. Predicting and quantifying sewage infiltration and transport in carbonate aquifers is complicated due to exchanges between slow flow in the matrix and fast flow in conduits (Geyer et al., 2007).

The micropollutants carbamazepine (CBZ) and caffeine (CAF) are both widely used as indicators of anthropogenic contamination in groundwater (Seiler et al., 1999). CBZ is generally accepted as a stable indicator of untreated/treated sewage (Clara et al., 2004; Dvory et al., 2018a; Fenz et al., 2005; Gasser et al., 2010). CAF concentration is often higher than CBZ concentration at the source (sewage) and consequently detected in soils, sediments (Bradley et al., 2007; Klabunde, 2016), surface water (Bueno et al., 2010; Chen et al., 2006; Ferreira, 2005; Kolpin et al., 2004), lakes, and seawater (Buerge et al., 2003; Gardinali and Zhao, 2002; Knee et al., 2010). Additionally, several studies have detected higher concentrations of CAF, as opposed to CBZ, in the groundwater as well (Godfrey et al., 2007; Lapworth et al., 2015; Manamsa et al., 2015; Metcalfe et al., 2011). CAF, like other micropollutants, is subject to degradation and sorption. These processes reflect major mechanisms for CAF attenuation in the environment (Hillebrand et al., 2015) and raise doubts as to its efficacy as a tracer for detecting and quantifying contamination from wastewater (Ahmed et al., 2008; Liu et al., 2014). The above implies that CAF concentration is often higher and easier to detect than CBZ close to the pollution source briefly after contamination occurs. As a result, CAF was found to be a possible indicator of sewage in rapid flow systems, such as karst aquifers (Hillebrand et al., 2012b).

CAF attenuates in the subsurface by biodegradation and sorption processes (Martínez-Hernández et al., 2016). Those processes are affected by local environmental conditions. Recent studies reveal fast biodegradation rates in carbonate aquifer conduits (Hillebrand et al., 2012b, 2015). In such aquifers, the flow is also affected by the connections between the matrix and conduit flow paths, which can influence the CAF attenuation. Lab (Arye et al., 2011; Conn and Siegrist, 2009; Hebig et al., 2017; Martínez-Hernández et al., 2016) and field (Hillebrand et al., 2012b; Zhang et al., 2013) experiments were carried out by many researchers in order to assess CAF sorption. Sorption parameters can vary as a result of the groundwater hosting media. Hebig et al. (2017) showed that full removal of CAF was observed in the presence of organic carbon and no sorption was detected in iron-coated sand. The attenuation of CAF in the unsaturated zone can be very high, as shown by Martínez-Hernández et al. (2017) by column test experiments and simulations.

Until the current study, simulations of micropollutant transport in the subsurface were mostly performed with either single (Bertelkamp et al., 2014) or dual (Geyer et al., 2007; Hillebrand et al., 2012a, b) porosity 1-D models for saturated or unsaturated (Martínez-Hernández et al., 2017) conditions. However, under field-scale conditions, the fate of contaminants is affected by both the transport through the vadose zone and lateral spread in the groundwater.

The objective of this study was to assess differences in CAF and CBZ attenuation by simulating their transport in the karst/fractured-porous unsaturated zone and groundwater system with a dual permeability mathematical model as described by Dvory et al. (2018a). The sorption and degradation parameters for CAF were estimated using observed concentrations in a single well (Fig. 1). Simulation results allowed for the characterization of CAF transportation and natural attenuation processes from the initial release of sewage at the Sorek creek bed, its infiltration into the unsaturated zone and transport into the aquifer. The study employed field observations and simulations of CAF and CBZ transport to reveal the effect of factors leading to differences in the attenuation of each micropollutant.

The full description of the field experiment, the mathematical model development, and the simulation results of CBZ transport are presented in Dvory et al. (2018a). In this paper we provide a short description of those processes, for the convenience of readers, and include only details which are essential for understanding of the presented material.

2 Materials and methods

2.1 Sewage release event

The Sorek creek watershed drains approximately 88 km² in the study area and is located west of the city of Jerusalem, Israel (Fig. 1). A local reservoir (Beit Zait), located 2.05 km upstream from the study site, collects surface flow and limits the natural runoff downstream. Periodically, the reservoir discharges its reserves; once every few rainy seasons regular controlled releases from the dam downstream may occur (Dvory et al., 2018b).

The geology of the area is comprised of a carbonate section of the Judea group (Dvory et al., 2016). The unsaturated zone is thick, spanning tens of meters to 200 m. The groundwater primarily flows in the southwesterly direction (Dvory et al., 2016).

The current study examined a discharge event, which took place between 2 and 19 April 2013, when wastewater was released from a main sewage pipeline on five separate occasions into the Sorek creek (Figs. 1 and 3c) (Dvory et al., 2018a). CBZ served as an indicator for the identification and quantification of sewage water migration into the aquifer (Dvory et al., 2018a) and CAF was used to assess its attenuation and suitability as a tracer for wastewater contamination characterization.

2.2 Field work

Field work, including water sampling and hydrological monitoring, was done to monitor the distribution of the discharged sewage (Dvory et al., 2018a). Over the course of 310 days, 23 groundwater samples were taken from a depth of 100 m below the ground surface. The intervals between sampling events ranged from 1 to 56 days, where the interval was shorter during the expected tracer breakthrough time in order to provide a higher temporal resolution of the CBZ and CAF breakthrough tracer curves (BTCs). Additionally, a pressure and temperature probe with data logging capabilities was installed in observation well EK11 (Fig. 1) in order to take groundwater level and temperature data measurements every 30 min (Solinst Levelogger). Hourly measured values for precipitation and evaporation rates were acquired from the local Israel Meteorological Service (IMS) weather station (“Tzuba Station”). Data on sewage and surface runoff discharge rates were obtained from gauging stations upstream and downstream from the well head (Fig. 1).

Figure 1 The upper Sorek Basin monitoring sites and flow and transport simulation domains (following Dvory et al., 2018a; aquifer boundaries from Dafny et al., 2010).

2.3 Analytical methods

A methodology using the Agilent G6410A Triple Quadrupole mass spectrometer (QQQ) with an electrospray ionization ion source (ESI) was used in order to identify and verify CBZ and CAF presence in collected water samples. The LC/MS/MS method's main characteristics are presented in Table 1. The method was developed from the EPA Method 1694 guidelines for acid and basic compound elution from solutions containing less than 1 % solids.

Table 1 LC/MS/MS method main characteristics for CBZ and CAF.

Measured amounts of labeled compounds (CBZ-d10 and CAF 13C3) were added to the water samples. The compounds were then extracted by solid-phase extraction (SPE) with Oasis HLB 60 mg cartridges (Waters, Milford, MA) using 1000 mL of each sample. Analytes were subsequently eluted with methanol and formic acid solutions, and the mixed extracts were concentrated to a final volume of 5 mL by nitrogen flow. Analytes were separated with an Agilent ZORBAX Eclipse Plus C18 (2.1 mm ID, 100 mm length, 3.5 µm particle size). Column temperature was set at 25 ^∘C. The mobile phase consisted of 10 % Acetonitrile, 90 % H₂O, and 0.1 % formic acid. The eluent composition included initial conditions, and 10 % Acetonitrile fed at 0.2 mL min⁻¹ for 5 min. After 6 min, the flow rate was increased to 0.3 mL min⁻¹, and eluent composition was changed to 60 % Acetonitrile, fed at a flow rate of 0.3 mL min⁻¹ until 24 min. In the last stage, eluent composition was gradually increased until it reached 100 % Acetonitrile at 30 min. Injection volume was 15 µL (Ferrer and Thurmann, 2008).

Quantifications were carried out by isotope-labeled internal standards for CBZ and CAF, by multipoint calibrations. Limits of quantification (LOQs) are shown in Table 1 and were calculated at 10 times the background levels, along with the recovery at 1000 ng L⁻¹. The linearity of the response of 3 orders of magnitude was demonstrated (R²>0.99) for both the micropollutants studied.

2.4 Modeling

2.4.1 Mathematical model

Variable saturation flow and contaminant transport in the vadose zone and the groundwater were simulated using a quasi-3-D model. The model conceptual sketch is shown in Fig. 2. Two overlapping continua representing highly conductive karst/fractures (conduits) (c) and the low-permeability matrix (m) were used to simulate the karst/fractured-porous medium in both the unsaturated and saturated zones. The quasi-3-D approach (Levy et al., 2017; Twarakavi et al., 2008; Yakirevich et al., 1998) was used as the basis for the numerical model. The model uses a series of 1-D equations in a variably saturated zone to simulate the “vadose zone–aquifer” system and 3-D equations for saturated flow and transport to simulate groundwater. At the dynamic phreatic surface, the 1-D and 3-D equations are coupled (Kuznetsov et al., 2012). Unsaturated flow in high- and low-permeability regions is described by two Richards equations accounting for the linear exchange kinetics between them (Gerke and van Genuchten, 1993; Dvory et al., 2016). Horizontal flow in the vadose zone was neglected. The pollutant transport equations in the conduits and matrix, respectively, are

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{r}_c\left[{\displaystyle \frac{\partial \left({\mathit{\theta }}_c{C}_c\right)}{\partial t}}+{\mathit{\rho }}_{\mathrm{b}}{\displaystyle \frac{\partial F_c}{\partial t}}\right]=\mathrm{\nabla }\cdot \left[r_c\left({\mathit{D}}_c{\mathit{\theta }}_c\mathrm{\nabla }{C}_c-{\mathit{q}}_c{C}_c\right)\right]\\ \text{(1)}& {\displaystyle }& {\displaystyle }-M_{cm}-\mathit{\lambda }{r}_c\left({\mathit{\theta }}_c{C}_c+{\mathit{\rho }}_{\mathrm{b}}{F}_c\right){\displaystyle }& {\displaystyle }\left(\mathrm{1}-r_c\right)\left[{\displaystyle \frac{\partial \left({\mathit{\theta }}_m{C}_m\right)}{\partial t}}+{\mathit{\rho }}_{\mathrm{b}}{\displaystyle \frac{\partial F_m}{\partial t}}\right]\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\partial }{\partial z}}\left[\left(\mathrm{1}-r_c\right)\left(D_m{\mathit{\theta }}_m{\displaystyle \frac{\partial C_m}{\partial z}}-q_{zm}{C}_m\right)\right]\\ \text{(2)}& {\displaystyle }& {\displaystyle }+M_{cm}-\mathit{\lambda }\left(\mathrm{1}-r_f\right)\left({\mathit{\theta }}_m{C}_m+{\mathit{\rho }}_{\mathrm{b}}{F}_m\right),\end{array}$$

where C_i and F_i are the concentrations in liquid and sorbed phases, respectively (M L⁻³ and M M⁻¹); i=c and m for the conduit and porous matrix; θ_i is water content (L³ L⁻³); ρ_b is the rock bulk density (M L⁻³); D_i is the hydrodynamic dispersion tensor (L² T⁻¹); q_c is Darcy's flux of water in conduits (L T⁻¹); q_zm is the vertical water flux in blocks (L T⁻¹); M_cm is a term which accounts for the solute exchange between the conduits and porous matrix, respectively (M L⁻³ T⁻¹); λ is the degradation rate (T⁻¹); r_c is the relative conduit volume; z is the vertical coordinate (L); and t is time (T).

Under field-scale conditions, when larger-scale hydrological models are applied, there is not usually detailed quantitative knowledge of these controlling factors, and simple first-order degradation is often assumed (Bradley et al., 2007). Therefore, in the present investigation the degradation rate coefficient for conduits and matrix is assumed to be the same.

Sorption of a solute is described by the following linear isotherm:

$$\begin{array}{}\text{(3)}& F_i=K_{\mathrm{D}}{C}_i,\end{array}$$

where K_D is the distribution coefficient (L³ M⁻¹).

At time zero, the initial flow condition prescribes the pressure head or water content distribution along the simulation profile. Boundary conditions are the type of transient head or flux (Fig. 2). The changes in concentration are defined at the inflow boundaries, while a zero concentration gradient was prescribed at the outflow boundaries. The distribution of CBZ initial concentrations was discussed in Dvory et al. (2018a), and based on several measurements. A zero-level concentration for CAF was used for the entire domain.

2.4.2 Numerical model setup

The aforementioned equations were solved using a method of finite differences. The MODFLOW model (Harbaugh et al., 2000) was modified to incorporate the 1-D Richards equations (Kuznetsov et al., 2012) and to account for a double permeability approach. In order to simulate solute transport in the vadose zone and groundwater, the MT3D (Zheng and Wang, 1999) numerical code was also modified. Pre- and post-processing data were conducted with the Groundwater Modelling Software (GMS 6.0, 2002).

A three-step approach was used to address flow and transport problems (Fig. 2): (1) flow was simulated in a large domain with a coarse grid and well-defined hydrogeological boundaries; (2) flow and CBZ transport were simulated in a small domain, using a more refined grid (Fig. 1); and (3) flow and CAF transport were simulated in the small domain. This sequential process was used in order to minimize processing time and increase the accuracy of the solution. Both the large and small domains had a uniform grid on the horizontal plane and a variable size grid in the vertical plane, which varied in respect to the ground level and aquifer base altitudes. The large domain was 1600×1400 m on the horizontal plane and varied from 170 to 445 m on the vertical plane. The grid size was 20×20 m on the horizontal plane. The vertical plane of the grid was composed of 38 layers which increased from 0.0002 T at the ground surface to 0.1 T at the aquifer bottom (T(x,y) is the aquifer thickness including the unsaturated zone). Simulations in the large domain were used to calibrate hydraulic parameters and the western lateral boundary condition for the flow model.

The size of the small domain was 590×460 m on the horizontal plane and varied from 185 to 280 m on the vertical plane. The grid was a uniform 10×10 m on the horizontal plane and the vertical plain was the same as for the grid in the larger domain. The transient boundary conditions for the small domain were obtained from the solution from the larger domain (GMS 6.0 Tutorials, vol. 2).

Figure 2 Model conceptual sketch.

2.4.3 Calibration and sensitivity analysis

The model flow component was calibrated to fit the observed aquifer water levels in EK11. As a result, a set of hydraulic parameters was estimated for both the vadose zone and the aquifer (Dvory et al., 2016, 2018a). The PEST software (Doherty, 2004) was used to calibrate the transport component of the model by minimizing the least squared errors between simulated and observed concentrations of CBZ and CAF in EK11. First, the CBZ breakthrough curve was used to define longitudinal dispersivity (a_L), ratios of transverse to longitudinal dispersivities and ratios of vertical to longitudinal dispersivities (a_T∕a_L and a_Z∕a_L), exchange rate parameter (η_c), and the distribution coefficient (K_D) for CBZ sorption, assuming a zero degradation rate of CBZ. Then, the distribution coefficient for CAF sorption (K_D) and its degradation rate (λ) were found using CAF concentration measurements.

An additional estimate of the CAF degradation rate in the karst aquifer was performed analytically using both normalized CBZ and CAF concentrations. By neglecting the differences in sorption of both contaminants, the CAF degradation rate was found by minimizing the following expression:

$$\begin{array}{}\text{(4)}& \begin{array}{c}min\\ \mathit{\lambda }\end{array}\sum _{j=\mathrm{1}}^n{\left[\ln \left({\displaystyle \frac{{\stackrel{\mathrm{‾}}{C}}_{\mathrm{CBZ},j}}{{\stackrel{\mathrm{‾}}{C}}_{\mathrm{CBZ},j}}}\right)+\mathit{\lambda }{t}_j\right]}^{\mathrm{2}},\end{array}$$

where ${\stackrel{\mathrm{‾}}{C}}_{\mathrm{CAF},j}$ and ${\stackrel{\mathrm{‾}}{C}}_{\mathrm{CBZ},j}$ are relative to the observed concentrations of CAF and CBZ in the sewage at time t_j, respectively. Expression (4) reflects the relative decrease in degradable CAF concentration compared to that of non-degradable CBZ.

The sensitivity analysis was performed in order to evaluate the impact of solute transport parameters on simulated concentrations. The analysis involved varying parameters one at a time from their best fit value in the calibrated model and comparing the root mean square error (RMSE) of the obtained simulation results to the calibrated model's RMSE (Geyer et al., 2007).

3 Results and discussion

3.1 Aquifer water level and pollutant concentration fluctuation

Observations revealed that the wastewater and surface water infiltration had an immediate impact on the water level in the aquifer and on contaminant concentrations (Fig. 3). Groundwater level increased as a result of winter precipitation, runoff and four significant sewage release events. Groundwater level rose to 9.3 m during leakage events and water level decline was observed between events. Reservoir water infiltration triggered a further increase in the water level, reaching 16.6 m above the initial natural groundwater level. In the dry season there was a measured decrease in the groundwater level, and the following winter, the aquifer water level rose again as a result of precipitation and a runoff infiltration caused by Beit Zait dam overflow.

CBZ concentration variation mirrored the aquifer water level fluctuation. On 7 April 2013, CBZ concentration in sewage was measured at 995 ng L⁻¹. During the sewage release events, CBZ concentrations in the aquifer rose to 21.0–25.9 ng L⁻¹ from background levels of 3.7–7.0 ng L⁻¹ in EK11 (Fig. 3). Later, due to releases from the Beit Zait dam into the creek, CBZ concentration rose to 47.0 ng L⁻¹. During the dry season CBZ concentration returned to background levels and rose again, to 11.0 ng L⁻¹, by the end of the year due to runoff infiltration.

CAF was not initially detected in EK11 (4 July 2012), but was detected later on, likely because of the subsequent sewage spillages on 15 February and 2 April 2013. Prior to the dam opening on 15 February and 23 April 2013 CAF concentrations in the reservoir were 100 and 245 ng L⁻¹, respectively. In the same period (7 April 2013), CAF concentration in sewage was 58 500 ng L⁻¹. Three main CAF concentration peaks were observed in groundwater samples during the sewage discharge event (100.0, 300.0 and 240.0 ng L⁻¹, respectively). The first peak was related to sewage discharge and the last two CAF peaks were connected to the controlled releases from the Beit Zait dam. CAF concentration level declined to 40.0 ng L⁻¹ between the sewage leakage episodes and the reservoir water discharge event. In addition, the last two CAF concentration peaks were higher than the first, in both scale and duration. This indicates that the reservoir water effectively flushed pollution from the thick vadose zone. In the dry season, CAF concentration declined to 0.01 ng L⁻¹ as detected on 26 May 2013 and was not detected afterward.

Figure 3 Time series data observation and calculation (following Dvory et al., 2018a). (a) Tzuba meteorological station daily precipitation rate; (b) dam runoff flow; (c) sewage surface flow; (d) measured and simulated aquifer water level at EK11. (e) measured CBZ and CAF concentration.

3.2 Modeling outcome

Simulated and measured concentrations of CBZ and CAF in EK11, actual and relative, are shown in Fig. 4a and b, respectively. Fair agreement between observed and simulated concentration was obtained. The model correctly predicted the major peaks of both CBZ and CAF concentration, with RMSEs of 7.4 and 53.9 ng L⁻¹, respectively. The larger RMSE of CAF (compared to that of CBZ) is due to higher CAF concentration oscillation. The model was not able to mimic a few significant fluctuations in CAF concentration. Despite the fact that CAF concentration (58 500 ng L⁻¹) in discharged sewage was 59 times higher than CBZ (995 ng L⁻¹), the observed maximum concentration of CAF was only 6.4 times greater (Fig. 4a). Thus, there is a more significant reduction of CAF's relative concentration (Fig. 4b). Assuming no significant differences between conduit and matrix solute exchange, a drastic reduction in relative CAF concentration can be attributed to sorption and degradation.

CBZ and CAF calibrated transport model parameters are presented in Table 2. These values, calculated in this study, represent combined vadose zone–groundwater model characteristics. Even though the presence of an air phase can influence the physico-chemical processes of contaminant transport and transformation, given the quality of the dataset available (breakthrough curves in one observation well) we can only obtain lumped parameters for both the unsaturated and saturated zones. The effect of variable water saturation on pollution dispersion and degradation is accounted for by multiplying these parameters by the water content (Eqs. 1 and 2). A longitudinal dispersivity value of a_L=6.44 m lines up with its expected field-scale order of magnitude (Neuman, 1990). As discussed in Dvory et al. (2018a), the transverse and vertical-to-longitudinal dispersivity ratios cannot be efficiently found under the considered conditions. CBZ and CAF distribution coefficients (K_D) values were 0.011 and 0.1 L kg⁻¹, respectively. This indicates a higher sorption of CAF than CBZ. In laboratory batch studies with riverbed sediments, Lin et al. (2010) showed that the sorption of CAF better fit the Freundlich isotherm, and its sorption was highest in comparison to other pharmaceutical compounds studied. A batch test by Martínez-Hernández et al. (2016) revealed that sorption played a key role during the first 48 h of contact with the soil, and gave way to biodegradation afterwards, with the fastest initial sorption velocities of CAF. Nevertheless, our results indicate that using linear sorption in simulating transport of micropollutants on a large field scale can also be appropriate.

Table 2 Transport model parameters.

¹ Calculated using molecular weight of CBZ (Schwarzenbach et al., 1993). ² Dvory et al. (2018a).

The degradation rate of CAF, found by inverse solution, was 0.091 d⁻¹. The value of this parameter, which was estimated analytically using expression (4), was 0.082 d⁻¹, which is a bit smaller than the result found by inverse solution. However, the analytical procedure does not account for the effect of CAF sorption on concentration. The porous aquifer biodegradation rate for CAF was evaluated by Swartz et al. (2006) as 0.014–0.07 d⁻¹. Bradley et al. (2007) showed that the rate of CAF biotransformation in stream water is sensitive to in situ redox conditions (0.72–3.14 d⁻¹) and Martínez-Hernández et al. (2016) calculated a similar value range in laboratory batch experiments (0.37–4.18 d⁻¹). The role of different redox conditions and biodegradable dissolved organic carbon in CAF biodegradation was examined in soil column experiments by Regnery et al. (2015). They estimated that CAF biodegradation rate was 0.37 d⁻¹ for oxic redox conditions when biodegradable dissolved organic carbon was present. Hillebrand et al. (2012b, 2015) evaluated CAF biodegradation rates between 0.16 and 0.19 d⁻¹ based on tracer experiments in karst conduits, which is similar within an order of magnitude to the value obtained in this study. The above studies showed that CAF degradation rate is expected to be higher in oxic carbonate aquifer rapid flow systems.

Figure 4 (a) Observed and simulated BTCs of CBZ and CAF in EK11; (b) relative concentration variations of CBZ and CAF in EK11 (CBZ data from Dvory et al., 2018a).

3.3 Sensitivity analysis

The sensitivity of the modeled CBZ concentrations with respect to the parameters of dispersivity and the solute exchange rate between conduits and matrix was discussed in Dvory et al. (2018a). This analysis is also valid for CAF. The sensitivity of the model to the CAF degradation rate (λ) and the distribution coefficient (K_D) was assessed by changing each parameter by 10 % and 25 % of their calibrated values (Fig. 5). The model was found to be sensitive to variations in both parameters. A decrease in λ by 10 % and 25 % increases peak CAF concentration by around 18 % and 54 %, respectively. Similar increases in λ decrease the maximum concentration by 15 % and 34 %, respectively. Increasing or decreasing the λ fitted value results in a narrower or wider breakthrough curve, respectively (Fig. 5a).

One of the model assumptions was to prescribe the same pollutant degradation rate in both conduits and the matrix. In order to evaluate the effect of CAF degradation in the matrix on the breakthrough curve, we performed an additional simulation by assigning λ=0 in the matrix (i.e., CAF degradation in the matrix was neglected) while keeping its calibrated value in conduits. The result shows that the simulated breakthrough is similar to the calibrated one, except that the former has a non-zero concentration (around 9 ng L⁻¹) tail which is caused by solute exchange between the matrix and conduits. Observations show that only a few months after sewage was discharged into the creek, CAF concentration was no longer detected in EK11; therefore, we hypothesize that degradation of CAF occurs in both the matrix and conduits. However, based on the existing dataset, we cannot evaluate the differences between them.

A similar sensitivity trend is observed with respect to the distribution coefficient K_D. A decrease in K_D by 10 % and 25 % leads to consequent increase in peak concentration by 16 % and 43 %, respectively. An increase in K_D by 10 % and 25 % leads to a 14 % and 31 % decrease, respectively. This also causes a delay in peak concentration arrival time by around 1 and 3 days, respectively, compared with the calibrated case.

Figure 5 Simulated CAF sensitivity to parameters changes (a) the degradation rate and (b) the distribution coefficient. The insets show the effect of parameters on RMSE.

3.4 CAF vs. CBZ as indicators of sewage contamination in a carbonate aquifer

Both CBZ and CAF can serve as tracers to monitor sewage contamination in carbonate aquifers. Each of these contaminants has specific advantages as a pollution indicator. CBZ is more stable than CAF, and this characteristic enables CBZ to migrate far from the pollution source into the aquifer. Thus, CBZ can be used to estimate the reach of sewage pollution in the aquifer. However, there are often low background levels of CBZ, which can make it hard to distinguish a new pollution event (Hillebrand et al., 2012a). The advantage of CAF as an indicator is its higher concentration levels proximate (in time and location) to the event. The simulation results presented in Fig. 4b elucidate the difference between CBZ and CAF transport and attenuation in the aquifer. Several weeks after sewage was discharged CAF was no longer detected in groundwater near the contaminant source. Dissimilarly from the CAF breakthrough curve, which shows a rapid decline in concentration, the CBZ breakthrough curve exhibits a long tail of low concentration which can increase during the rainy season due to enhanced exchange between porous matrix and conduits. The tail of the low CBZ concentration during the dry season is a result of low saturation in the vadose zone. This reduces the hydraulic conductivity and the exchange between matrix and conduits, resulting in low CBZ transport rates toward the aquifer. In line with this, the simulation results demonstrated that by the end of the year, a relatively high volume of CBZ (around 95 % of that discharged) remained in the vadose zone matrix close to the source, and tens of meters downgradient in the groundwater, while all the CAF was degraded soon after sewage discharge stopped.

In the present study, the half-life of CAF was calculated to be 7.6 days. This result was similar, within an order of magnitude, however almost twice as long as the half-life estimated by Hillebrand et al. (2012b, 2015) which ranged between 3.7 and 4.3 days. However, the present study accounts for simulation of flow and transport in a thick unsaturated zone in both the conduits and matrix. Other conditions affecting CAF attenuation were also different.

4 Conclusion

A quasi-3-D dual permeability mathematical model was used to simulate CAF transport and attenuation in the vadose and saturated zones of the studied carbonate aquifer. The CAF sorption and degradation rate parameters were found through calibration with monitoring data from a single well. Simulation results from the model were compared to previous simulations of CBZ transport on the site (Dvory et al., 2018a). Both pollutants were found to be effective tracers for characterization and assessment of aquifer contamination by wastewater. CBZ was found to be more suitable for assessing the contamination boundaries, while CAF can be used as a contaminant tracer only shortly after a contamination event occurs. In addition, CAF's relatively ephemeral nature could present an advantage when trying to detect new contamination events in aquifers that have background CBZ contamination.

The estimated degradation rate and the distribution coefficient of CAF were 0.091 d⁻¹ and 0.1 L kg⁻¹, respectively, which may explain its high attenuation potential in the studied aquifer. The calculated CAF degradation rate was slightly smaller compared with previous studies of CAF decay in karst aquifers. The sensitivity analysis results revealed that the model is highly sensitive to the CAF degradation rate and the distribution coefficient of linear sorption.

The simulations reveal that by the end of the year, about 95 % of CBZ mass remained in the vadose zone porous matrix near the pollution source and tens of meters downstream in the groundwater, while all the CAF was degraded a few months after leakage stopped.