1 Introduction

One of the most common hydrocarbons contamination detected over the years in urban environments is the leakage of gasoline from underground storage tanks, mainly located at gas stations [1]. Actions responsible for soil and groundwater contamination in urban areas can be large. However, oil leaks from gas station are a common fact in cities around the world [2]. Leaks produced over the years from underground hydrocarbon storage system (UHSS) and associated pipes generated a gradual accumulation of hydrocarbons in the unsaturated zone (UZ). Fuel can be mobilized as a non-aqueous liquid phase (NALP) by the action of gravity and capillary forces when the accumulated quantity is enough and the porous medium allows it, until they reach the capillary fringe [3, 4]. However, as the NALP progresses through the UZ, part of the hydrocarbon is retained in the soil media as a residual phase. NALP in the form of immobile residual contamination, or pools of mobile or potentially mobile NALP, can represent continuing source of groundwater contamination [5]. The immobile hydrocarbon, trapped in the pores, keeps its competition with soil water. Interaction of retained hydrocarbon with water generates the partial dissolution of its constituents, mainly a group of aromatic organic chemical compounds called BTEX (benzene, toluene, ethylbenzene and xylene) [6]. The dissolved phase is highly mobile, and it is responsible for the transport of contaminant over long distances from the source [7, 8].

Dissolution mass transfer from NALP and immobilized hydrocarbons is a very complex process and is controlled by many factors such as saturation and distribution of fuel, heterogeneity and geochemical properties of porous media, soil water content and fluid flow [9, 10]. This manuscript studies the fate and transport of BTEX dissolved from the residual phase and the role of several factors such as sediment heterogeneity, water table depth, water flux and biodegradation rate. The focus is to study the role of UZ properties on the transport of BTEX that have been retained in the site for the last five years. Numerical simulations of several field scenarios, given the combination of the factors described above, were used to assess the risk of groundwater contamination in each case. In addition, a sensitivity analysis was carried out to acknowledge the influence of input parameters over model results and risk as proposed by Saltelli [11]. The input parameters considered were: fraction of organic content of carbon in the soil, permeability, water table depth and biodegradation rate.

Risk analysis was carried out considering theoretical data from three lithologies that can be extrapolated to sites around the world. Lithologies chosen were silt loam, loamy sand and clay. It is well known that contaminants adsorption is largely regulated by the particle size distribution of the soil. As proposed by Capri et al. [12], the finest particle sizes have the largest surface area and, therefore, a high capacity to retain organic and inorganic compounds. Otherwise, sandy sediments have lower retention capacities. Under this conditions, chemical compounds can travel further in the UZ. The model also considers different water table depths (1.5, 2.5, 5 and 8 m below ground surface). These depths were chosen to represent phreatic aquifer conditions in sites contaminated by hydrocarbons. Many works evidenced the importance of considering the fluctuations of the water table [13,14,15,16,17,18]. When the upward movement of the water level occurs, part of the NALP follows the rise of the water table, and another part is trapped below it, due to a capillary hysteresis that reduces its mobility. In the opposite direction, when the water level drops, the water drains from the porous medium and the NALP agglutinate, increasing their saturation and mobility [19,20,21,22,23]. However, NALP presence as an independent phase was not fixed in the model, so groundwater level was considered steady.

The scenarios resulting from the combination of subsoil textures and water table depths were doubled, to analyse maximum and null biodegradation in transport. The real value of biodegradation rate of BTEX compounds is difficult to obtain so two extreme cases were considered. Aerobic and anaerobic hydrocarbon biodegradation rates have been determined in many studies using laboratory microcosms, in situ field microcosms, and observed hydrocarbon concentration distributions at sites undergoing natural attenuation [24]. The rate of dissolution from a non-aqueous phase has been studied in laboratory experiments [25] and in numerical studies at the pore scale [26, 27] and field scale [28,29,30].

Petrol is composed of various hydrocarbons, the group of monoaromatics BTEX makes up 18% of these fuels approximately. These compounds are found in gasoline, constituting a large part of its soluble fraction. BTEX are the most soluble part of gasoline, thus environmental problems associated with gasoline directly involves the study of these substances [31]. BTEX percentages found in gasolines may vary depending on the country and refining methods. Generally, is considered that a 2nd grade gasoline is composed of < 1.2% benzene, 1–5% ethylbenzene, 5–10% toluene and total xylene. Leusch and Bartkow [32] informed the volume of each BTEX in the gasoline as 1% benzene, 5% toluene, 1% ethylbenzene and 7% xylenes.

The novelty of this study is the assessment of BTEX transport from the residual phase in the UZ. As mentioned above, hydrocarbon leaks from UHSS are common, but NALP need time to move as an individual phase and be detected. If the contamination is already present and remediation is carried out, the residual phase is the most difficult to remove from the soil, as it represents a continuous source of contaminant. If the hydrocarbons are not yet detected, the model can be used as a method to predict the behavior of different organic compounds in the soil over time. In recent decades, many transport models have been developed to simulate specific environmental processes under certain environmental conditions, on a laboratory and field scale. This includes modeling vertical transport of contaminants in the UZ [33,34,35,36,37], assessment of migration risks to aquifers [38, 39] and transport in aquifers [37, 40]. Valsala and Govindarajan [41] carried out a numerical investigation into the interaction of various BTEX transport processes in groundwater. Farahani and Mahmoudi [42] used HYDRUS software to solve environmental problems associated with hydrocarbons contamination. Most of the existing works on the subject simulate BTEX transport through the direct and total dissolution from the pollutant source or NALP. However, a numerical model for BTEX transport derived from residual fuel in unsaturated conditions have not been applied yet.

The aim of this work is to assess the risk of groundwater contamination associated with BTEX dissolution from fuels as a residual phase. The application of the methodology will allows to verify the use of the HYDRUS 2D/3D software in the simulation of BTEX dissolved from the residual phase.

1.1 Study area

The city of Bahía Blanca is located in the southwestern area of Buenos Aires, Argentina, approximately at 38°43′′ South Latitude and 62°16′′ West Longitude. The general landscape of the city is comprising three main geomorphological units: plateau, colluvium-alluvium and tidal plain [43]. Figure 1 shows geomorphological units and the location of monitoring wells. These wells are located in gas stations that have reported leaks from their UHSS in the past. Behalf, remediation task have been applied to remove the NALP, groundwater still shows concentrations of BTEX. Thus, these gas stations are a good example to evaluate the mobility of BTEX dissolved from residual fuel phase.

Fig.1
figure 1

Ubication map with geomorphological zones and location of monitoring wells

Full size image

The plateau is composed of silt loam sediments enriched in calcium carbonate corresponding to the Pampeano Formation. The alluvial cones originated by fluvial action have great extension, and their coalescence makes them look topographically as alluvial plains. The colluvium-alluvium complex is made up of fine sand with an abundant clayey matrix. On top of this, medium sized sands with calcareous cement were deposited and together they form the Bahía Blanca Formation. Below the 20 m level, the alluvial plain of the Napostá stream is present and continues toward the southeast with an old tidal plain, now emerged mainly composed of clayey material from the Maldonado Formation [43].

In the plateau area, residential land use predominates, while the greatest urban settlement is on the alluvial plain. The port, industrial and residential activities of Ingeniero White are located on tidal plain [44].

The phreatic aquifer of the city is characterized by different hydrogeological conditions due to the important formational variations of the subsoil. In general, it is composed by medium to fine calcareous sands, friable and without stratification that alternate with coarse sands and gravels, corresponding to the Bahía Blanca Formation. Hydrologically, this unit behaves like a shallow aquifer with good permeability and specific variable flows between 0.5 and 5 m3/h per meter of depression. From the hydrochemical point of view, the water from the phreatic aquifer is classified as sodium chloride to sodium chloride sulfate [45]. Soluble salts are freely mobilized in the unsaturated profile, mostly linked to upward vertical water flows, where the position of the water table exerts an effective control on the hydrodynamics of the unsaturated system [46]. Depth of the water table in Bahía Blanca is between 3 and 5 m according to the season and the annual alternation of rainfall cycles [47].

2 Model description

Theoretical model was divided into three steps (Fig. 2). Initially, according to real data from fuel leaks, widely known formulas were applied to calculate BTEX concentrations in the dissolve, volatile and solid phase. Next, the HYDRUS 2D/3D [48] software was used to execute sixty simulations and model BTEX transport in each of these scenarios. Finally, the maximum concentrations leached to the aquifer in each scenario were highlighted and compared with an international guide to calculate the risk factor.

Fig. 2
figure 2

Theoretical model of work development

Full size image

2.1 Step 1: theoretical framework

The spatial arrangement of BTEX in the subsoil, can be expressed by the places that they can occupy. If we consider the presence of hydrocarbons in the soil volume, the total porosity (n) can be expressed as:

$$n = n_{{\text{w}}} + n_{{\text{a}}} + n_{{\text{H}}}$$
(1)

where nw is volumetric water content (-), na is volumetric air content (-), and nH is residual hydrocarbons in the soil (-).

The capacity of soil to contain hydrocarbons in the dissolve, volatile and solid phases is finite. The maximum capacity can be estimated by the soil saturation concentration (Csat):

$$C_{{{\text{sat}}}} \left( {\frac{{{\text{mg}}}}{{{\text{kg}}}}} \right) = \left( {S*\frac{{n_{w} }}{\rho }} \right) + \left( {S*H^{\prime } *\frac{{n_{a} }}{\rho }} \right) + \left( {K_{d} *S} \right)$$
(2)

where S is water effective solubility for each compound (mg/L), ρ is soil bulk density (kg/L), H’ is Henry’s Law constant (-), and Kd is partitioning coefficient for soil–water solution (L/kg).

When the hydrocarbon content is kept below Csat, they will distribute between the dissolve, volatile and solid phases. Thereby, mobile NALP will be absent. Otherwise, if the saturation content is exceeded, mobile NALP will be present within the three described phases above. Thus, if hydrocarbons are present as a residual phase without the presence of a mobile NALP, maximum concentration of compounds in each soil phase can be expressed as:

$$C_{w} = S*n_{w}$$
(3)
$$C_{a} = S*H*n_{a}$$
(4)
$$C_{S} = S*K_{d}$$
(5)

where Cw is soil water concentration (mg/L), Ca is soil air concentration (mg/L), and Cs is adsorbed concentration in soil media (mg/kg).

Effective solubility (S) of an organic substance in a mixture can be calculated using Raoult’s law [49]. Estimation consists of multiplying molar fraction (MF) of the compound in the mixture by solubility of pure substance (S0).

$$S = S_{O} *{\text{MF}}$$
(6)

Henry’s Law explained volatilization process [50]. Water–air partitioning coefficient (\(K_{h}\)) is represented by:

$$K_{h} = \frac{{K^{\prime }_{h} }}{R*T}$$
(7)

where \(K'_{h}\) is Henry’s Law constant for each compound (atm.m3/mol), R is the universal constant of gases (atm.m3/mol.K−1), and T is temperature (°K).

Soil–water partitioning coefficient (Kd) can be estimated from partition coefficient in organic carbon (Koc) (ml/g) and fraction of organic content of the soil (foc). If adsorption effect of pollutant on the finest fraction of the soil is discarded [51], Kd is normally expressed by the relationship:

$$K_{d} = K_{{{\text{oc}}}} *F_{{{\text{oc}}}}$$
(8)

Koc is related to the octanol partition coefficient (Kow), which can be estimated experimentally in laboratory for each contaminant [52].

In aerobic environments, fuel biodegradation can generate large decreases in concentrations. Biodegradation processes can be explained mathematically according to a first-order decay model as:

$$C_{2} = C_{1} *e^{ - tk}$$
(9)

where C1 is initial concentration in soil (mg/kg), C2 is actual concentration in soil (mg/kg), k is biodegradation first-order decay constant (1/d), and t is time (d).

The conceptual transport model for BTEX dissolved from a residual phase is illustrated in Fig. 3. Simulation of BTEX can be carried out under the following premises:

  1. (1)

    Hydrocarbons as a residual phase is constant in time and space. This condition is feasible at gas stations, where leaks of small proportions tend to be generated over long continuous periods.

  2. (2)

    Maximum BTEX concentration in soil water would correspond to the maximum soil water retention capacity. This volumetric water content is equal to field capacity. Hydrocarbons will occupy the largest pores while water will occupy the smaller ones.

  3. (3)

    Dissolved BTEX concentration in soil water is at equilibrium with the residual phase, being equal to the saturation concentration (Csat).

  4. (4)

    BTEX are transported by water flow that moves downwards the damping depth.

  5. (5)

    Soil media is vertically and laterally homogeneous.

Fig. 3
figure 3

Conceptual model of BTEX transport in the UZ from the dissolution of residual fuel phase

Full size image

2.2 Step 2: Numerical simulation

Transport simulation was performed by HYDRUS 2D/3D software [48]. Mathematical code is an executable software in Windows domain; widely use in water flow and solute transport simulation in saturated and unsaturated soil media. The software allows resolution of modified Richards Eq. (10) for water flow and convection–dispersion Eq. (11) for solute transport. Similarly, the procedure could be applied in the freeware 1D software version.

$$\frac{\partial \theta }{{\partial t}} = \frac{\partial }{{\partial x_{i} }}\left[ {K\left( {K_{ij}^{A} \frac{\partial h}{{\partial x_{j} }} + K_{iz}^{A} } \right)} \right]$$
(10)

where θ is volumetric water content (-), h is pressure head (cm), xi (i = 1,2,..,n) are spatial coordinates, t is time (d), \(K_{iz}^{A}\) are components of a dimensionless anisotropy tensor y K is unsaturated hydraulic conductivity (cm/d).

$$\frac{\partial \theta c}{{\partial t}} + \frac{\partial \rho S}{{\partial t}} + \frac{{\partial n_{a} g}}{\partial t} = \frac{\partial }{{\partial x_{i} }}\left( {\theta D_{ij}^{w} \frac{\partial c}{{\partial x}}} \right) + \frac{\partial }{{\partial x_{i} }}\left( {n_{a} D_{ij}^{g} \frac{\partial g}{{\partial x}}} \right) - \frac{\partial qC}{{\partial x}} - k_{w} \theta C - k_{g} n_{a} - k_{s} \rho S$$
(11)

where c, S and g are solute concentrations in liquid (mg/L), solid (mg/kg) and gaseous (mg/L) phases, respectively. Q is the volumetric flux density (cm/d), kw, kg y ks are first-order rate constants for solutes in liquid, solid and gas phases (1/d). ρ is bulk density (kg/L), na is volumetric air content (-), θ is volumetric water content (-), and \(D_{ij}^{w}\) y \(D_{ij}^{g}\) are dispersion coefficient tensor (cm2/d) for the liquid phase and diffusion coefficient tensor (cm2/d) for the gas phase, respectively.

To solve Richard’s equation, the following constitutive relationships are required. Van Genuchten [53] proposed such relations as given in Eqs. (12) and (13):

$${\text{Se}} = \frac{{\theta - \theta_{r} }}{{\theta_{s} - \theta_{r} }} = 1 + \left[ {\left( {\alpha h} \right)^{n} } \right]^{ - m}$$
(12)
$$K\left( {{\text{Se}}} \right) = Ks{\text{Se}}^{l} \left[ {1 - \left( {1 - {\text{Se}}^{\frac{1}{m}} } \right)^{m} } \right]^{2}$$
(13)

where Se is effective water saturation (-), θr is residual water content (-), θs is saturated water content (-). α = ha−1 is related to the inverse of the air-entry suction (cm−1), h is tension head (cm), n and m are empirical parameters m = 1–1/n (-). Ks is saturated hydraulic conductivity (cm/d), and l is pore connectivity parameter in the hydraulic conductivity function (-), estimated to be about 0.5 as an average for many soils [54].

A 2D numerical domain was created to simulate water flow and BTEX transport in the UZ. Vertical length was considered as the distance between the lower limit of soil affected by leaks and the position of the water table. The upper boundary was defined as a time-variable flux condition, and the lower boundary was considered as a constant pressure head condition (h = 0 cm). Lateral boundaries were considered as no flux condition (Fig. 4). The initial water content was estimated as field capacity for each lithology.

Fig. 4
figure 4

Spatial dimensions and boundary condition used in simulations

Full size image

To evaluate diverse scenarios, sixty (60) simulations were performed. Different combinations of water fluxes, sedimentary materials, water table depths and the occurrence of biodegradation during transport were evaluated. The time considered in the simulation was five years since it is an appropriate period between the occurrence of fuel leak in UHSS and the consequent detection of BTEX in groundwater. BTEX dissolve phase source was set at one-meter depth from the surface in all scenarios. Affected area was considered as 1m2. Water flow was considered as the flux that moves downwards from the damping depth. Rainfall values considered are from Bahía Blanca (Argentina) for the period January 2013 to December 2018. Daily precipitation was obtained from a meteorological station of the Argentine National Meteorological Service. The average rainfall value was 580 mm/year during the period 1908–2008 [55]. Carrica and Lexow [56] determined that the average recharge in the basin varies between 7.5% and 11% of rainfall (R). Since gas stations commonly have waterproofed surfaces, it is difficult to estimate the actual recharge rate. Thus, three potential recharge rates were taken as an approximation: 0.5% R, 1% R and 3% R.

Three types of sediments from Bahía Blanca subsoil were evaluated: silt loam, loamy sand and clay (Table 1). Samples were collected in field by a helical drill for textural analyses and by core method [57] for bulk density determinations. Sediments correspond to three subsoils that differ in textural and physicochemical characteristics. The textural analysis was done with the pipetting technique [58]. Hydraulic parameters of each material were obtained from Rosseta pedotransfer function [59], based on soil texture and bulk density. Rosseta was obtained from a large number of soil hydraulic data and soil properties from three worldwide databases. Total organic content of sediments was determined using a dry combustion technique [60].

Table 1 Hydraulic and physicochemical characteristics of the sedimentary materials evaluated
Full size table

Groundwater contamination by BTEX is expected to be greater in areas where the water table is closed to the surface and minor in basins where aquifers are deep. To verify the effects of vadose zone thickness over BTEX transport, four position of the water table were considered (1.5, 2.5, 5 and 8 m below ground surface).

Longitudinal dispersivity (DL) was calculated according to the following equation [61]:

$$D_{L} = 0.83L^{2.414}$$
(14)

where L is the vertical distance in the flow domain.

Transverse dispersivity (Dt) was considered as 10% of DL. The increase of dispersivity with increasing travel distance in soils agrees with reported dispersivities measurements in groundwater tracer studies [62]. The solute transport parameters used in the model for each BTEX are present in Table 2:

Table 2 Solute transport parameters used in HYDRUS
Full size table

2.3 Step 3: groundwater contamination risk assessment

To assess the risk of groundwater contamination, maximum concentration leached to the aquifer for each BTEX were compared with the limits established by the Environmental Protection Agency–United States [69] for drinking water. The potential risk factor (Rf) of groundwater contamination can be defined as:

$$R_{f} = \frac{{C_{\max } }}{{C_{{{\text{limit}}}} }}$$
(15)

where Cmax is the maximum concentration of each BTEX leached into the aquifer during simulations, and Climit is the concentration established by the international guideline [69]. Cmax value was established according to solute flux at the bottom boundary in HYDRUS 2D/3D software. Accordingly, groundwater contamination risk is calculated for the maximum solute mass leached into the aquifer, equal to:

$${\text{BTEX}}_{{{\text{mass}} }} = C_{{{\text{max}}}} * {\text{Recharge rate}}$$
(16)

Guideline has a limit value of 0.005 mg/L, 1 mg/L, 0.7 mg/L and 10 mg/L for benzene, toluene, ethylbenzene and xylene, respectively. Risk categories were defined as: very high when Rf exceeds 1, high when it is between 1 and 0.67, moderate if value is from 0.66 to 0.34 and low between 0.33 and 0.01 of limit concentration. When Rf value was less than 0.01, risk was considered null. Risk categories were established by dividing the EPA drinking water limit for each contaminant into thirds [39].

2.4 Model validation

To validate the model, comparisons were made between simulated and observed concentrations of BTEX in groundwater. Recorded data were obtained from environmental reports from Bahía Blanca Government. BTEX concentration was monitored in wells from gas stations in the year 2013. Wells were chosen according to the lithology and the depth of the water table. The Mean Absolute Error (MAE) was used to assess model validation, according to the formula proposed by Stauffer and Seaman [70]:

$${\text{MAE}} = \mathop \sum \limits_{i = 1}^{n} \frac{{\left| {X_{i} - X_{{i {\text{sim}}}} } \right|}}{n}$$

where Xi corresponds to real values, Xi sim corresponds to simulated values, and n is the number of dates.

2.5 Sensitivity analysis

A systematic sensitivity analysis (SA) was carried out to evaluate the parameters contribution to modeled results [11]. The SA was based on several simulation run over the same temporal period by increasing and decreasing individual parameters while fixing other input parameters constant. The magnitude of change of maximum concentration of BTEX leached into the aquifer was used to determine the sensibility of output to variations in each parameter. The input parameters analyzed were biodegradation rate, soil organic carbon fraction, saturated hydraulic conductivity, and water table depth.

3 Results

3.1 Dissolved phase behavior

Distribution of BTEX in the dissolved phase varies depending on lithology as well as the UZ thickness considered. Figures 5 and 6 show the behavior in depth of dissolved contaminants for two types of sediments: silt loam and loamy sand. Results are shown for vertical transport distances of 0.50, 1.50 and 3.50 m, respectively (t = 1278 days).

Fig. 5
figure 5

Concentration profiles for BTEX dissolved phase in silt loam sediment. Simulations are shown with active biodegradation and non-active biodegradation for three vertical transport distances: a 0.50, b 1.50 and c 3.50 m below the polluting source

Full size image
Fig. 6
figure 6

Concentration profiles for BTEX dissolved phase in loamy sand sediment. Simulations are shown with active biodegradation and non-active biodegradation for three vertical transport distances: a 0.50, b 1.50 and c 3.50 m below the polluting source

Full size image

Silt loam has low saturated hydraulic conductivity (0.1 m/d) and medium to high organic content (0.76% TOC). Saturation concentrations, calculated according to step 1, are 33.5 mg/L benzene, 33.4 mg/L toluene, 95.9 mg/L ethylbenzene and 167.9 mg/L xylene. These concentrations represent the maximum capacity that soil water could dissolve, without generation of NALP as an independent phase. Figure 5 shows concentration profiles for dissolved phase in silt loam sediment. Compounds migrate to greater depths when biodegradation was not considered, increasing the risk of groundwater pollution. The most mobile substances are benzene and toluene. Ethylbenzene and xylene are mainly concentrated near the source zone. These compounds have Kd values up to 10 times greater than benzene and toluene. When biodegradation process was considered, dissolved concentrations are notably reduced. After three years simulation, BTEX did not reach phreatic aquifer.

In contrast, loamy sand sediments have high hydraulic conductivity (1.21 m/d) and very low organic content (0.072% TOC). Saturation concentrations are 6.8 mg/L benzene, 10.5 mg/L toluene, 10.9 mg/L ethylbenzene and 19 mg/L xylene. Since water retention capacity of loamy sand is low, capacity of BTEX to dissolve from fuel decrease. Figure 6 shows concentration profiles for the dissolved phase in the loamy sand sediment. When biodegradation was not considered, all compounds reached the aquifer regardless of the UZ thickness. The maximum concentrations leached were generated in shallow aquifers, decreasing with an unsaturated thickness increase. When biodegradation was considered, dissolved concentrations were reduced. Anyway, toluene continued leaching into the aquifer in considerable quantities.

The last sediment analyzed corresponds to a clayey soil, which has low hydraulic conductivity (0.14 m/d) and high organic content (1% TOC). Like silt loam, BTEX mobility was poor. Figure 7 shows the variation of BTEX concentration over time for different depths from the source zone. The analysis carried out in this lithology was different due to the low mobility of the compounds. Observation nodes were placed at 0.30 m below the source (Fig. 7a), 0.90 m below the source, depth that correspond to the upper limit of the capillary fringe (Fig. 7b) and 1.50 m below the source, corresponding to the water table (Fig. 7c). Sediment saturation concentrations are: 42 mg/L benzene, 45 mg/L toluene, 110 mg/L ethylbenzene and 193 mg/L xylene.

Fig. 7
figure 7

source: a 30 cm, b 90 cm: upper limit of the capillary fringe c 150 cm: water table position

Dissolved BTEX concentration vs time in clay sediments. Observation nodes located below

Full size image

Only toluene reached the water table, when maximum biodegradation was considered. However, concentrations leached are very low. When biodegradation was not considered, ethylbenzene and xylene did not reach the water table, being adsorbed close to the source. Benzene and toluene leached into the aquifer with concentrations of 7 ug/L and 77 ug/L, respectively.

3.2 Solid phase behavior

Behavior of BTEX was also analyzed in the solid phase for three lithologies. Figure 8 shows adsorbed concentrations for two sub-sections of UZ for the five years simulation. The sub-sections considered correspond to depths from 0 to 0.30 cm and 0.31 to 0.60 m below the source. To assess the influence of unsaturated thickness, two water table positions were considered, 1.5 m (Fig. 8a) and 4 m (Fig. 8b) below ground surface. Results corresponded to scenarios when water flux was equal to 3% Rainfall. BTEX adsorption is similar for all lithologies (Fig. 8). Compounds are mainly concentrated in the first 30 cm from the source, equally for both water table conditions. The affinity for the solid phase is: xylene > ethylbenzene > toluene > benzene. Order of magnitude of concentrations adsorbed in sands is up to 100 times lesser than fine-grained sediments. Besides, when water table was shallow, adsorbed concentrations are minimal to null at depths greater than 30 cm (Fig. 8a); but when water table was located at greater depths (Fig. 8b), adsorbed concentrations increase in the lower sub-section.

Fig. 8
figure 8

Simulated concentration of absorbed phase. Results are presented for the three lithologies and two water table depths a 1.5 m and b 4 m below the surface

Full size image

3.3 Contamination risk assessment

Risk assessment will allow recognize the factors that enhance groundwater protection. Also, it is essential to evaluate the effects of biodegradation process during BTEX transport through the unsaturated medium. The Rf values for five years simulation with non-active biodegradation are listed in Table 3.

Table 3 Risk factor for five years simulations without considering biodegradation. Risk classification (Rf) > 1 very high (red), 1–0.67 high (orange), 0.66–0.34 medium (light orange), 0.33–0.01 low (yellow), < 0.01 null (green)
Full size table

In silt loam sediments, when leak is generated near the water table, risk of benzene contamination is very high. Maximum leachate concentrations are approximately 13, 28 and 117 orders of magnitude greater than the limit established by EPA, considering a water flux of 0.5% R, 1% R and 3% R, respectively. However, as vertical transport distance increased, mass of BTEX that reached groundwater decrease. When the water table was placed at 2.5 m depth and the highest rate of water flux was considered (3% R), risk is cataloged as low.

Loamy sand sediments showed a different behavior. In this case, Rf for benzene is very high in all scenarios considered. For the five years simulation, benzene leachate concentrations greatly exceed the limit established by EPA. This condition happened even when the hypothetically water table was 8 m deep. Toluene, have a very high risk only when a shallow water table and the highest water flux were considered. Toluene risk decreased to moderate or low in opposite cases. Similarly, ethylbenzene risk was estimated from moderate to low. For deep water tables, risk is present only when the highest water flux was considered. Xylene is the compound with the most limited mobility, and Rf is always null to low.

Clay sediments have the highest sorption capacity. For this reason, only shallow aquifers (1.5–2.5 m) could be cataloged at high risk of contamination when long-term leakages are present.

When degradation was considered, results differ greatly from the previous ones. In silt loam and clay soils, risk became null for all the scenarios simulated. These results were obtained even for shallow aquifers with the highest recharge rate. Table 4 listed Rf values calculated for 5 year simulation in loamy sand sediments.

Table 4 Risk Factor for five years simulation considering degradation in loamy sand material. Risk classification (Rf) > 1 very high, 1–0.67 high, 0.66–0.34 medium, 0.33–0.01 low, < 0.01 null
Full size table

Maximum concentration of benzene leached into the aquifer was 14.95 μg/L. This concentration is equal as an Rf of 2.99, which stand for very high risk of contamination. The worst scenario is generated when combining the shortest vertical transport length and the highest water flux rate. For the same water flux (3% R), calculated risk was moderate when the water table depth was between 2.5 and 4 m deep. When water table was 8 m deep, risk was low to null. Under degradation scenarios, the risk calculated for toluene ranged from low to null. Similarly, xylene and ethylbenzene did not represent a risk in any case.

Since benzene is the most dangerous compound for human health, concentrations that would reach to the aquifer were analyzed. Figure 9 shows concentration of benzene leaching into the aquifer for three water flux rates: 0.5% R (Fig. 9a), 1% R (Fig. 9b) and 3% R (Fig. 9c). Results correspond to loamy sand with active biodegradation. Under this condition, it is expected to have the greatest attenuation. The blue dotted line indicates the maximum concentration allowed by EPA [61] for drinking water. When the lowest water flux (0.5% R) was considered, concentrations are below the prohibition limit. However, when water flux increased to 3% R, leached concentrations increased too. For the shallowest water table, benzene concentration leachate reached as far as 15 ug/L. When increasing the water table depth, leached concentrations are negligent.

Fig. 9
figure 9

source and water table

Benzene concentration leaching to the phreatic aquifer for the 2014–2018 simulated period. Results correspond to water fluxes of: a 0.5% R b 1% R c 3% R. Color lines indicate vertical distance between pollutant

Full size image

3.4 Model validation

Model validation was done by comparing simulated concentrations of BTEX leachate in the aquifer with existing information from monitoring wells at different gas stations in Bahia Blanca city. Three gas stations were chosen (Fig. 1), due to their similarity in subsoil and hydrologic conditions with the simulated scenarios. The PP site presents a shallow water table (1.3 m deep), and the lithology is clay. The other two sites, named CH and AM, present a sandy subsoil and water tables at 2.6 m and 4.21 m deep, respectively. Figure 10 shows simulated and observed concentrations of BTEX in the aquifer. The simulated values corresponded to leached concentrations into groundwater, according to three recharge values considered by the model (0.5, 1 and 3% R).

Fig.10
figure 10

Simulated leachate for three position of water table (1.5, 2.5 and 4 m) and three recharge flux (0.5, 1 and 3%R) was compared with real values recorded from three gas stations with similar properties (PP, CH and AM)

Full size image

The highest benzene simulated leachate was recorded when phreatic level was at 1.5 m. The predicted value for the scenario of highest water flux (3% R) (0.829 mg/L) duplicated the highest observed value in the aquifer (0.434 mg/L). However, the other simulated concentrations (0.232 and 0.016 mg/L for 1 and 0.5% R, respectively) were close to the observed values (0.341 and 0.050 mg/L). When water table was at 2.5 m, simulated concentrations for the highest recharge (0.220 mg/L–3% R) was higher than the highest observed value recorded (0.084 mg/L). The other simulated values (0.072 and 0.035 mg/L for 1 and 0.5% R, respectively) were close to observed values (0.021 and 0.022 mg/L). When groundwater level was at 4 m, the measured concentrations in the aquifer (0.160, 0.130 and 0.060 mg/L) and simulated values (0.171, 0.058 and 0.029 mg/L) were similar for all recharge rates. MAE for benzene was calculated for three position of water table, being higher when water table was shallower (17.9%), intermediate when phreatic level was at 2.5 m (6.7%) and low for the deepest water table position (3.80%).

Toluene has the greatest differences between simulated leachate and real concentrations. Simulated values are almost double that observed concentrations when phreatic level was at 1.5 m. Predicted values were 1, 0.702 and 0.07 mg/L for the scenario with 3, 1 and 0.5% R recharge rate, respectively; while recorded values were 0.70, 0.33 and 0.30 mg/L. When water table was at 2.5 m, differences between values are smaller. Simulated values were 0.428, 0.140 and 0.069 mg/L for the different scenarios with recharge rates 3, 1 and 0.5% R; while real values were 0.036, 0.020 and 0.016 mg/L. Finally, when the water table was at 4 m, values recorded (0.390, 0.2 and 0.05 mg/L) are even more similar to simulated values (0.339, 0.116 and 0.058 mg/L) for the three recharge rates. MAE for measurements made at 1.5 m was 59%, while for measurements made at 2.5 and 4 m was 11.8 and 4.8%, respectively.

Ethylbenzene showed a good overall fit between simulated leachate and real data. When water table was positioned at 1.5 m, simulated values are close to zero for all the recharge values; however, recorded values are higher being 0.184, 0.084 and 0.052 mg/L. When water table was at 2.5 and 4 m, simulated and real values were close. In the first case, predicted values were 0.052, 0.016 and 0.008 mg/L for the recharge values of 3, 1 and 0.5% R, respectively, and measured values were 0.026, 0.017 and 0.008 mg/L. When water table was at 4 m, simulated values were 0.041, 0.013 and 0.007 mg/L for the recharge values of 3, 1 and 0.5% R, respectively, and recorded values were 0.060; 0.020 and 0.010 mg/L. MAE for ethylbenzene was 10.70% for the most superficial water table, 0.90% and 1% for the deepest positions (2.5 and 4 m, respectively).

Xylene compound also showed a good overall fit between predicted and real data. When water table was near to the surface (1.5 m), simulated values were close to zero for all recharge values. However, observed concentrations were 0.168, 0.080 and 0.050 mg/L. For the deepest levels (2.5 and 4 m), simulated values were close to the detected values. The highest simulated values were 0.001, 0.123 and 0.091 mg/L for 3% R recharge value, for 1% R recharge value were 0.037 and 0.030 mg/L and for the lowest percentage of recharge were 0.018 and 0.015 mg/L. Real values for water table at 2.5 m were 0.184, 0.029 and 0.016 mg/L, and for the deeper water table were 0.1, 0.08 and 0.02 mg/L. MAE for xylene was 9.90% when water table was near to the surface, 2.40% and 2.10 when water table was deeper (2.5 and 4 m, respectively).

3.5 Sensitivity analysis

SA was carried out to recognize the most important variables in BTEX fate and risk calculation. Parameters considered for this analysis were: fraction of organic carbon (Foc) (0,07% and 1%), which directly affects compounds Kd, saturated hydraulic conductivity (10 and 100 cm/day), water level depth (2.5 and 4 m) and active or null biodegradation. Values considered in each case were chosen according to extreme values recorded during the simulation of the previous scenarios and were listed in Table 5.

Table 5 Parameters considered in the sensitive analysis and % of BTEX leachate reduction
Full size table

To run the simulations, the lithology was considered as sand, and the water flux was equal to 3% of rainfall. The comparisons were made based on the maximum simulated concentration reaching the water table. Table shows the parameters that were considered in the sensitivity analysis and percentage of reduction of BTEX leached concentration in each case. The maximum concentrations leached into groundwater were considered as 100% in case 2 for the four compounds.

The greatest reduction in BTEX compounds was generated by biodegradation (Table 5). Comparing case 2 with case 6, only difference between them is maximum or null biodegradation, toluene reduction was 94%, benzene 99% and 100% for ethylbenzene and xylene. Parameter Foc affects directly Kd values of compounds and generated an important reduction too. Comparing case 2 with case 4 (Table 5), toluene reduction was 68%, benzene 85% and 99.5% for ethylbenzene and xylene. However, when phreatic level was at 4 m, variation of Foc produced a lower compounds reduction (case 10 and 12), 37% for toluene, 40% for benzene and 12% and 14% for ethylbenzene and xylene. The depth of the water table is a parameter that greatly varies, in this case two extremes were considered (Table 5). When comparing case 2 with case 10, toluene reduction was 55.5%, 58% for benzene, 88% and 86% for ethylbenzene and xylene, respectively. The last parameter considered was the saturated hydraulic conductivity (Ks). This parameter showed a minor important in model sensitivity. If we compare case 1 with case 2, the greatest reduction is presented by toluene (3%), followed by benzene (2.7%) and finally, ethylbenzene and xylene with lower values (0.2 and 0.4%, respectively). Ks was directly related to the variation of Foc. When Foc was higher, K variation produced a minor reduction of leachate. Comparing cases 3 and 4, leachate reduction percentages were always lower than 1%.

4 Discussion

The behavior of BTEX dissolved from the residual fuel soil phase has been investigated using the software HYDRUS 2D/3D. Different theoretical scenarios were simulated to consider lithological and hydrological variability and to evaluate the conditions that enhance or decrease the mobility of BTEX in the UZ. The method presented in this work can be applied to any type of soil and/or aquifer in a simple way. The simulations are based on soil hydraulic properties for the application of a numerical model. Estimation can be made based on soil physical properties if laboratory determinations for the hydraulic parameters are not available, e.g., pedotransfer functions. Although HYDRUS 2D/3D can represent multi-layered systems, pure lithologies were considered as allows a simplification of the complex reality of Bahía Blanca city subsoil. Anyway, vertical variations in the soil profile can be incorporated into the model if necessary.

Simulations were carried out considering different: lithologies (silt loam, sandy loam and clay), water fluxes (0.5%, 1% or 3% Rainfall), water table depths (1.5 m, 2.5 m, 5 m or 8 m below ground surface) and biodegradation rate (active or null). The model was validated by comparing observed concentration of BTEX and simulated concentrations leached into the aquifer. This software is a powerful tool for modeling the behavior of BTEX compounds in the UZ when the pollution source was at a safe distance from the aquifer (1.5 m). This was demonstrated by the low MAE values (≤ 11%) obtained for the deepest water table positions. These error values are consistent with the relative error reported by Farahani and Mahmoudi [42]. They used HYDRUS 1D to model the changes of PAHs concentration in soil columns and verify its correlation with field data. They reported a good conformity between experimental and laboratory data with relative errors less than 10%. Gupta and Yadav [70] also reported a good agreement between simulated data and laboratory results. They had done laboratory transport test in a 3D sand tank to assess LNAPL behavior.

The mobility of each BTEX compound is closely related to its physicochemical characteristics and properties of each sediment. However, vertical transport velocity was remarkably reduced when BTEX reached the capillary fringe. Water diffusivity of BTEX is several orders lower than air diffusivity. Thus, BTEX mobility notable reduced under humidity conditions close to saturation. For the three lithologies analyzed, benzene and toluene were the compounds with the highest mobility in the dissolved phase. Similarly, Shores et al. [39] established that benzene and toluene were the compounds with greatest mobility when studied BTEX transport for two sandy loam and clay loam. Furthermore, they determine that these compounds represent the greatest risk to human health and the environment, because their limit for drinking water established by the EPA is the lowest. Likewise, same results were obtained when simulating BTEX transport in silt sediments for a fifty years’ period [31]. Toluene demonstrated to be the most mobile BTEX in the dissolved phase. However, calculated risk for toluene varies from low to null depending on the scenario. This is mainly because concentration limit for drinking water is 1000 μg/L, and risk is present only when great amounts of toluene are transported into the aquifer. Despite this, benzene is the most dangerous compound for health and the environment from BTEX group. Benzene has greater mobility and toxicity in relation to others BTEX, in addition to being a carcinogenic compound for humans [71]. Its greater mobility is given by a greater effective solubility and a lower sorption potential, so benzene tends to generate more extensive dissolving plumes in relation to other hydrocarbons [72, 73]. In fine-grained sediments when biodegradation was not considered, benzene concentrations leached into the aquifer represented a very high risk for shallow aquifers (> 1.5 m depth) and moderate to low risk for deeper aquifers. In loamy sand sediment, risk is classified as very high for aquifers up to 8 m deep.

Xylene and ethylbenzene tended to concentrate close to the source zone. These two compounds predominantly concentrated in the solid phase. Similar results were obtained by Shores et al. [39] in numerical simulations with HYDRUS. Furthermore, according to authors ethylbenzene and xylene leaching never reached sufficiently high concentrations to be consider as a risk, when analyzing oil surface spills. Soil type and its composition play an important role for pollutants retention. In general, coarse-grained soils exhibit lower tendency for contaminants adsorption than fine-grained soils [74]. In this work, it was demonstrated that, in general, concentrations of ethylbenzene and xylene leaching to the aquifer are low or null, representing a low to null risk to human health.

Leached concentrations for different scenarios allowed calculating the risk of groundwater contamination. Risk was proportional to water flux rate and inversely proportional to the unsaturated thickness. Water percolation through soil affected the surface reactions of soil particles and the attraction forces between water and other fluids in the soil matrix, resulting in enhanced drainage or leaching of contaminants from the soil [75]. As well, Thomé et al. [76] concluded that rainfall infiltration conditions can decrease the retention and biodegradation of contaminants and allow a greater leaching. In silt loam and clay sediments, characterized by low hydraulic conductivities and moderate to high organic content, BTEX mobility was very limited. Otherwise, in loamy sand sediments BTEX concentrations leached into the aquifer can be large. Transit time of dissolved compounds is greatly reduced in this soil, due rapid percolation of water and low sorption capacity of contaminants.

The SA verified that biodegradation and sorption processes highly affect BTEX transport. Same results were exposed by Valsala and Govindarajan [41] for highly soluble BTEX constituents. Biodegradation generates the highest reduction of BTEX compounds. In fine grain soils, risk became null for all BTEX, even for shallow aquifers and the highest water flux rate. In loamy sand sediments, risk was very high only for the worst scenario considered (shallow water tables and high recharge rate). When the water table depth was increased, risk decreased from moderate to low, depending on the water flux rate. It is widely documented that BTEX are biodegraded by a wide variety of microorganisms [77,78,79]. However, biodegradation process could be inhibited when compounds concentrations exceed toxic level for microorganisms [17]. Thus, it is advisable to assess the risk when biodegradation is null since BTEX mass transfer into the aquifer will be much higher under these conditions.

5 Conclusions

  • The proposed theorical model proved to be a good simulation tool when pollution source was located at distances greater than 2 m from the water table. However, when the source was near to the aquifer, the error between simulated and real concentrations increased.

  • Distribution of BTEX in the dissolved phase varies depending on lithology and biodegradation as well as the UZ thickness considered. In silt loam sediment, compounds migrate to greater depths when biodegradation was null but, when biodegradation was maximum, dissolved concentrations were notable reduced. In loamy sand sediment, when biodegradation was null, all BTEX reach the aquifer regardless of the UZ thickness. In spite of dissolved concentration were reduced when biodegradation was considered, high concentration of toluene still reaches groundwater. In clay sediment, BTEX mobility was poor. Only toluene reached the aquifer with low values when biodegradation was considered, while without biodegradation benzene and toluene reached the phreatic level.

  • The BTEX affinity for solid phase is similar for all lithologies xylene > ethylbenzene > toluene > benzene. However, concentration adsorbed in sand is up to 100 times lesser than fine grain sediments.

  • For the three lithologies analyzed, benzene and toluene were the compounds with the highest mobility in the dissolved phase, and toluene demonstrated to be the most mobile BTEX. Contrary, xylene and ethylbenzene tended to concentrate close to the source zone. These two compounds predominantly concentrated in the solid phase.

  • Risk was proportional to water flux rate and inversely proportional to the unsaturated thickness. In fine-grained sediments when biodegradation was not considered, benzene concentrations leached into the aquifer represented a very high risk for shallow aquifers (< 1.5 m depth) and moderate to low risk for deeper aquifers. In loamy sand sediment, risk is classified as very high for aquifers up to 8 m deep. When biodegradation was considered, leached concentrations were greatly reduced for all BTEX. In fine grain soils, risk became null for all BTEX, even for shallow aquifers and the highest water flux rate. In loamy sand sediments, risk was very high only for the worst scenario considered (shallow water tables and high recharge rate). When the water table depth was increased, risk decreased from moderate to low, depending on the water flux rate.

  • Toluene is the aromatic compound with less fit between recorded and simulated data and ethylbenzene and xylene are the compounds with a better fit. In all cases, the highest values of the MAE are given when water table is shallow (1.5 m). However, values of MAE are acceptable (in general less than 11%) when water table is deeper (2.5 and 4 m).

  • The SA allowed to recognize the parameters most sensitive to changes such as: biodegradation > Foc > water table depth > permeability.