Hydrogeological conceptual model of andesitic watersheds revealed by high-resolution heliborne geophysics

We conducted a multidisciplinary study at the watershed scale of an andesitic-type volcanic island in order to better characterize the hydrogeological functioning of aquifers and to better evaluate groundwater resources. A heliborne Time Domain ElectroMagnetic (TDEM) survey was conducted over Martinique in order to investigate underground volcanic structures and lithology, characterized by high lateral and vertical geological variability and resulting in a very high heterogeneity of their hydrogeological characteristics. Correlations were made on three adjacent watersheds between resistivity data along flight lines and geological and hydrogeological data from 51 boreholes and 24 springs, showing that the younger the formations, the higher their resistivity. Correlation between resistivity, geology and transmissivity data of three aquifers is attested: within the interval 10–100 ohm m and within a range of 1 to 5.5 Ma, the older the formation, the lower its resistivity, and the older the formation, the higher its transmissivity. Moreover, we demonstrate that the main geological structures lead to preferential flow circulations and that hydrogeological watersheds can differ from topographical watersheds. The consequence is that, even if the topographical watershed is small, underground flows from an adjacent watershed can add significant amounts of water to such a catchment. This effect is amplified when lava domes and their roots are situated upstream, as they present very high hydraulic conductivity leading to deep preferential groundwater flow circulations. We also reveal, unlike basaltic-type volcanic islands, that hydraulic conductivity increases with age in this andesitic-type volcanic island. This trend is interpreted as the consequence of tectonic fracturing associated with earthquakes in this subduction zone, related to andesitic volcanic islands. Finally, our approach allows characterization in detail of the hydrogeological functioning and identification of the properties of the main aquifer and aquitard units, leading to the proposition of a hydrogeological conceptual model at the watershed scale. This working scale seems particularly suitable due to the complexity of edifices, with heterogeneous geological formations presenting high lateral and vertical variability. Moreover, our study offers new guidelines for accurate correlations between resistivity, geology and hydraulic conductivity for volcanic islands. Finally, our results will also help stakeholders toward a better management of water resources.

1 Introduction

Water resource management on volcanic islands is challenging as these territories are often densely populated, subject to several natural hazards (volcanism, earthquakes, tsunamis, landslides, erosion and sea-level rise, etc.), and with increasing water demands (for irrigation, drinking water, etc.) or overexploitation of rivers or aquifers. Understanding the hydrogeological functioning of these islands is thus a major issue to achieve a sustainable management of their water resources. Hydrogeology of volcanic islands is challenging, taking into consideration the complexity of these edifices and the difficulties encountered when acquiring accurate in situ data (such as steep slopes, tropical vegetation, few access tracks, distance from laboratories, extreme climatic and hydrometric conditions for equipment). Indeed, as shown by Ingebritsen et al. (2006), volcanic formations exhibit extreme spatial variability or heterogeneity, both among geologic units and within particular units, with large variation from core scale to regional scale, permeability being, especially in volcanic environments, a scale-dependent property.

Historically, basaltic islands have been widely studied (e.g. Hawaii: Peterson, 1972; Macdonald et al., 1983; Canary Islands: Ecker, 1976; Custodio et al., 1988; Custodio, 2005; Custodio and Cabrera, 2008; Cruz-Fuentes et al., 2014; Izquierdo, 2014; Iceland: Sigurðsson and Einarsson, 1988; Réunion: Violette et al., 1997; Join et al., 2005; Azores: Cruz and Silva, 2001; Cruz, 2003; Galapagos Islands: d'Ozouville et al., 2008; Pryet et al., 2012; Violette et al., 2014; Jeju: Hamm et al., 2005; Won et al., 2005, 2006; Hagedorn et al., 2011; or Mayotte: Vittecoq et al., 2014), leading to several hydrogeological conceptual models, essentially at the island scale, each model being intrinsically dependent on the dynamic of volcanism activity, on the number and history of volcanoes and on their effusive and rest phases, generating a more or less complex geometry within which water infiltrates and circulates in a complex pattern, according to the recharge conditions.

Andesitic islands in subduction zones, and especially the Caribbean ones, are less known, and a limited number of hydrogeological studies have been conducted and published in these archipelagos, mainly at the island scale (e.g. Unesco, 1986; Falkland and Custodio, 1991; Davies and Peart, 2003; Gourcy et al., 2009; Vittecoq et al., 2010; Robins, 2013; Hemmings et al., 2015). Charlier et al. (2011) showed the interest in working at the watershed scale to define a hydrogeological scheme of a tiny site (45 ha) in Guadeloupe. Hydrogeological analyses of volcanic formations at several scales are indeed essential, especially for andesitic volcanism, characterized by heterogeneous geological formations, with alternation between intense eruptive phases marked by andesitic lava flows, pyroclastic flows, lahars, etc., interspersed with quieter phases marked by the dismantling of the volcano with debris avalanches and meteoric and alluvial erosion (Westercamp et al., 1989, 1990). Furthermore, andesitic stratovolcanoes display volcanic facies trends with variation and lateral distribution between central, proximal, medial and distal zones, depending on the valley and interfluve dynamics (Vessell and Davies, 1981; Bogie and Mackenzie, 1998; Selles et al., 2015). Finally, meteorological and hydrothermal weathering processes are superimposed on these lithological heterogeneities. This high lateral and vertical geological variability thus induces a very high heterogeneity of their hydrogeological characteristics. As shown by most of these studies, without in-depth data, it is not possible to understand relevant geological structures and consequently to understand the hydrogeological functioning.

Recently, heliborne geophysical surveys (e.g. Sorensen and Auken, 2004) started providing new regional in-depth data, which contribute to solving this scientific and technical challenge. High-resolution heliborne EM (ElectroMagnetic) resistivity data provide information down to the first 100 m along flight lines, and allow a continuous imagery of resistivity variations. Geological structures and hydrogeological properties can then be interpreted from these geophysical data to determine and constrain accurate conceptual models. To be relevant, and because resistivity is not a univocal parameter, this dataset analysis must be constrained with as much direct observation data (outcrop, borehole geological log, hydraulic conductivity data, etc.) as possible (see for instance Vittecoq et al., 2014) for Mayotte basaltic island.

Vittecoq et al. (2015), in studying an andesitic coastal aquifer in Martinique, demonstrate the relevance of working with and analysing heliborne EM data at the aquifer scale to characterize geological and hydrogeological heterogeneities of a 15 Ma old geological formation. At this scale, this approach is corroborated thanks to a very long-term pumping experiment. The working scale should indeed be sufficiently fine to be relevant to the structural specificities of these andesitic volcanic islands. However, working scale should also include surface and hydrogeological watersheds to integrate water balance estimation, interaction between groundwater and surface water, and potential contribution of different aquifers and vertical downward transfers, for a comprehensive view of the water cycle, so that stakeholders can use the results for sustainable management of water and energy resources.

Considering these issues, we conducted a multi-disciplinary approach at a watershed scale, based on the correlation of geological, hydrological, hydrogeological and heliborne Time Domain ElectroMagnetic (TDEM) data. We focus on a few strategic watersheds situated in Martinique, a predominantly andesitic volcanic island (Westercamp et al., 1989) located in the Lesser Antilles volcanic arc, in the subduction zone between the Atlantic plate and the Caribbean plate. The goals of our study are thus to (i) characterize the structure and hydrogeological functioning of Martinique andesitic aquifers at the watershed scale, (ii) show the influence of geological structures on groundwater flows and the consequence for the interactions between rivers and aquifers, (iii) assess the adequacy and difference between hydrological watersheds and hydrogeological watersheds, (iv) propose a conceptual model at the watershed scale, and (v) strengthen the hypothesis of Vittecoq et al. (2015) that, in contrast to the basaltic islands, hydraulic conductivity may increase with age in andesitic-type volcanic islands.

2 Martinique and studied watersheds

2.1 Site location and climate

Martinique (Fig. 1) is the largest volcanic island (1080 km²) of the Lesser Antilles archipelago. Its relief is mountainous in the north (highest volcano at 1397 m) and gentler in the south (highest hill at 504 m). Rainfall is characteristic of a humid tropical climate controlled by trade winds and orographic effects (Guiscafre, 1976; Vittecoq et al., 2010), with the rainy season between July and November and the dry season between January and April, interspersed with fluctuating transition periods. Annual temperatures vary between 18 and 32 ^∘C at Fort-de-France and an easterly trade wind regime ensures relatively constant ventilation. Average annual precipitation (Fig. 1c) is high in the northern part, reaching 5000 to 6500 mm yr⁻¹ at the summits and between 1200 and 1500 mm yr⁻¹ in the south.

Figure 1 The location of the island of Martinique (a) on the scale of the Americas and (b) on the scale of the Lesser Antilles. Locations of the watersheds (c) on the scale of the northern part of Martinique with annual rainfall map. (d) Location of rivers, water supply dams, gauging stations, watersheds, thermal springs, cold springs, boreholes and HTEM soundings along flight lines.

The three studied watersheds (Fig. 1) are located just near the capital city of Fort-de-France whose urban area includes half the population of the island (376 500 inhabitants on the island in 2016). Three dams are located on the Case Navire River, and provide an average of 5.9×10⁶ m³ yr⁻¹ to the urban area. During the driest seasons, the river is often dry over several hundred metres downstream of the dams, causing strong environmental impacts. Consequently, scientific studies are expected to understand the hydrological and hydrogeological functioning of this area, in order to propose alternative water resource management.

2.2 Geology

The volcanic activity of Martinique (Westercamp et al., 1989; Germa et al., 2010, 2011), which began more than 25 Ma ago, is characterized by a succession of many volcanic formations, mainly andesitic, set up from a dozen principal volcanic edifices, active during successive phases, with alternating periods of construction and erosion, sometimes contemporary.

The geology of the study area (Figs. 2a and S1 in the Supplement) is concerned with two distinct phases and volcanic edifices (Westercamp et al., 1989): the Morne Jacob shield volcano and the Carbets volcanic complex (Fig. 2b). The Morne Jacob shield volcano is the largest edifice on the island and lasted 3.3 Ma. Given its position, offset from the pre-existing reliefs, the first phase is first submarine and then aerial. First-phase formations are mostly weathered, because of a long period of rest and erosion of at least 1 Ma before the next phase. The strong aerial effusive volcanic activity of the second phase of the Morne Jacob volcano is shown on the field by massive flows (²α) up to 200 m thick. The Carbets volcanic complex developed on the western flank of the Morne Jacob shield volcano and lasted 1.8 Ma with four main aerial phases.

Figure 2 (a) Geological map (adapted from Westercamp et al., 1990) of the studied watersheds. Location of the piezometers (piezometric chronicles in Fig. 4): (1) Case Navire, National number 1177ZZ0165, (2) Fond Lahaye, National number 1177ZZ0161, (3) Fond Lahaye, National number 1177ZZ0177 and (4) Case Pilote, National number 1177ZZ0173. Cross-sectional location in white (cross sections in Figs. 6 and 7). (b) Litho-stratigraphic scale (adapted from Westercamp et al., 1990, and Germa et al., 2011).

Despite this detailed knowledge of the nature and location of the geological formations constituting the watersheds, and their lateral extension at the 1:50 000 scale, it remains difficult to have a precise and 3-D vision of their geometries and relationship at depth.

2.3 Hydrogeology

The position of the springs and the available drilling data (Figs. 1 and 2, Tables S1 and S2) suggest that aquifers could be associated with almost every volcanic phase of each edifice.

2.3.1 Cold springs

Springs (Figs. 1 and 2 and Table S1) are located mainly in the upper parts of the watersheds. Spring water discharges are most of the time a few litres per second. They are associated with four main geological formations (Fig. 3a and b). Seven springs, situated between 440 and 580 m a.m.s.l., emerge from andesitic domes and lava flows ⁹αbi (0.3 to 0.9 Ma). This geological formation is the last main event of the Pitons du Carbet Complex, strongly marking the landscape with several monolithic domes. In addition to observed springs, many perennial rivers flow from these peaks, so the aquifer that feeds these springs and rivers can be considered an important perched aquifer. Three springs, situated between 473 and 505 m a.m.s.l., emerge from andesitic and dacitic lavas ⁸ρα (0.9–1.2 Ma). Nine springs, emerging between 135 and 631 m a.m.s.l., are associated with andesitic lavas ²α (2.2–2.8 Ma), and four springs, emerging between 296 and 350 m a.m.s.l., are associated with basaltic lavas ¹βol (4–5.5 Ma). These springs are mostly situated at the slope foot, at slope breaks or at the top of gullies. Andesitic lavas ²α and ⁸ρα and basaltic lavas ¹βol are thus permeable and considered aquifer formations. Finally, one spring emerges from debris flow (⁶B) associated with the first phase of construction of the old Carbet (2 Ma).

Figure 3 (a) Comparison between borehole and spring elevation and associated piezometric level (for 26 boreholes). The piezometric level is on average 7 m below ground level, following this linear relationship: zw $=\text{elev}-\mathrm{6.94}$ (R²=0.99), where “zw” is the piezometric level (m) and “elev” the elevation (m). (b) Topographic profiles of the Case Navire, Dumauzé, Duclos and Fond Lahaye rivers, the piezometric levels of boreholes in the associated watersheds with the lithology of the aquifer, and the elevation and lithology of the aquifer of springs.

2.3.2 Thermal springs

Two thermal springs, Didier (210 m a.m.s.l. – 32 ^∘C – 1850 µS cm⁻¹) and Absalon (350 m a.m.s.l. – 36 ^∘C – 1730 µS cm⁻¹) are situated in the middle of the Case Navire watershed (Figs. 1 and 2), at a distance of 1700 m from each other. Their waters are mainly bicarbonated Ca–Na–Mg and are associated with high emissions of magmatic CO₂ and precipitation of iron hydroxides (Gadalia et al., 2014). The geochemical model (Gadalia et al., 2014) proposes an evolution in three stages: (1) deep mixing between water of meteoric origin and marine water (around 0.1 %), during a first partial chemical and isotope equilibrium; (2) water–rock and magmatic CO₂ interaction at medium temperature (90–140 ^∘C) in a residual geothermal system; and (3) mixing with fresh waters during the ascent, at a lower temperature. The Absalon spring emerges within fissured and fractured andesitic lavas ²α. The geological context of the Didier spring is poorly known because the bottling plant masks the outcrops. A borehole drilled 200 m from the spring shows, under a thickness of 16 m of pyroclastic flow, andesitic lavas ²α over 80 m thick. Waters of those two springs are thus mixed with the waters of the aquifer of andesitic lavas ²α.

2.3.3 Boreholes

Fifty-one boreholes (Table S2) were drilled on these watersheds or in the close vicinity (Figs. 1d and 2a). Transmissivity data are available for 19 boreholes (¹βol, ¹α, ²α, and ⁶B) and vary by 2 orders of magnitude between $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}$ and $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{3}}$ m² s⁻¹, with an average value of $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$ m² s⁻¹ (standard deviation: $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{4}}$ m² s⁻¹). Hydraulic conductivity varies between $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{7}}$ and $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹, with an average value of $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹ (standard deviation: $\mathrm{9}\times {\mathrm{10}}^{-\mathrm{6}}$ m s⁻¹). As aquifers are fissured or fractured with heterogeneities along the screen height, and as data were calculated by dividing transmissivity by the height of the screened saturated aquifer, calculated hydraulic conductivities have to be considered minimum values.

Piezometric-level measurements (Fig. 3a) show that the piezometric level is on average 7 m below ground level and shows that the hydrogeological functioning is not marked by a basal groundwater body with a low hydraulic gradient. In addition, two main typologies of aquifer are distinguished in Fig. 3b: on the one hand perched aquifers with springs located at altitudes above 400 m a.m.s.l. and on the other hand aquifers crossed by boreholes in the valleys with water levels close to the ground level.

Piezometric-level monitoring (Fig. 4) presents evidence of unconfined aquifers (piezometers 1, 2 and 4), with annual dynamics and a well-defined seasonal cycle (with fluctuations between 1 and 2 m): low groundwater levels occur during dry seasons (April to July) and high ones during rainy seasons (August to December). In contrast, piezometer 3 (situated 1 km above piezometer 2) characterizes a confined aquifer with multiannual dynamics, with a minor influence of the seasonal cycle.

Figure 4 Piezometric levels (monitored in the framework of the European piezometric network) and effective rainfall monitoring between 2005 and 2015. The first piezometer is located on the Case Navire watershed (1 in Fig. 2a), the next two are on the Fond Lahaye watershed (2 and 3 in Fig. 2a) and the last one is in Case Pilote (4 in Fig. 2a), 3 km west of Fond Lahaye. Piezometric levels of piezometers 2, 3 and 4 have been modified to fit on the same graph (−25, −53 and +10 m compared to their initial value, respectively).

2.3.4 Water balance

The water budget hydrological terms of the studied watersheds have been computed in Fig. 5 in order to show a synthetized view of the annual water balance and the contribution of each hydrological term. Rainfall and potential evapotranspiration are provided by the national meteorological agency for the period 1991–2015 (the annual rainfall map is shown in Fig. 1c). River discharge is monitored by the Ministry of the Environment. Real evapotranspiration and effective rainfall are 1 km² spatialized data calculated by Arnaud and Lanini (2014) (over the period 1991–2010), following a methodology detailed in Vittecoq et al. (2010) and based on the Thornthwaite model. The ratio runoff ∕ infiltration and groundwater contribution to river discharge have been calculated (1) for the Case Navire River by Vittecoq et al. (2007) (over the period 1987–1990) based on inverse modelling (e.g. Pinault and Schomburgk, 2006) with Tempo software (Pinault, 2001) and (2) for the Alma River by Stollsteiner and Taïlamé (2017) (over the period 2010–2015) based on lumped hydrologic modelling (e.g. Thiery, 2010) using Gardenia software (Thiery, 2014). Both methods used daily meteorological data series (rainfall, potential evapotranspiration) and river flow rates.

Figure 5 Annual water balance of the Case-Navire and Alma watersheds. Rainfall (R), real evapotranspiration (RET) and effective rainfall (Reff) are from Vittecoq et al. (2010) and Arnaud and Lanini (2014). Average annual river discharge (AARivD) and river-specific discharge (RSD) are calculated from gauging stations. Runoff: Reff contribution to river discharge. WD: volume of water for the drinking water distribution system. RivD: river discharge. RivTD: river total discharge ($=\text{RivD}+\text{WD}$). GwR: groundwater contribution to river discharge. Deepinf: deep groundwater flow.

The Alma watershed is the highest and smallest one, located upland, and is exclusively covered with tropical forest. This watershed is equipped with a gauging station with valid data since July 2010 (specific discharge of about 112 L s⁻¹ km⁻²). Water balance calculation (Fig. 5) shows that the difference between total effective rainfall and average annual flow in the Alma River is about 2.3×10⁶ m³ yr⁻¹ (18 % of effective rainfall volume). This volume of water (1) infiltrates in depth and/or (2) joins another stream/nearby hydrological watershed, if the hydrogeological catchment area differs from the topographic catchment. This volume infiltrated in depth or flowing towards an adjacent catchment area is therefore to be considered a minimum value, as measured rainfall gauges are situated at elevations not exceeding 600 m, whereas the watershed peak culminates at 1197 m. The national climatic agency (Météo-France) considers values up to 7000 mm yr⁻¹ of rain to be quite possible on the summits. Considering this highest value, and the various uncertainties in the water budget parameters, the deep infiltrated volume could reach a maximum of 8×10⁶ m³ yr⁻¹.

The Fond Lahaye watershed culminates at 532 m of elevation and its stream joins the sea 4 km downstream. Since there is no gauging station on this river, it is difficult to define a water balance. In the maximalist hypothesis where 100 % of the effective rainfall returns to the river, its maximum specific discharge would be about 11 L s⁻¹ km⁻² (corresponding to 10 % of the nearby Alma watershed-specific discharge).

The Case Navire watershed culminates at 1197 m of elevation and its stream joins the sea 10 km downstream. Its upstream part is divided into two sub-basins (Duclos and Dumauzé rivers) that meet 5 km before reaching the sea (Fig. 3b). Three dams are located in the upstream part of the Case Navire River (one on the Dumauzé River and two on the Duclos River; cf. Fig. 1). The annual volume of the three dams on arrival at the main distribution tank is 5.9×10⁶ m³ yr⁻¹ (over the period 2009–2012), corresponding to an average of 16 300 m³ d⁻¹ (and corresponding to 19 % of the annual effective rainfall). During the driest seasons, the river is often dry downstream of the dams, causing strong environmental impacts. The gauging station is situated on the Case Navire River a few hundred metres before reaching the sea (Figs. 1 and 2a), 5 to 6.5 km downstream of the dams, which allows calculation of water balances (Fig. 5). The supposed natural flow rate of the Case Navire River is about 18.7×10⁶ m³ yr⁻¹, corresponding to 60 % of the annual effective rainfall, by adding water abstraction volume by dams. Consequently, the volume of groundwater circulating in this watershed that does not return to the river is about 12.4×10⁶ m³ yr⁻¹. This volume may infiltrate in depth and circulates in the aquifers or to another watershed, or flows into the sea.

These water balance calculations show the main key component of the hydrological cycle of each watershed and provide first evaluations of the groundwater budget. In particular, they reveal significant quantities of deep infiltrated water (14.7 to 20.4×10⁶m ³ yr⁻¹) equal to 2 or 3 times the surface water intakes in the Case Navire River. There is thus a necessity to better understand aquifer nature and hydrodynamic characteristics, extension, thickness and groundwater preferential flows and interactions with rivers, and to locate recharge areas, in order to propose appropriate hydrogeological conceptual models necessary for sustainable management of water resources.

Figure 6 Internal resistivity and hydrogeological structure along three cross sections: C1, C2 and C3.

3 Heliborne TDEM method

Our methodology is based on a multidisciplinary approach combining geology, hydrogeology and a heliborne TDEM geophysical survey, in order to identify relationships between ground-based punctual geological and hydrogeological data on the one hand and in-depth geophysical information derived all over the area on the other.

3.1 The survey

A heliborne TDEM survey was conducted from February to March 2013 with the SkyTEM 304 system (Sørensen and Auken, 2004) over the entire Martinique. This survey, fully described by Deparis et al. (2014) and Vittecoq et al. (2015), was supervised by BRGM (the French geological survey) for geological and hydrogeological purposes. Over the studied watersheds, the survey was flown mainly along the N–S direction with 400 m line spacing and along the W–E direction with 4000 m line spacing. The spacing between each EM sounding along flight lines is approximately 30 m. In the lower part of the watershed fewer to no data have been acquired because of the urbanization. Finally, 13 596 TDEM soundings were processed in the study area. The TDEM method allows imaging of the conductivity–resistivity contrasts of the subsurface, inducing eddy currents in the ground (Ward and Hohmann, 1988). Locally, the depth of investigation (DOI) of the method depends on the emitted magnetic moment, the bandwidth used, the subsurface conductivity and the signal ∕ noise ratio (Spies, 1989). In this study, the average depth of investigation is around 150 to 200 m.

3.2 TDEM data processing

The ground clearance of the loop was obtained by degrading an available 1 m digital elevation model to a 25 m grid (consistently with the Airborne ElectroMagnetic (AEM) footprint) and subtracting it to the Differential Global Positioning System (DGPS) altitude of the frame; we did not use the data from the laser, which proved to be noisy in such a rough relief environment. Tilt measurements were processed by taking into account the local topography in order to consider an effective tilt at each TDEM data location (Reninger et al., 2015). As part of an environmental study in an entropized area, particular attention was paid to properly removing noise from the TDEM data. They were processed with a singular value decomposition (SVD) filter (Reninger et al., 2011). The SVD allows explanation of a dataset with only a few components, each datum being a linear combination of these components. Thanks to this decomposition we are able to identify and remove several types of noise, making the processing less time-consuming and subjective and reducing the amount of careful editing. In addition, a trapezoidal stack (Auken et al., 2009) was applied to the data. The trapezoid shape is consistent with the increase in the footprint of the EM method with time. The stack size was adapted to the noise level along flight lines. Thanks to this filter we manage to recover some noisy windows, which are unusable otherwise (Reninger et al., 2018). The aim of the applied processing was to keep as much resolution as possible (Reninger et al., 2018). Finally, manual editing was performed, mainly to remove remaining inductive–galvanic coupling noises. In order to improve the coverage of the dataset, good-quality portions of ferry lines were also considered during the processing (Reninger et al., 2018). Figures 1 and 2 display the position of the TDEM dataset after processing. Data were then inverted using the spatially constrained inversion (SCI) algorithm (Viezzoli et al., 2008). Each TDEM datum was interpreted as a 1-D earth model (EM sounding) divided into n layers, each one being defined by a thickness and a resistivity. During the inversion, constraints were applied vertically and horizontally to nearby soundings (independently of the flight lines and the ground clearance); weak constraints were applied for this study in order to limit the smoothing of the inversion procedure. Results were obtained with a smooth inversion (consisting of 23 layers from 0 to 200 m depth). This inversion method is effective for imaging complex geological structures with the lowest dependency on the starting model. In addition, the altitude of the transmitter was inverted for, and the DOI was evaluated as, a final step of the inversion (Christiansen and Auken, 2012).

Figure 7 Internal resistivity and hydrogeological structure along two cross sections: C4 and C5.

4 Resistivity profiles and correlations between resistivity, geological and hydrogeological data

Five resistivity profiles obtained by inverting TDEM data are provided in Figs. 6 and 7 (localization in Fig. 2). Confronting these profiles with geological and hydrogeological data (springs, boreholes, observations and outcrops, geological maps, etc.), the TDEM data can be interpreted in terms of geological or hydrogeological contrasts, and show the main internal geological structures and associated aquifers, at depths up to around 200 m. Thus, as shown by Vittecoq et al. (2014, 2015), geological and transmissivity data of each borehole can be compared to the closest TDEM sounding in order to get information on the resistivity of the aquifers and aquitards and better constrain their extension and thickness. However, in such a particularly rugged and contrasted environment, attention must be paid to how this comparison is achieved, mainly in terms of distance and elevation. This was done on 18 boreholes. They are located at an average distance of 35 m (with a maximum distance of 90 m) to the closest EM sounding, with a difference in elevation of less than 10 m. At each of these TDEM soundings, we looked at the average of the resistivity falling in each associated borehole geological formation. To complete aquifer characterization, a specific analysis was conducted on the springs. Resistivity values of the cells located upstream of the 24 springs (Table S1), corresponding to supposed aquifer formations, were manually extracted from the 3-D resistivity models.

Figure 8a displays borehole (BR) and spring (SR) aquifer resistivity ranges. Alluvial deposits display a relatively large resistivity range (12–74 ohm m) because of the heterogeneity of alluvial materials (in terms of granulometry, nature, ages, etc.). Except for alluvial deposits, a good correlation appears between resistivity and the age of the geological formations, showing the relationship between weathering process and resistivity: the older the formation, the lower its resistivity. Correlation between geology and hydrodynamic properties (Fig. 8b) also displays a trend: the older the formation, the higher its transmissivity or its hydraulic conductivity. In particular, the relatively large resistivity and hydraulic conductivity range for andesitic lavas ²α could be related to their intrinsic heterogeneity. These correlations are relevant for the three studied aquifers, within the interval 10–100 ohm m and within a range of 1 to 5.5 Ma. Below 10 ohm m, several authors (e.g. d'Ozouville et al., 2008; Pryet et al., 2012; Vittecoq et al., 2014) present evidence that very low resistivity layers can correspond to high-permeability formations saturated with saltwater (old confined water or seawater intrusion) or to impermeable clays resulting from meteorological or hydrothermal weathering processes. Beyond 100 ohm m, there are no boreholes on the studied watershed, with transmissivity or hydraulic conductivity values crossing formations with resistivity values higher than 100 ohm m. The correlation between age and hydraulic properties is valid for the same kinds of rocks (i.e. andesite and basalt lava flows in the context of a subduction zone volcanic arc island), but cannot be considered for domes. Indeed, eruptive mechanisms of andesite and basalt flows on the one hand, and intrusive domes on the other hand, are different, and domes are only observed in this area between 0.3 and 0.9 Ma.

Figure 8 (a) Borehole (BR) and spring (SR) resistivity ranges according to their lithological facies and age (Fig. 2). The younger the formation, the higher its resistivity. (b) Comparison between transmissivity, hydraulic conductivity and resistivity for three aquifer formations considering borehole values.

These correlations demonstrate the necessity and advantage of coupling hydrogeological data (springs, boreholes, etc.) and geological and geophysical data for an advanced interpretation of resistivity data, as such information is scarce in volcanic island environments and because resistivity data alone does not allow differentiation between age, nature of geological formations or aquifer identification.

5 Hydrogeological conceptual model

Our methodology and associated correlations allow identification and characterization of the main aquifer and aquitard formations (synthetized in Table 1) as well as their lateral extent and thickness, enabling the construction of a hydrogeological conceptual model at the hydrogeological watershed scale. This conceptual model, synthetized in Fig. 9, characterizes the structure and hydrogeological functioning of andesitic aquifers at the watershed scale and highlights the influence of geological structures on groundwater flows and the consequence for the interactions between rivers and aquifers. Joint analysis of water balance and geological structure also presents evidence of the differences between hydrological watersheds and hydrogeological watersheds.

Figure 9 Hydrogeological conceptual model of an andesitic complex in the subduction zone at a watershed scale.

5.1 The upper major perched aquifer of andesitic domes

The conceptual model is marked by the presence of andesitic domes and lava flows (⁹αbi), occupying the upper part corresponding to half of the Case Navire watershed and the entire Alma watershed. Water balance calculated on the Alma River suggests that 85 % of effective rainfall (Reff) infiltrates in these andesitic domes. Considering the high resistivity values of the domes (cf. Fig. 8a, spring resistivity analysis: 150–300 ohm m) and in comparison with other volcanic islands (d'Ozouville et al., 2008; Pryet et al., 2012; Vittecoq et al., 2014), it is assumed that these andesitic domes (⁹αbi^D) are highly fissured and fractured, conferring a high hydraulic conductivity to this aquifer (with an order of magnitude of $\mathrm{7}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹ by similarity to a borehole drilled in a dacitic dome 6 km north of the studied watershed). Given the rooting of endogenous domes within the volcano, and as shown thanks to the water balance calculation (Fig. 5), up to 40 % of Reff seeps in depth within dome roots, through fissures and fractures, and recharges underlying aquifers.

In addition, the unsaturated zone should present significant thickness, and since some springs have relatively low flow rates, we consider that they could emerge thanks to small and low hydraulic conductivity horizons, such as paleo-soils, geological heterogeneities (for instance between ⁹αbi and underlying formations), or structural discontinuities. The main rivers have their sources in this important perched aquifer, with significant flow rates (as shown in Fig. 5, for instance, the Alma River specific discharge is 112 L s⁻¹ km⁻²). On the western and eastern topographic ridges of the Case Navire River watershed, andesitic lava flows (⁹αbi^C) also constitute the first aquifer receiving rainfall and from which flow some non-perennial springs during the rainy season and few perennial springs during the dry season.

Table 1 Geological, geophysical and hydrogeological characteristics of the main aquifer and aquitard formations.

5.2 The lower aquifer of andesitic lavas

In this conceptual model the upper major perched aquifer, described above, underlies the second main aquifer of thick andesitic lavas ²α, marked by a relatively “smooth” morphology or paleo-topography of their top, consistent with the structure of lava cooling along a shield volcano. Hydraulic conductivity data dispersion over 2 orders of magnitude is in agreement with the heterogeneity of theses andesitic lavas. The various facies that were observed at the outcrop (Supplement S1) are (1) auto-brechified breccias and lavas, (2) massive facies more or less fractured according to the cooling rate of the lava, (3) facies with flow structures showing significant horizontal cracking parallel to the substratum, and (4) breccias and scorias associated with the base of the lava flow. Tectonic fracturing superimposes on these heterogeneities and can contribute to maintaining and developing the hydraulic conductivity of volcanic formations, as shown by Vittecoq et al. (2015). In this type of andesitic formation, boreholes can also be dry if no fissured or permeable zone is intersected.

The recharge of this aquifer is quite atypical as in the upper part of the Case Navire and Dumauzé watersheds, effective rainfall is high, and permeable andesitic domes and lava flows ⁹αbi overlay andesitic lavas ²α. As suggested by numerous springs in the limit of extension of ⁹αbi, and the low flow rates observed in the two Fond Baron boreholes screened into andesitic lavas ²α (Senergues, 2014), effective rainfall infiltration into ²α should be limited by paleo-soils and/or the hydraulic conductivity contrast between the two formations, acting as semi-permeable hydraulic obstacles. In the lower part of the watersheds, effective rainfall is limited (200 to 800 mm yr⁻¹) compared to the upper part, and furthermore the plateaus located on both sides of the rivers are overlain by low hydraulic conductivity breccias. Effective rainfall infiltration towards andesitic lavas ²α is thus also small in the lower part of the watersheds. Then, the recharge of this aquifer should follow four main steps. Firstly, a part of effective rainfall (18 % to 40 %, depending the watershed, as shown in Fig. 5) deeply infiltrates through the fractures and in the rooting of andesitic domes ⁹αbi. Secondly, as andesitic lavas ²α were crossed through faulting by ⁹αbi lavas, this deeply infiltrated water then flows deeper towards andesitic lavas ²α, thanks to geological heterogeneity inside the old volcanic chimney. Thirdly, groundwater flows into andesitic lavas ²α and, lastly, the ²α aquifer, incised by the river, allows this deeply infiltrated water to be drained by the river and the sea.

5.3 The regional aquitard

Hyaloclastites ¹H, mainly observable in Fig. 6 (C3) and Fig. 7 (C4 and C5) at altitudes below 100 m a.m.s.l., are the lower boundary of the watersheds and, more generally, of a major northern part of the island. In Fig. 6 (C2), they are suspected between 200 and 300 m a.m.s.l. to the east of the cross section, probably due to the displacement generated by major faults: this topographical limit is interpreted by Boudon et al. (2007) as the eastern limit of a large flank collapse with a horseshoe-shaped structure opened westward. The weathering grade observed on the outcrop in the Case Navire River, associated with their very low resistivity, leads us to consider the hyaloclastites mainly as a very low permeable formation and are then considered the regional aquitard.

5.4 Difference between hydrological watershed and hydrogeological watershed

The continuity of andesitic lava flows ²α along the resistivity cross sections (Fig. 7), from north to south and especially under the “Morne Jeanette” (C5), clearly suggests a continuity of groundwater flows, through andesitic lava-permeable facies, beyond the Duclos River watershed and in the direction of the Fond Lahaye watershed. This hypothesis of a clear difference between hydrological watershed and hydrogeological watershed is supported by (1) the piezometric fluctuations (Fig. 4), showing that the Fond Lahaye upper borehole is in a captive aquifer with multiannual dynamic fluctuations, (2) groundwater mineralization and long-duration time transfers (>50 yr by CFC groundwater dating, Gourcy et al., 2009) and (3) the high flow rates of the Fond Lahaye and Case Navire boreholes (more than 1.2×10⁶ m³ yr⁻¹ have been calculated by Ollagnier et al., 2007; Vittecoq et al., 2008; Vittecoq and Arnaud, 2014).

5.5 Geothermal insights

The very low resistivity (6–10 ohm m) of hyaloclastites ¹H cannot correspond to actual saltwater intrusion, as they are situated higher above sea level. Their very low resistivity could rather result from weathering during the 1 Ma period of rest before the next volcanic phase and from hydrothermal weathering. This low resistivity layer (<10 ohm m) should indeed be evidence of a smectite-bearing hydrothermally altered caprock (e.g. Browne, 1970; Simmons and Browne, 1990) of an underlying geothermal system. The two thermal springs (Didier – 32 ^∘C, 1850 µS cm⁻¹ and Absalon – 36 ^∘C, 1730 µS cm⁻¹) could be leaks of this geothermal system, through faults allowing the rise of mineralized gaseous waters (it must be noted that the supposed fault interpreted in Fig. 6 (C2 – 2200 m) is aligned with Didier springs, Absalon springs and the Alma and Dumauzé domes).

Then, geothermal fluid circulations could follow five steps: (1) deep infiltration of effective rainfall through andesitic domes (⁹αbi) and associated deep rooting, (2) deep mixing at temperatures between 100 and 140 ^∘C, according to Gadalia et al. (2014), (3) interaction with CO₂ and ascent along faults, (4) mixing with andesitic aquifer ²α, and (5) emergence in thermal springs. The flow rate of these springs being relatively low, we can state that a part of the ascending enriched fluids do not emerge at the surface and diffuse in the andesitic aquifer ²α. The higher groundwater mineralization downstream (1000 µS cm⁻¹ in Fond Lahaye boreholes), compared to the range of water electrical conductivity of cold springs (50–350 µS cm⁻¹) emerging from the aquifers upstream (cf. Table 1), clearly support this hypothesis.

6 Discussion

Heliborne TDEM data reveal in-depth resistivity contrasts. Their interpretation with borehole and spring data allowed constraint of a detailed hydrogeological conceptual model. Working at the watershed scale brings new elements of hydrogeological functioning of andesitic volcanic complexes. Vessell and Davies (1981), Bogie and Mackenzie (1998) and Selles et al. (2015) proposed a geological conceptual model of andesitic stratovolcanoes presenting evidence of central (0–2 km from the vent), proximal (5–10 km), medial (10–15 km) and distal (15–40 km) facies variations. The originality of our work is to focus on improving the hydrogeological functioning of the central and proximal parts of such an andesitic system. Indeed, medial and distal parts, on which hydrogeological studies are generally focused on continental volcanoes (Selles, 2014), corresponding to lower and accessible areas, are in our case under the sea.

On the scale of the island of Martinique, the proposed hydrogeological functioning conceptual model (and also our methodology) could likely be extended to the other watersheds situated on the Carbet volcanic complex and on the Morne Jacob shield volcano. Extrapolation to the entire Martinique is nevertheless not considered, as a specific hydrogeological functioning has been demonstrated for the centre of the island (Vittecoq et al., 2015), as effective rainfall is significantly lower (<1500 mm) in the central–southern half of the island (Vittecoq et al., 2010) and because our conceptual model, concerning mainly fissured and fractured lava, could not fit with the Mount Pelée stratovolcano located in the north of Martinique (covering 15 % of the Martinique area) and constituted by pyroclastic flows (Traineau et al., 1989).

Secondly, our conceptual model could also enhance, with new insights, existing characterizations of the hydrogeology of small volcanic islands and especially West Indies and Caribbean volcanic islands (Unesco, 1986; Falkland and Custodio, 1991; Davies and Peart, 2003; Robins, 2013). However, Hemmings et al. (2015), studying Montserrat, an andesitic island located in the Lesser Antilles, concluded that they did not know which model of the Hawaiian and Canarian models could fit Montserrat. As both concerned only basaltic islands, with different geological structure and recharge conditions, other models have to be proposed. Our conceptual model, thanks to the high-resolution heliborne geophysical survey and correlations with geological and hydrogeological data, could then help better understand the hydrogeological functioning of other Lesser Antilles andesitic islands. For instance, the 3-fold division of the West Indies hydrogeological classification by Robins et al. (1990) could be updated with a fourth category considering groundwater in permeable perched high-rise volcanic domes and in underlying fractured volcanic rocks.

The main geological structures highlighted lead to preferential flow circulations and to a non-adequacy between hydrogeological and topographical watersheds, as supposed by Charlier et al. (2011) at a smaller scale (45 ha) in Guadeloupe. The consequence is that even if the topographical watershed is small, underground flow circulations can add significant amounts of water to the river watershed's water balance, if aquifers are situated above (in elevation or upstream). We thus support the necessity to include and characterize neighbouring watersheds to extend our methodology and results to other areas or islands. This can be even emphasized if lava domes and associated roots are situated upstream, as they present very high hydraulic conductivity and preferential flow circulations vertical in depth.

Thanks to the interpretation of the geological, geophysical and hydrogeological data, we highlight, for the present study (i.e. the watersheds and the three studied aquifers, within the interval 10–100 ohm m and within a range of 1 to 5.5 Ma), that (1) the older the formation, the lower its resistivity and (2) the older the formation, the higher its transmissivity or hydraulic conductivity. This last result is also consistent considering the results of Vittecoq et al. (2015) obtained on an older aquifer (15 Ma) on Martinique, with higher hydraulic conductivity and lower resistivity than the ones observed in the present study. Consequently, unlike hotspot basaltic islands (Custodio, 2005; Vittecoq et al., 2014), hydraulic conductivity of the studied aquifers of subduction zone andesitic volcanism does not decrease with age. On the contrary, our results show an increase with age. Nevertheless, time itself is not the activating factor, and only a few geological processes can cause an enhancement in hydraulic conductivity. Given (1) the tectonic and seismic context of the subduction zone, (2) the fact that earthquakes are known for increasing hydraulic conductivity (e.g. Rojstaczer et al., 1995; Ingebritsen et al., 2006) and (3) the fact that earthquake-induced modification of hydraulic conductivity have been observed in Martinique (Lachassagne et al., 2011), we interpret the observed hydraulic conductivity increase as the consequence of earthquake tectonic fracturing.

The accuracy of correlations between boreholes and TDEM soundings is highly dependent on the distance to the nearest TDEM flight line. Accordingly, particular attention must be paid to the way this comparison is achieved, mainly in terms of distance and elevation difference. This being said, heliborne geophysical survey is certainly the best cost-efficiency method, and probably the only method providing this density of data down to 200 m depth, allowing a detailed geological and hydrogeological characterization at this working scale. Nevertheless, a minimum of ground-based geological and hydrogeological data are necessary, thanks to borehole data.

7 Conclusions

From an operational point of view, our data and results should be very helpful for local stakeholders facing environmental impacts and overexploitation of the Case Navire River. We show that large volumes of water infiltrate and flow in several aquifers. Sustainable management of water resources will require a better repartition between rivers and aquifers. Aquifers, and especially downstream the watersheds, could be exploited in order to decrease the use of the dams, especially in dry seasons. Future drilling programmes could be launched considering our conceptual model. We also provide some insights into potential geothermal resources such as the pathway of deep infiltrated water through the roots of the andesitic dome, the presence of a low-resistivity regional aquitard and the link with the thermal springs.

In conclusion, our multidisciplinary approach and results allow characterization in detail of the hydrogeological functioning and characteristics of the main aquifer and aquitard units, leading to the proposition of a hydrogeological conceptual model of an andesitic island at the watershed scale, presenting evidence of the key role of geological structures and volcanic domes in groundwater flows. We also demonstrate, for the studied geological formations, that hydraulic conductivity increases with age in this andesitic-type volcanic island. Moreover, the working scale seems particularly suitable due to the complexity of edifices, with heterogeneous geological formations presenting high lateral and vertical variability. Andesitic-type volcanic islands being little known and studied, our work offers, in addition to the proposed conceptual model and thanks to the high-resolution heliborne geophysical survey, new guidelines for accurate correlations between resistivity, geology and hydraulic conductivity for other volcanic islands.