Etienne Marti

Sarah Leray
email: saleray@ing.puc.cl

Daniela Villela

José Maringue

Gonzalo Yáñez

Esteban Salazar

Fernando Poblete

José Jimenez

Gabriela Reyes

Guillermo Poblete

Zeidy Huamán

Ronny Figueroa

Jaime Araya Vargas

Jorge Sanhueza

Marjorie Muñoz

Reynaldo Charrier

Gabriel Fernández

Unravelling geological controls on groundwater flow and surface water-groundwater interaction in mountain systems: A multi-disciplinary approach

Keywords: Structural geology, UAV-based magnetic, hydrogeophysics, fault zone, mountain aquifer, groundwater-surface water interaction

come   L'archive ouverte pluridisciplinaire

Introduction

Interest in mountain catchments has been growing recently (Markovich et al., 2019 and references within). Mountains serve as water reservoirs for downstream alluvial basins due to the heavy precipitation they receive [START_REF] Immerzeel | Importance and vulnerability of the world's water towers[END_REF][START_REF] Somers | A review of groundwater in high mountain environments[END_REF]. More than half of humanity depends on mountain hydrosystems for their water supply [START_REF] Grêt-Regamey | Mountain Ecosystem Services: Who Cares? Mountain Research and Development[END_REF].

Snow-and glacier-dominated hydrosystems supply more than one-sixth of the world's population Unravelling geological controls on groundwater flow and surface water-groundwater interaction in mountain systems: a multi-disciplinary approach 4/42 [START_REF] Barnett | Potential impacts of a warming climate on water availability in snow-dominated regions[END_REF]. Human dependence on these systems is all the more critical as they are particularly vulnerable to climate change [START_REF] Beniston | Climatic change in mountain regions: A review of possible impacts[END_REF][START_REF] Huss | Toward mountains without permanent snow and ice[END_REF][START_REF] Viviroli | Climate change and mountain water resources: overview and recommendations for research, management and policy[END_REF].

Understanding the functioning of mountain blocks and their hydraulic connectivity with anthropized basins downstream is thus critical for water resource management [START_REF] Kamenik | The influence of catchment characteristics on the water chemistry of mountain lakes[END_REF][START_REF] Mcdonnell | Where does water go when it rains? Moving beyond the variable source area concept of rainfall-runoff response[END_REF][START_REF] Soulsby | Inferring groundwater influences on surface water in montane catchments from hydrochemical surveys of springs and streamwaters[END_REF], as well as for society's resilience to climate change [START_REF] Hayashi | Alpine Hydrogeology: The Critical Role of Groundwater in Sourcing the Headwaters of the World[END_REF][START_REF] Viviroli | Increasing dependence of lowland populations on mountain water resources[END_REF]. Furthermore, mountain catchments constitute preserved spaces and biosphere reserves, in some cases harbouring endemic or native species [START_REF] Clarke | Macroinvertebrate diversity in headwater streams: a review[END_REF][START_REF] Freeman | Hydrologic connectivity and the contribution of stream headwaters to ecological integrity at regional scales[END_REF][START_REF] Meyer | The contribution of headwater streams to biodiversity in river networks[END_REF] and providing various ecosystem services [START_REF] Egan | Mountain ecosystem services and climate change: A global overview of potential threats and strategies for adaptation[END_REF].

However, despite their importance, mountain hydrogeological systems remain poorly understood [START_REF] Bishop | Aqua Incognita: the unknown headwaters[END_REF] owing to their inherent complexity stemming from significant topographic, climatic, and geological gradients. In particular, the interplay of geology, geomorphology, and climate on groundwater flow and related variables is barely understood [START_REF] Condon | Evaluating the relationship between topography and groundwater using outputs from a continental-scale integrated hydrology model[END_REF][START_REF] Gleeson | Regional groundwater flow in mountainous terrain: Threedimensional simulations of topographic and hydrogeologic controls[END_REF][START_REF] Winter | The concept of hydrological landscapes[END_REF]: for example, groundwater-flow scale, depth, distribution [START_REF] Markovich | Mountain-Block Recharge: A Review of Current Understanding[END_REF][START_REF] Roques | Hydrological behavior of a deep sub-vertical fault in crystalline basement and relationships with surrounding reservoirs[END_REF], and recharge [START_REF] Mccord | Impact of geologic heterogeneity on recharge estimation using environmental tracers: Numerical modeling investigation[END_REF][START_REF] Nicolas | Managed Aquifer Recharge in fractured crystalline rock aquifers: Impact of horizontal preferential flow on recharge dynamics[END_REF][START_REF] Sanford | Recharge and groundwater models: an overview[END_REF][START_REF] White | Evidence for structurally controlled recharge in the Blue Ridge Province, Virginia, USA[END_REF], groundwater-surface water interaction [START_REF] Brunner | Hydrogeologic controls on disconnection between surface water and groundwater[END_REF][START_REF] Fleckenstein | Groundwater-surface water interactions: New methods and models to improve understanding of processes and dynamics[END_REF], groundwater timescales [START_REF] Jasechko | Substantial proportion of global streamflow less than three months old[END_REF], and river ecology and biogeochemistry [START_REF] Kolbe | Coupling 3D groundwater modeling with CFC-based age dating to classify local groundwater circulation in an unconfined crystalline aquifer[END_REF][START_REF] Weyer | Catchments as heterogeneous and multi-species reactors: An integral approach for identifying biogeochemical hot-spots at the catchment scale[END_REF]. The geological heterogeneity of mountain blocks related to multiscale fracturing, lithological discontinuities, weathering, etc., is likely to substantially control groundwater flow as significant gradients of the hydraulic conductivity induce flow anisotropy [START_REF] Leray | Numerical modeling of the productivity of vertical to shallowly dipping fractured zones in crystalline rocks[END_REF][START_REF] Leray | Contribution of age data to the characterization of complex aquifers[END_REF], flow channelling, intercatchment flow [START_REF] Fan | Are catchments leaky?[END_REF], and possibly groundwater exfiltration [START_REF] Bense | Fault zone hydrogeology[END_REF][START_REF] Caine | Fault zone architecture and permeability structure[END_REF][START_REF] Scibek | The biases and trends in fault zone hydrogeology conceptual models: global compilation and categorical data analysis[END_REF].

Although studies indicate that geological structures substantially impact groundwater-surface water interaction [START_REF] Schaller | River basins as groundwater exporters and importers: Implications for water cycle and climate modeling[END_REF] and stream flowrate and regime, most studies have focused on the role of climate and geomorphology [START_REF] Condon | Evaluating the relationship between topography and groundwater using outputs from a continental-scale integrated hydrology model[END_REF][START_REF] Cuthbert | Global patterns and dynamics of climate-groundwater interactions[END_REF][START_REF] Gleeson | Regional groundwater flow in mountainous terrain: Threedimensional simulations of topographic and hydrogeologic controls[END_REF][START_REF] Welch | Consistency of groundwater flow patterns in mountainous topography: Implications for valley bottom water replenishment and for defining groundwater flow boundaries[END_REF]. The geological complexity of hydrogeological systems Unravelling geological controls on groundwater flow and surface water-groundwater interaction in mountain systems: a multi-disciplinary approach 5/42 strongly impedes the generalization of in situ observations and site-specific results and accordingly limits the development of unified hydrogeological frameworks [START_REF] Markovich | Mountain-Block Recharge: A Review of Current Understanding[END_REF][START_REF] Wilson | Mountain-block hydrology and mountain-front recharge[END_REF]. Worse, mountain hydrogeological systems remain mainly uncharacterized as their steep topographic gradients and remoteness strongly limit access and site instrumentation. The scientific bottleneck regarding the understanding of the hydrogeological functioning of mountain groundwater systems is therefore above all a methodological gap, the filling of which requires the improved characterization of their subsurface geologic framework [START_REF] Bense | Fault zone hydrogeology[END_REF][START_REF] Markovich | Mountain-Block Recharge: A Review of Current Understanding[END_REF].

Groundwater flow paths and origins in mountain catchments have traditionally been inferred using hydrogeochemistry, water stable isotopes, and environmental tracers such as noble gases and anthropogenic tracers [START_REF] Alvarez-Campos | Evidence for high-elevation salar recharge and interbasin groundwater flow in the Western Cordillera of the Peruvian Andes[END_REF][START_REF] Baraer | Contribution of groundwater to the outflow from ungauged glacierized catchments: a multi-site study in the tropical Cordillera Blanca, Peru[END_REF][START_REF] Hilberg | Review: Natural tracers in fractured hard-rock aquifers in the Austrian part of the Eastern Alps-previous approaches and future perspectives for hydrogeology in mountain regions[END_REF][START_REF] Jung | Assessing seasonal variations in water sources of streamflow in a temperate mesoscale catchment with granitic bedrocks using hydrochemistry and stable isotopes[END_REF][START_REF] Ma | Groundwater mixing and mineralization processes in a mountainoasis-desert basin, northwest China: hydrogeochemistry and environmental tracer indicators[END_REF][START_REF] Manning | An integrated environmental tracer approach to characterizing groundwater circulation in a mountain block[END_REF]. As they directly derive from flow processes, these methods can theoretically compensate for the lack of direct observations of hydraulic conductivity and hydraulic gradients from standard hydraulic methods (boreholes, pumping tests) [START_REF] Ginn | Notes on Groundwater Age in Forward and Inverse Modeling[END_REF].

However, because of water mixing processes, they constitute integrative methods that remain poorly effective to differentiate between hydrogeological models at the local scale [START_REF] Troldborg | The importance of alternative conceptual models for simulation of concentrations in a multi-aquifer system[END_REF] and to identify hydrogeologically pertinent structures in complex systems [START_REF] Leray | Contribution of age data to the characterization of complex aquifers[END_REF]. Hence, in the absence of end-member waters with significantly different chemical signatures, the deconvolution of tracers and hydrogeochemical data is extremely challenging.

Water balance analysis is a straightforward method that has been mainly used for the estimation of water balance components and for closing the water balance [START_REF] Andermann | Impact of transient groundwater storage on the discharge of Himalayan rivers[END_REF][START_REF] Cochand | Groundwater Storage in High Alpine Catchments and Its Contribution to Streamflow[END_REF][START_REF] Hood | Characterization of snowmelt flux and groundwater storage in an alpine headwater basin[END_REF][START_REF] Mcclymont | Groundwater flow and storage within an alpine meadow-talus complex[END_REF]. However, closing the water balance may be strongly impaired in small ungauged catchments where the large amount of data required is not satisfied as ground data are rarely available and remote sensing data are often too coarse [START_REF] Acharya | Unmanned Aerial Vehicles in Hydrology and Water Management: Applications, Challenges, and Perspectives[END_REF][START_REF] Hobouchian | Assessment of satellite precipitation estimates over the slopes of the subtropical Andes[END_REF][START_REF] Longuevergne | GRACE Hydrological estimates for small basins: Evaluating processing approaches on the High Plains Aquifer, USA[END_REF][START_REF] Zambrano-Bigiarini | Temporal and spatial evaluation of satellite-based rainfall estimates across the complex topographical and climatic gradients of Chile[END_REF]. Besides, addressing the role of geological structures on groundwater flow and surface water-groundwater Unravelling geological controls on groundwater flow and surface water-groundwater interaction in mountain systems: a multi-disciplinary approach 6/42 interaction implicitly requires the selection of characterization methods that are able to image geological discontinuities, that is, those with a relatively fine resolution (tens of meters). In this context, geophysics has emerged as a non-destructive method to survey the subsurface at various scales and depths and in some occasions hydrogeological processes [START_REF] Binley | The emergence of hydrogeophysics for improved understanding of subsurface processes over multiple scales[END_REF][START_REF] Robinson | Advancing processbased watershed hydrological research using near-surface geophysics: a vision for, and review of, electrical and magnetic geophysical methods[END_REF]. Electrical Resistivity Tomography (ERT) has been widely used [START_REF] Samouelian | Electrical resistivity survey in soil science: a review[END_REF] for soil moisture studies [START_REF] Vargas | Soil electrical resistivity monitoring as a practical tool for evaluating irrigation systems efficiency at the orchard scale: Unravelling geological controls on groundwater flow and surface water-groundwater interaction in mountain systems: a multi-disciplinary approach 27/42 a case study in a vineyard in Central Chile[END_REF], the identification of groundwater aquifer geometry, and the study of geological structures [START_REF] Ball | Constraining fault-zone hydrogeology through integrated hydrological and geoelectrical analysis[END_REF][START_REF] Martos-Rosillo | The oldest managed aquifer recharge system in Europe: New insights from the Espino recharge channel (Sierra Nevada, southern Spain)[END_REF][START_REF] Robert | Assessing the contribution of electrical resistivity tomography (ERT) and self-potential (SP) methods for a water well drilling program in fractured/karstified limestones[END_REF].

ERT cross-sections can be aligned to reveal resistivity contrasts resulting from geological structures and to construct a 3D representation based on previous geological knowledge [START_REF] Figueroa | Deciphering groundwater flow-paths in fault-controlled semiarid mountain front zones (Central Chile)[END_REF].

Gravity surveys can be performed at different spatial scales, from local to regional, and can identify tiny changes in gravity to obtain structural information [START_REF] Telford | Applied Geophysics[END_REF]. Gravity surveys are generally used to delineate the basement depth of sedimentary basins [START_REF] Abbott | Depth to bedrock using gravimetry in the Reno and Carson City, Nevada, area basins[END_REF][START_REF] Figueroa | Deciphering groundwater flow-paths in fault-controlled semiarid mountain front zones (Central Chile)[END_REF][START_REF] Murty | The gravity method in groundwater exploration in crystalline rocks: a study in the peninsular granitic region of Hyderabad, India[END_REF][START_REF] Yáñez | Gravity derived depth to basement in Santiago Basin, Chile: implications for its geological evolution, hydrogeology, low enthalpy geothermal, soil characterization and geo-hazards[END_REF]; however, they can also be used to identify geological contacts and discontinuities that do not outcrop [START_REF] Ruelleu | Impact of gently dipping discontinuities on basement aquifer recharge: An example from Ploemeur (Brittany, France)[END_REF]. While they can be used to image geological structures, the abovementioned methods and many more geophysical methods are unable to image hydrogeological features without an additional hydraulic-based constraint such as groundwater level, water resurgences, etc. Besides, ground geophysical surveys of mountain catchments are limited by their remoteness, which has raised interest in the use of aerial surveys. So far, aerial geophysical surveys have been conducted by helicopter or airplane for largescale studies (>1,000 km 2 ) focusing on major regional geological structures (Dumont et al., 2021;[START_REF] Hussien | An integrated approach for identification of potential aquifer zones in structurally controlled terrain: Wadi Qena basin, Egypt[END_REF][START_REF] Lopez Loera | Geophysical Exploration of Fractured-Media Aquifers at the Mexican Mesa Central: Satellite City, San Luis Potosi, Mexico[END_REF][START_REF] Reeves | Aeromagnetic Surveys: Principles, Practice and Interpretation[END_REF][START_REF] Vittecoq | Hydrogeological conceptual model of andesitic watersheds revealed by high-resolution heliborne geophysics[END_REF]). Yet, the recent emergence of drones (UAV) has offered a potential means to explore the subsurface framework at much smaller scales and more systematically, thanks to their smaller clearance distances (less than 50 m) and closer profiles.

Unravelling geological controls on groundwater flow and surface water-groundwater interaction in mountain systems: a multi-disciplinary approach 7/42

Hydrogeology is fundamentally at the cross-roads of various disciplines and multidisciplinary approaches have naturally shown their value for characterizing, understanding and managing groundwater as the joint use of various methods makes it possible to more widely cover the parameters space [START_REF] Kreamer | Chapter 37 -The future of groundwater science and research[END_REF]. Regarding mountain catchments, studies based on various methods became essential to develop a more in-depth knowledge of hydrogeological systems and to propose new conceptual models of their hydrogeological functioning [START_REF] Christensen | Hydrogeological characterization of an alpine aquifer system in the Canadian Rocky Mountains[END_REF][START_REF] Manning | Direct Observation of the Depth of Active Groundwater Circulation in an Alpine Watershed[END_REF][START_REF] Markovich | A Mountain-Front Recharge Component Characterization Approach Combining Groundwater Age Distributions, Noble Gas Thermometry, and Fluid and Energy Transport Modeling[END_REF][START_REF] Ruiz-Pereira | Hydrological connections in a glaciated Andean catchment under permafrost conditions (33°S)[END_REF][START_REF] Somers | Quantifying groundwater-surface water interactions in a proglacial valley, Cordillera Blanca, Peru[END_REF]. However, the choice of data that are informative and complementary enough and the interaction between disciplines that may not share the same approaches and perceptions are in itself challenging. The main objective of this study is to develop a novel multi-disciplinary approach to assess the hydrogeological role of geological features at the catchment scale in ungauged mountain catchments. To do so, we perform a thorough geological characterization through a multiscale geological mapping and ground and aerial geophysical instrumentation to accurately outline both the surface and the subsurface geological framework. To clearly establish the impact of geological features and heterogeneity on hydrogeological structures, hydrogeological information is retrieved from various sources and various data types, including spatially-distributed baseflow measurements, the Normalized Difference Vegetation Index (NDVI), and a thermal spring location and temperature. The multidisciplinary methodology is designed to get a comprehensive understanding of the system and establish a causal link between geological and hydrogeological features which cannot be achieved by each of the methods separately. The scope of this study is functional and pragmatical aiming to identify the key geological structures for groundwater flow but not the entire geological framework and focusing on a key hydrological variable which is the streamflow.

In the following section, we present the application of the methodology to the Parque Nacional del Río Clarillo, an ungauged mountain catchment in Chile whose large geological heterogeneity is illustrative of both the scientific and methodological challenges to the analysis of mountain groundwater systems.

Unravelling geological controls on groundwater flow and surface water-groundwater interaction in mountain systems: a multi-disciplinary approach 8/42 After presenting the results in section 3, we discuss the capacity of various hydrogeological conceptual models to explain the observed streamflow variations and hydric singularities as well as the differentiating aspects and novelty of the proposed approach. The Parque Nacional del Río Clarillo is an ungauged catchment with an area of 130 km² located 35 km southeast of Santiago (-33.72°, -70.49°, Figure 1 a andb). Situated in a semi-arid climate, the park is suffering from the uninterrupted drought that has been affecting Central Chile since 2010 and had led to a decrease in precipitation of about 30% in this region [START_REF] Garreaud | The Central Chile Mega Drought (2010-2018): A climate dynamics perspective[END_REF]. Though presenting significant interannual variability, the average annual precipitation downstream of the park historically amounted to approximately 700 mm during the period [START_REF] Telford | Applied Geophysics[END_REF][START_REF] Kaminsky | ZondRes2D. Software for two-dimensional interpretation of DC resistivity and IP data[END_REF][START_REF] Dmc | Hydrogeophysical Characterization in a Volcanic Context From Local to Regional Scales Combining Airborne Electromagnetism and Magnetism[END_REF]. However, since the beginning of the drought, it very often does not exceed 300 mm and reached 400 mm on very few occasions [START_REF] Dmc | Hydrogeophysical Characterization in a Volcanic Context From Local to Regional Scales Combining Airborne Electromagnetism and Magnetism[END_REF]. Accordingly, observations by park rangers and technical studies [START_REF] Niemeyer | Historia natural de la Reserva Nacional Río Clarillo: un espacio para aprender ecología[END_REF] point at a severe reduction of the Río Clarillo baseflow of close to 75%. Still, these observations should be interpreted with caution as they remain very sparse. The Río Clarillo regime is pluvionival with streamflow maintained by snowmelt in the upper part of the catchment. While snowmelt historically lasted up to the end of March, the past few years have shown an earlier melting, which rarely persists after January according to observations by park ranger.

Study site and methodology

Geological context

The present-day configuration of the Andes of Central Chile is marked by the development of two margin-parallel ranges, the Coastal Cordillera to the west and the Principal Cordillera to the east, separated by the Central Depression alluvial plain. The study area is located on the western margin of approach 9/42 the Principal Cordillera (Figure 1a), and is uplifted by more than 2,000 m over the Central Depression by the cumulative displacement of the west-vergent San Ramón-Pocuro Fault System (SRPFS) [START_REF] Armijo | The West Andean Thrust, the San Ramón Fault, and the seismic hazard for Santiago, Chile[END_REF][START_REF] Rauld | Análisis morfoestructural del frente cordillerano Santiago Oriente entre el río Mapocho y la Quebrada de Macul[END_REF].

The geology of the Río Clarillo catchment records the evolution of the western margin of the Eocene-Miocene Abanico Basin (Figure 1c). This includes the volcanic and sedimentary rocks of the extensional phase of the Abanico Basin, the Abanico Formation (Oligocene-Miocene), and the volcanic rocks that represent the tectonic inversion of the basin, that is, the Miocene Farellones Formation [START_REF] Charrier | Evidence for Cenozoic extensional basin development and tectonic inversion south of the flat-Unravelling geological controls on groundwater flow and surface water-groundwater interaction in mountain systems: a multi-disciplinary approach 28/42 slab segment, southern Central Andes, Chile[END_REF]. These stratified units are intruded by granitic rocks of the La Obra pluton (c.a. 21 Ma), (Figure 1c). However, at the scale of this work, there is no consensus on whether the rock bodies outcropping in the Río Clarillo catchment should be assigned to these formal units [START_REF] Fock | Cronología y tectónica de la exhumación en el neógeno de Los Andes[END_REF][START_REF] Piquer | The Piuquencillo fault system: a long-lived, Andeantransverse fault system and its relationship with magmatic and hydrothermal activity[END_REF][START_REF] Piquer | Along-strike segmentation of the Abanico Basin, central Chile: New chronological, geochemical and structural constraints[END_REF][START_REF] Thiele | Hoja Santiago: región metropolitana: carta geológica de Chile escala 1: 250.000[END_REF]. At the regional scale, the Abanico Formation shows intense deformation of tight west-vergent folds in the westernmost outcrops that gradually decreases towards the east to wide and open folds. The Farellones Formation lies unconformably over the tight folds affecting the Abanico Formation in the westernmost outcrops, showing less intense folding, and becomes gradually parallel to the underlying Abanico formation towards the east [START_REF] Charrier | Evidence for Cenozoic extensional basin development and tectonic inversion south of the flat-Unravelling geological controls on groundwater flow and surface water-groundwater interaction in mountain systems: a multi-disciplinary approach 28/42 slab segment, southern Central Andes, Chile[END_REF][START_REF] Fock | Cronología y tectónica de la exhumación en el neógeno de Los Andes[END_REF].

Scarce Pliocene to Holocene alluvial and fluvial deposits are preserved in the main fluvial valley and in tributary ravines. From the mouth of the Río Clarillo valley, alluvial deposits show a strong increase in thickness towards the Central Depression, reaching a thickness of up to 500 m [START_REF] González | Characterization of the depocenters and the basement structure, below the central Chile Andean Forearc: A 3D geophysical modelling in Santiago Basin area[END_REF][START_REF] Yáñez | Gravity derived depth to basement in Santiago Basin, Chile: implications for its geological evolution, hydrogeology, low enthalpy geothermal, soil characterization and geo-hazards[END_REF].

Geological methodology

The geological survey consisted of a regional geological mapping (1:100,000 scale) of the whole Río Clarillo catchment (Figure 1c), refining the previous geological maps in the area. To assess controls on primary permeability, since there is no consensus on the assignment of these rock units to the formal Unravelling geological controls on groundwater flow and surface water-groundwater interaction in mountain systems: a multi-disciplinary approach 10/42 lithostratigraphic units in the area [START_REF] Fock | Cronología y tectónica de la exhumación en el neógeno de Los Andes[END_REF][START_REF] Piquer | The Piuquencillo fault system: a long-lived, Andeantransverse fault system and its relationship with magmatic and hydrothermal activity[END_REF][START_REF] Thiele | Hoja Santiago: región metropolitana: carta geológica de Chile escala 1: 250.000[END_REF], strictly lithological units were defined and mapped. As for the controls on structural features on secondary permeability, faults, folds, and fracture zones were surveyed at a finer scale, focusing on fracture systems along the riverbed and on hilltops where rock outcrops are observed.

Geophysical methods

The geophysical study involves the joint use of a UAV-based aeromagnetic system (with only limited previous applications in the Andes [START_REF] Sanhueza | Rheological, petrophysical and geometrical constraints of a subduction channel from a numerical model perspective: Insights from La Cabana Paleozoic peridotites, Coastal Cordillera of south-central Chile[END_REF][START_REF] Yáñez | Shallow Anatomy of the San Ramon Fault (Chile) Constrained by Geophysical Methods: Implications for its Role in the Andean Deformation[END_REF]) and on-ground transect data obtained from ERT and gravimetry. As shown in Figure 1c, transect (A-A') and the drone coverage are oriented east-west to perpendicularly intersect the main structures present in the study site as the SRFPS is N-S oriented (see section 2.1.2 for a description of the geological context). As shown in Figure 1c, a second roughly E-W transect (B-B') is also performed. This transect is located 10 km northeast of the Río Clarillo catchment, in the Río Coipo valley, and was performed to gain insights on the geological structures on the eastern flank of the La Obra pluton, where the rough topography prevents the use of ground-based geophysical methods along the primary transect A-A'. Furthermore, having two parallel profiles allows a 3D view of the geology to be obtained (N-S along strike variations). We use the same methodological approach in both transects. The results of this second transect are presented in the supplementary material; nevertheless, their contribution to our geological interpretation of the study area is included in the discussion.

Electrical Resistivity

Electrical Resistivity Tomography is an exploration method that measures the subsurface distribution of electrical resistivity with high resolution [START_REF] Dahlin | The development of DC resistivity imaging techniques[END_REF]. It can be used to detect water-saturated areas and the depth of the contact between the sediment surface layer and the basement unit, as well as to identify fracture and fault zones in the basement unit, which show lower resistivity compared to Unravelling geological controls on groundwater flow and surface water-groundwater interaction in mountain systems: a multi-disciplinary approach 11/42 the surrounding unfractured rock. As the ERT method has optimal definition in the shallow subsurface, it is limited to depths of around 250-300 m, depending on the dipole separation. Data were acquired using an ABEM Terrameter LS2 (Guideline Geo, Stockholm, Sweden). The measurement protocol involved a Dipole-Dipole array with 48 electrodes and a 20 m inter-electrode spacing. Each profile was 1 km long with an investigation depth of about 250-350 m. The Dipole-Dipole array was chosen for its particularly good definition of the vertical structures typically present in faults and geological layers [START_REF] Loke | Tutorial: 2-D and 3-D electrical imaging surveys[END_REF][START_REF] Samouelian | Electrical resistivity survey in soil science: a review[END_REF].

The Río Clarillo and Río Coipo transects are respectively composed of six and seven profiles (the former having fewer profiles due to field conditions). Both profiles present a 500 m overlap between each profile to obtain a better characterization of the structures. The apparent resistivity datasets obtained in the field for each profile were inverted separately to obtain 2D resistivity cross-sections. The data were inverted using the ZonRes2D software [START_REF] Kaminsky | ZondRes2D. Software for two-dimensional interpretation of DC resistivity and IP data[END_REF]. Each profile was then "stitched" to the adjacent one to obtain the full transect using the stitching tool in the Oasis Montaj software (Seequent, Christchurch, New Zealand). This stitching method is necessary to obtain the full-length transect without noticeable edge effects at the extremes of the profiles.

Gravimetry

Gravimetry measures the gravity field caused by differences in mass distribution below the surface.

Variations in the residual gravity anomaly i.e. the difference between properly corrected gravity observations and the regional field (see a description of this process later in this section) can be indirectly associated with the subsurface density distribution, which can in turn be used to infer geological information along the measured profile. We used a CG-5 Autograv relative gravimeter (Scintrex, Concord, ON, Canada) with a resolution of 5 μGal and a Trimble R4 GNSS (Trimble, Westminster, CO, USA) receiver that gives the position and altitude with a precision of less than 30 cm.

The Río Clarillo transect (A-A') has 56 gravimetric stations spaced ≈50 m apart along a 2.5 km long Unravelling geological controls on groundwater flow and surface water-groundwater interaction in mountain systems: a multi-disciplinary approach 12/42 profile (Figure 1c). The Río Coipo transect (B-B') presents 72 gravimetric stations spaced ≈100 m apart along a 7 km long profile (Figure 1c). Gravity data were corrected using standard procedures [START_REF] Telford | Applied Geophysics[END_REF]: (1) The tidal correction is computed by the instrument during all the measurements and the drift correction is estimated as a linear function during the day with measurements at the base station at the beginning and end of the day; (2) The latitude correction was computed using the WGS84 reference ellipsoid; (3) The absolute gravity value was obtained using the base station at the Department of Geophysics of the University of Chile (Universidad de Chile) (Gravity DGF: 979,416.07 [mGal]); (4) The Free-air and Bouguer corrections were computed using the altitude obtained from the differential GNSS and assuming a mean subsurface density of 2.67 g/cm^3; and (5) The Terrain correction was calculated with a highresolution digital terrain model with a 5 m resolution (AW3D Standard product, acquired from NTT DATA Corporation). Finally, the regional effect was removed to obtain the residual anomaly along the profile due to the effect of the low-density bodies under the Principal Cordillera at deeper crustal levels: for this purpose, the Bouguer anomaly of Santiago obtained by [START_REF] Yáñez | Gravity derived depth to basement in Santiago Basin, Chile: implications for its geological evolution, hydrogeology, low enthalpy geothermal, soil characterization and geo-hazards[END_REF] was removed from our data. Then, the result was processed using the rock measurement performed at a point of zero-anomaly as a constraint. After correcting the gravity observations, we obtained a residual gravity anomaly suitable for geologically constrained gravity modelling. The modelling process consisted of a forward modelling using the ModelVision software (Encom), with the model output being compared to the residual gravity anomaly.

UAV-based magnetic survey

In the magnetic method, the magnetic properties of rocks are measured using a magnetometer. It provides a valuable tool to map geology and structural networks if they involve magnetic contrasts.

The operational simplicity of total field magnetometers, in which the magnetic field is recorded almost instantaneously regardless of the sensor or platform orientation, make them very suitable for use in mobile platforms, such as airplanes, ships, or UAVs. We used a Bartington MAG03 fluxgate Unravelling geological controls on groundwater flow and surface water-groundwater interaction in mountain systems: a multi-disciplinary approach 13/42 magnetometer (Bartington Instruments, Oxon, United Kingdom) with three orthogonal components mounted on a UAV system. We surveyed one block (see spatial distribution in Figure 1c, purple polygon). The polygon is 6,600 m long and 1,800 m wide, covering large parts of the reserve and both side contacts of the intrusive unit. Lines were flown in the N-S direction with a line separation of 200 m and a flight altitude of 50-80 m above terrain level. The location was recorded with a standard GPS with an accuracy of 5m.

The following corrections were applied to the raw magnetic field observations: (1) removal of diurnal variations; (2) removal of the Earth's magnetic field (IGRF, International Geomagnetic Reference Field); and (3) recovery of the total magnetic field from the three-component fluxgate magnetometer (for details see [START_REF] Sanhueza | Rheological, petrophysical and geometrical constraints of a subduction channel from a numerical model perspective: Insights from La Cabana Paleozoic peridotites, Coastal Cordillera of south-central Chile[END_REF]). After data correction, we obtained a residual magnetic field dominated by magnetic sources in the shallow crust. The data were then reduced to the magnetic pole (RTP), an approach which aims to centre anomalies over their source [START_REF] Blakely | Potential Theory in Gravity and Magnetic Applications[END_REF]. The magnetic data were then interpreted qualitatively by comparing the RTP map with geological, hydraulic, and additional geophysical information.

Streamflow measurements, remote sensing data, and hydric markers

The streamflow rate was measured from the outlet of the national park up to the confluence of the two main tributaries of the Río Clarillo, 9 km further upstream (Figure 1b). The measurements were carried out at least every kilometre and when possible at smaller intervals in the sectors of interest, that is, where the flow exhibited significant variation. Data were acquired with an Acoustic Doppler Velocimeter (FlowTracker, SonTek). Three campaigns were completed during the 2019 recession period (austral autumn) to better evidence the contribution of groundwater to the baseflow. The first and third campaigns are complete campaign from the outlet P0 to the confluence P8 while the second campaign is simply a small-scale campaign between points P4 and P5. If the streamflow measurements can be quantitatively compared between each campaign, the three campaigns exhibit similar trends Unravelling geological controls on groundwater flow and surface water-groundwater interaction in mountain systems: a multi-disciplinary approach 14/42 and features and can serve as cross-validation. Some data points are missing due to technical, logistical or staff complications on site (e.g., point P5 during the third campaign). We also considered different types of soft hydric data (i.e., qualitative data), which are essential in ungauged catchments [START_REF] Seibert | Gauging the Ungauged Basin: Relative Value of Soft and Hard Data[END_REF]. The soft data were mainly obtained from the knowledge of the nearby community of Pirque, the testimony of park rangers, and the authors' observations during field campaigns.

To complement the streamflow analysis and to obtain a long-term assessment of the streamflow dynamics and the persistence of possible resurgence zones, we analysed NDVI time series. The NDVI measures the difference between infrared and red light to quantify vegetation presence and density using reflectance bands. It is the most common index used in remote sensing to understand drought processes [START_REF] Gu | Evaluation of MODIS NDVI and NDWI for vegetation drought monitoring using Oklahoma Mesonet soil moisture data[END_REF][START_REF] Gu | A five-year analysis of MODIS NDVI and NDWI for grassland drought assessment over the central Great Plains of the United States[END_REF], related plant water stress [START_REF] Ballester | Monitoring the Effects of Water Stress in Cotton using the Green Red Vegetation Index and Red Edge Ratio[END_REF][START_REF] Katsoulas | Crop reflectance monitoring as a tool for water stress detection in greenhouses: A review[END_REF], and to evaluate groundwater fluctuation and vegetation cover [START_REF] Pan | Groundwater level monitoring model using multitemporal images in arid region of northwest China[END_REF][START_REF] Seeyan | Investigation of the Relationship between Groundwater Level Fluctuation and Vegetation Cover by using NDVI for Shaqlawa Basin, Kurdistan Region -Iraq[END_REF]. The NDVI was obtained from the MODIS/MCD43A4 product [START_REF] Schaaf | MCD43A4 MODIS/Terra+Aqua BRDF/Albedo Nadir BRDF Adjusted Ref Daily L3 Global -500m V006[END_REF] through the Google Earth Engine [START_REF] Gorelick | Google Earth Engine: Planetary-scale geospatial analysis for everyone[END_REF] from 24 Feb 2000 to 31 Dec 2019, that is, covering the drought period and the previous decade.

Conceptual hydrogeological models

We develop various conceptual hydrogeological models and discuss hypotheses in relation to the streamflow evolution through the catchment. As argued later in section 3.4, the groundwater-surface water interaction is highly likely to be related to geological discontinuities. Geological controls on groundwater-surface water interaction are further reinforced by a parallel sensitivity study (not shown here) of baseflow evolution through the catchment to varying uniform recharge rates on homogeneous numerical models in which topography-driven groundwater flows do not reproduce the streamflow measurements carried out in our catchment. Therefore, the conceptual models differ mostly in terms of the drainage surfaces and geological structures involved in the groundwater-surface Unravelling geological controls on groundwater flow and surface water-groundwater interaction in mountain systems: a multi-disciplinary approach 15/42 water interaction. The conceptual models are presented in detail in section 4.1 following the joint interpretation of data presented in section 3.4.

To offer a more quantitative point of view, the proposed conceptual models can be segregated by simple mass balance. Assuming steady-state conditions with no storage in the soil nor in the aquifer, a relatively low direct runoff coefficient considering the forest cover and relatively constant baseflow over the year, the effective precipitation in each catchment (defined as total rainfall minus actual evapotranspiration) can be directly compared to the baseflow per drainage area, i.e., what goes in goes out. A mean yearly value for effective precipitation is estimated from ground-based meteorological data for the period 1985-2013, as acquired by the park administration in the vicinity of the catchment outlet P0 (Figure 1). The estimated effective precipitation over the catchment in this period is close to 120 mm.

Results

We successively present the results of the geological mapping for hydrogeological considerations, the geophysical results, and the results related to streamflow measurements and remote sensing analysis.

Geology

Strictly based on lithological characteristics (Figure 1c), four rock units were defined and mapped in the studied area, some of which were recognized in more than one stratigraphic position. The first is an "Andesitic unit" composed mainly of andesitic lava flows and hypabysal porphyritic intrusives, with minor andesitic tuffs and sandstones. This unit was recognized west of the newly identified Papagayo fault presented in the next paragraph and at the base and in the middle of this fault hanging wall. The second is a "Sedimentary unit", which is a sandstone-dominated unit that also includes conglomerate, lithic tuff, and andesite intercalations and crops out discreetly in the footwall of the Papagayo faultextensively in the middle section of its hanging wall and as a small wedge within its upper section. The Unravelling geological controls on groundwater flow and surface water-groundwater interaction in mountain systems: a multi-disciplinary approach 16/42 third is a stratified "Dacitic unit", which is composed of brightly coloured dacitic lavas and tuffs recognized in the upper section of the Papagayo fault hanging wall. The fourth unit is a "Granitic unit" composed of monzogranites, syenogranites, and granodiorites, which outcrops in the western flank of the Principal Cordillera, intruding the stratified volcanic and sedimentary rocks.

The structural features of the study area allow the recognition of two distinct structural domains. The first is a western domain limited by the main trace of the west-vergent San Ramón fault and a newly identified N-S reverse fault, which here is called the Papagayo fault (Figure 1c). This domain shows tight short-wavelength (2-4 km) folds and minor inverse faults, generating highly dipping strata. The second domain is located to the east of the Papagayo fault and shows low dips on the flanks of longwavelength (6-8 km) folds. However, 2 km east of the Papagayo fault, a west-vergent anticline is recognized in the hilltops affecting the andesitic and sedimentary units. Directly beneath the anticline, faulting and fracturing is recognized affecting the granitic rocks outcropping along the riverbanks. A detailed structural survey in this area showed a zone with a high fracture density with vertical and closely spaced NS-and EW-striking sets and minor inclined NW-and NE-striking sets (criss-crossed zone in Figure 1c).

Geophysical study

Figure 2 presents the magnetic field reduced to the magnetic pole. The advantage of this representation is that it tends to show a direct correlation between the magnetic anomaly and the causative magnetic source (i.e., a positive correlation when the magnetic source is more magnetic than the surroundings, and a negative correlation for the opposite case). On this map, we identify four zones, mainly oriented north-south. Going from west to east, the first 2,000 m of the profile (zone 1) presents a strong positive magnetic anomaly that is relatively homogenous in the NS direction. Inside this zone, a short-wavelength (less than 200 m) but persistent negative anomaly is observed at longitude 361,000E. Secondly, in zone 2, between 362,000E and 364,000E, we observe a negative Unravelling geological controls on groundwater flow and surface water-groundwater interaction in mountain systems: a multi-disciplinary approach 17/42 anomaly with its larger amplitude in the northern half, and tapering off southward. This broad negative domain is followed to the east by another negative domain (zone 3) of the same length and signal amplitude (-60 to -80 nT) but decreasing southward and westward. The limit between zone 2 and zone 3 is a fairly linear positive NS anomaly around 364,000E. Finally, zone 4 extends east of 366,000E, and is characterized by a broad positive anomaly that seems to extend beyond the study area.

In agreement with the geological map, the mineralogical composition of geological units and susceptibility samples taken from outcrops (Table 1), the magnetic sources associated with zones 1 and 4 represent volcano-sedimentary units, with mid-to-high contents of ferromagnetic minerals. In contrast, the negative magnetic anomaly of zones 2 and 3 can be associated with the intrusive body mapped in the geology. This intrusive has a consistently more felsic mineralogy and low magnetic susceptibility with a clear negative contrast compared with the volcanic units in the flanking domains (Table 1). The sharp transition between magnetic zones 3 and 4 can be associated with the mapped contact between the intrusive unit and the volcanic unit at the eastern border of the intrusion. Meanwhile, at the western limit of the study area, we can identify a similarly sharp transition between zones 1 and 2, in good agreement with the contact between the intrusive and volcanic rocks. We postulate that the positive anomaly at 364,000 E (which separates zones 2 and 3) represents a fault that provided room for the emplacement of more mafic magmatism (the andesites of the flanking zones 1 and 4) and/or a path for hydrothermal circulation and secondary magnetic enhancement due to contact metamorphism. This anomaly is hence associated to the newly identified Papagayo fault which appears to be a large-scale regional fault according to the gravimetric survey carried out in the Río Coipo catchment. Looking at the path of the river entering zone 2 going west, we see a change in direction from E-W to S-N, which may be another sign of the influence of this fault on the river network.

Finally, in zone 1, the relative negative anomaly at 361,000E may be related to a demagnetization of the volcanic rocks caused by the hydrothermal fluid circulation at the piedmont fault zone.

Unravelling geological controls on groundwater flow and surface water-groundwater interaction in mountain systems: a multi-disciplinary approach 18/42

The Río Clarillo ERT profile (Figure 3) can be separated into two well-defined zones, identified as zones 1 and 2 on the aeromagnetic map. Zone 1 is defined at distances of between 0-1,300 m along the profile, where the environment is much more conductive than the other part of the profile (100-500 ohm.m versus more than 1,000 ohm.m eastward). In this zone there are some small-wavelength heterogeneities (100-200 m) that are relatively positive and tabular (along-profile distances of 200-400 m and 800-1,000 m) and a shallow sub-horizontal negative anomaly (along-profile distance of 100-500 m). Magnetic zone 2 is partially represented between 1,300 -3,600 m on the ERT profile. It contains a shallow (~ 50 m depth), almost continuous and relatively conductive layer (<300 ohm-m) above a resistant body (> 1,000 ohm-m).

The geological interpretation is based on the locations of the two aeromagnetic zones. The higher resolution of the ERT profile compared to the magnetic data suggest the presence of a shallow layer of sediments overlying the volcanic and intrusive units. Zone 1 evidences a more conductive layer which likely represents the alluvions that cover the volcanoclastic rocks of the western flank of the study area. The magnetic signal in this zone cannot be used to discriminate this sedimentary unit, probably due to its low magnetic content. The vertical conductive structure between 400-500 m and the horizontal thin layer between 0-400 m are likely to be associated with the location of the thermal spring, as this really low resistivity can be associated with saturated sediments. Zone 2, between distances of 1,300 and 3,600 m the eastern portion of the ERT profile, occurs within the granitic intrusion and is covered by a thin layer of alluvial material. The high and rather homogeneous resistivity underneath the thin conductive cover suggests that this intrusive is a relatively homogeneous unit that is probably cut by structures at distances of 2,000 m and 2,800 m along the ERT profile. Even though no regional large-scale structures are defined in this area, [START_REF] Piquer | The Piuquencillo fault system: a long-lived, Andeantransverse fault system and its relationship with magmatic and hydrothermal activity[END_REF] observed and mapped secondary structures inside the granitic intrusion at similar locations, which supports our working hypothesis on the discontinuity of the granitic body. These results indicate that Unravelling geological controls on groundwater flow and surface water-groundwater interaction in mountain systems: a multi-disciplinary approach 19/42 these heterogeneities represent faults or fold axes in the geology associated with the domain of high electrical resistivity.

The gravity anomaly along the Río Clarillo valley transect (Figure 4 upper panel) shows two distinctive sectors, a western sector characterized by a relative negative anomaly (-5 mGal) between distances of 0-500 m along the profile. To the east is a wide zone with a relatively steep gradient (2 mGal/km) between distances of 500-1,500 m along the profile. The remainder of the profile, between distances of 1,500-3,500 m, shows a relatively flat gravimetric response at -2.5 mGal. The proposed density model considers a background density of 2.67 gr/cc, as well as the structures and layers defined at the subsurface by ERT and magnetic interpretations. Densities for intrusive and volcanic rocks are based on in-situ rock sample measurements (Table 1), while sediment density is obtained from [START_REF] Yáñez | Gravity derived depth to basement in Santiago Basin, Chile: implications for its geological evolution, hydrogeology, low enthalpy geothermal, soil characterization and geo-hazards[END_REF].

As shown in Figure 4 upper panel, the model response is in good agreement with observations, and is thus a plausible solution that is consistent with geological observations. Low density and thick sedimentary cover above the volcanic basement explain the lower gravity of the western side. The thickness of the sedimentary unit decreases eastward in line with the decrease of the amplitude in the negative gravity signal. This sedimentary body represents the alluvial fan at the piedmont. Considering the interpretation of the ERT profile, the low-resistivity zone at a distance of around 500 m is here represented as a fractured area in the volcanic unit and a low-density area in the sediment layer (1.85 gr/cc). As seen on the ERT profile as a low-resistivity zone (below 50 ohm.m), this low-density zone extends between distances of 0 and 400 m, in the sediment layer, and probably represents a claydominated body associated with the thermal spring zone recognized in the area. Towards the east (zone 2), an almost flat signal is modelled as the gravimetric response of the more felsic intrusive body (2.63 gr/cc) covered by a relatively thin layer of sediments. In this zone, the sediment cover is approximately 50 m, which agrees well with the ERT transect, which shows a covering of intermediateresistivity material of a similar thickness. At the contact between the two zones, at a distance of around Unravelling geological controls on groundwater flow and surface water-groundwater interaction in mountain systems: a multi-disciplinary approach 20/42 1,300 m, there is a small decrease in the gravity anomaly. This negative anomaly is associated with a low-density fractured zone. This contact zone is associated with the fractures induced by the intrusive emplacement. The fracture zone identified at a distance of 500 m could favour groundwater exfiltration, creating a preferential upwelling flow, according to the location of the thermal spring over this fault. The proximity of the two fracture zones (which are located at along-profile distances of 500 m and 1,300 m) and their similar vertical structure indicate that a joint fault system may be connected at depth and act as listric faults.

Streamflow measurements and remote sensing analysis

Streamflow measurements acquired during the 2019 recession are displayed with regards to drainage (or upslope) area (Figure 5). A few interesting features are observed. First, we observe a significant increase of the streamflow at point P1 (grey triangles in Figure 5). However, this increase is not systematically observed within the same recession period (orange squares in Figure 5). Second, and most notably, it exhibits a significant increase of approximately 50% within a 1km distance (between P5 and P4-P3). We interpret this variation as groundwater exfiltration as there is no significant increase of the drainage area and there is no affluent observed in the area or a superficial supply of any type (snowmelt, water bodies). Groundwater exfiltration has been observed in the Carrizal area (Figure 1c) by park rangers, who have noted that "water is usually seen disappearing in the river" while they are going back up the Río Clarillo which is consistent with the lower streamflow rate from P4 to P5.

Additionally, some hydric features are observed downstream of point P5, such as a permanent puddle about 50 m uphill and a curtain of dripping water on the east at the riverbank.

The persistency of the groundwater exfiltration is confirmed by the difference in NDVI between the Carrizal area and the entire national park. The Carrizal area is displayed here as a strip 1,000 m in length and 500 m in width along the east side of the river where geological and hydric features have been observed. Over the period 2000-2019, the mean NDVI in the Carrizal area was 0.57, while it was 0.36

Unravelling geological controls on groundwater flow and surface water-groundwater interaction in mountain systems: a multi-disciplinary approach 21/42 in the entire national park. The difference in NDVI between the Carrizal area and the entire national park is almost systematically positive (Figure 6), indicating that the Carrizal area has denser vegetation due to a higher water availability. In conclusion, both qualitative observations and quantitative data indicate that the groundwater exfiltration between P5 and P4 has persisted for at least the last two decades.

Joint interpretation of geophysical and geological results and hydraulic markers

The interpretation of the magnetic map allows us to describe four north-south zones. The ERT and the gravity transects overlap zones 1 and 2 and indicate a similar spatial extent of these zones, which reinforces our results and interpretations. Moreover, the thermal spring is identified using each of the geophysical methods, indicating a relationship between this feature and the presence of a fault inside the volcanic unit. The temperature of the thermal spring is about 23°C. As the geothermal gradient is about 3°C/100 m in this region [START_REF] Muñoz | Heat flow and temperature gradients in Chile[END_REF][START_REF] Valdenegro | A model for thermal gradient and heat flow in central Chile: The role of thermal properties[END_REF], it is possible to approximate at first order that this thermal source is fed along a 300 m deep structure at least, which is supported by the ERT interpretation and the gravity modelling. Both methods indicate that the thermal spring is likely associated to upwelling flow along a fault or fractures zone as seen on other locations along the Andes in Central Chile [START_REF] Figueroa | Deciphering groundwater flow-paths in fault-controlled semiarid mountain front zones (Central Chile)[END_REF]. Our second hydraulic indicator is the increase of baseflow that was measured between points 5 and 4. This area was only monitored with a UAV-based magnetic study. The contact between the volcanic unit at the eastern end of the reserve and the intrusive unit might be of significance for the increase of the river baseflow. Nevertheless, the location of the increase in the river flow is not exactly located over the contact defined on the aeromagnetic map (Figure 2). According to the geological map, the resurgence area close to the contact, presents a pattern of local east-west fractures (criss-crossed zone in Figure 1c). As seen at the Unravelling geological controls on groundwater flow and surface water-groundwater interaction in mountain systems: a multi-disciplinary approach 22/42 western edge of the contact between the intrusive unit and the volcanic unit (Figure 4), a highly damaged zone engenders the structural control of water flows.

Discussion

Based on the various observations described in the previous section, we first discuss the capacity of various conceptual hydrogeological models to explain the observed streamflow variations and hydric singularities. We then discuss the differentiating aspects and novelty of the proposed multi-disciplinary approach for characterizing groundwater flow patterns in mountain catchments and for identifying structural controls on groundwater flow patterns and surface water-groundwater interaction.

Conceptual hydrogeological models

We propose three different conceptual models based on the observed geological features. First, groundwater exfiltration between P5 and P4-P3 could simply originate from the local drainage of lateral hillslopes in the Carrizal area, which are partly composed of alluvial deposits (Model 1 of Table 2, Figure 7). However, the Carrizal area offers a very small drainage area (3 km 2 from P4 to P5 and 9 km 2 from P3 to P5), which we argue is likely insufficient to explain the full extent of the resurgence.

Second, as the exfiltration occurs in the vicinity of the contact between the intrusive and volcanic units, it could be related to a reduction of the groundwater seepage area, that is, a reduction of the active depth of the system due to a much less permeable underlying unit outcropping (Model 2 of Table 2, Figure 7), consistent with the results from section 3.2. In that case, the drainage area could be the entire catchment area upstream of point P5, which is close to 100 km 2 . Finally, as shown on the geological map (Figure 1c) and described in section 3.1, the Carrizal area shows a large number of heterogeneities inside the intrusive unit. It is described as highly fractured in relation to a local E-W steeply dipping fault zone (criss-crossed zone in Figure 1c) and close to the N-S major Papagayo fault Unravelling geological controls on groundwater flow and surface water-groundwater interaction in mountain systems: a multi-disciplinary approach 23/42 This local fault zone and the associated water resurgence may be a local effect of the E-W faults and the fracture network in this area with no influence of regional structures. Still, it is likely that the principal N-S structure is directly related to these E-W faults forming a fault network with secondary structures perpendicular to the major one. Supporting evidence for this last model is provided by the very low streamflow from the nearest downstream tributary, Los Bueyes (green catchment in Figure 7), with no significant streamflow increase measured at P2. The catchment contributing to this tributary is directly stroke by the N-S fault identified by both the geological mapping and UAV-based magnetic results. This N-S structure likely exports groundwater of the contributing subcatchment to the identified resurgence area (Model 3 of Table 2, Figure 7).

As presented in Table 2 and above, each conceptual model corresponds to significantly different drainage areas. Model 2 can be discarded as the estimated flow from simple mass balance (as presented in section 2.5) is seven times larger than the resurgence flow. In model 1 and model 3, the estimated flow is 0.22 and 0.90 times the resurgence flow, respectively. Thus, model 3 is the most likely scenario considering the drainage area of each sub catchment. Hypothetically, a combination of models 1 and 3 might as well explain the resurgence processes considering the uncertainty in the simple hydrological balance presented in section 2.5.

Differentiating aspects and novelty of the approach

This work presents a methodology to identify structural controls on groundwater flow patterns and surface water-groundwater interaction. The scope of the methodology is above all functional in the sense that it aims to identify the key geological structures for groundwater flow but in no way to characterize the entire geological framework. This functional scope necessarily impacts the scale of the study which is neither at the fault scale (for example, there is no topological analysis of the faults zone), nor at the catchment scale as the elementary bulk unit. This functional scope besides impacts the methodological choices made, following the informational content of well-known data and Unravelling geological controls on groundwater flow and surface water-groundwater interaction in mountain systems: a multi-disciplinary approach 24/42 technics. As the current approach aims at linking geological structures to hydrogeological structures, it is inevitable to acquire hydraulic-related data. For that aspect, widely-used methods based on dissolved elements (environmental tracers, water isotopes, hydrogeochemistry) are directly related to transport processes and hence to groundwater flow and hydraulic gradient. However, their inability to clearly identify local scale processes mainly due to mixing, limits their use when addressing problems related to groundwater flow architecture as it is the case here. On the contrary, baseflow measurement is an easily accessible data that can be related to local flow processes as the resolution of spatially-distributed measurements can be arbitrary fixed. Beyond this consideration, streamflow is a variable of interest in various disciplines including hydroecology. In that sense, the comprehensive methodology proposed here could be easily extended to hydroecology or hydrobiology.

The novelty of the proposed methodology lies mainly in the innovative implementation of classical methods from various disciplines that are conjointly designed and iteratively carried out. The geological study has been indeed conducted at multiple scales: at a traditional main scale (1:100,000) and at a finer scale to incorporate some elements (e.g., the N-E secondary fracture zone) that would confirm hydraulic observations. It results in a hybrid geological map partly fed by hydrogeological considerations. Regarding geophysical surveys, they significantly compliment the multiscale geological mapping. The inclusion of a second geophysical transect in the north of the study area (transect B-B' n Figure 1c) confirms the general N-S orientation and regional scale of some structures (results presented in supplementary material). Meanwhile, the UAV-based magnetic study, performed here at high-resolution in comparison to previous studies, does not only provide valuable information in remote areas but also constitutes a unifying method allowing the creation of a coherent geological model thanks to its fast and large covering.

Finally, while spatially-distributed streamflow campaigns are not commonly used in hydrology, they have shown here to be valuable to identify and confirm critical features of the baseflow evolution within the catchment. If spatially-distributed flow measurements may seem tedious, low-cost Unravelling geological controls on groundwater flow and surface water-groundwater interaction in mountain systems: a multi-disciplinary approach 25/42 instrumentation and recently-developed optical methods could allow a systematization of measurements.

Conclusions

We propose a comprehensive methodology to unravel geological controls on groundwater flow and surface water-groundwater interaction in mountain catchments. Geological mapping, on-ground and aerial geophysics, and hydraulic data are acquired to simultaneously image geological structures and their hydrogeological manifestations.

Using the novel methodology, we identify peculiar hydric behaviours in an ungauged Andean catchment, including a significant and localized groundwater exfiltration that seems to be imported from a neighbouring subcatchment and related to the presence of a regional N-S fault and a secondary set of E-W local faults.

The novelty of the methodology mainly lies in the innovative implementation of classical methods: multi-scale geological mapping, spatially distributed streamflow measurements, high-resolution aerial geophysics. This novel methodology is distinguished by its pragmatism in that it focuses on the geological structures that are the most significant for flow without aiming to characterize the entire geological framework.

The findings have significant implications for the management and monitoring of the study site as flow singularities such as the abrupt streamflow change that may strongly affect the biogeochemical conditions in the river, like biogeochemical hotspots.

Unravelling geological controls on groundwater flow and surface water-groundwater interaction in mountain systems: a multi-disciplinary approach 26/42

Acknowledgments

This study received financial support from the project ANID-Conicyt/Fondecyt de Iniciación n°11170380. The authors would like to thank the Pirque community for their incredible kindness, their constant support and their interest in the HydroPirque project. Miles y miles de gracias a la CONAF: Don Carlos, Don Luis P., Don Luis C., Doña Catalina y todos los guardaparques de la reserva. Miles de gracias también al comité de Agua Potable Rural El Principal: a la señora Nora y a Erik, a la municipalidad de Pirque y a Nicolas. Thanks to Mr. Carlos Uribe del Río Coipo: sometimes research hangs by a thread, or rather by a padlock key. The authors are thankful to the two anonymous reviewers and the associated editor for their fruitful comments that significantly improved the manuscript. understanding, it includes projections of the hydraulic markers on the profile (red dot for the thermal spring and white dots for the river flow measurements) as well as the magnetic zones limits (Figure 2) presented above the transect.

Figure 1

 1 Figure 1: a) Satellite view of the Santiago area, the white rectangle outlines the study area. b)

Figure 4 :

 4 Figure 4: Representation of the gravity study for the Río Clarillo transect (A-A', Figure 1) and the

Model

Main process Surface and river flow measurements (black-circled white dots). Numbers 1 to 3 correspond to the possible drainage areas contributing to the resurgence between points P4 and P5, based on the three models mentioned in sections 2.5 and 4.1. Grey area within the catchment is not part of any of the proposed models.


Carole Bessette
email: carole.bessette@uha.fr

Mathieu Decrette

Michel Tourlonias

Jean-François Osselin

Marie-Ange Bueno

Adrienne Leca

Dominique Coupé

Warp yarn damages prediction regarding to the evolution of yarn tension during 3D weaving process

Keywords: yarn damages, yarn friction, tension, weaving process

Interlock weaving process is still difficult to manage and lead to significant damages on warp yarns, causing a serious loss of mechanical properties on the final composite, and affecting the good execution of the process. The aim of this study is to have a contribution to 3D weaving process optimization, regarding to yarn-toyarn friction. To predict damages due to this friction, evolution of warp yarn tension can be a good indicator. Firstly, evolution of the yarn tension during a weaving cycle has been precisely analysed. Then, great differences of tension have been highlighted, function of types of crossing. Specific yarn crossings generate very high peaks of tension on the studied yarn, leading to severe damages until yarn breakage. A significant difference has been noticed between the polyester and glass investigated yarns. The study shows the evolution of tension can be helpful for process monitoring and optimization. Nevertheless, in order to anticipate yarn damages, a direct measurement of yarn interactions has to be developed.

Introduction

Woven reinforcements for high performance composite applications are already widespread in several domains mostly as replacements of former pieces made of steel or aluminium alloys. Woven composites are particularly used in aeronautic industries, like military and commercial aviation, thanks to their great mechanical properties and their ability to make pieces of complex shapes [START_REF] Marsh | Composites consolidation in commercial aviation[END_REF][START_REF] Coupé | 3D textile structures for aeronautic applications[END_REF]. In order to stay competitive with other materials, composites have to be robustly manufactured. Textiles processes such as weaving, give a large diversity of 2D or 3D dry preforms which are perfectly compatible with composite production. 2D fabrics are the first kind of fabrics used for composite reinforcements, different layers of fabrics are assembled in a stratified structure to form a composite piece. Stratified composites are quite easy to manufacture but have a major drawback due to their structure, which is delamination between layers. Delamination is a serious threat leading to composite breakage, therefore solutions were found to reduce this weakness. 3D woven composite is a good option to avoid this issue, which explains their fast growth. Many studies have presented a comparison of mechanical properties between 2D and 3D. In almost all cases 3D woven reinforcement have showed better properties. Mouritz et al. presented that 3D woven composites have better mechanical properties than 2D reinforcement composites such as ballistic damage resistance, impact tolerance resistance, tensile strain, and interlaminar fracture toughness [START_REF] Coupé | 3D textile structures for aeronautic applications[END_REF]. Ding et al found also that 3D composites have greater performance in fatigue behaviours than 2D structures [START_REF] Mouritz | Review of applications for advanced threedimensional fibre textile composites[END_REF].

It exists many ways to weave 3D fabrics, one of the most common is 3D interlock process, which have the advantage to be employed on a conventional weaving loom. In this study, we are focusing on 3D warp interlock weaving, in which some warp yarns link several or all layers together [START_REF] Ding | Comparison of the fatigue behavior of 2-D and 3-D woven fabric reinforced composites[END_REF]. Interlock weaving process is not fully managed yet and leads to significant damages on warp yarns, causing a serious loss of mechanical properties on the final composite and impacting on the good execution of the process [START_REF] Wendling | Consistent geometrical modelling of interlock fabrics[END_REF][START_REF] Decrette | Jacquard UNIVAL 100 parameters study for high-density weaving optimization[END_REF]. Many studies have been performed about the effects of friction between yarns and weaving loom elements on mechanical properties at different scales. L. Lee et al. worked on the effect of weaving process on the tensile properties of glass yarns woven reinforcements composites. No significant effect have been noticed on tensile stiffness of yarns, but a 30 % decrease of the tensile strength value of dry woven yarns has been observed [START_REF] Archer | The effect of 3D weaving and consolidation on carbon fiber tows, fabrics and composites[END_REF]. In the same document, they observed warp yarns breakage due to friction against loom elements. Rudov-Clark et al. preconized to cover loom machinery with wear-resistant coating to reduce abrasion damage [START_REF] Lee | Effect of weaving damage on the tensile properties of three-dimensional woven composites[END_REF]. Lefebvre et al. investigated the loss of mechanical properties on para-aramid warp and weft yarns [START_REF] Rudov-Clark | Fibre damage in the manufacture of advanced three-dimensional woven composites[END_REF]. The study showed that warp yarns are more impacted than weft yarns, and interlock pattern also has an influence on the loss of properties.

It was shown that Warp yarn damages during weaving process impact the reinforcement quality, but only few studies had investigated the causes, and even less presented results on yarn-to-yarn friction. Decrette et al., presented an optimization strategy based on warp yarns degradations count during a weaving process, and proposed to reduce yarns damages by using Jacquard UNIVAL 100 parameters [START_REF] Decrette | Jacquard UNIVAL 100 parameters study for high-density weaving optimization[END_REF]. The aim of this study, is to make a contribution to 3D weaving process optimization and monitoring, regarding evolution of warp yarns tension during weaving process. Yarn tension is linked to several parameters such as yarn-to-yarn interactions, yarn wear or fibres breakage, moreover, this parameter is quite simple to measure during the process. A specific pattern has been designed and woven to induce different kinds of yarn-to-yarn interactions, and yarn tension has been analysed during this process to confirm if it is relevant as a good indicator for damage prediction.

Materials and methods

3D weaving process

Experiments have been performed on a sequential weaving loom fitted out with a Jacquard machine shedding device. The shedding system is a UNIVAL 100 working with several individual motors where each one drives a Jacquard string. Therefore, each motor is linked to one yarn and controls its movement and position. This mechanism offering lots of possibility in pattern design, appears to be really appropriate for 3D weaving. To let each yarn evolves independently, a 512 bobbins creel is used in our case. The weaving loom has other characteristics like:

• Weft insertion with unilateral rapier

• Rectilinear reed movement,

• Horizontal take-up,

• 16 yarns per reed dent,

• 1.5 mm of space between reed dent The aim of this part is not to weave a specific interlock fabric but to investigate different types of crossing interactions between yarns during that kind of process. A test pattern has been done to scan a sample group of nine crossings and determine their effects on the yarn tension value and evolution in time. This weave pattern is presented in figure 2, it has been designed regarding 16 yarns of the same reed dent corresponding to a single harness column, to be sure that yarn to yarn interactions are generated during a measurement.

The red frame, in figure 2, indicates the evolution of the studied yarn, which is connected to a tensiometer. In order to contribute to process optimization and monitoring, crossings could be classified function of their impact on the measured yarn tension. Crossings can be depicted by other graphical representations, as it is illustrated on figure 3. In this representation, red circles and red arrows represent respectively the studied yarn and its movement, bold grey arrows show the movement of the whole other yarns of the column and small grey arrows the movement of a single yarn which is not the studied one. Experiments has been performed on two different yarns:

• 170 tex,PET yarns

• 600 tex E-glass yarns, PET and glass yarns have very different mechanical behaviours, like brittleness which is interesting to compare on tension evolution.

High speed camera

The sequential organisation of the weaving loom implies that only one machine element is working at a time, and waiting times isolate each sequence. Waiting times were unfortunately unknown and not editable, but are essential to understand the global machine time line. High speed camera has been used to define precisely all machine movement and the sequences duration. This study is imperative to know which stage happens at what time. This information is used for tension evolution analysis. Weft insertion stages do not affect warp yarns interactions, no friction happens at that time, this is why those stages are not taken into account in table 1. Figure 4 presents a screen shot of a camera recording.

Yarn tension measurement

Warp yarn tension is an easily measurable parameter on a weaving loom, and possibly useful for process management. In this study, a 0-10N range tensiometer has been placed on the machine, between harness back and in front of the condensation slab, which makes a geometrical link between creel dimension and harness. Figure 5 is an illustration of the set-up; the arrow shows the tensiometer. Preliminary studies have been hold to make sure that tensiometer position does not interfere in the weaving process, and does not create new friction surfaces between yarns. This way, tension is measured continuously during the process, and can be correlated with high speed camera images. Thereby, every tension variation or peak is associated to a machine stage. Tension value is determined for every crossing and its evolution in time is measured. Tension function of time curves give some information:

• Which kind of stage impacts the most yarn tension?

• Which kind of pattern crossing affects the most yarn tension? Afterwards, a spectral analysis has been made on tension curves during hundreds of cycles, to identify a correlation between yarn damages and yarn tension. All the results are showed and analysed in results and discussion part below. First of all, yarn tension function of time has been analysed during several weaving cycles. Preliminary measures have been hold to make sure that yarn-to-yarn friction impacts on yarn tension and can be detected by the tensiometer. Crossing pattern has been woven, on 16 PET yarns and then on 16 glass yarns.

Results and discussion

Yarn tension function of time

PET yarns

Figure 6a illustrates a significant example of yarn tension as a function of time for PET yarns during a complete weaving pattern cycle. Two areas can be distinguished: i) a huge range of tension variation, one of them being highlighted by a blue frame on the figure and ii) a low amplitude of variation, like a steady state, represented by a black frame. Thanks to the high speed camera, machine stages have been linked to the tension evolution. The high range variation is corresponding to beat-up, shedding, and back reed movement steps. At the opposite, low variation occurs during the weft insertion stages. Those results are consistent with the hypothesis that no friction occurs during weft insertion, and is in accordance with Boussu et al. study on weaving process monitoring [START_REF] Lefebvre | Influence of high-performance yarns degradation inside three-dimensional warp interlock fabric[END_REF]. Moreover, focusing on the high variation area, difference between crossings is highlighted. Some tension peaks due to specific crossings are higher than others. Series of five identically experiments have been hold to confirm these tendencies. Results are similar to the one presented in figure 6a. Figure 6b is an enlargement of blue frame, i.e. tension variation where peaks are more perceptible.

The high speed camera results are overlaid on curve of tension evolution of a one crossing in figure 7. This crossing corresponds to maximum value of tension obtained during experiments. This curve is the Fig. 7: Yarn tension evolution with machine sequence Fig. 6: Warp yarn tension evolution for PET yarns average of 5 different tests for the same crossing. Several peaks are observed, due to different events.

At 100 ms after the beginning of reed movement, a peak appears, corresponding of reed beat-up. A second huge peak, around 270 ms is a consequence of shedding, and a third peak at 430 ms is linked to reed back movement. Small peaks are also observable, at 130 ms and 520 ms, they must result from yarn vibration on the tensiometer sensor. Thanks to the analysis of a high camera recordings, it has been found that the third peak is also due to yarn vibration. The back movement of the reed seems to have an insignificant effect.

This kind of curve clearly indicates that interactions occurring during shedding strongly impact yarn tension. It can be noticed that it is during that specific sequence. That yarn tension changes the most The maximum peak value is reached during shedding for any type of crossings. This maximum value is the chosen parameter to classify crossings with each other. Several crossings emerge as "critical" ones, given their high tension generation. They are represented on figure 8 (crossings 6, 7 and 8). The same colour code as employed in figure 3 is used and the graphical representation has been made considering a reed dent front view.

For every critical crossing mentioned, the majority of yarns are in movement, the tension is severely impacted even if the studied yarn is motionless.

Glass Yarns

Glass yarns are used for technical applications, and their mechanical behaviours is quite different compared to PET. For example, glass yarns are more sensitive to fibrillation, because of material brittleness. Same tests as PET yarns have been conducted. Results for the complete weave pattern (as figure 6a), are not presented here, as they are very close to PET results. The maximum values of tension appear at different crossings, but the same shape of curve is obtained. Tension peaks occur at the same time and weft insertion has similarly no impact on tension. Figure 9 shows the evolution of the mean glass yarn tension (5 repetitions) during the same crossing presented in figure 7 for PET (No. 6 of the pattern), overlaid with machine sequences. The first peak appears at the same time, during beat-up stage, and is higher than PET results. The second peak during shedding is comparable to the one obtained with PET yarns. Contrary to PET results, any third peak is noticed and yarn vibrations are lower maybe due to higher yarn count.

Other differences are noticeable between PET and glass yarns results; critical interactions occur for different crossings. For glass yarns, three critical crossings are found, which are No. 2, 4 and 6. Those critical crossings are depicted in figure 10 with a reed dent front view. In the glass yarns case, critical interactions occur only when the studied yarn is in movement with an ascending motion. Descending motion have obviously less effect on tension. However, crossing No. 6 is critical for PET and glass, interactions occurring during this shedding strongly impact the warp yarn tension, due to an important friction mechanism.

Spectral analysis

Tension can be used to distinguish some crossings from the others, and seems to be sensitive to yarn-toyarn friction. To determine if tension can be used as a damage state indicator for a warp yarn, a spectral analysis has been done for hundreds insertions in a row. Time-frequency graphs are obtained and the corresponding spectral power has been calculated. From a weft insertion to the following, a FTT spectrum is obtained. So for every crossing, a FFT spectrum is related. Each spectrum is successively plotted in a timefrequency graph, showing if some interactions seems to be more critical than others. Moreover, for each FTT spectrum, corresponding to a specific crossing, spectral power is calculated using the following relationship: f: frequency (Hz) Δf: spectre frequency range used (Hz), RMS(f): Root Mean Square of yarn tension at frequency f (N), C: coefficient corresponding at Hanning window function used in this case, equal to 1.5.

Spectral power function of pick number is plotted for three hundred cycles in a raw for PET and plotted until yarn breakage for glass yarns.

PET Yarns

Only one crossing is really notable in the graph presented on figure 11, corresponding to insertion No. 6. Peaks are present at 12, 24 and 32 Hz. This result is matching with tension evolution analysis, that kind on graphical representation can be used as an indicator. The advantage of this time-frequency graph is to visualize more cycles in a raw, which is more convenient for process monitoring. Spectral power evolution is shown on figure 12. Measurement is made on 300 cycles in a raw, making a global evolution noticeable. From 1-50 cycles, spectral power increases linearly, then from 51-end, power is stabilized around 0.7 N. Moreover, power value is dependant of the type of crossing. Once more, crossing No. 6 gives the higher value of spectral power. No yarn breakage happened during those tests, but a significant wear appears on yarns surface.

Glass Yarns

As it has been done for PET yarns, spectral analysis has been hold with glass yarns. Figure 11 presents a time frequency graph and Figure 13, the spectral power evolution is illustrated. Like for PET yarns, crossing No. 6 is highlighted by the representation: it is the only crossing with high peaks visible at 24 Hz and 36 Hz frequencies. Spectral power evolution is really relevant in this case, because studied yarn has broken during the test. This breakage occurred around 180 th cycles, and corresponds at power drop on the curve. Before tension fall, three areas can be distinguished. From 1 to 40 weft insertions, power increases linearly, in the same way as PET result. From 41 to 140 weaving cycles, spectral power is stabilized around the same value as PET spectral power threshold, which is 0.7N. From the 140 th insertion to the yarn breakage a sharp increase occurs leading to breakage. This curve is really interesting because a clear correlation between PET and glass results can be made. PET studied yarn did not break during experiment, but has been progressively damaged, whereas the monitored glass yarn has been hardly damaged and then broken. These different damage states seems to match with spectral power evolution. A spectral power evolution theoretical profile composed of three areas, can be proposed in figure 15. The first is a linearly increase corresponding at an initialisation leading to a stabilization of the value. The second part is a steady state of spectral power, when yarns are not too hardly damaged. The third part is consisting in a huge tension increase, corresponding to a critical damaging zone for yarns, leading quickly to breakage. The duration of each area depends on the material considered and certainly other parameters which need to be investigated like yarn sizing.

Conclusion and Perspectives

Those tension tests respond to several issues, such as identifying which kinds of crossing and which events of the weaving cycle lead to very strong yarn-to-yarn interactions. The results obtained in this study prove that tension evolution measurements can capture yarn damages, which can be helpful for process monitoring and optimization. A correlation between yarn damages and spectral power of tension signal has been made, which represent a high potential of use for process monitoring. Nevertheless, yarn-to-yarn interactions need to be more precisely characterized in order to anticipate, and at long term, to avoid yarn damages.

Fig. 1 :

 1 Fig. 1: Weaving machine loom.

Fig. 3 :

 3 Fig. 3: Few crossings from test pattern on lateral harness view.

Fig. 5 :

 5 Fig. 5: Tensiometer position.

Fig. 4 :

 4 Fig. 4: Screenshot from camera recording.

Fig. 9 :Fig. 8 :

 98 Fig. 9: Glass yarn tension evolution with machine sequence

𝐶Fig. 10 :

 10 Fig. 10: Critical crossings for Glass yarns

Fig. 11 :

 11 Fig. 11: Time frequency graph during 100 crossing cycles for PET yarns

Fig. 14 :

 14 Fig. 14: Time frequency graph during 100 crossing cycles for glass yarns

Table 1 :

 1 Weaving machine time line

	Stages	Mean time (ms)	Standard deviation (ms)	Variation coefficient
	Beat up	107.5	1,7	1,5%
	Dead time	57.37	6,1	10,5%
	Shedding	153.71	8,7	5,6%
	Dead time	135.02	6,8	5,0%
	Back reed movement	234.80	8,3	3,5%
	Whole 3 stages	688.40
	Remaining part of cycle	1252.	7,0	0,6%