Groundwater recharge assessment in a rural, arid, mid-mountain basin in North-Central Chile

ABSTRACT

A basic component of any hydrogeological study is the magnitude and temporal variation of groundwater recharge. This can be difficult to assess accurately, particularly in arid and semi-arid rainfed mid-mountain zones, as is the situation in the rural, low population density zones of North-Central Chile. In this study, recharge in the Punitaqui Basin, North-Central Chile, was characterized, contrasting the results of two methods: a modified Thornthwaite-Mather (MTM) and discharge recession analysis (DRA). We found a recharge rate of between 1 and 4% of average annual precipitation. Average recharge estimated by the MTM method is consistently higher than that estimated by DRA. Also, DRA tends to smooth the recharge values, resulting in a lower inter-annual variation coefficient. Both methods identified a threshold value of total annual precipitation, above which recharge can be expected to occur, of the order of 180 mm year⁻¹, consistent with values reported in similar areas.

KEYWORDS:

• mountain block recharge
• Punitaqui Basin
• Chile
• dryland
• water balance
• recession analysis

Editor A. CastellarinAssociate editor L. Bouchaou

1 Introduction

Arid and semi-arid zones, where potential evapotranspiration exceeds precipitation, cover more than 30% of the Earth’s surface. Although access to water resources is invariably restricted in these areas, they are increasingly subject to population pressure, and present population growth rates higher than for any other ecological zone (Scanlon et al. 2006, Touhami et al. 2014, UN 2017).

In addition to their socio-economic and political importance, the mountainous portions of arid and semi-arid zones have clear hydrological and environmental consequences for large numbers of people. Mountainous zones, which are areas at altitudes greater than 1000 m a.s.l., or between 200 and 1000 m a.s.l., but with slopes greater than 20%, make up more than 39% of the global continental surface, excluding Antarctica and Greenland (Viviroli et al. 2007). Worldwide, around 67% of surface runoff originates in mountainous basins (Ajami et al. 2011). In this context, a good understanding of hydrological aspects in mid-mountain areas, particularly those located in arid and semi-arid regions, takes on special importance.

Although arid and semi-arid mountainous regions are significant contributors to the availability of water resources for areas at lower elevations, groundwater is commonly the only reliable and permanent source of water in arid and semi-arid regions themselves (Herczeg and Leaney 2011). The evaluation and characterization of these resources is especially difficult in large and/or remote areas, given the costs associated with such investigations, and this is particularly true for basin-scale studies. As a result, few hydrogeological studies are carried out at this scale, especially in rural areas with low population densities where the need for such studies is likely the greatest (Van Camp et al. 2013, Oyarzún et al. 2014a).

A key component of any hydrogeological study is the correct estimation of recharge. Recharge is defined as the fraction of the water infiltrating at the land surface that eventually reaches the water table, “…adding to groundwater storage…” (Healy 2010). Very frequently the estimation of precipitation-derived recharge is a difficult task, due to the spatial and temporal variability of rainfall, and the limited ability to observe and analyse aquifer dynamics. These difficulties are particularly evident in arid and semi-arid areas, as precipitation is affected highly by climatic phenomena such as El Niño Southern Oscillation (ENSO) or the Pacific Decadal Oscillation (PDO) (Núñez et al. 2013). This in turn determines that recharge actually occurs in a sporadic fashion (UNEP/DEWA 2002). In addition, in arid zones precipitation amounts are usually low in comparison to the potential evapotranspiration (PET), and the fraction of precipitation that eventually becomes recharge is therefore comparable in magnitude to the range of uncertainty of the recharge quantification methods (Scanlon et al. 2002).

Despite the fact that groundwater recharge studies and the characterization of recharge events in arid zones have been identified in Chile as key variables for the sustainable use of water resources in the country, especially for rural or small urban communities (CNID 2016), there is in general a lack of such studies in rainfed mid-mountain areas in North-Central Chile. In a context of declining rainfall amounts (Garreaud et al. 2017), there is an urgent need for such knowledge, particularly on “comunidades agrícolas” (i.e. rural communities) areas, which depend on rain as the water source, either directly through rainwater capture or indirectly through groundwater recharge (Oyarzún et al. 2014a, 2014b, 2015). This type of social organization in North-Central Chile, especially important in the Coquimbo Region, dates back to colonial times (Livenais and Aranda 2003, Wilkins and Green 2014), and thus possesses enormous cultural significance. More than 16 000 people, or “comuneros” (members) participate in comunidades agrícolas, and land managed under the comunidades agrícola tradition occupies nearly 25% of the Coquimbo Region (MBN 2017).

In this context, the present work aims mainly to characterize and quantify recharge processes, at the mountain block scale (hereafter referred to as mountain block recharge – MBR), in a representative basin of North-Central Chile, considering the Punitaqui Basin, in the Coquimbo Region, as a case study. Whereas a previous study in this basin estimated recharge at around 5% of mean annual precipitation (DGA 2008), it is important to recognize that groundwater recharge in arid zones such as the Punitaqui Basin is usually dominated by episodic events related to extreme (and thus infrequent) precipitation events (Ajami et al. 2011), and influenced by phenomena such as ENSO and the PDO (Scanlon et al. 2006). In spite of the importance of episodic, high-impact precipitation events that are likely to result in high inter-annual variability in recharge rates, the DGA (2008) study cited above did not consider such issues. Thus, this study was undertaken with the expectation of using the knowledge gained from the Punitaqui Basin to better characterize recharge occurrence and dynamics in arid mid-mountain areas.

2 Methods and materials

2.1 Study area

The study area is the Punitaqui Basin, located in the Limarí Province, Coquimbo Region, North-Central Chile (Fig. 1). The criteria by which this basin was selected include: (a) its climatic and geomorphological representativeness of middle-mountain dry-land areas (Oyarzún et al. 2014a); (b) the existence of flow-gauging information at the basin outlet, a rather uncommon situation in such basins in North-Central Chile (Nauditt et al. 2017); and (c) the existence of previous work on the determination of hydrogeological properties at the basin scale from recession flow analysis (Oyarzún et al. 2014a).

Figure 1. Punitaqui Basin area.

The Punitaqui Basin has a surface area of approx. 1240 km². The main stem draining the basin is the Punitaqui Stream, with an approximate length of 60 km (Oyarzún et al. 2014a). The basin is roughly north–south to northwest–southeast in orientation, and the majority of the area ranges between 200 and 800 m a.s.l.; numerous sub-basins located in the surrounding hills to the east and west contribute runoff that converges to the main stem (the Punitaqui Stream) (Huerta 2009).

The climate is predominantly arid, with a precipitation to potential evapotranspiration ratio in the range of 0.05 to 0.2 (Scanlon et al. 2006, UNESCO 2010). The average annual temperature is of the order of 16°C, ranging between a minimum of 9°C and a maximum of 24°C (Cade-Idepe 2004). Long-term precipitation averages around 160 mm year⁻¹, with significant inter-annual variations associated with the ENSO phenomenon (Veit 1996), but concentrated in the winter months (May–August). The greatest streamflows occur between June and July, and generally attain values of less than 1000 L s⁻¹ (50% exceedence probability), while the lowest streamflows occur between December and March, of the order of 50 L s⁻¹ (also at a 50% exceedence probability) (Cade-Idepe 2004).

Igneous rocks, such as thick-grained, dark grey granodiorites and tonalites of Jurassic and Cretaceous ages, with moderate weathering, predominate in the Punitaqui Basin. Also important are andesitic rock outcroppings, intercalated with sedimentary rocks of the lower Cretaceous (Sernageomin 2003, DGA 2008). Partially consolidated sediments of Pliocene-Pleistocene age are distributed in the northwest of the study area, horizontally arranged over the batholithic rocks, and composed of gravels, sands and clays of alluvial and fluvial origin (Huerta 2009). Finally, a narrow band of recent alluvial unconsolidated gravels and sandy gravels, distributed in the central valley and neighbouring gullies in the northern part of the study area, make up the current Punitaqui Stream floodplain (Huerta 2009). In geomorphological terms, four formations predominate: (a) alluvial terraces, present in the northern part of the basin, (b) low coastal ridges in the northwestern sector, (c) coastal mountains that cross the territory diagonally from the northeast to southwest, and (d) low interior ridges in the southwestern sector (Pouget et al. 1996).

Finally, it is possible to identify hydrogeological units that exist primarily within surficial sediments, and can be subdivided into zones of high- and medium-permeability (DGA 2008). The high-permeability hydrogeological unit is located in the Punitaqui Stream fluvial deposits; it is composed of cobbles, gravels and sands, with finer sediments in lesser quantity, and demonstrates hydraulic conductivities ranging from 10⁻² to 10⁻⁴ m s⁻¹. Locally saturated aquifers in this unit are mainly unconfined, with water tables located at depths of between 5 and 30 m, and specific flows on the order of 1.2 L s⁻¹ m⁻¹. Meanwhile, the alluvial terrace deposits, in the northern part of the Punitaqui Basin, are a semi-permeable unit composed of sands, gravels and silts. The hydraulic potential of these sediments is moderate, with hydraulic conductivities between 10⁻⁴ and 10⁻⁸ m s⁻¹, and specific flows of 0.4 L s⁻¹ m⁻¹. Saturated aquifers in this semi-permeable hydrofacies are also typically unconfined, with water tables located between depths of 20 and 30 m (DGA 2008). Shallow groundwater in the Punitaqui Basin also occurs in weathered and/or fractured granitic rocks in various proportions, with moderate hydrogeological potential, at depths of 3–30 m. When these aquifers occur, they may merge or interact with, or be composed of in part by, colluvial deposits. In an integrative approach, Oyarzún et al. (2014a) determined a hydraulic conductivity on the order of 3 × 10⁻⁵ (m s⁻¹) and a drainable porosity of 7 × 10⁻⁵, for effective parameter values at basin scale.

Within the described physiographic context, this work is focused specifically on MBR evaluated at the basin level. In agreement with Ajami et al. (2011), we define MBR as precipitation that infiltrates through the mountain block on its way to recharge groundwater at lower elevations as well as to feed surface water discharge. For the area of study, this is schematically represented in Figure 2. Therefore, MBR influences both the storage dynamics of groundwater at the basin level, and the behaviour of the surface flow.

Figure 2. Conceptual representation of MBR (curved blue dashed arrows) for the Punitaqui Basin based on three transverse profiles of the study area (depicted in Fig. 1). The vertical green arrow indicates the location of the Punitaqui Stream.

2.2 Basic hydrological characterization

An analysis of the historical hydrological behaviour of the basin and an isotope characterization of surface water components were carried out.

The hydrological analysis addressed the historical dynamics (both monthly and annually) of precipitation and surface discharge patterns. For this, publicly available data obtained from the Chilean Water Authority (DGA, Dirección General de Aguas) were considered.

Regarding the isotope characterization, Punitaqui Stream samples were taken between May 2016 and April 2017, every 1–2 months, at two sites located in the middle part and the end (i.e. confluence with the Limarí River) of the stream, at Chalinga and Desembocadura, respectively (Fig. 1). These samples were analysed for ²H, ¹⁸O and ²²²Rn isotopes, which have been used in earlier hydrogeological investigation in the area to understand surface–groundwater interaction (Oyarzún et al. 2014a, 2014b, 2015). In addition, precipitation samples were obtained during 2015 and 2016 from a rain collector located in the town of Punitaqui (Fig. 1). Isotope analyses of the samples collected in 2015 (rainfall) were performed at the Environmental Isotope Laboratory of the Chilean Nuclear Energy Commission (CCHEN), while those collected in 2016/17 (rainfall, stream water) were analysed at the University of Waterloo-Environmental Isotope Laboratory (uwEILAB). In both laboratories a Liquid Water Stable Isotope Analyzer (Los Gatos Research, CA, USA) was used. Isotope ratios are expressed in per mil (‰) using the usual δ notation:

(1)

where R refers to the ²H/¹H or ¹⁸O/¹⁶O ratio of a sample and the standard (VSMOW), respectively. Informed analytic errors for ¹⁸O and ²H at CCHEN were 0.08 and 1‰ (2015 samples), whereas at uwEILAB (2016/17) these were 0.2 and 0.8‰ for ¹⁸O and ²H, respectively. The ²²²Rn measurements (stream water) were made using alpha spectrometry on RadH₂O equipment (Durridge Co., MA, USA) at the Environment Laboratory of the University of La Serena, the day after each sampling date. Corrections were made to account for radioactive decay between time of sampling and time of analysis in the laboratory.

2.3 MBR quantification methods

2.3.1 Hydro-meteorological information

We employed two methods to estimate recharge in the Punitaqui Basin: the modified Thornthwaite-Mather method, and the discharge recession method. Both methods require the availability of historical series of basic hydro-meteorological data, i.e. precipitation and runoff. There was one raingauge and one stream discharge gauge in the Punitaqui Basin with suitable data (in terms of both temporal extension and continuity). The names and locations of the gauges, and the time period for which data are available for each of the stations, are shown in Figure 1 and Table 1. Note that, for the purposes of this analysis, we expressed discharge in terms of equivalent water depth, i.e. the discharge divided by the surface area of the basin.

Table 1. Rain and flow gauging stations in the Punitaqui Stream basin and data availability.

Station type	Station name	Availability
Raingauge	Punitaqui	1962–2014
Flow gauge	Punitaqui Stream before Limarí River confluence	1962–1977 and 1993–2015

2.3.2 Modified Thornthwaite-Mather method

The first method used for MBR quantification was a modified Thornthwaite-Mather (MTM) method. The Thornthwaite-Mather approach is a well-known technique for estimating recharge, and is of proven utility in hydrogeological studies (e.g. Scozzafava and Tallini 2001, Walraevens et al. 2009, Westenbroek et al. 2010, Uribe et al. 2013, Touhami et al. 2014). In particular, the scheme employed in the present study is based on the work of Uribe et al. (2013). The method relies on a mass balance, taken over a monthly time interval, as expressed by:

(2)

where θ is the unsaturated zone moisture; P is precipitation; ETr is real evapotranspiration; Qr is surface runoff; and R is groundwater recharge, assuming the studied basin does not receive external inputs from other sub-basins (Uribe et al. 2013). The subscript m indicates the time step month (m = 1, 2, 3, …, where m − 1 = 0 is the initial condition at the start of the calculation). Within the balance, it must be taken into account that the unsaturated zone, as a container, has a limited storage capacity. This is known as the maximum storage capacity (θ_max), which is a function of the texture and depth (thickness) of the unsaturated zone (Coral et al. 2015). If the calculated quantity of water in the system by the end of the month (θ_cal m) is greater than θ_max, the excess is assumed to be partitioned to groundwater recharge (Westenbroek et al. 2010, Uribe et al. 2013, Coral et al. 2015). The monthly balance embodied in Equation (2) was implemented using the algorithm presented in Figure 3. In order to apply the algorithm, a value must be established for θ_max; furthermore, an initial moisture content (θ_m−1) must be supplied for each iteration. Unlike the study of Uribe et al. (2013), in the present work we calculated monthly potential evapotranspiration (ETP_m) and determined the crop evapotranspiration (ETc_m) using the crop coefficient (kc), and then defined the value of ETr_m for months in which calculations indicated the existence of moisture in the soil (i.e. (θ_m−1 + P_m) − Qr_m > 0).

Figure 3. Flow diagram of the algorithm used for the MTM method.

Monthly precipitation (P_m) and runoff (Qr_m) values were obtained from the daily records at the selected rain and stream discharge stations (Table 1). The ETP_m was estimated using the Thornthwaite method according to average monthly temperature, with a correction for astronomical day length and the number of days in the month (Xu and Singh 2001, Almorox 2007). Due to the lack of historical temperature data in the study zone, global-scale gridded (GSG) temperature data were used (Matsuura and Willmott 2015). We checked that our estimates of ETP were reasonable by comparing the calculated values against short time series (May 2011–December 2014), available at the Camarico station of the CEAZA-Met network (www.ceazamet.cl; Fig. 1).

To determine the parameters θ_max and kc, we used information from several sources. Land use was simplified to two classes: cultivated surface area and non-cultivated surface area (natural vegetation), and regions of the study area were classified as belonging to one category or the other as determined from (a) soil-use GIS coverage from the CONAF (National Forestry Office) Land Information System (CONAF 2014), (b) irrigated surface GIS coverage of the Limarí River basin (CAZALAC-RHODOS 2006) and (c) soil characteristic GIS coverage (CIREN 2012). For the cultivated area, and to distinguish the crop types according to their evapotranspiration capacities, a representative kc parameter was determined as a function of vegetation type. Also, seasonal variability was considered due to the variable demand of water required during different stages of the crops present, obtained from CAZALAC-RHODOS (2006) and FAO 56 (Allen et al. 2006). An examination of the crop types existing in the study area and respective kc values suggested the use of a basin average value of kc equal to 0.75 (CAZALAC-RHODOS 2006). However, because various species present kc close to 1 (and sometimes greater than that) in the summer period, but much lower in winter (particularly in the case of deciduous fruit trees), we decided to correct kc based on the time of year. Thus, the preliminary basin average value of 0.75 was multiplied by 0.8 for months with lower water demand (April–September) and 1.2 for months with greater requirements (October–March). For the non-cultivated surface, a fixed kc of 0.065, i.e. approximately 1/10 of the average kc of the cultivated surface, was used. The chosen value is based on the work of Kalthoff et al. (2006). These authors determined real evapotranspiration (ETr) rates of 750 mm year⁻¹ in cultivated areas and 65 mm year⁻¹ in rainfed areas with natural vegetation (i.e. approx. 1/10 of the figure of the cultivated area) in the middle (east–west) portion of the Elqui Basin, approx. 80 km to the north of the Punitaqui Basin, with similar physiographic and environmental characteristics.

With respect to the unsaturated zone hydrological properties, we used the field capacity (FC) as a reasonable approximation of the maximum storage capacity (θ_max), and assumed the permanent wilting point (PWP) as the initial soil moisture content (θ₀). The use of the PWP as a starting value for the iterative calculation of moisture content is valid when the initial time (m = 0) corresponds to a dry month (summer), such as January for the area of study; in any case, the calculated soil moisture contents become asymptotically independent of the initial condition after a few time steps. The FC and PWP values were assigned on the basis of the dominant unsaturated zone texture for the area and literature values from Israelsen and Hansen (1975). In the case of cultivated surfaces the soil texture, which is predominantly fine (clayey to clay) loam in the study area, was obtained from CIREN (2012). For non-cultivated surfaces, the unsaturated zone texture was taken from Casanova et al. (2010) and was characterized as loam or silty loam on the surface (0–32 cm) and clay loam or silty clay loam at deeper horizons (32–100 cm). Thus, FC values of 36 and 27% and PMP values of 17 and 13% were used for cultivated and non-cultivated zones, respectively. The unsaturated zone depth was derived from the work of Casanova et al. (2010) and values obtained from CIREN (2012) which, when taken together, suggested using a depth value of 70 cm as an average for the study area. The monthly water balance given by Equation (2), calculation of the model parameters, and processing of the meteorological time-series data (i.e. precipitation) were all implemented in MS Excel.

2.3.3 Discharge recession analysis method

The second method used to estimate recharge was discharge recession analysis (DRA), which is based on the initial work of Kirchner (2009) and the subsequent work of Ajami et al. (2011). This method estimates recharge based on changes in groundwater storage dynamics, which in turn are obtained from local storage–surface runoff discharge (S–Q) relationships acquired in a basin from recession analysis of runoff discharge. Conceptually, the DRA method attempts to infer a functional relationship between groundwater storage and discharge processes from observations of discharge (baseflow) recession. Thus, the storage–discharge function methodology provides an indirect method to quantify mountain block recharge by estimating changes in mountain bedrock storage as a function of changes in stream discharge. This is somewhat different than the commonly used concept of mountain block recharge (e.g. Wilson and Guan 2004). While this does not mean that lateral flow does not exist, we acknowledge that we lack of enough information to quantify it directly. The general approach and specific stages of the method applied in this study are described below; for a more detailed explanation of the theoretical foundations, the reader is referred to Kirchner (2009) and Ajami et al. (2011).

The analysis requires a record of time-series discharge observations when source/sink terms are not contributing significantly to streamflow. Data cannot be influenced by precipitation, and must take place when evapotranspiration is at a minimum. In our case, daily streamflow data for the May–September period of each year of record, i.e. the autumn–winter period (given lower ET, as well as lower potential influence of irrigation return flows from agricultural areas) were used. However, because the May–September period coincides with the occurrence of precipitation events in the area, we used the criteria of Oyarzún et al. (2014a), eliminating from the analysis any discharge data that demonstrated any of the following three characteristics: (a) days with rainfall greater than 10 mm, (b) two days following a day with rainfall greater than 10 mm, and (c) consecutive days with increasing discharge. From this, we plotted Q*, the average discharge of two consecutive days:

(3)

versus the corresponding rate of change in discharge, −dQ/dt, approximated as:

(4)

where the time step, ∆t, is taken as one day.

The Q* versus −dQ/dt values were plotted on logarithmic axes, and the point cloud grouped into bins or sets. To accomplish this binning, two methods were considered: (a) the equal interval method (Oyarzún et al. 2014a) and (b) the quantile method (Ajami et al. 2011). The equal interval method divides the range of the logarithm of discharges into a selected number of bins of equal log-length:

(5)

In the general case, when intervals are chosen in this fashion the number of points per bin is likely to vary widely. For the present analysis we chose 15 as a reasonable number of bins, following heuristic guidelines from previous investigators (Mendoza et al. 2003, Oyarzún et al. 2014a).

For the quantile method, bins are chosen such that each interval contains the same number of data points; thus, the total number of data pairs (Q* and −dQ/dt) are divided by the desired number of bins (in this case, 15 intervals) to obtain the number of data points per bin:

(6)

For each bin (in both criteria) we calculated the average values of both Q* and −dQ/dt, as well as the standard error:

(7)

where is the standard error of the mean, s is the sample standard deviation, and n is the number of points in the bin.

Following the criteria of Ajami et al. (2011), only bins in which the standard error of −dQ/dt was lower than or equal to the average value of (−dQ/dt) were ultimately retained. Thus, with the selected bins, the central trend of the discharge recession data was determined by linear regression. Based on this, the groundwater storage–surface discharge (S–Q) relationship for each year was obtained using the following equation (Ajami et al. 2011):

(8)

where a is the intercept and b is the slope of the corresponding linear regression between the binned −dQ/dt vs Q*. From this, MBR was obtained as:

(9)

The discharge (Q) used in Equations (8) and (9) (Q_max, Q_min) corresponds to the maximum and minimum daily values for each year, being careful that Q_max did not correspond to a day with more than 10 mm of precipitation or the subsequent two days.

3 Results and discussion

3.1 Hydrological behaviour

In Figure 4 we present monthly total precipitation and discharge. Along with the marked variability, it is observed that discharges are rather temporally delayed, by approx. one month, with respect to precipitation. In addition, although precipitation tends to decrease in late winter–early spring (August–September), significant discharge is still recorded in the stream then. Finally, even during the dry period of the late spring–summer (November–January), a “baseflow” is recorded. Thus, it can be inferred that the underlying hydrogeology of the study area exerts important control of surface discharge through the baseflow discharge. Thus, some fraction of the precipitation in the basin infiltrates the mountain block, recharging groundwater and later feeding surface discharge during the dry (non-rainfall) months.

Figure 4. Monthly total precipitation and discharge in the study area.

Figure 5 presents the historical precipitation and discharge in the study area. From inspection of the data, it is clear that both rainfall and discharge are highly variable from year to year, and there are frequent periods of low annual precipitation, which is a common occurrence in rural, arid and semi-arid rainfed areas both in Chile and abroad (Oyarzún et al. 2014a).

Figure 5. Temporal evolution of monthly precipitation and average monthly discharge in the Punitaqui Basin.

There is also a period from 1978 to 1992 for which discharge data are unfortunately not available. However, the comparison of years in which there were similar amounts of precipitation during the periods of available data, i.e. 1962–1977 and 1993–2015, shows that there are still differences in reported discharge. From this it can be inferred that other factors (e.g. rainfall spatial or temporal distribution for the year) influence precipitation–runoff processes in the study area. Furthermore, since 1998 the Punitaqui area (and North-Central Chile in general) has been affected by a period of water scarcity (low and decreasing precipitation), which has translated into historical discharges that are also low and decreasing.

Figure 6 presents the temporal evolution of the isotope signals (²H, ¹⁸O, ²²²Rn) in the 2016/17 period in the Punitaqui Stream at two stations, located at Chalinga and Desembocadura. This particular period was a year with relatively low precipitation, which favoured the use of isotopic compositions for understanding the contributions of the groundwater reservoir to streamflow. Indeed, between 1 April 2016 and 30 April 2017, only 4 days were recorded with precipitation greater than 10 mm, and the total recorded rainfall for that period was approx. 66 mm. The sampling of 26 July 2016 was carried out a day after the greatest precipitation event of the year (approx. 27 mm between 24 and 25 July). Thus, it is possible to compare the isotopic ratios in the baseflow contribution to surface discharge (i.e. streamflow that is maintained by groundwater discharge) to those associated with surficial and shorter-term rainfall–runoff processes. With the exception of the 26 July sampling, ²H and ¹⁸O values in the Punitaqui Stream for 2016/17 are rather similar to those described by Strauch et al. (2006) and Oyarzún et al. (2014b, 2015) for wells and springs associated with hills and low granite ridges in the first third (southern portion) of the basin. Likewise, these values are consistent with the values obtained from precipitation samples in 2016 (n = 3) and 2017 (n = 1) from the collector located in the town of Punitaqui (averages of −5.7‰ and −32‰ for ¹⁸O and ²H, respectively). Although the δ¹⁸O and δ²H ratios generally vary around the average values, the samples of surface water runoff obtained on 26 July 2016 (−9.0/−8.0‰ for δ¹⁸O, and −59.1/−52.1‰ for δ²H, for Chalinga/Desembocadura, respectively) stand out from the samples collected on other dates. Given that the rainfall isotopic signal of the 24–25 July event was −5.9‰ and −37.0‰ for δ¹⁸O and δ²H, respectively, the depleted values in stream water of 26 July could be attributable to a significant contribution to surface water of precipitation falling in the basin headwater, upstream of the town of Punitaqui (where the rainfall collector was installed).

Figure 6. Values of ²H, ¹⁸O and ²²²Rn for the Punitaqui Stream at (a) Chalinga and (b) Desembocadura in the 2016–2017 period.

The interpretations for the oxygen and hydrogen isotope ratios are confirmed by the ²²²Rn data. Green and Stewart (2008) indicated that values of ²²²Rn activity above 1000 Bq m⁻³ in surface water can be considered indicative of some degree of exfiltration of groundwater (i.e. the stream is in a gaining, or influent, condition), and the values of ²²²Rn activity, generally above 1400 Bq m⁻³, appear to support the supposition that discharge in the Punitaqui Stream contains a non-negligible fraction of baseflow at all sampled times during the year. As with the δ¹⁸O and δ²H ratios, the exception to this general behaviour is provided by the 26 July 2016 sampling event, when direct surface runoff from the antecedent rainfall represents a major part of stream discharge. Taken together, both the analysis of the historical precipitation and discharge data and the isotopic evidence suggest that surface water discharge in the Punitaqui Stream contains a variable amount of surface runoff, along with an identifiable fraction of groundwater discharge (baseflow) that is the result of MBR.

3.2 MBR estimates

The results of the MBR estimations with the MTM and DRA methods are presented below. Based on the final availability of data, the following assessment sub-periods (for both methods) were selected: 1963–1966, 1973–1975, 1994–1998 and 2003–2010.

3.2.1 MTM method

Table 2 shows the estimates of cultivated surface area in the basin, obtained from the review and processing of various GIS layers, as discussed previously. The percentage of the total area that is cultivated fluctuates, depending on the GIS layer reviewed, between 8 and 11% of the basin; for the present study, we chose to use a cultivated surface area value of 10%, with the non-cultivated surface area (with natural vegetation) set at 90%.

Table 2. Cultivated surface area of each basin according to various GIS layers.

Source	Surface area (ha)	Surface area (%)
Visual inspection using Google Earth	12 510	10
CIREN (2012)	11 377	9
CONAF (2014)	14 469	11
CAZALAC-RHODOS (2006)	10 394	18

With the proportion of each type of plant coverage defined, a weighted crop coefficient (kc) value was obtained for the entire basin. In addition, as mentioned in the methods section, weighted values of FC and PWP were determined, as shown in Table 3.

Table 3. Weighted parameter values used in the modified Thornthwaite-Mather water balance.

Parameter	Value
Weighted kc (Apr.–Sep.)	0.12
Weighted kc (Oct.–Mar.)	0.15
FC (mm)	258
PWP (mm)	124

Based on the indicated values and in accord with the described procedures, the results of the water balance are presented at annual scale in Table 4, for each of the four sub-periods of assessment. It is clear that, while precipitation and discharge present significant inter-annual variability, the same does not occur in the case of evapotranspiration, which is more stable over time. In fact, with the exception of the years of its maximum and minimum values (1997 and 1966, respectively), ETP values are close to the average for all the years. In addition, it is observed that in most of the assessed years, the value of ETr, as a percentage of precipitation, exceeds the threshold of 100, an expected situation in an arid climate.

Table 4. Summary of hydrological variables considered in the water balance. P: precipitation, Qr: runoff, ETP: potential evapotranspiration, ETr: real evapotranspiration, MBR: recharge. SD: standard deviation; CV: coefficient of variation.

		Qr			ETr		MBR
Year	P (mm)	(mm)	(% of P)	ETP (mm)	(mm)	(% of P)	(mm)	(% of P)
1963	300.0	20	7%	1012	142	47%	51	17%
1964	75.0	7	9%	1014	142	190%	0	0%
1965	289.2	45	16%	1019	143	49%	45	16%
1966	188.6	19	10%	942	132	70%	32	17%
1973	70.5	12	17%	1028	144	204%	0	0%
1974	106.5	8	8%	993	116	108%	0	0%
1975	110.0	3	2%	1000	116	106%	0	0%
1994	74.7	8	10%	1029	144	193%	0	0%
1995	39.0	6	15%	1018	92	236%	0	0%
1996	76.0	3	4%	1014	143	188%	0	0%
1997	510.5	19	4%	1068	149	29%	191	37%
1998	10.2	5	49%	1032	144	1416%	0	0%
2003	140.5	16	11%	1020	143	102%	0	0%
2004	182.5	15	8%	1024	143	79%	0	0%
2005	71.3	4	6%	1024	144	201%	0	0%
2006	110.0	3	2%	1041	137	124%	0	0%
2007	27.2	2	6%	1008	77	282%	0	0%
2008	151.5	3	2%	1034	79	52%	0	0%
2009	102.5	1	1%	1022	143	140%	0	0%
2010	149.5	2	1%	1013	117	78%	0	0%
Average	139.3	10	9.4%	1018	129	195%	16	4%
SD	116.0	10.5	10.5%	24	23	296%	44	10%
CV	83%	105%	111%	2%	17%	152%	276%	226%

In contrast to the estimated values of ETP, estimated MBR appears to be quite variable and rather episodic, only attaining an observable magnitude during the few years with elevated precipitation (above approx. 188 mm year⁻¹). The long-term (historical average) estimated recharge rate is on the order of 4% of precipitation; however, this average value is somewhat misleading. The rate of recharge is actually highly variable (coefficient of variation of 200%, much greater than that obtained for discharges and evapotranspiration). In fact, the mode of the distribution of the annual recharge is zero, and no recharge events are indicated for the most recent time period, 2003–2010.

3.2.2 DRA method

We plotted the data as Q* vs −dQ/dt relationships for both log-equal intervals and intervals containing equal numbers of data points. The results are presented in Figure 7.

Figure 7. Q* vs −dQ/dt relationship (blue dots), and bins (orange dots) determined using (a) the equal interval method and (b) the quantile method.

Figure 8. Temporal evolution (± standard deviation) of the average water table levels (lowest annual value) from the six DGA monitoring wells in the Punitaqui Stream basin.

In the case of the equal interval method, only 12 of the 15 originally defined bins were kept (given the standard error criterion, Section 3.2.3). The bins defined by this method tend to be concentrated towards the highest discharge values, which also coincide with somewhat lower data quantity. In contrast, the bins defined by the quantile method were more evenly distributed across all values of −dQ/dt. All the bins defined by the quantile method met the standard error criterion and were retained in the analysis. We applied (log) linear regression to the binned data values shown in Figure 7, and obtained slope and intercept estimates, as well as estimates of the significance of the values (Table 5).

Table 5. Linear and polynomial regression coefficients for the quantile and equal interval methods.

			Coefficients
Grouping method	R^2	p value	a	b
Equal intervals	0.907	<0.01	0.208	1.204
Quantile	0.934	<0.01	0.104	0.926

It is interesting to note the high R² values obtained in the analysis, which are even greater than some of the values presented in Ajami et al. (2011). Although not a quantitative assessment of correlation (as developed, for example, from a Pearson correlation test), the calculated R² is an indication of the fraction of variability in the dependent variable that can be predicted by the independent variable (in this case, >90%). The quality of the fit implies that the S–Q (storage–discharge) relationship in the study area can be described adequately by the central trend of the recession analysis (Ajami et al. 2011), thus fulfilling a basic assumption of the method. Of the two clustering schemes, the data grouped by quantile yield slightly greater regression coefficients. The inferred storage–discharge relationships were MBR = 6.04Q^0.796 and MBR = 8.59Q^1.074, based on the equal interval and quantile grouping methods, respectively. In spite of the different relationships and slightly different R² values obtained from the separate binning methods, the annual MBR values obtained (Table 6) are ultimately not significantly different. Also, the long-term MBR to precipitation ratio was observed to be 1% (for both methods of bin grouping).

Table 6. Hydrological year period MBR estimates, with maximum and minimum discharge each year, using the equal interval (EI) and quantile (QU) bin grouping methods. Dates (columns 4 and 6) in the form dd/mm.

Year	P (mm)	Qmax (mm d^−1)	Date Qmax	Qmin (mm d^−1)	Date Qmin	MBR [EI] (mm year^−1)	MBR/P (%)	MBR [QU] (mm year^−1)	MBR/P (%)
1963	299	0.6055	07/09	0.0058	23/03	3.9	1.3	5.2	1.7
1964	75	0.1080	21/08	0.0058	20/03	0.9	1.2	0.8	1.1
1965	289	0.7848	18/08	0.0039	19/04	4.9	1.7	6.9	2.4
1966	189	0.5242	14/07	0.0008	07/07	3.6	1.9	4.5	2.4
1973	71	0.1642	06/05	0.0006	24/02	1.4	2.0	1.3	1.8
1974	107	0.1462	10/06	0.0036	22/03	1.2	1.1	1.1	1.0
1975	110	0.0340	05/07	0.0001	29/02	0.4	0.4	0.2	0.2
1994	54	0.0828	28/07	0.0014	04/02	0.8	1.5	0.6	1.1
1995	32	0.0660	22/09	0.0013	31/03	0.7	2.2	0.5	1.6
1996	76	0.1195	09/07	0.0006	29/03	1.1	1.4	0.9	1.2
1997	511	0.1592	14/10	0.0006	07/04	6.8	1.3	10.4	2.0
1998	12	0.0547	29/04	0.0003	05/01	0.6	4.9	0.4	3.3
2003	146	0.2462	24/05	0.0024	19/02	1.9	1.3	2.0	1.4
2004	178	0.3514	05/08	0.0025	23/03	2.6	1.5	2.9	1.6
2005	71	0.0616	16/09	0.0030	05/04	0.6	0.8	0.4	0.6
2006	110	0.0263	29/07	0.0029	03/02	0.3	0.3	0.2	0.2
2007	27	0.0085	31/08	0.0022	27/02	0.1	0.4	0.0	0.0
2008	152	0.0325	19/08	0.0008	10/01	0.4	0.3	0.2	0.1
2009	103	0.0158	21/08	0.0006	31/01	0.2	0.2	0.1	0.1
2010	150	0.0242	09/07	0.0006	28/05	0.3	0.2	0.2	0.1
Average	138	0.1808		0.0020		2	1	2	1
SD	117	0.2188		0.0017		2	1	3	1
CV (%)	85	121		87		112	82	143	78

3.2.3 Integrated analysis

The average annual recharge estimated for the Punitaqui Basin, expressed as a percentage of average precipitation, is 4 and 1% for the MTM and DRA methods, respectively. These values are of the same order of magnitude as the average figure described by DGA (2008). However, the current study also allowed us to estimate the temporal recharge dynamics. Indeed, annual recharge values obtained by MTM and DRA are very variable throughout the period of study.

The average recharge from DRA is lower than that obtained with the MTM. This could be explained by the fact that the “dynamic storage” (or MBR) is lower than the “total storage” of a basin, since basins “can retain significant volumes of residual water, even under drought conditions” (Kirchner 2009). Similarly, Ajami et al. (2011) indicated that “stream gauging stations may only capture part of the MBR flow path”. Thus, recharge determined by the DRA should be regarded as a “lower limit” of recharge for a given area.

Both approaches yield no recharge from 2005 onward, consistent with the period of high water scarcity (low rainfall) that has affected the Coquimbo Region in recent years (Garreaud et al. 2017, Nauditt et al. 2017). This result is consistent with water table dynamics recorded in six wells of the Chilean Water Authority (DGA) monitoring network in the Punitaqui Basin (Figure 8). It should be acknowledged that these wells are primarily located in the northern part of the basin, in the alluvial terrace area (Fig. 1), so they may be influenced by agricultural activity in their vicinity (e.g. groundwater extraction, irrigation return flows). However, they actually show decreasing groundwater levels for the last 10 years, which could be related with, among other factors, extremely low (or non-existing) recharge, as identified in this work.

Also, it is interesting to note that the MTM method predicts many years of zero recharge. This outcome contrasts with the DRA results, for which a zero-recharge situation is not expected to occur. Because some infiltration is retained within the mountain block, rather than being discharged during the year it is recharged, the DRA tends to distribute or “smooth” inter-annual variations in recharge. Thus, years that yield a zero-recharge value by the MTM are predicted to have a finite positive recharge rate by DRA, reflecting MBR released from storage. Indeed, while 1998 was the driest year (least rain) of the analysed data record, with an annual precipitation of only 12 mm, DRA nonetheless registers a recharge value, albeit a low one. In other words, part of the recharge that occurred in previous years, especially in 1997 (which registered the greatest precipitation value, 511 mm) must have remained stored in the mountain system and been released gradually. This explains why the DRA method indicates the existence of recharge in years in which the water balance method predicts negligible recharge as a consequence of “insufficient rainfall”.

It should be pointed out that, in general, balance methods such as Thornthwaite-Mather are usually used for small- or medium-sized basins, with surface areas from dozens to a few hundred square kilometres (Westenbroek et al. 2010, Gong et al. 2012). Similarly, in the work of both Kirchner (2009) and Ajami et al. (2011), the DRA method was applied to small- and medium-sized basins (i.e. 20 and 100 km², respectively). In fact, Kirchner (2009) highlighted the need for additional studies to assess the method in basins larger than those assessed by that author. Given that the results obtained in the present work with both MTM and DRA are generally consistent and plausible, it can be inferred that these approaches may be applicable in larger basins, e.g. on the order of 1000 km², at least in general terms with regard to the determination of average recharge rates over a period of several years. The satisfactory results obtained from the DRA method in the present study are also interesting in terms of application of the method within an arid to semi-arid climate context. In the study of Kirchner (2009), precipitation levels in the study area exceeded 2500 mm year⁻¹, while in the case of Ajami et al. (2011), they fluctuated between 300 and 900 mm year⁻¹. Thus, reasonably consistent results obtained from the DRA and MTM methods for the Punitaqui Basin suggest that these methods may produce acceptable results over a wider range of conditions (both size and climatic) than those to which they have previously been applied.

Finally, both methods identify a minimum total annual precipitation threshold, above which recharge could be expected to occur. This is around 190 mm for the MTM, as for years with precipitation lower than this value no recharge is determined. In the case of the DRA method, and by taking 2 mm as a “reasonable” minimum recharge value (considering the inherent uncertainties of the precipitation and discharge records required by the method), the precipitation threshold is around 180 mm. These values are consistent with a result obtained by Scanlon et al. (2006), whose observations indicate that “rainfall below 200 mm usually results in negligible recharge”. The identification and validation of a minimum precipitation value, below which recharge would be unlikely to occur, has great practical importance in arid and semi-arid regions such as the Punitaqui Basin and other similar situations elsewhere. For example, a recent study financed by the Coquimbo Regional Government advocated the use of “rain harvesting” initiatives as a method for artificial recharge to aquifers in arid rainfed areas, even under conditions of very low precipitation (Romero 2015). The results of the present study indicate that practical application of such tactics for increasing recharge would not be effective unless considerable attention was paid to the nature, amount and variability of the rainfall itself, as well as to methods for circumventing the impacts of evapotranspiration and near-surface storage. Without considering the implications of the results presented here, it is likely that many attempted mitigations would be ineffective and a wasteful use of public resources.

Although the good correspondence between the results achieved by the two methods (MTM and DRA) provides confidence, it is reasonable to expect that some uncertainties are associated with the study findings. In particular, we note that the data analysed in this study for the MTM method begin in the early 1960s and continue until the present; however, over this timeframe it is likely that crop characteristics, such as crop types and areas of coverage, have changed and evolved. Since it was not possible to make corrections for such changes, it is not clear to what extent the results may have been affected by land-use trends. It should also be noted that the availability of data from only one raingauge in the basis introduces uncertainty into the analysis, particularly considering the size of the basin and the spatial variability in precipitation that could reasonably be expected at the basin scale. In spite of these uncertainties, the results obtained are supported not only by the consistency between the two methods, but by the fact that they are of the same order of magnitude as found by other investigators (especially for studies of arid zones), particularly in terms of the occurrence (timing, magnitude) of recharge (Scanlon et al. 2006, Healy 2010). While we do not present these results as exact values, they can be interpreted as providing a robust order of magnitude assessment of the expected recharge rates and timing in arid zones, such as the Punitaqui Basin and others elsewhere, suitable for planning purposes and as a guide to the design of future water management measures.

4 Conclusions

The work used a lumped approach to characterize the mountain block recharge (MBR) process in the Punitaqui Basin, considered representative of rainfed mid-mountain systems in North-Central Chile. First, using isotope and hydrological data allowed us to verify that a major source of surface flow corresponds to MBR. Then, we applied modified Thornthwaite-Mather (MTM) and discharge recession analysis (DRA) methods to examine the magnitude and temporal variation of groundwater recharge in the basin. On the basis of the results, it appears possible to use either the water balance or the discharge recession analysis method to quantify average recharge rates and temporal recharge dynamics, and we expect the general trends inferred for the Punitaqui Basin to be representative of many of the rainfed, mid-mountain systems in the arid North-Central region of Chile, and possibly of other, similar arid and semi-arid basins elsewhere.

Although there are year-to-year differences between the MTM and DRA methods, overall they converge toward a determination of recharge rates at around 1–4% of the annual precipitation value as a long-term average, consistent with values described in the literature for arid zones. However, the two methods predict slightly different behaviours when examined in detail. In particular, inter-annual variability in precipitation and crop type/cover observed in the data is more clearly reflected in the recharge values determined by the MTM method. In contrast, we observe that the DRA approach smoothed recharge on an annual basis by capturing annual carry-over storage within the mountain block. As a result of this smoothing, the method registers the existence of recharge in years in which the MTM method predicts no recharge due to “insufficient precipitation”. In any case, our analyses imply that significant recharge is unlikely except in years for which precipitation exceeds the annual average by around 20–25%, i.e. when annual precipitation reaches (or surpasses) approximately 180 mm.

The methods of analysis presented here are relatively simple in terms of their use, and our study shows they can be applied to good effect in arid to semi-arid basins on a medium scale (of the order of hundreds to a few thousand square kilometres), provided the basic information required for their implementation is available. In this sense, the DRA method is simpler to implement than the MTM method, as it is less demanding in terms of information requirements (i.e. “only” precipitation and daily discharge values).

Acknowledgements

The authors would like to thank Professor E. Muñoz (UCSC) for providing the gridded temperature data. This work was developed within the framework of Conicyt/Programa Fondecyt Project 1150587, the Water Resources and Environment Programme (PRHIMA, for its acronym in Spanish) of the Department of Mining Engineering of the Universidad de La Serena, and Conicyt/Fondap/15130015 Center CRHIAM. The paper benefited from the comments of two anonymous reviewers.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by the Comisión Nacional de Investigación Científica y Tecnológica [Fondap 15130015]; Fondo de Fomento al Desarrollo Científico y Tecnológico [Fondecyt 1150587].