Arti ﬁ cial recharge by means of careo channels versus natural aquifer recharge in a semi-arid, high-mountain watershed (Sierra Nevada, Spain)

The acequias de careo are ancestral water channels excavated during the early Al-Andalus period (8th – 10th centuries), which are used to recharge aquifers in the watersheds of the Sierra Nevada mountain range (Southeastern Spain). The water channels are maintained by local communities, and their main function is collecting snowmelt, but also runoff fromrainfallfromtheheadwatersofriverbasinsanddistributingitthroughouttheupperpartsoftheslopes.Thismethod of aquifer arti ﬁ cial recharge extends the availability of water resources in the lowlands of the river basins during the dry season when there is almost no precipitation and water demand is higher. This study investigates the contribution of the careo channels in the watershed of Bérchules concerning the total aquifer recharge during the 2014 – 2015 hydrological year.Several channels were gauged, andthe runoffdata werecompared with thoseobtained froma semi-distributed hy-drologicalmodelappliedtothesamehydrologicalbasin.Thenaturalin ﬁ ltrationofmeteoricwatersaccountedfor52%of the total recharge, while the remaining 48% corresponded to water transported and in ﬁ ltrated by the careo channels. In otherwords,the careo rechargesystemenhancesby92%thenaturalrechargetotheaquifer.Ourresultsdemonstratethe importance of this ancestral and ef ﬁ cient channel system for recharging slope aquifers developed in hard rocks. The ace- quiasdecareo arenature-basedsolutionsforincreasingwaterresourcesavailabilitythathavecontributedtoaprosperous life in the Sierra Nevada. Its long history ( > 1200 years) suggests that the system has remarkable resilience properties, which have allowed adaptation and permance for centuries in drastically changing climatic and socioeconomic conditions. This recharge system could also be applied to — or inspire similar adaptation measures in — semi-arid mountain areas around the world where it may help in mitigating climate change effects.

The acequias de careo are ancestral water channels excavated during the early Al-Andalus period (8th-10th centuries), which are used to recharge aquifers in the watersheds of the Sierra Nevada mountain range (Southeastern Spain). The water channels are maintained by local communities, and their main function is collecting snowmelt, but also runoff from rainfall from the headwaters of river basins and distributing it throughout the upper parts of the slopes. This method of aquifer artificial recharge extends the availability of water resources in the lowlands of the river basins during the dry season when there is almost no precipitation and water demand is higher. This study investigates the contribution of the careo channels in the watershed of Bérchules concerning the total aquifer recharge during the 2014-2015 hydrological year. Several channels were gauged, and the runoff data were compared with those obtained from a semi-distributed hydrological model applied to the same hydrological basin. The natural infiltration of meteoric waters accounted for 52% of the total recharge, while the remaining 48% corresponded to water transported and infiltrated by the careo channels. In other words, the careo recharge system enhances by 92% the natural recharge to the aquifer. Our results demonstrate the importance of this ancestral and efficient channel system for recharging slope aquifers developed in hard rocks. The acequias de careo are nature-based solutions for increasing water resources availability that have contributed to a prosperous life in the Sierra Nevada. Its long history (>1200 years) suggests that the system has remarkable resilience properties, which have allowed adaptation and permance for centuries in drastically changing climatic and socioeconomic conditions. This recharge system could also be applied to -or inspire similar adaptation measures in-semi-arid mountain areas around the world where it may help in mitigating climate change effects.

Keywords:
Careo channel Slope aquifer Nature-based solution Managed aquifer recharge Science of the Total Environment 825 (2022) 153937

Introduction
Sierra Nevada is a high mountain range located in the southeastern Iberian Peninsula. It formed as a consequence of the convergence between the Eurasian and African plates during the Alpine orogeny. Its most prominent peaks are Mulhacén (3482 m a.s.l.) and Veleta (3398 m a.s.l.), the highest elevations of the peninsula. Its glaciers retreated, almost disappearing, from the massif 13-15 ka BP, largely due to the proximity of the African continent, which imposes semi-arid climatic conditions on this massif and its surroundings (Gómez-Ortiz et al., 2012.
The past and present hydro-climatic conditions of the Sierra Nevada massif have never been ideal for the development of permanent communities because of the high average elevation and the prevailing semi-arid climate conditions. Although the first evidence of human occupation can be traced to Late Neolithic times, around 3000 BCE (Calatrava and Sayadi, 2019), it was not until the Al-Andalus period -under Muslim rule, between the 8th and 15th centuries-that human activity became evident. With the construction of extensive terrace cultivations all around the mountain range, early settlers began to change the natural landscape drastically, especially on and around the southern slopes, where both a warmer climate and a gentler topography can be found. The oldest permanent human settlements of the area (e.g., Lanjarón, Trevélez, Bérchules) ( Fig. 1) still exist today (Martín Civantos, 2007, 2008. To solve the problem of growing water demand due to agricultural activity and associated settlements, people began excavating vast networks of water channels, locally known as acequias de careo. These channels gently descend along the contour lines of the slopes, have no impermeable lining, range in size from 0.5 to 2.5 m in width, and can be several kilometers long (Pulido-Bosch and Sbih, 1995). The snowmelt, which mostly flows through the mountain streams in spring, is diverted into such channels with the objectives of (1) transporting and distributing water towards irrigated areas, (2) transferring surplus water among neighboring hydrological basins, and (3) regulating basin resources through the recharge of slope aquifers with snowmelt water (Martos-Rosillo et al., 2019a). This system of water use and management, based on the principle of water sowing and harvesting using acequias de careo, has been documented since the Early Middle Ages (8th-10th centuries) (Martín Civantos, 2010), and is still in operation today, exemplifying the new paradigm of Integrated Water Management (Vivas et al., 2009). Acequias de careo are solutions inspired by and supported by nature that use or mimic natural processes to help improve water management (WWAP, 2018) and can, thus, be regarded as naturebased solutions for water management .
Ancestral water management systems based on indigenous people's wisdom exist in various other places in the world. In South America, particularly in Andean arid and semiarid regions (Yapa, 2013;Ochoa-Tocachi et al., 2019;Martos-Rosillo et al., 2020), there is a variety of ancient water sowing/harvesting systems (WSHS) to regulate aquifer recharge . They include the pre-Hispanic amunas of Peru (Ochoa-Tocachi et al., 2019), which developed independently on a different continent but resemble careo channels of Sierra Nevada in design and operation (Martos-Rosillo et al., 2019a, 2019b. In Perú, hundreds of kilometers of abandoned amunas are currently being restored (Cárdenas-Panduro, 2020), all of them in river basins that provide water to Lima. Another good example of WSHS is the cocha or albarrada in Ecuador, Peru and Bolivia. These are ponds of low crest heights built for water infiltration, which is enhanced by using earth construction materials (MINAGRI, 2016;Martos-Rosillo et al., 2020;Albarracín et al., 2021). Similar to cochas, but smaller, are Peru's cuchacuchas, infiltrating ponds with diameters ranging between 2 and 15 m (Yapa, 2016). The tape, another type of WSHS that can be found in Ecuador, consists of tiny walls set up along the main channels of intermittent streams and rivers with the aim to collect water during the rainy months. The infiltrated water can then be tapped via wells or drainage galleries (Carrión et al., 2018;Martos-Rosillo et al., 2020). Surprisingly, this WSHS can be also found in Kenya (Lasage et al., 2008). Apart from these WSHS, numerous Andean high-altitude wetlands known as bofedales are watered and expanded by constructing ditch networks to irrigate pastures and infiltrate water . Downstream of all these WSHS, biodiversity is boosted; even species suited for more humid conditions flourish through the operation of these systems, allowing a diverse range of floral species with widely varying water requirements to coexist (Yapa, 2013;Martos-Rosillo et al., 2020;Albarracín et al., 2021).
All the aforementioned ancient WSHS approaches can be thought of being Managed Aquifer Recharge systems combining "blue" and "green infrastructure" and/or being a generator of them (Benedict and McMahon, 2012). All of them are clear examples of Nature Based Solution for Water Managemenet (NbSWM), representing an antagonistic concept to water management with grey infrastructure (i.e., large concrete hydraulic works) on which current water management models rely. Dams or reservoirs are a good example of such grey infrastructure. They are believed to mitigate natural calamities such as floods, landslides, and debris flows while facilitating groundwater recharge. Nevertheless, the option of building more reservoirs is socially, environmentally, and economically controversial (Ribeiro, 2021). In this line, Spain is one of the countries in the Mediterranean zone with the largest number of big dams per inhabitant. Despite of that, Spain is suffering increasing intensity of long-lasting droughts, which are challenging the capacity of such "grey infrastructures", alone, for satisfying human water demands without serious threats to water-dependent ecosystems. The ancient WSHS as NbSWMs are more efficient and ecosystem-friendly forms of water storage as groundwater recharge than grey infrastructures, which in turn are less sustainable and cost-effective than NbSWMs (WWAP, 2018).
This study focuses on the Bérchules river basin, located in the central sector of the southern slopes of the Sierra Nevada ( Fig. 1), where a WSHS system of acequias de careo still operates. The study area has a great demand for water resources for two reasons: (1) the environmental policy of afforestation with monoculture of coniferous forest, which began in the mid 20th century and continues today (Jiménez-Olivenza et al., 2015), and (2) the transformation of agricultural land from traditionally rainfed to intensively irrigated land for horticultural production . Although the hydrological balance of the basin has been documented Jódar et al., 2017Jódar et al., , 2018, the contribution of the careo channels to the total water resources has not yet been quantified at the hydrological basin scale, and the contribution of the careo channels to the total aquifer recharge is still unknown.
The objective of this work is twofold: 1) to present a methodological approach for estimating the infiltration capacity of the careo channels, and 2) to evaluate the careo channels' contribution to total groundwater recharge in the Bérchules basin. This will help understand the resilience of such ancestral water management systems to past drastic social (Martín Civantos, 2007, 2010 and climate changes (García-Alix et al., 2020;Ramos-Román et al., 2016) that occurred in Sierra Nevada (Spain) from the Middle Ages to date, while giving sound arguments to maintain these aquifer recharge channel systems as adaptation measures to minimize the impact of the forecasted climate change.

Study area
The study area lies on the southern face of Sierra Nevada, in southeastern Spain (Fig. 1). The altitude of the Bérchules basin varies between 2913 m a.s.l. at Cerro del Gallo peak and 979 m a.s.l. in the Narila gauging station, located at the outlet of the basin, where the average discharge flow for the period 1970-2015 was 12.6 hm 3 /year. This basin has a surface area of 68 km 2 . Its main drainage network is formed by the Grande de Bérchules River, with a total length of 17.3 km from the headwaters of the basin (to the NE), at 2600 m a.s.l., to Narila's gauging station. The Chico River is the largest tributary of the Bérchules River. It originates, over 2900 m a.s.l., to the northwest of the basin's headwaters. This tributary joins the Grande de Bérchules River in the central part of the basin.
From a meteorological point of view, following the classification of Köppen-Geiger (Peel et al., 2007), the area has a cold climate, with mild and cool dry summers, together with significant altitudinal and thermal variations (Rigueiro-Rodríguez et al., 2008;Gómez-Zotano et al., 2015). The mean temperature (T) at the meteorological station of Bérchules (1319 m a.s.l.) (Fig. 1D) is 13.3°C. Mean precipitation (P) and potential evapotranspiration (PET) are 677 and 1012 mm/year, respectively. P and T show a clear altitudinal dependence on the southern slope of Sierra Nevada (Fig. 2). In the basin, the presence and persistence of snow is restricted to elevations higher than 2400 m a.s.l. and to a period between November and May. Scrubs cover >50% of the total basin area, and conifers, which have been introduced for reforestation purposes, especially during the second half of the 20th century, cover 15% of the basin. In the lower part of the watershed, the climate is temperate Mediterranean. Irrigated horticultural cropland extends over 1.3% of the total basin surface, up to heights above 2200 m a.s.l. .
The Bérchules watershed is situated over poorly permeable rocks, mostly schists, of the Nevado-Filábride Complex (Fig. 1). Weathering processes have contributed to the development of a relatively permeable layer over most of its surface, with a fissured and fractured zone topping the unaltered rock ( Fig. 3). Quaternary formations containing glacial and periglacial sediments, including colluvium and slope debris, are mainly found at heights above 2000 m a.s.l., overlying and complementing the alteration zone. From the hydrogeological point of view, the altered schists and Quaternary materials together constitute a 30 to 40 m thick formation that holds a shallow unconfined aquifer with transmissivity ranging between 0.05 and 0.5 m 2 /s (IGME, 2015), with the most transmissive zone located between 20 and 30 m depth. The permeable outcrops cover 59 km 2 of the total 69 km 2 basin area. Unaltered hard rocks can be found throughout the rest of the watershed. The boundaries of this upper aquifer coincide with those of the surface watershed.
In the Sierra Nevada and surrounding areas, a long-term shortage of water supplies has been a major impediment to the establishment of permanent communities. It was not until the Al-Andalus period in the 8th century that people of Arabic-Berber origin started to tackle this problem efficiently. They developed a water management system that makes use of an extensive network of water infiltration channels (Delaigue, 1995;Martos-Rosillo et al., 2019b). These hydraulic structures, locally known as acequias de careo (Fig. 1C), are relatively shallow and narrow channels dug with simple tools and strengthened with local materials, but not sealed (Espín et al., 2010). The channels' main function is to collect water and facilitate its infiltration. The study area features three types of infiltration channels (Fig. 3): 1. Careo channels: Channels that are solely used for aquifer recharge, a practice known locally as careo. Their main function is to collect runoff from the headwaters of the rivers and to infiltrate it along its way and into pre-determined infiltration spots, locally called simas, where the soil infiltration capacity is known to be high. They usually operate from the beginning of autumn (October) until the end of the snowmelt season (June), with some variation depending on the meteorological characteristics of each hydrological year. 2. Recharge and irrigation channels: Channels that are used for both groundwater recharge (i.e., careo; October-June) and crop irrigation of agricultural areas situated below the channels (July-September).
Land that is irrigated with the traditional flooding method, or more recently using localized irrigation methods, comprises approximately 0.88 km 2 of the area (~1.3% of the surface area of the Bérchules watershed; Carpintero et al., 2018). 3. Recharge and water transfer channels: Channels that are used for groundwater recharge in a watershed while simultaneously transporting water into a neighboring one. The water resources are shared equally between both neighboring basins during autumn and winter. From April until the end of the snowmelt season (June), the distribution of such water resources changes: the total channel flowrate is divided between the original (i.e., where the snowmelt is generated) and neighboring basins, 3/7 and 4/7 respectively. From June onwards, usually starting at the summer solstice (June 21st), the water is kept entirely within the original watershed. The Bérchules watershed's channel system consists of 19 major channels with a total length of 57.45 km (Fig. 1D). Around 21 km of these channels are exclusively dedicated to aquifer recharge. The most important careo channel of the area, the Acequia del Espino, originates in the upper part of the watershed, where a small dyke is located at 1998 m a.s.l. diverts part of the Chico de Bérchules River flow into the recharge channel (Fig. 3A). This very simple hydraulic intervention, known locally as toma, is maintained, and periodically controlled by local irrigation communities (Fig. 3B), who not only ensure that enough water is entering the channel but also have to guarantee the ecological flow of the river (CIS, 2015). This channel has a length of around 7 km, an average width of 1.5 m, and an average slope of 6.8%. Along its path, the channels cross several infiltration zones or simas, characterized by a high infiltration capacity, where water is partially released from the channels for watering the topsoil and recharging the aquifer. At an altitude of 1820 m a.s.l., the channel water reaches its main destination, entering the sima de Bérchules, which consists of a gently sloping, grassy, and highly permeable surface covering an area of 32,300 m 2 (Fig. 3C). All the water that reaches this zone ends up infiltrating entirely into the sima, including peak flows that may reach 390 L/s, which implies infiltration rates of over 1 m/day, as measured during the hydrological year 2014-2015.
Aquifer discharge plays a noteworthy role in the hydrodynamic response of the basin. Groundwater amounts to 95% of the total discharge of the basin , which is reflected by the flat shape of the average hydrograph ( Fig. 2A). An inventory of the existing springs in the Bérchules watershed listed 609 springs (9 springs/km 2 ), 95% of them discharging <0.2 L/s, and most drying up in the dry period (June to August) of the year . Close to 30% of the total basin springs are found in the higher parts of the watershed, at elevations between 2165 and 2790 m a.s.l. (Fig. 4A; González-Ramón et al., 2015). They are associated with changes in the terrain's slope or decreasing thickness of the periglacial deposits (Fig. 4B). For the rest of the springs, there is a clear relationship between their position and that of the careo channels , as they regulate both the spatial distribution and the amount of discharge of the springs in the watershed. Consequently, their use affects the whole system's hydrogeological functioning. Recharging the aquifer using these channels modifies the distribution of discharge zones (e.g., springs), increasing both the inertia of the hydrogeological system and the number of springs, which in turn increases discharge flow on the slopes of the watershed.

Hydrometeorological data
The Narila gauging station, located at the exit of the Bérchules watershed (Fig. 1), provides a daily series of the basin outflow for the period 1970-2015. For the same period, daily precipitation (P) and temperature (T) time series were measured in the meteorological station located in the village of Bérchules. These time series are used to estimate the potential evapotranspiration (PET) time series by applying Hargreaves's method (Hargreaves and Samani, 1982).

Evaluation of the total aquifer recharge by application of the HBV model
The HBV model (Bergström, 1976;Seibert and Bergström, 2021) is a conceptual rainfall-runoff model used for watershed modeling. It is based on the following water balance: where P [LT −1 ] is precipitation in the form of rain and snow, E [LT −1 ] is evapotranspiration, Q calc [LT −1 ] is the basin total runoff, SP [L] is water storage as snow, SM [L] is the volume of water stored as soil moisture, and UZ [L] and LZ [L] represent water storage in the upper and lower aquifer layers, respectively. The upper aquifer layer generates the interflow discharge, whereas the lower aquifer layer produces the groundwater discharge of the system (Fig. 5).
The HBV model provides the daily discharge of a hydrological basin, based on daily precipitation and temperature time series, and calculates the recharge to the saturated zone. Precipitation enters as rain or snow depending on a predetermined threshold temperature. The first processing unit is the snow routine, which simulates the snow accretion and melting processes. It works as intermittent hydrological storage, accumulating snow and letting it enter the system as snowmelt when daily temperatures exceed the freezing point. Rain and snowmelt together enter the soil moisture routine, which distinguishes between soil recharge by direct infiltration (RT) and the contribution to soil moisture potential (Seibert, 2005). Furthermore, the soil moisture routine computes the direct superficial runoff Q S [LT −1 ] and the potential evapotranspiration PET [LT −1 ] as a linear function of the soil moisture (Seibert, 2005). If water flow exceeds the soil moisture capacity, a third processing unit -the groundwater response routine-starts to distribute it between two interconnected reservoirs. The upper one acts like a saturated subsurface layer generating surface and subsurface runoff, Q 0 [LT −1 ] and Q 1 [LT −1 ], respectively (Fig. 5). The lower reservoir imitates an aquifer that is generating groundwater runoff Q 2 [LT −1 ] (Fig. 5). A predefined and constant percolation rate, which coincides in steady-state with the value of Q 2 , limits the amount of water that recharges the aquifer by entering the lower reservoir. Finally, the runoff estimation routine calculates the total runoff at the outflow point of the catchment (Q calc ) by summing up all the runoff components (i.e., Q S , Q 0 , Q 1, and Q 2 ). This work uses the HBV model adapted by Jódar et al. (2018) to simulate the hydrodynamic response of the Bérchules watershed. To this end, the availability of the meteorological data (i.e., P, T, and PET) and the Berchules River runoff time series is required. All the time series of the applied model cover the period 1970-2015.
The HBV model provides an output time series for several variables, including water infiltration, which is assumed here as total groundwater recharge (RT). As can be seen in Fig. 5, water infiltration is the water flux added to the upper groundwater box. A detailed description of the HBV model can be found in Bergström (1992Bergström ( , 1995 and Seibert (1999).

Evaluation of careo channels' contribution to total groundwater recharge
Flow rates in five channels -including the Trevélez and Mecina interbasin transfer channels, Acequia Real and Acequia Nueva irrigation channels, and Acequia del Espino careo channel-were measured during the 2014-2015 hydrological year to assess the contribution of the recharge channel system to the total aquifer recharge in the Bérchules watershed.
In addition to periodical manual flow measurements in these five channels, the total infiltration capacity along the Acequia del Espino was characterized on the 28th of April 2015, and the careo recharge campaign of this channel was monitored during its operation in the hydrological year 2014-2015. For this purpose, two gauging stations were installed (Fig. 6A) at the beginning (i.e., toma; control point 2) and the endpoint (i.e., Sima de Bérchules; control point 13) of the channel. Both stations were equipped with an Odyssey™ automatic capacitive probe that measured the water column height in each section with an hourly frequency. In addition, the water flow rate was measured at both channel points with an approximately monthly frequency over one year, using a C2-OTT® small current meter, while applying the area-velocity method to measure the flow rate through each channel section. Monthly flow rate data were later correlated with the hourly water column height to derive hourly flow measures for each section.
The infiltration capacity of the El Espino channel was determined by gauging the flow rate in thirteen sections along the channel (control points in Fig. 6A). Fig. 6B represents the flow rate measured for these sections on the 28th of April 2015. The most permeable intervals of the channel, in which the infiltration capacity is higher than average, are marked in red in Fig. 6A and C. They correspond to the channel sections where the flow rate variation between two consecutive control points is greatest. Fig. 6B shows the hydrographs of control points 2 and 13 for the hydrological year 2014-2015.
In this work, the flow rate of the careo channel network in the Bérchules basin has been monitored during the hydrological year 2014-2015. This information provided the basis for estimating the volume of water that was recharged by the careo channel network into the basin during one specific hydrological year. With this data, however, it is hard to estimate the interannual variability of the careo recharge and its corresponding uncertainty. To estimate such variability, it is assumed that the annually diverted volume of surface runoff from the river into the careo channel network is proportional to the annual precipitation in the basin. With this assumption, and being known both the annual recharged careo volume (RC ref [L]) and the annual precipitation volume (P Ref [L]) for a given hydrological year (e.g., 2014-2015), it is possible to estimate the annual recharged careo volume time series RC(t) [L] in terms of the available annual precipitation time series P(t) [L] as: where t is the hydrological year to which RC is estimated. Once RC(t) is obtained, the annual natural recharge time series RN (t) [L] can be evaluated as where RT(t) [L] is the aquifer annual recharge volume time series, which is calculated by postprocessing the HBV results.

HBV modeling results
The HBV model adapted by Jódar et al. (2018) for the Bérchules basin allowed us to simulate the hydrodynamic response of the watershed for the time interval 1970-2014 and thus evaluate the average basin discharge for this period as 13.25 hm 3 /year. The seasonal hydrograph at the Narila gauging station (Fig. 2) shows a temporal delay and a higher baseflow than would be expected for a hydrogeological basin developed in materials of low permeability and subjected to an assumingly pluvio-nival hydroclimatic regime, where marked maximums and minimums in river flow rates are expected in spring and summer, respectively (Beckinsale, 2021). The observed hydrological behavior may be due to a high percentage of groundwater discharge in the total discharge of the watershed. The average total recharge of the aquifer resulting from the model application for the period 1970-2014 is 12.59 hm 3 /year (Fig. 7A), which represents 28% of the average rainfall in the basin (667 mm/year) and 95% of the total discharge of the hydrological system . Similar fractions of  groundwater in the total discharge, only slightly lower than those measured here, have been observed in high mountain watersheds in the Andes, where the substrate consists almost entirely of highly permeable, glacially altered materials, with ratios of total groundwater discharge to total stream discharge of 72% (Somers et al., 2019), 77% (Saberi et al., 2019) and 80% (Baraer et al., 2015). Likewise, in the Himalayas, maximum values of 77% are reached (Williams et al., 2016). These high groundwater contributions to the total streamflow have been also observed in high mountain karst aquifers, which are known to store and transmit large volumes of water (Somers and Mckenzie, 2020). In this line, Chen et al. (2018) reported a contribution of 87% in the northern Alps on the Germany/ Austria border, whereas Jódar et al. (2020) reported a contribution of 75% in the Ordesa and Monte Perdido National Park, which constitutes the highest karst system in Western Europe. Contrary to what happens in these karst systems, the permeable outcrops in the Bérchules watershed correspond almost exclusively to weathered materials (Fig. 4B), as most of the area consists of schists. It seems clear that the elevated groundwater component of the mean hydrograph (Fig. 2) is partially driven by the artificial recharge conducted with the careo channels, as pointed out by Barberá et al. (2018), who analyzed the isotope content of precipitation, snowmelt, groundwater, and surface water in the Bérchules Basin.
Rainfall was 509 mm in the hydrological year 2014-2015. This value represents a reduction of 24% relative to the mean annual rainfall in the Bérchules watershed for the period 1970-2014 while having an exceedance probability of 64% for the annual precipitation for the same period. Similarly, the total basin discharge for this period is 5.3 hm 3 , which is a 60% reduction of the average annual discharge for the 1970-2014 period. The lower discharge value highlights the hydro-climatic stress experienced in the watershed during the 2014/15 hydrological year. Still, the HBV model calculates a recharge for this hydrological year of 112 mm (7.62 hm 3 ), which is a value below the average recharge value (Fig. 8) and corresponds to 22% of the precipitation recorded for this period. Deserving mention is the fact that practically all the meteoric water generated in the headwaters of the Bérchules River was captured and infiltrated by the careo activities during the 2014-2015 hydrological year.

Artificial recharge via careo channels
During its operation time between April 1st and May 15th, 2015, the El Espino channel captured a total volume of 1.99 hm 3 of water from the Grande de Bérchules River. Of this volume, 82% (1.64 hm 3 ) infiltrated along the course of the channel, while the remaining 18% (0.34 hm 3 ) infiltrated at the end of the channel into the permeable materials of the sima de Bérchules (Fig. 6A).
The total volume of water transferred from the Trevélez to the Bérchules watersheds via the Trevélez-Juviles interbasin transfer channel (Fig. 4A) was 0.06 hm 3 during the 2014-2015 hydrological year. This volume was infiltrated completely through the transfer channel into the Bérchules basin. While this relatively low volume was due to leakages in several of its sections, the Mecina channel (Fig. 4A) transferred 1.2 hm 3 of water from the Bérchules River to the Mecina watershed during the same period. According to local agreements between the Bérchules and the Mecina watershed communities, not all the water diverted into the Mecina interbasin transfer channel is transferred to the neighboring Mecina watershed. A part of the diverted water is recharged in the Bérchules watershed. This volume amounted to 0.57 hm 3 during the 2014-2015 hydrological year. For the same period, the total water volume diverted from the Bérchules River into the two main irrigation channels Acequia Real and Acequia Nueva (Fig. 4A) amounted to 1.44 hm 3 . Unlike the water channels exclusively dedicated to aquifer recharge, part of the water flow in these two irrigation channels is for orchards, which cover an area of 0.88 km 2 , thus representing 1.3% of the basin surface . Considering that the irrigation season runs from July to September, and assuming a mean evapotranspiration rate of the orchard crops of 5 mm/day , the corresponding evapotranspiration loss during the irrigation season would be 0.40 hm 3 , and the estimate of total aquifer recharge with these channels would be 1.04 hm 3 in the period 2014-2015.
On a whole, the contribution of the five studied recharge channels -El Espino (1.99 hm 3 ), Trevélez (0.06 hm 3 ), Mecina (0.57 hm 3 ), Acequia Real, and Acequia Nueva (1.04 hm 3 )-to the total aquifer recharge in the Bérchules watershed during the hydrological year 2014-2015 was 3.66 hm 3 , which is equivalent to 70% of the river water flow discharge at the Narila gauging station (5.3 hm 3 ) and constitutes 48% of the total aquifer recharge estimated for this period (7.62 hm 3 ). The natural (i.e., diffuse) recharge of the aquifer thus accounts for 52% of the total recharge (Fig. 7). Furthermore, it should be emphasized that this ancestral form of water "sowing" increases the total groundwater storage of the aquifer, which in turn enhances water "harvest" for the following hydrological year(s).
The mean careo and natural recharge for the period 1970-2015 are 4.8 ± 2.0 hm 3 /yr and 7.4 ± 10.6 hm 3 /yr, respectively. In both cases, the variation interval is expressed as the standard deviation of their corresponding time series. As can be shown, the mean annual recharge and annual variability associated with the careo activities are 1.5 and 5.3 times lower, respectively than those associated with the natural annual recharge (Fig. 7A). Nevertheless, the median (i.e., percentile 50%) of the annual recharge associated with both the careo and natural recharges is 4.1 hm 3 /yr. The disparity between the average and median values associated with the natural recharge is conditioned by the relationship between the natural and total recharge (Eq. (3); Fig. 1-Suppl. Mat), which in turn is driven by the linear dependence between the careo recharge and precipitation (Eq. (2)). As a result, the larger the annual precipitation the larger the difference between the natural and careo recharges (Fig. 8). In any case, the aquifer recharge driven by the careo channels seems to essentially change the proportions of focused and diffuse recharge. In the case of low to mean annual precipitations, both the careo and natural recharges play a similar role concerning the total annual recharge, as can be shown in Fig. 7B. As the careo channels mostly work during the thawing season, this result reveals the high efficiency of the careo channels driving the aquifer recharge. However, these results should be considered as a first approach only, given that they rely on an assumption (Eq. (2)), to be confirmed in future research. Measuring groundwater recharge is not easy and very much determined by local conditions. There are different techniques, including lysimeters, environmental tracers, and water table fluctuations (Custodio, 2019), which are not easily applied to a high mountain basin context. Moreover, recharge may include different processes such as diffuse and focused recharge or mountain system recharge (Meixner et al., 2016), which occur at different scales in time and space, and may not be adequately delineated by point measurements. The role of such hydrological processes in groundwater recharge is usually estimated by numerical modeling. In this line, Somers et al. (2018) and LaFevor and Ramos-Scharrón (2021) analyzed the effects of hillslope trenching on surface water infiltration, the former in the nonglacierized basin of the Shullcas River Watershed in the Cordillera Central of Perú, and the latter in several subalpine forested catchments of Mexico, where the mean precipitation is 800 and 827 mm/yr, respectively. In both cases, the origin of the infiltrated water through the dug trenches was runoff from rainfall, which is produced in several rainfall events distributed over the rainy season. The increases in total infiltration driven by the channels were modest, varying between 1.2 and 2.5% of the total recharge. The discrepancy with the results obtained in the Bérchules basin, where the mean precipitation is 677 mm/yr, and the careo recharge represents an increment of 65% concerning the natural recharge, may rely on the origin of the meteoric water being infiltrated. In the case of the Bérchules basin, most of the infiltrated water through the careo channels corresponds to snowmelt, which provides a continuous source of water to be infiltrated along the thawing season. As a result, during the careo period, the concentrated recharge along the channel network allows less time for evapotranspiration than that associated with intermittent precipitation generating diffuse natural recharge. Besides, the water flowing through the careo channels generate a saturated bulb under the channels that maximizes the infiltration, and hence groundwater recharge (Appels et al., 2015). Moreover, during the careo period, a daily inspection and maintenance work of the channel network is carried out to guarantee its operation without downtime, preventing the bottom of the ditches from becoming clogged and losing infiltration capacity, as often happens in ditches and infiltration ponds.
The data presented in this work confirms that using the ancient careo channels may lead to improvements in water resources management. Similar aquifer recharge systems can be found in the Andes (South America) where channels (amunas) are in use since before the arrival of the Spanish in the 15th century (Escolero et al., 2017;Ochoa-Tocachi et al., 2019). Reusing, and replicating these WSHS would allow quality water recharge in many aquifers sustainably, at least complementing and even replacing water management grey infrastructures. This will improve the availability of local water resources from both quantitative and qualitative perspectives. Moreover, these ancient WSHS could be of interest for many Mediterranean high mountain watersheds, such as those found in the Pyrenees (Spain, France), the Atlas Mountains (Morocco), Mount Etna (Italy), the Dinaric Alps (Croatia), the Taurus Mountains (Turkey) or Mount Lebanon (Lebanon), especially as efficient adaptation measures to mitigate the forecasted impact of climate change on water resources (Abd-Elmabod et al., 2020;Haro-Monteagudo et al., 2020) and on ecosystems (Peñuelas et al., 2018).
The ancient WSHSs were developed in hillslope areas with a moderately permeable substrate, where exist both semiarid climate and high interannual variability of precipitation. To replicate the ancient WSHS successfully in other areas in the world more research is needed, especially in characterizing the hydro-geo-morphological variables (e.g., basin concavity profile and the nature and variability of the soil bedrock interface) that control the efficiency of the careo recharge at the hillslope scale. This knowledge is still lacking and might be part of future research endeavors in different basins of Sierra Nevada where the practice of groundwater recharge with careo channels is still alive (Fig. 1). At this point, it is worth stressing that to ensure the successful replication of these ancient WSHS, the technical details are as important as the social aspects, since there must be a community that needs such water resources surplus and is willing to build and maintain the channel network. (Fig. 1).

Conclusions
The response of the mean annual hydrograph of the Bérchules River is not typical of a high mountain river, which should show a nival or at least pluvio-nival hydrological regime. Instead, it exhibits the characteristics of a river connected to an inertial (slow response time) aquifer. The Bérchules watershed sits on metamorphic rocks, mainly schists. Nevertheless, these materials became permeable as a result of weathering and glacial and periglacial alteration processes, which allowed the formation of a shallow aquifer. The water from artificial aquifer recharge plus irrigation returns, both achieved by capturing headwaters from the river with careo channels, end up slowly discharging into the river. As a result of the recharge process, the snow flow peak is eliminated from the hydrograph of the Bérchules or other watersheds in the Sierra Nevada basins where this form of water management is practiced.
The obtained results demonstrate that aquifer recharge occurs through two different processes in the Bérchules watershed: (1) a natural and spatially distributed recharge associated with infiltration of meteoric waters (rainwater and snowmelt), and (2) an artificial and concentrated recharge along the course of the careo channels, with irrigation and inter boundary transfer channels that capture and transport river waters for this purpose. For the hydrological year 2014-2015, distributed natural recharge and artificial recharge respectively accounted for 52% and 48% of the total recharge of the hydrogeological system.
Although this study is undertaken in a particular watershed in the southern edge of Sierra Nevada (Spain), the methods are generically applicable to different basins in other geographical settings where dug channels are used to recharge groundwater. Further study is required to characterize the behavior of such recharge systems in other geographical conditions as they may be an efficient adaptation measure against the impact of climate change.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.