Hydrodynamic conceptual model of groundwater in the headwater of the Rio de Oro, Santander (Colombia) by geochemical and isotope tools

The study of groundwater in high mountain regions is an important challenge for the local management of water resources. Numerous shallow aquifers of small dimension are interconnected with surface water, providing their base flow, and representing the main source of water available during these periods (Cartwright & Morgenstern 2016; Somers et al. 2018; Ghimire et al. 2019). For this reason, it is necessary to present a conceptual hydrogeological model allowing a better knowledge of the whole water system, including the physical and geological aspects of the area and water paths through schematic diagram blocks that display the main aspects of the system (Banks et al. 2009; Sappa et al. 2012).

In high mountain systems, numerous fractured hard-rock aquifers have been identified (Abbott et al. 2000; Ghimire et al. 2019), as well as karstic aquifers, where recharge and discharge processes are related to fractures, fault zones and weathering zones of the rocky massif (Hilberg & Riepler 2016). Water chemistry and environmental tracers have been widely used to better understand groundwater systems in mountain systems (Azzaz et al. 2008; Chu et al. 2017). Hydrochemistry allows the determination of the geochemical evolution of groundwater and water–rock interaction processes (Appelo & Postma 2005). With ²H, ³H and ¹⁸O isotopic compositions, it is possible to determine the signature of the water origin, the recharge mode and discharge processes, the groundwater flows dynamics, the interconnections between aquifers (Clark & Fritz 1997), the estimation of water age (Cook & Böhlke 2000), etc.

The present study was carried out at the headwater of the Rio de Oro (north-east Colombia), framed by the La Moza micro-basin stream, in a high mountain area adjacent to the Berlin Paramo, Santander (Figure 1). The aim of the present study is to propose a conceptual model from hydrodynamic, geochemical and isotopic data identifying the flow circulation, water origin, recharge zone, and residence time, as well as their relationship with surface water streams that maintain the downstream base flow that supply water to approximately 150,000 people.

The study area is located in the Eastern Cordillera of Colombia, at the Santander Massif, 46 km from the city of Bucaramanga (Figure 1). The area is demarcated by the hydrographic micro-basin of La Moza stream and its main tributary the Picacho stream, covering an area of 9 km². They are part of the Rio de Oro basin that extends over 540 km². The topography of the area is abrupt, varying from 2,500 to 3,500 masl, with elongated and steep strips.

The area shows an annual average temperature of 6–12 °C (Ideam 2005). A bimodal precipitation regime has been established in the Santander Massif, with two rainy seasons in March–April and October–November, and two dry seasons in June–July and December–January, controlled mainly by the passage of the Intertropical Convergence Zone (ITCZ). The rainfall station (El Picacho station, years 1967–2017, 3,455 masl) reports a mean annual rainfall of 1,544 mm. The estimation of weighted average rainfall over the basin is 1,300 mm/year and the real evapotranspiration (ETR) is 700 mm/year, resulting in an available water supply of 600 mm/year with an approximate potential recharge of 340 mm/year (Cetina 2019).

The study area is located in the Santander Massif, bounded to the south-west by the Bucaramanga Fault and to the south-east by the Boconó Fault System, which generated a strong fracture inside the rocky massif (Cediel et al. 2003). In this physical environment, an igneous-metamorphic basement emerges. Quartz-feldespatic schists from the Silgará Formation (pDs), with a Lower Paleozoic metamorphism age, intruded by granites from the La Corcova Cuarzomonzonite (JTRcl) and granodiorites from the Santa Bárbara Cuarzomonzonite (Jcs) units, which report an age of 195 and 192 Ma (Ward et al. 1973). A sedimentary succession rests on these bodies. Sandy limolites and conglomerate sandstones of the Jurassic Girón Formation (Jg), quartzarenites belonging to the Valangenian–Hauterivian Los Santos Formation (Kils), limestones with levels of claystones and immature sandstones from the Hauterivian Rosablanca Formation (Kir) and quaternary deposits (Qc) (Ward et al. 1973) (Figure 1). The area is affected by dextral and thrust faults, which present two main trends, NNE and NW (Figure 1), generating high fracturing in the rocks, creating spaces where groundwater can flow.

Most of the rock units have low to no primary porosity, so the hydrogeological potential is limited to the development of secondary porosity generated by the fracturing. The secondary porosity developed in the limestones of the Rosablanca Formation is increased by carbonate dissolution processes (Hartmann et al. 2014). On the other hand, the igneous units acquired a sandy texture due to the meteorization phenomenon that also increases the porosity of the rocks (Sukhija et al. 2006).

Nine springs were identified (Figure 1). M1 and M9 springs are located between El Picacho and Las Cruces faults, to the north of the study area. The M1 spring emerges in a strong change of slope, in carbonate rocks of the Rosablanca Formation, which present cracks and sinkholes by dissolution processes, configuring in this zone a karstic system. The M9 spring also appears in a change of slope, in fractured quartzarenites of the Los Santos Formation.

M6 and M7 springs are located in the hanging block of the Santana Fault, to the east of the study area (Figure 1). The springs emerge mid-slope in weathered coarse-grained granodiorites from the Santa Barbara Cuarzomonzonite unit. M2, M3, M4, M5 and M8 springs are located south-east of the study area (Figure 1). The M2 spring emerges at the base of an escarpment (change of slope) of carbonate rocks of the Rosablanca Formation, where there are cracks and sinkholes by dissolution processes, configuring (like the M1 spring) a karstic system. The M3 spring emerges in quartzarenites from the Los Santos Formation, in a slight change of slope near to the watershed of the micro-basin. The M4 spring emerges on a level of immature fine-grained sandstones from the Rosablanca Formation. The M5 spring emerges in the middle of a moderate slope on a residual clay soil of the Los Santos Formation. The M8 spring outcrops mid-slope on carbonate rocks from the Rosablanca Formation, meters below the M2 spring, where part of its water comes from this spring.

The present study is based mainly on geochemical and isotopic data acquired during 2017–2018. Temperature, pH, EC and discharge were measured in the field. Sixteen samples from springs (seven in the rainy season and nine in the dry season) and one sample from surface water (in the dry season) were analyzed for major ion chemistry. Seventy-nine samples were measured for δ¹⁸O and δ²H, including monthly rainwater from two stations (n = 24), springs (n = 42) and surface water (n = 13). To complete the sampling, Tritium was measured on six springs during the dry season to carry out a first evaluation of water dating.

The rainwater isotopic composition was determined from composite samples taken at two rainwater collecting stations called Picacho (3,550 masl) and Trucheria El Horno (2,440 masl). The sampling of the nine springs was made during the rainy season in October 2017 and the dry season in February 2018, complemented with samplings during May 2017, March, May and June 2018. The surface water was sampled along the main La Moza and Picacho streams, in February and August 2018, complemented with a monthly sampling at the outlet of La Moza stream (January–September 2018). The hydrochemical analyses were determined at the Industrial Consulting Chemical Laboratory of the Santander Industrial University (UIS), using titrimetric, turbidimetric, spectrophotometric and atomic absorption methods. The isotopic analyses were carried out at the LAMA Stable Isotope Laboratory of the University of Montpelier in France, using an ISOPRIME isotope ratio mass spectrometer (IRMS) with an analytical error of ±0.8‰ for deuterium and ±0.08‰ for oxygen-18; values are expressed in ‰ vs VSMOW. Tritium was measured at the LHA-Avignon (analytical error 0.2–0.4 UT). The water residence time was estimated from a simple piston-flow model and a natural annual rainfall value of 2.08 UT (based on the data reported by the IAEA-GNIP at the Bogotá station over the period available (2000–2009), where no residual fallout of Tritium from historical bomb atmospheric tests (1950–1963) was detected) (https://nucleus.iaea.org/wiser).

Groundwaters showed predominantly calcium bicarbonate facies, with low geochemical temporal evolution (Table 1 and Figure 2). The resulting chemical composition can be attributed mainly to the hydrolysis process (Appelo & Postma 2005). The springs that outcrop in carbonate rocks presented higher concentrations of HCO₃^– and Ca²⁺, dissolution products of carbonate minerals. The river sample from the middle part of La Moza Stream presented the same hydrochemical facies as the groundwater, showing the link between the groundwater and the river (Figure 2). The saturation indexes (Is) of calcite and dolomite indicate that all spring waters are undersaturated with respect to calcite and dolomite, and there is no potential for mineral precipitation (Table 1).

The springs temperature registered an average of 10.4 °C, showing values close to the ambient annual temperature of the higher parts of the study area, and the pH showed neutral diluted waters (7.7) with changes to slightly acid (6.5) and slightly basic (8.5) (Table 2). The spring water's EC showed a strong link with the rock formation crossed, presenting low values in the non-carbonated rock units (20–79 μS/cm, Jcs and Kils) and medium values in the carbonated units (93–525 μS/cm, Kir). EC showed a temporal variability, with highest values reported in the dry season and lowest values in the rainy season (Table 2).

The isotopic composition of rainwater (δ¹⁸O and δ²H) is shown in Table 3. The variability of δ¹⁸O is relatively consistent with the amount effect predominant in the tropical zone (Clark & Fritz 1997), and in our case linked to seasonality (ITCZ), where a greater amount of rain causes the impoverishment in heavy isotopes. However, the correlation between precipitation and δ¹⁸O in the Picacho and Truchería El Horno stations is moderate (R² = 0.2953 and R² = 0.2511, respectively, Figure 3), because other meteorological factors can influence the isotopic composition, such as cloud type (stratiforme or convective), origin and direction of air mass, vapor volume remaining in the cloud, state of continental vapor recycling, etc.

Picacho Station reported δ¹⁸O between −15.38 and −4.49‰ and a weighted annual average of −11.09‰ (δ¹⁸O_w = Σ(δ¹⁸O_m*P_m)/P_annual where m = monthly value), while Truchería El Horno station presented values between −14.40 and −2.75‰ and a weighted average of −9.86‰, showing that the highest rainfall collector (Picacho Station) presents more depleted values in accordance with the altitude effect (Clark & Fritz 1997). The deuterium excess ‘d’ defined as d = δ²H −8 × δ¹⁸O indicates evaporation processes and water vapor recycling in the atmosphere (Aemisegger et al. 2014). Most monthly ‘d’ values stood between +8 and +12‰, linked to global air masses circulation (except one value at +6.6‰ possibly due to a sample conservation problem). Some months showed a greater deuterium excess (between +12.7 and +18.4‰) during the rainy season, when recycling of continental water vapor in the Amazonian zone mostly occurs (Villacís et al. 2008). However, over the Santander Massif, the origin of atmospheric humidity reaching mountain tops greater than 3,500 masl has not been studied yet, and a study of air mass trajectories is needed to confirm this hypothesis.

The isotopic local meteoric line is represented by the equation δ²H = 8.05 × δ¹⁸O + 12.8, not dissimilar to the Global Meteoric Water Line defined by Craig (1961) (Figure 4).

The isotopic composition of springs and surface water are indicated in Table 2. The spring samples showed δ¹⁸O and δ²H values aggregated about the local meteoric line, close to values of the annual weighted mean rainfall of the study zone, not marked by evaporation (Figure 4). For the surface flow water samples, the same isotopic response like groundwater was recorded, which allows assumption of a same origin between both flows (at least at the moment of sampling).

In the northern part, M1, M9, M6, M7 spring samples showed oxgen-18 values ranging from −11.08 to −10.53‰, −10.56 to −10.11‰, −10.58 to −10.20‰ and −10.65 to −10.64‰, respectively, showing low temporal variability, evidencing the groundwater mixing time at least at the annual scale and the same recharge altitude zone (Table 2). The M1 spring (karstic phreatic aquifer) presented an approximately constant discharge throughout the sampled period (0.21–0.37 L/s), with EC medium values of 225–365 μS/cm (Table 2). The ³H determined in M1 (1.4 UT) indicates a residence time of 7 years, which supports the hydrodynamic and mineralization evolution observed during the year.

The M6 and M7 springs are outlet stream flows from a phreatic aquifer of fractured granodiorites showing an important development of the weathering phenomenon. They showed low annual variation of EC (24–75 μS/cm) but a higher flow variation during the year (0.15–1.89 L/s) compared with the M1 spring (Table 2). Calculated residence time for M6 is 4.7 years (1.6 UT). This shorter residence time compared to M1 seems consistent with a higher discharge variability, however it would be necessary to confirm the value of tritium in the wet season, to test a possible mixing between the aquifer base flow and a rapid infiltration from current rainfall close to the spring.

M9 springs out from a fractured well-cemented quartzarenite phreatic aquifer, with weak mineralization due to their silicate nature (40–77 μS/cm), with a strong annual variability of the discharge (0.07–1.05 L/s) like M6 and M7, and probably a low water residence time.

Two groundwater behaviors are observed in the south-east part. The oxygen-18 values ranged between −10.98 and −10.51‰ for M2, and a slightly higher variability, −10.81 to −10.12‰ for M3, and also showed an important well mixing time and a recharge altitude comparable to the northern part (Table 2). The tritium contents measured indicates different estimated residence times, 7 years in M2 (1.4 UT) and 2.5 years in M3 (1.8 UT), showing faster transit flows in M3, related to a greater degree of fracturing, possibly due to its proximity to the Santana Fault (Figure 1). Discharge annual variability is also greater in M3, indicating a direct contribution of infiltrated rainfall. The M2 spring pertains to a second system of karstic phreatic aquifer in the south-east part, without connection with the northern karstic part, and shows the highest EC (240–525 μS/cm). The M3 is a well-cemented quartzarenite fractured phreatic aquifer like M9, with low EC (20–79 μS/cm) (Table 2).

On the other hand, M4 (fractured phreatic aquifer in a level of immature sandstone part of the south-east karst system), M5 (colluvial phreatic aquifer) and M8 (carbonate rocks) springs showed a larger temporal variability of oxygen-18 over the period, from −11.51 to −8.84‰, −10.48 to −9.25‰ and −10.85 to −9.43‰, respectively, the most enriched values being registered in February 2018 (dry season). Only these three springs showed this behavior. Monthly rainfall in January and February showed enriched values (δ¹⁸O = − 5.75 and −4.49‰ at Picacho Station), compared to −10.92 to −11.36‰ from October to December. This suggests a direct contribution of January–February rainwaters to the spring flow. The high intra annual flow variability observed at these springs confirms the direct rainfall input to the discharge. Tritium is in good agreement with a low residence time for M4 (1.7 UT – 3.6 years). For M8 (1.4 UT – 7 years), the residence time is similar to the other carbonate rock springs M1 and M2, but with a different behavior for oxygen-18 and a higher flow variability. This can result from differing fracturing and karstification levels enabling two types of water circulation, slower inside the karst through different reservoirs with efficient water mixing, and faster through the fractured system in direct connection with springs.

The isotopic compositions of surface waters in the whole watershed of La Moza micro-basin presented low spatial variability with δ¹⁸O values between −10.71 and −10.27‰ (Table 2 – sampling of 11/08/2018 after 3 days without rainfall). The outlet of the catchment at La Moza showed a smoothed annual variability, between −10.68 and −9.47‰, which highlights the important contribution of aquifers to the river with respect to direct rainfall (fast surface and subsurface flows). However, one sampling showed a different behavior in the dry season (14/02/2018) with a very depleted value (−14.56‰). Very little rainfall (11 mm) was detected at Picacho station during the 23 days prior. A plausible hypothesis for this anomalous value could be an important but localized convective rainfall event in the middle or lower part of the catchment.

Gómez et al. (2015) determined an altitudinal oxygen-18 gradient of −0.16‰ per 100 m in the Santander Massif. In this case study, with only two rainfall stations with an altitude difference of 1,000 m, the gradient value is close to −0.12‰, consistent with (but less representative than) the Gomez et al. (2015) study. The weighted average of the two rainfall stations are compared with the altitudinal oxygen-18 gradient in Figure 5. Most sampled springs present a local recharge in terms of altitude, with isotopic values close to the mean annual rainfall of their altitudinal range. The enriched February 2018 values for M4, M5 and M8 correspond to the period of the year of enriched isotope rainfall contents. This may indicate a shorter water residence time and the direct contribution from rainfalls to the springs, impacting the isotope signal in these springs during the dry season. (Figure 5).

A basic mass balance estimation (Rainfall, Runoff, Infiltration) shows an excess of groundwater flow with respect to rainwater falling over the study zone, which could indicate that a part of groundwater flow comes from the adjacent basin Jordan River (Berlin Paramo), which constitutes a relatively flat surface where the recharge process is important. For example, in the base flow period from 27/11/2017 to 20/03/2018, a total rainfall of 166 mm was observed, corresponding to a volume of rainwater of 3,320 m³ over the M2 and M6 watersheds (SAA = 0.02 km²) and 4,980 m³ over the M7 watershed (SAA = 0.03 km²). For the same period the measured base flow (minimum value registered on the period 13/02–24/03/18) gave respective discharges of 2,659 m³ (M2), 2,955 m³ (M6) and 4,137 m³ (M7), corresponding to an infiltration between 80 and 90% while minimizing the spring flow. In steep slope systems, runoff prevails over infiltration. For instance, runoff coefficients have been estimated to be around 56% in a high slope Mediterranean catchment (Gallart et al. 2002). The exceeding low runoff value estimated on the period (10–20%) supports the hypothesis of a contribution coming from another connected watershed feeding the discharge of these springs.

A schematic block diagram (Figure 6) is used to represent the relationship between the groundwaters present in the fractured and karstic rock massifs. The amount of water available in the area depends on an average annual rainfall (P) of 1,300 mm, the real evapotranspiration (RET) value of 600 mm, the estimated potential recharge (P.R) of 440 mm and surface runoff (Sur. Runoff) of 260 mm. The conceptual model shows the variability in the isotopic composition of δ¹⁸O of rainwater (−9.5 to −11.1‰), of each spring and surface water (middle part of La Moza stream −10.5 to −10.7‰ and Outlet La Moza stream −14.5 to −9.4 ‰) and the chemical facies of the groundwater (Calcium bicarbonate-HCO₃Ca). The spring isotopic composition shows circulating patterns that surpass the division line of surface waters with a transit at monthly and annual temporal scales. The north part is not connected with the south part. In the south part, the flows inside the granodioritic rock (Jcs) are connected with the sedimentary rocks (Kir and Kils). All the flows circulate through secondary porosity in the whole geological formations. Finally, the diagram shows the recharge zone of the springs that surpasses the surface water division line, covering mainly parts of the adjacent Jordán basin.

The groundwaters in the La Moza micro-basin are associated with different phreatic aquifers. All geological formations are fractured by tectonic movement and water flows mainly through the secondary porosity (Los Santos Formation) which may increase by dissolution processes impacting the carbonate formation (Rosablanca Formation) producing karstic systems, and by the weathering processes mainly affecting granodioritic rocks (Santa Bárbara Cuarzomonzonite). Water geochemistry in springs shows low to medium electrical conductivity (EC) linked to the rock types (carbonate or silicate rocks). The geochemical facies of the waters is calcium bicarbonate in the whole zone. The spring water isotope content (oxygen-18 and deuterium) indicates circulation flow mode (also supported by the temporal variability of EC and discharge) with a meteoric origin and local infiltration. The tritium dating shows young waters with low residence times (2.5–7 years estimated from a Piston Flow Model) with a more or less important contribution from monthly or annual rainfall.

The surface water isotopic composition shows that a main part of its flow is fed by subterranean systems, due to the fact that the isotopic signatures present low temporal and spatial variations. Finally, a conceptual model of the recharge zone is proposed that crosses the surface watershed and covers part of the adjacent Rio Jordán basin, where the potentially important recharge zone Berlin Paramo is located.