M Rajaobelison
email: miora_mirah.rajaobelison@inrs.ca

J Raymond

M Malo

C Dezayes

S Larmagnat

Understanding heat transfer along extensional faults: the case of the Ambilobe and Ambanja geothermal systems of Madagascar

Keywords: Energy, Geothermal, Heat transfer, Faults, Numerical modeling, Madagascar Energy, Geothermal, Heat transfer, Faults, Numerical modeling, Madagascar

Understanding the role of faults where forced convective heat transfer is the dominant mechanism giving rise to hot springs is critical in geothermal exploration in extensional environments. This study uses two-dimensional models of coupled fluid flow and heat transfer along cross-sections perpendicular to faults and the regional topography to identify favorable fault conditions for geothermal system development in northern Madagascar.

Structural data collected at surface were used to define fault scenarios and simulate the ascension of hot fluids to reproduce hot spring temperatures in the Ambilobe normal fault zone area and the Ambanja graben structure. Fault dips facing topography-driven groundwater flow was shown to be favorable, and hot spring temperatures could be reproduced when the fault permeability was >10 -14 m 2 . Faults located in a discharge zone near a river were the most favorable for fluid ascension, regardless of their dip.

Constraining the model with a basal heat flow between 90 and 148 mWm -2 at a depth of 10 km allowed the reservoir temperature to reach 150-200 °C at depths of 2 km or shallower along favorable faults.

Introduction

Economic activity in the northern region of Madagascar relies on export products such as vanilla, cocoa and coffee that are grown and transformed locally. Madagascar is one of the largest vanilla exporters in the world (UNComtrade, 2018). The energy used to transform these products relies heavily on fossil fuels and wood. Off-grid communities, mainly in remote areas, use traditional biomass such as firewood and charcoal as their main source of energy, inducing deforestation and water resource degradation, along with biodiversity and soil loss, as has already happened in many Madagascar communities. Renewable sources, such as geothermal energy, are expected to play a more prominent role in future energy production.

Ambilobe and Ambanja, two areas in northern Madagascar, have been characterized for their geothermal potential and respectively classified as a graben border-fault liquiddominated moderate-temperature system and a fossil magmatic liquid-dominated moderate-temperature system (Rajaobelison et al. 2020). In this classification scheme based on the extensional type, the fault and lithologic-diagenetic characteristics affecting porosity are thought to control formation permeability. As reported by Rajaobelison et al. (2020), the average matrix permeability of fracture-free rocks in these two areas was evaluated between 10 -13 and 10 -18 m 2 , which is low, while the average porosity is ~5.7 %, which is low or poor (5-10 %) according to Levorsen and Berry (1967); Tiab and Donaldson (2016). However, this matrix permeability has to be considered the lower bound of formation permeability as fractures are common in such an extensional tectonic setting.

The occurrence of thermal springs along a major fault zone in Ambilobe and a graben structure in Ambanja confirms that advective heat transfer can take place within fractured rocks (Rajaobelison et al. 2020). The equilibrium formation temperature at depth, estimated using chemical geothermometers (Gunnlaugsson et al. 1981), ranges from 140 to 200 °C, high enough to be used for scalding vanilla beans (~60 °C), heating greenhouses (> 40 °C;Lindal 1973), and generating electricity (> 80 °C; Grasby et al. 2011;Tomarov and Shipkov 2017). Fluids from geothermal sources with temperatures below 120 °C can be used to generate electricity using an optimized Organic Rankine Cycle (e.g. [START_REF] Chagnon-Lessard | Optimal design of geothermal power plants: A comparison of single-pressure and dual-pressure organic Rankine cycles[END_REF]Liu et al. 2017;Shi and Pan 2019).

No exploration boreholes have been drilled in northern Madagascar, so the optimal depth for finding geothermal sources remains unknown, as do the geothermal gradient and heat flux density. Previous work classified geothermal play types and evaluated rock thermohydraulic properties from surface outcrops to better understand heat transfer mechanisms at depth (Rajaobelison et al. 2020;Rajaobelison et al. 2021). Crustal-scale faults have been identified as potential targets for geothermal exploration, but further work is needed to better understand their role in the subsurface heat transfer mechanisms that control the depth and extent of geothermal sources. Fault influence must be taken into account when exploring geothermal energy sources, especially the influence of the large-scale faults identified by surface mapping and shown on the geological maps of Madagascar by Roig et al. (2012) and Tucker et al. (2014). Therefore, numerical modeling was used to simulate groundwater flow and heat transfer according to various fault scenarios affecting the geothermal energy sources associated with petrothermal systems in this extensional domain of Madagascar.

Numerical modeling has been previously used to simulate the depth distribution of temperature to better understand fractured geothermal reservoirs (e.g. Bakhsh et al. 2016; Quaternary direction of extension (N65°E to N85°E) observed in the volcanic region of Montagne d'Ambre, in the extreme north of Madagascar [START_REF] Chorowicz | An approach using SAR ERS images to relate extension fractures to volcanic vents: examples from iceland and Madagascar[END_REF]) enclosing the Ambilobe and Ambanja study areas.

The Ambilobe area is located in the Andavakoera fault zone, a major NE-SW normal fault delineating the contact between sedimentary sequences and basement rocks of the northern Bemarivo domain. Numerous undifferentiated NW-SE faults and NW-SE strike-slip faults with dextral or sinistral movements can be seen on the geological map (Figure 2b) and in the field (Figure 3a). The strike-slip faults include the faults near the Betsiaka hot springs and undifferentiated normal faults of Ranomafana and Ankatakata (Figure 2b). The main fault and fracture systems are dominated by NNW-SSE, NE-SW and N-S orientations.

Fracture orientations from field measurements and fault orientations from geological map are shown on the rose diagram of Figure 2b.

The Ambanja area is located in the graben structure in the Proterozoic basement units of the southern terrane of the Bemarivo domain. The graben is oriented NW-SE, and its average width is 4 km. The main faults and fractures dominantly strike NNW-SSE. The rose diagrams of fracture orientations from field measurements as well as the faults orientations from the geological map are shown on Figure 2c. Some metric fractures were observed at the outcrop scale (Figure 3b)

Favorable structural patterns

In a comprehensive inventory of the structural settings of geothermal systems in the Great Basin region of western USA, [START_REF] Faulds | Favorable tectonic and structural settings of geothermal systems in the Great Basin region, western USA: Proxies for discovering blind geothermal systems[END_REF] identified the most favorable dominant faulting patterns to host geothermal activity in the extensional to transtensional terranes, and more specifically in areas where the bulk of the geothermal resources are likely hidden. Applying the same principles to the Ambilobe and Ambanja areas, the identified favorable fault structures are: 1) major normal fault segments, i.e., near the displacement maxima of the NE-SW Andavakoera fault and the NW-SE fault-borders of the Ambanja graben structure intersections between two normal faults (Figure 2b and2c),

and 2) or between normal faults and transverse oblique-slip faults (Figure 2c). At the surface, hot springs reach a temperature varying from 50 to 72 °C in both areas [START_REF] Besairie | Contribution à l'étude des sources minérales et des eaux souterraines de Madagascar[END_REF]Gunnlaugsson et al. 1981). Waters from the hot springs have high Mg concentrations and fall into the immature field on the ternary Na-K-Mg diagram of Giggenbach (1988), indicating dilution of the geothermal fluid along its path to the surface (Figure 4). The geothermal fluid thus mixed with meteoric and shallow groundwater.

Using chemical geothermometers (Gunnlaugsson et al. 1981) and the Giggenbah diagram (Figure 4), the reservoir temperatures in the Ambilobe and Ambanja areas are expected to be 140 and 200 °C, respectively.

Figure 4:

Chemical composition of waters from the Ambilobe (AMBL) and Ambanja (AMBJ) hot springs plotted on a ternary Na-K-Mg diagram (Giggenbach 1988). The Na, K, Mg concentrations were reported in Gunnlaugsson et al. (1981).

Lithosphere structure and geological formations

In the Bemarivo domain, the Moho is at a depth of 33 km (Andriampenomanana et al. 2017). The thickness of the upper crust is up to 20.5 km, and its composition is assumed to be primarily paragneiss and sedimentary rocks which have an average thickness of 4 km.

The lower crust is assumed to be composed of mafic rocks with an average thickness of 12.5 km.

The geological formations of the Ambilobe and Ambanja areas and their thermo-hydraulic properties have been described in Rajaobelison et al. (2020);and Rajaobelison et al. (2021).

The Permo-Triassic sedimentary rocks of the Ambilobe Basin, overlying the Proterozoic basement of the Bemarivo domain, are the main formations in the Ambilobe area. The Proterozoic basement rocks of the Bemarivo domain are the only geological formations exposed in the Ambanja area (Figure 2).

Development numerical heat transfer model

Conceptual model and assumptions

The geometry of the models (Figure 5) accounts for the topography on NNW-SSE (AA')

and NE-SW (BB') cross-sections (Figures 2b and2c) that pass through, respectively, the Andavakoera normal fault zone at the contact between sedimentary sequences and basement rocks in the Ambilobe area, and the graben structure in the Ambanja area. The dominant sedimentary and basement units have been considered in the models. The crosssections are perpendicular to the main undifferentiated faults observed near hot spring occurrences. The faults represented in the models are N010° Ankatakata post-Triassic normal oblique-slip faults in the Ambilobe area and N140° border faults in the Ambanja graben (Figures 2b and2c). These faults appear on the geological map of Madagascar produced by Roig et al. (2012).

Different modeling scenarios were defined to investigate the influence of fault dip on groundwater flow and heat transfer processes (Table 1). The dips in the two models range from 60° to 85° since dips from 60° to vertical have been shown to cause variations in the pressure and temperature of subsurface fluids [START_REF] Cherubini | Impact of single inclined faults on the fluid flow and heat transport: results from 3-D finite element simulations[END_REF]). These modeling scenarios were based on the geological maps of Ambanja and Ambilobe, although many uncertainties remain about fault dip and down-dip continuity at depth, which justifies the use of various scenarios. The models simulate the potential temperature distribution for the Ambilobe and Ambanja geothermal areas according to fault dip influence.

Governing equations

The models were developed using the subsurface flow module with a finite element method in the COMSOL Multiphysics software. The module couples groundwater flow simulated with Darcy's law and heat transfer in porous media. This made it possible to simulate forced convective heat transfer, which occurs along major faults where the host rock is considered an equivalent porous medium. Steady-state flow and heat transfer were numerically resolved according to the partial differential equation 1 and the input parameters summarized in Table 2. Sedimentary rocks and granitic and paragneiss basement rocks were considered low-permeability 2D units using the equivalent porous medium, where fluids can circulate. Faults are discrete 1D linear elements superposed on the 2D porous medium, where flow and heat transfer are solved simultaneously. Coupling between heat transfer and fluid motion occurs through the velocity field.

𝜆 b ( 𝜕 2 𝑇 𝜕𝑥 2 + 𝜕 2 𝑇 𝜕𝑧 2 ) -𝜌 f c f ( 𝑢 𝑥 𝜕𝑇 𝜕𝑥 + 𝑢 𝑧 𝜕𝑇 𝜕𝑧 ) + 𝐴 = 0 (1)
where λb (Wm -1 K -1 ) is the bulk thermal conductivity; T (°C) is absolute temperature; 𝑢 (m s -1 ) is the velocity vector determined from simulations of fluid flow quantified with Darcy's law, ρf (kg m -3 ) is the density of the fluid, cf (J kg -1 K -1 ) is the heat capacity of the fluid, and A (10 -6 W m -3 ) is the radiogenic heat production.

Geometry and subsurface properties

Geometry

 Ambilobe

The model has a width of 10 km and a depth of 10 km (Figure 5a). The sedimentary units dip 20-30° to the NE. Table 2 details the material properties of the 5 lithological domains.

The reverse fault F marks the contact between the mica schist and orthogneiss formations, and F1 and F2 are the mains faults represented in the models.

 Ambanja

The model has a rectangular geometry 20 km wide and 10 km deep (Figure 5b). Table 2 details the material properties of the 3 lithological domains. The properties of the border faults of the graben, represented by F1 and F2, are also included in the table.

Thermohydraulic properties and internal heat generation

The thermohydraulic properties of the model units were based on average properties measured in the laboratory for each geological formation reported by Rajaobelison et al. (2021), as summarized in Table 2. They are assumed to be constant with depth. Porosity ∅ and permeability k (m 2 ) were evaluated with transient measurements using the combined AP-608 gas permeameter-porosimeter from Core Test Systems with nitrogen as the gas source. The average bulk density ρ (kg m -3 ) was obtained from porosity assessments. Transient thermal conductivity λ (Wm -1 K -1 ) and diffusivity α (m 2 s -1 ) were measured with the optical scanning method using an infrared scanner developed by LGM Lippmann. The heat capacity, c (J kg -1 K -1 ), was deduced from the thermal conductivity and thermal diffusivity assessment (Popov et al. 2016;Popov et al. 1999).

The internal radiogenic heat generation (A; 10 -6 W m -3 ) for each rock presented in Table 2 was estimated from the concentration of the radiogenic elements U, Th and K, as reported in Rajaobelison et al. (2021).

Fault properties

This study assumed fault aperture and permeability from 250 µm to 0.2 cm and from 10 -14 to 10 -7 m 2 , respectively (Table 2). These ranges are based on the literature, taken from various in-situ hydraulic tests performed at geothermal sites (e.g. Holzbecher et al. 2011;Jeong et al. 2001) and examples from numerical simulation studies of groundwater flow in discrete fracture models (e.g. [START_REF] Bisdom | The impact of in-situ stress and outcropbased fracture geometry on hydraulic aperture and upscaled permeability in fractured reservoirs[END_REF]Guillou-Frottier et al. 2013). The minimum fault aperture range corresponds to that observed under a normal stress of 20 Mpa in a laboratory experiment and which validated the cubic law for laminar flow of fluids through open fractures (Witherspoon et al. 1980)

Mesh

Refinement of the mesh increases from bottom to top and close to the faults to avoid convergence problems. These locations require more accuracy due to large flow and temperature contrasts (Figure 6). A non-structured mesh with triangular elements was used, and the mesh statistics are summarized in Table 3.

Figure 6

Finite element mesh of the 2D numerical models for a) the Andavakoera normal fault zone at the basement-sediment contact in the Ambilobe area, and b) the Ambanja graben. A mesh convergence study was undertaken to verify that the temperature solution is independent of the element size; i.e., a chosen cutline was sampled uniformly from the models at different mesh sizes. The parameter σ was then computed according to equation 2:

𝜎 = √ 1 𝑁 ∑ (∆𝑇 𝑖 ) 2 𝑁 𝑖=1
(2

)
where N is the sample number, ∆Ti is the temperature variation at depth xi=(i-1)D/N-1 along the chosen cutline where the mesh size of the model is refined. Here, D is the depth of the model. It demonstrated that a stable solution was obtained above 10 4 elements with a difference of about 1 °C or less (Figure 7).

Figure 7

Mesh convergence study for the a) Ambilobe and b) Ambanja models.

Boundary conditions

Groundwater flow conditions

Topography-driven groundwater flow affecting forced convection was reproduced by imposing lateral flow boundary conditions (López and Smith 1995;Smith and Chapman 1983;Tóth 2009;Woodbury and Smith 1985). Without any information about the water table level, the lateral boundaries of the model are assigned a constant hydraulic head equal to the topographic elevation (Figure 5). Each segment representing a river in the model valley is also assigned a constant hydraulic head value equal to its nodal elevation. An outflow boundary was assigned at the upper-end of the faults to allow fluids to exit at the upper-end of the faults. The base of the models was assumed impermeable.

Heat transfer conditions

A constant temperature was specified at the surface of the models and assumed to be a linear function of the altitude. Equation 3 establishes the relationship between the surface temperature (Ts) and altitude (z): The lower boundary of the models, at a depth of 10 km, was assigned a second type condition of constant heat flow (Table 5), defined according to the Curie point depth (CPD; Tanaka et al. 1999) and the internal heat generation in the crust. Heat flow at the CPD (QCDP, mWm 2 ) was initially calculated using equation 4, which takes into account the surface (T0 =20 °C) and Curie (TCDP=580 °C) temperatures:

Ts= -0.0046𝑧 + 26.85 (
𝑄 𝑐𝑝𝑑 = 𝜆 𝑇 CPD-𝑇 0 𝑍 CPD -𝐴 𝑍 CPD 2 (4)
where A (10 -6 W m -3 ) is the average heat generation due to the radioactive decay of elements in the crust, which is the weighted average of the lower and upper crust heat generation. The internal heat generation in the upper and lower crusts was estimated at 1.5 × 10 -6 W m -3 and 0.72 × 10 -6 W m -3 , respectively (Table 2; Horton et al. 2016).

λ (W m -1 K -1 ) is the average thermal conductivity of the crust and ZCPD (m) is the CPD.

The range of the CPD was estimated at 10 km to 33 km, corresponding to the average Moho depth in northern Madagascar, using the centroid method of Tanaka et al. (1999) and according to the available world map of the CPD (Li et al. 2017).

The heat flow at 10 km (Q10, mWm -2 ) was then calculated with:

𝑄 10 = 𝑄 cpd + 𝐴 LC 𝐻 LC + 𝐴 UC (𝐻 UC -𝐻 model ) (5)
where HLC and HUC are the thickness of the lower and upper crusts, respectively, Hmodel is the thickness of the model, and ALC and AUC are the heat generation of the lower and upper crusts, respectively (Table 2).

Adiabatic heat transfer conditions were applied to the lateral boundaries.

Parameter sensitivity

The sensitivity of the model to the uncertain parameters was evaluated using the parametric sweep function implemented in COMSOL Multiphysics. The value of the parameters was set from the minimum to maximum. These parameters are basal heat flow, fault aperture and fault permeability (Tables 2 and5). Combining these three parameters with the different fault dip scenarios (Table 1) yielded 96 simulations for the Ambilobe 2D model and 81 for Ambanja. Results and analysis

2D subsurface temperature distribution and flow velocity

The different scenarios shown in Figures 8 and 10 are for base conditions using minimum, median and maximum heat flow values and average fault properties (Tables 1 and2). The simulation results confirm that the highest temperatures occur when heat flow is maximum, and that fluids and heat ascend along faults compared the fault-free base case (scenario 0).

Temperature profiles as a function of depth along fault lines (Figures 9 and 11

Figure 11

Temperature profiles as a function of depth along fault lines for the Ambanja 2D model. The temperature profiles in scenario 0 were obtained from an average cutline drawn in the fault-free model.

The magnitude of the Darcy velocity varies between 10 -9 and 10 -2 m s -1 along the faults (e.g. Figure 12). The simulated outflow velocity at the upper-ends of favorable faults (Table 6 and Figure 12) shows no difference or a 1-2 magnitude difference compared to the measured flow velocity at some hot springs, which varies between 10 -9 and 10 -7 m s -1 (Table 6; Gunnlaugsson et al. 1981). Simulated outflow velocity at upperends of fault (m s -1 )

6.47  10 -7 to 1.21  10 -6

7.07  10 -7 to 1.11  10 -6

9.63 10 -7 and 1.24  10 -6 to 2. 6 10 -7 2.14 10 -7 and 2.72  10 -7 to 1.24  10 -6

2.32  10 -8 to 1.19  10 -3 9.37  10 -5 to 9.74  10 -4

Figure 13

Variation of simulated flow velocity and temperature as a function of depth along the fault. Ambilobe 2D model, scenario 3, fault F2 and median heat flow value.

Influence of faults

Faults play a fundamental role in the movement of fluids and pore pressure distribution and therefore exert considerable influence over the temperature distribution. This is evident when comparing scenario 0 (without faults) to faulted scenarios (Figures 8 and10). The faults serve as conduits for fluid upflow, with Darcy velocity vectors converging upward along the faults toward the hot springs location. The temperature along faults can be up to 150 °C higher than adjacent areas (Figure 8 and 10). Heat transfer becomes conductive in the surrounding host rock, where bulk permeability is lower than the high-permeability zones because fluids drain towards the fault where convection occurs. Faults create permeability heterogeneities in the crust, which affect fluid flow and convective heat transfer, thereby influencing the location of hot springs [START_REF] Faulds | Characterizing structural controls of geothermal reservoirs in the Great Basin, USA, and Western Turkey: Developing successful exploration strategies in extended terranes[END_REF][START_REF] Faulds | Favorable tectonic and structural settings of geothermal systems in the Great Basin region, western USA: Proxies for discovering blind geothermal systems[END_REF], especially at their intersections (Belgrano et al. 2016;[START_REF] Person | Analytical and numerical models of hydrothermal fluid flow at fault intersections[END_REF]. This is consistent with the numerical modeling of deep fluid flow associated with crustal-scale normal faults in the south-central Canadian Cordillera [START_REF] Ferguson | What do aqueous geothermometers really tell us?[END_REF]). The latter reported that the permeability contrasts between host-rocks and faults acting as conduits for preferential fluid flow, have the greatest influence on the subsurface temperature distributions consistent with observed thermal springs. However, the presence of a fault does not necessarily imply the formation of a hot spring unless the favorable conditions are present that allow fluids and heat to ascend and reach the surface (e.g. López and Smith 1995;Thiébaud 2008) Simulation scenario 3 for the Ambilobe 2D model clearly illustrates favorable conditions for fluid ascension, with higher simulated temperatures near the upper-end of both faults (Figure 8; Table 1) and fault dips of < 90° facing the fluid flow. As shown in Figure 13, the fluid at the surface, reaches the highest Darcy velocity (10 -6 m s -1 ) which match the temperature of the hot spring (50 °C). A flow velocity of 10 -8 m s -1 is found at the depth having the reservoir temperature of 180 °C.

On the other hand, if the fault is facing away from the fluid flow direction, like F1 in scenarios 1 and 3 for the Ambanja graben model (Figure 10; Table 1), then the hydraulic gradient between both ends of the fault is smaller and the simulated temperature near the upper-end of the faults is lower. Similar results were reported by [START_REF] Moreno | Fault characterization and heat-transfer modeling to the Northwest of Nevado del Ruiz Volcano[END_REF], who simulated the impact of faults on the circulation of fluids originating from the basement underlying the Nevado del Ruiz volcanic complex of Colombia, where hot springs are present. However, scenarios 1 and 2 for the Ambilobe 2D model showed that fluids are able to rise to the surface along a fault with a dip facing away from the fluid flow direction (F2 and F1, respectively) if it is close to a river that acts as a discharge zone (Figure 8 and9). This is also seen in scenario 4. It suggests that both the fault dip and its distance from a discharge zone where the hydraulic head can be low are important characteristics in favor of fluid ascension.

Additionally, convective heat transfer along faults implies that ascending fluids are diluted as they mix with groundwater of shallower origin, to the degree that half of the ascending fluids will reach the upper 75 % of the fault length [START_REF] Ferguson | What do aqueous geothermometers really tell us?[END_REF]. This can explain the finding of significant mixing of geothermal fluids with shallower groundwater (Giggenbah diagram of Figure 5) and the absence of boiling hot springs in the Ambilobe and Ambanja geothermal fields due to cooling along the fault between the reservoir depth and the surface.

Parameter sensitivity

For each scenario, we investigated the variation in depth at which reservoir temperature is attained as a function of basal heat flow, fault aperture and permeability (Figures 15 and16). The results reveal that the depth of reservoir temperature decreases with increasing heat flow, fault aperture and fault permeability. For example, in scenario 1 (Figure 14a), with the minimum fault aperture (Df) value and maximum fault permeability (kf), the depth of reservoir temperature decreases from >2500 to 200 m along fault F2 when heat flow increases from the minimum to the maximum. For the same scenario, using maximum heat flow and mean fault permeability, the depth of reservoir temperature decreases from 1400 to 1100 m along fault F2 by increasing the value of fault aperture. In scenario 3 (Figure 14b), using minimum fault permeability and maximum fault aperture, the depth of reservoir temperature decreases from >2500 to 750 m along fault F2 as function of the heat flow value. For the same scenario, using maximum heat flow and mean fault permeability, the depth of reservoir temperature decreases from 800 to 400 m along fault F2 by increasing the fault aperture. The least favorable cases in scenarios 1 and 2 are when minimum heat flow fails to produce any points on the plots because the reservoir temperature is always deeper than 2500 m despite increases in fault aperture and permeability. The parameter sensitivity trend described for different fault dip scenarios in the Ambilobe model is the same for the Ambanja model (Figure 15). The least favorable cases correspond to scenarios 1 and 2 using the minimum heat flow value because the reservoir temperature is deeper than 2500 m despite increasing fault aperture and permeability.

We assessed the influence of basal heat flow, fault aperture and fault permeability given the uncertainty associated with those inputs. Varying any of these parameters results in significant differences in pressure, velocity and field temperature, enhancing or diminishing fluid pathways along the fault.

The basal heat flow used in our models provides a first estimation of the heat flow range for the geothermal region in northern Madagascar, which is comparable to the known heat flows of the proposed analogue geothermal systems for Ambilobe and Ambanja discussed below (Rajaobelison et al. 2020). In zones of extension and continental rifting, regional heat flow estimates are high, generally from 75 to 125 mWm -2 (Jaupart et al. (2015). The range of heat flow for the Great Basin (USA) is from 50 to ~120 mWm -2 , averaging roughly 90 mWm -2 [START_REF] Blackwell | Geothermal Map of North America[END_REF]. In Western Turkey, particularly in the Gediz graben, the heat flow varies from ~50 to 140 mWm -2 , with values exceeding ~100 mWm -2 (Yilmazer et al. 2010). In the present study, the basal heat flow affects the temperature distribution at depth and thus influences the geothermal gradient along permeable faults (Figure 9,8 and 10,[START_REF] Faulds | Favorable tectonic and structural settings of geothermal systems in the Great Basin region, western USA: Proxies for discovering blind geothermal systems[END_REF]. Moreover, a higher basal heat flow amplifies the convective regime when the fault permeability also increases (López and Smith 1995). In turn, the reservoir temperature is reached at shallower depths. However, the thermal spring distribution does not appear to be a function of heat flow. Indeed, high heat flow is associated with higher temperatures near the upper-ends of the fault, while hot spring temperatures can be found near the upper-ends of the fault in any of the favorable fault scenarios regardless of heat flow (Figures 2).

Fault aperture impacts heat transfer along the fault, which is amplified when fault permeability is high (kf > 10 -14 m 2 ), allowing for fluid ascension as shown in scenarios 1 and 2 for the Ambilobe 2D model of Figure 14a. Indeed, a small fault aperture (Df > 250 μm) combined with a moderate to high fault permeability (10 -10 -10 -7 m 2 ) allows hot fluids to be found at shallow depths, ˂ 2 km from the surface. Increasing fault permeability shifts the system from petrothermal to transitional, thereby making this parameter a critical factor in numerical models. By varying fault permeability in the model, the transmissivity to flow along faults becomes significant when fault permeability is higher than 10 -14 m 2 , attracting hot upwelling fluids at depth (Figure 14b and15b). An optimum combination of country-rock (host rock) and fault permeabilities is also needed to obtain subsurface temperature distributions consistent with observed thermal springs (Forster and Smith 1989). This means, the permeability difference between the country rock and the fault can increase or decrease springs temperatures [START_REF] Ferguson | What do aqueous geothermometers really tell us?[END_REF]Forster and Smith 1989;López and Smith 1995). In the present study, a permeability difference of 4 to 7 orders of magnitude between the country rock and the fault depresses the thermal field. The decrease is due to the downward flux of recharge water, while hightemperature discharge associated with upward flow occurs along the fault plane. In fact, the pressure that forces the water to the surface is mainly supplied by the difference in elevation between the recharge and discharge areas.

However, it is well known that fractures and faults can seal and re-open over time through multiple tectonic constraints (e.g. Ingebritsen and Manning 2010). The current seismicity (magnitude Mw > 5 in the study areas; (Bertil and Regnoult 1998) suggests that faults may maintain efficient permeability over time.

Discussion

Groundwater flow boundary conditions

The hydraulic head applied to the lateral boundaries of the model reproduces topographydriven groundwater flow from NNE to SSE in the Ambilobe model and from NW to SE in the Ambanja model. Surface recharge is neglected since the goal is to evaluate the impact of faults intercepting the regional flow path. At the river segment, the flow vectors are directed toward the surface and the river constitutes a discharge zone (Figures 8 and10).

The magnitude of ascending flow vectors increases from 10 -10 to 10 -2 m s -1 when faults intercept ascending flow vectors. The impact of the regional flow path was previously highlighted by Taillefer et al. (2018), who suggested that topography (i.e., an elevated water table) is the key factor controlling fluid pathways and the hot spring distribution in the Têt hydrothermal system of the Eastern Pyrénées in France. Similar conclusions can be made for Ambilobe and Ambanja, where the topography-driven groundwater flow affects forced convective heat transfer along the faults. The flow lines within the models show that the fluids originating from higher elevations reach considerable depth where they can be heated and redirected towards the surface through permeable faults (Figures 8 and 10). This also confirms that hydraulic heads and hydrological budgets control the regional groundwater flow systems in an extensional domain geothermal play, as mentioned by [START_REF] Moeck | Catalog of geothermal play types based on geologic controls[END_REF].

Model limitations

As a first step, a simplified model is necessary to understand the flow and heat transfer mechanisms influencing the temperature distribution in a geothermal system where crustal-scale faults control the fluid migration pathway. Fault and host rock permeabilities are assumed constant with depth, although they can decrease with increasing confining pressure (Rajaobelison et al. 2021). With this assumption, forced convection versus conduction is clearly favored. Indeed, permeability probably decreases with depth in the basement due to the contrasting mechanical properties of the minerals present in the rocks, such as mica and feldspar, and in faults due to sealing by the mineral precipitation (Stober and Bucher 2015).

The scarcity of structural information in the study area is an important limitation to the models. More structural field data, including dip angle and fracture density measurements, would better constrain the fault dip, aperture and permeability values in the numerical models. Moreover, statistical features summarized from field data are needed for the stochastic processes in discrete fracture models (e.g. Ivanova et al. 2014). It is currently difficult to anticipate fracture behavior at depth. Ggeophysical surveys (e.g. gravity and magnetotellurics) would help to image deep structures, build more accurate 2D models, and eventually develop 3D models around faults.

Even though these aspects limited our modeling approach, the present work successfully reproduces hot spring temperatures observed at surface, with consistent reservoir temperatures at reasonable depths.

Geothermal exploration

This study identified favorable fault characteristics that trigger fluid ascension and favor the formation of geothermal energy sources, knowledge that can assist geothermal exploration efforts to better identify exploration targets. More specifically, the target should be normal faults with ˂ 90° dips facing the regional fluid flow path driven by topography. Fault permeability must be higher than the host rock, corresponding to the critical range of permeability under which forced convection occurs, combined with a minimum aperture of 250 μm.

The simulation results also revealed that an average drilling depth of 1 to 2 km along favorable faults should reach the temperature interval of 150-200 °C (Figures 8 and10).

This finding can positively impact an exploration project since drilling cost increases nonlinearly with depth (e.g. Augustine et al. 2006 ).

Conclusions

This study presents, for the first time, 2D numerical models simulating groundwater flow and heat transfer according to a realistic topography and fault hypothesis in Madagascar.

The aim is to better understand heat transfer mechanisms involved with the geothermal systems of the Ambilobe normal fault zone and Ambanja graben in northern Madagascar.

Introduction

Economic activity in the northern region of Madagascar relies on export products such as vanilla, cocoa and coffee that are grown and transformed locally. Madagascar is one of the largest vanilla exporters in the world (UNComtrade, 2018). The energy used to transform these products relies heavily on fossil fuels and wood. Off-grid communities, mainly in remote areas, use traditional biomass such as firewood and charcoal as their main source of energy, inducing deforestation and water resource degradation, along with biodiversity and soil loss, as has already happened in many Madagascar communities. Renewable sources, such as geothermal energy, are expected to play a more prominent role in future energy production.

Ambilobe and Ambanja, two areas in northern Madagascar, have been characterized for their geothermal potential and respectively classified as a graben border-fault liquiddominated moderate-temperature system and a fossil magmatic liquid-dominated moderate-temperature system (Rajaobelison et al. 2020). In this classification scheme based on the extensional type, the fault and lithologic-diagenetic characteristics affecting porosity are thought to control formation permeability. As reported by Rajaobelison et al.

(2020), the average matrix permeability of fracture-free rocks in these two areas was evaluated between 10 -13 and 10 -18 m 2 , which is low, while the average porosity is ~5.7 %, which is low or poor (5-10 %) according to Levorsen and Berry (1967); Tiab and Donaldson (2016). However, this matrix permeability has to be considered the lower bound of formation permeability as fractures are common in such an extensional tectonic setting.

The occurrence of thermal springs along a major fault zone in Ambilobe and a graben structure in Ambanja confirms that advective heat transfer can take place within fractured rocks (Rajaobelison et al. 2020). The equilibrium formation temperature at depth, estimated using chemical geothermometers (Gunnlaugsson et al. 1981) Furthermore, the prediction and modeling of coupled fluid flow and heat transfer processes in naturally fractured rock systems are critical components in energy recovery analysis (Hao et al. 2012). In addition, numerical modeling of hydrothermal systems has often been used to explain hot spring occurrences and explore their controlling factors (e.g. This study aims to understand the role of faults on forced convective heat transfer that gave rise to hot springs in northern Madagascar, thereby providing fundamental knowledge that can be used to improve geothermal exploration concepts. Our working hypothesis is that fault dip at depth, which is currently unknown, can affect the rise of hot fluids. We developed numerical groundwater flow and heat transfer models to reproduce the nearsurface temperatures of hot springs since this is the only information available in the area to validate simulations. We present innovative concepts supported by 2D numerical simulations of deep subsurface temperature in northern Madagascar that test seven basecase fault scenarios or hypotheses. We identified fault characteristics that favored the formation of geothermal energy sources to help advance geothermal exploration.

Geological setting and geothermal context 2.1 Tectonic framework and regional fault system

The Bemarivo domain, one of the six geodynamic domains forming the Madagascar with dextral or sinistral movements can be seen on the geological map (Figure 2b) and in the field (Figure 3a). The strike-slip faults include the faults near the Betsiaka hot springs and undifferentiated normal faults of Ranomafana and Ankatakata (Figure 2b). The main fault and fracture systems are dominated by NNW-SSE, NE-SW and N-S orientations.

Fracture orientations from field measurements and fault orientations from geological map are shown on the rose diagram of Figure 2b.

The Ambanja area is located in the graben structure in the Proterozoic basement units of the southern terrane of the Bemarivo domain. The graben is oriented NW-SE, and its average width is 4 km. The main faults and fractures dominantly strike NNW-SSE. The rose diagrams of fracture orientations from field measurements as well as the faults orientations from the geological map are shown on Figure 2c. Some metric fractures were observed at the outcrop scale (Figure 3b)

Favorable structural patterns

In a comprehensive inventory of the structural settings of geothermal systems in the Great Basin region of western USA, [START_REF] Faulds | Favorable tectonic and structural settings of geothermal systems in the Great Basin region, western USA: Proxies for discovering blind geothermal systems[END_REF] identified the most favorable dominant faulting patterns to host geothermal activity in the extensional to transtensional terranes, and more specifically in areas where the bulk of the geothermal resources are likely hidden. Applying the same principles to the Ambilobe and Ambanja areas, the identified favorable fault structures are: 1) major normal fault segments, i.e., near the displacement maxima of the NE-SW Andavakoera fault and the NW-SE fault-borders of the Ambanja graben structure intersections between two normal faults (Figure 2b and2c),

and 2) or between normal faults and transverse oblique-slip faults (Figure 2c). At the surface, hot springs reach a temperature varying from 50 to 72 °C in both areas [START_REF] Besairie | Contribution à l'étude des sources minérales et des eaux souterraines de Madagascar[END_REF]Gunnlaugsson et al. 1981). Waters from the hot springs have high Mg concentrations and fall into the immature field on the ternary Na-K-Mg diagram of Giggenbach (1988), indicating dilution of the geothermal fluid along its path to the surface (Figure 4). The geothermal fluid thus mixed with meteoric and shallow groundwater.

Using chemical geothermometers (Gunnlaugsson et al. 1981) and the Giggenbah diagram (Figure 4), the reservoir temperatures in the Ambilobe and Ambanja areas are expected to be 140 and 200 °C, respectively.

Figure 4:

Chemical composition of waters from the Ambilobe (AMBL) and Ambanja (AMBJ) hot springs plotted on a ternary Na-K-Mg diagram (Giggenbach 1988). The Na, K, Mg concentrations were reported in Gunnlaugsson et al. (1981).

Lithosphere structure and geological formations

In the Bemarivo domain, the Moho is at a depth of 33 km (Andriampenomanana et al. 2017). The thickness of the upper crust is up to 20.5 km, and its composition is assumed to be primarily paragneiss and sedimentary rocks which have an average thickness of 4 km.

The lower crust is assumed to be composed of mafic rocks with an average thickness of 12.5 km.

The geological formations of the Ambilobe and Ambanja areas and their thermo-hydraulic properties have been described in Rajaobelison et al. (2020);and Rajaobelison et al. (2021).

The Permo-Triassic sedimentary rocks of the Ambilobe Basin, overlying the Proterozoic basement of the Bemarivo domain, are the main formations in the Ambilobe area. The Proterozoic basement rocks of the Bemarivo domain are the only geological formations exposed in the Ambanja area (Figure 2).

Development numerical heat transfer model

Conceptual model and assumptions

The geometry of the models (Figure 5) accounts for the topography on NNW-SSE (AA')

and NE-SW (BB') cross-sections (Figures 2b and2c) that pass through, respectively, the Andavakoera normal fault zone at the contact between sedimentary sequences and basement rocks in the Ambilobe area, and the graben structure in the Ambanja area. The dominant sedimentary and basement units have been considered in the models. The crosssections are perpendicular to the main undifferentiated faults observed near hot spring occurrences. The faults represented in the models are N010° Ankatakata post-Triassic normal oblique-slip faults in the Ambilobe area and N140° border faults in the Ambanja graben (Figures 2b and2c). These faults appear on the geological map of Madagascar produced by Roig et al. (2012).

Different modeling scenarios were defined to investigate the influence of fault dip on groundwater flow and heat transfer processes (Table 1). The dips in the two models range from 60° to 85° since dips from 60° to vertical have been shown to cause variations in the pressure and temperature of subsurface fluids [START_REF] Cherubini | Impact of single inclined faults on the fluid flow and heat transport: results from 3-D finite element simulations[END_REF]). These modeling scenarios were based on the geological maps of Ambanja and Ambilobe, although many uncertainties remain about fault dip and down-dip continuity at depth, which justifies the use of various scenarios. The models simulate the potential temperature distribution for the Ambilobe and Ambanja geothermal areas according to fault dip influence.

Governing equations

The models were developed using the subsurface flow module with a finite element method in the COMSOL Multiphysics software. The module couples groundwater flow simulated with Darcy's law and heat transfer in porous media. This made it possible to simulate forced convective heat transfer, which occurs along major faults where the host rock is considered an equivalent porous medium. Steady-state flow and heat transfer were numerically resolved according to the partial differential equation 1 and the input parameters summarized in Table 2. Sedimentary rocks and granitic and paragneiss basement rocks were considered low-permeability 2D units using the equivalent porous medium, where fluids can circulate. Faults are discrete 1D linear elements superposed on the 2D porous medium, where flow and heat transfer are solved simultaneously. Coupling between heat transfer and fluid motion occurs through the velocity field.

𝜆 b ( 𝜕 2 𝑇 𝜕𝑥 2 + 𝜕 2 𝑇 𝜕𝑧 2 ) -𝜌 f c f ( 𝑢 𝑥 𝜕𝑇 𝜕𝑥 + 𝑢 𝑧 𝜕𝑇 𝜕𝑧 ) + 𝐴 = 0 (1)
where λb (Wm -1 K -1 ) is the bulk thermal conductivity; T (°C) is absolute temperature; 𝑢 (m s -1 ) is the velocity vector determined from simulations of fluid flow quantified with Darcy's law, ρf (kg m -3 ) is the density of the fluid, cf (J kg -1 K -1 ) is the heat capacity of the fluid, and A (10 -6 W m -3 ) is the radiogenic heat production.

Geometry and subsurface properties

Geometry

 Ambilobe

The model has a width of 10 km and a depth of 10 km (Figure 5a). The sedimentary units dip 20-30° to the NE. Table 2 details the material properties of the 5 lithological domains.

The reverse fault F marks the contact between the mica schist and orthogneiss formations, and F1 and F2 are the mains faults represented in the models.

 Ambanja

The model has a rectangular geometry 20 km wide and 10 km deep (Figure 5b). Table 2 details the material properties of the 3 lithological domains. The properties of the border faults of the graben, represented by F1 and F2, are also included in the table.

Thermohydraulic properties and internal heat generation

The thermohydraulic properties of the model units were based on average properties measured in the laboratory for each geological formation reported by Rajaobelison et al. (2021), as summarized in Table 2. They are assumed to be constant with depth. Porosity ∅ and permeability k (m 2 ) were evaluated with transient measurements using the combined AP-608 gas permeameter-porosimeter from Core Test Systems with nitrogen as the gas source. The average bulk density ρ (kg m -3 ) was obtained from porosity assessments. Transient thermal conductivity λ (Wm -1 K -1 ) and diffusivity α (m 2 s -1 ) were measured with the optical scanning method using an infrared scanner developed by LGM Lippmann. The heat capacity, c (J kg -1 K -1 ), was deduced from the thermal conductivity and thermal diffusivity assessment (Popov et al. 2016;Popov et al. 1999).

The internal radiogenic heat generation (A; 10 -6 W m -3 ) for each rock presented in Table 2 was estimated from the concentration of the radiogenic elements U, Th and K, as reported in Rajaobelison et al. (2021).

Fault properties

This study assumed fault aperture and permeability from 250 µm to 0.2 cm and from 10 -14 to 10 -7 m 2 , respectively (Table 2). These ranges are based on the literature, taken from various in-situ hydraulic tests performed at geothermal sites (e.g. Holzbecher et al. 2011;Jeong et al. 2001) and examples from numerical simulation studies of groundwater flow in discrete fracture models (e.g. [START_REF] Bisdom | The impact of in-situ stress and outcropbased fracture geometry on hydraulic aperture and upscaled permeability in fractured reservoirs[END_REF]Guillou-Frottier et al. 2013). The minimum fault aperture range corresponds to that observed under a normal stress of 20 Mpa in a laboratory experiment and which validated the cubic law for laminar flow of fluids through open fractures (Witherspoon et al. 1980)

Mesh

Refinement of the mesh increases from bottom to top and close to the faults to avoid convergence problems. These locations require more accuracy due to large flow and temperature contrasts (Figure 6). A non-structured mesh with triangular elements was used, and the mesh statistics are summarized in Table 3.

Figure 6

Finite element mesh of the 2D numerical models for a) the Andavakoera normal fault zone at the basement-sediment contact in the Ambilobe area, and b) the Ambanja graben. A mesh convergence study was undertaken to verify that the temperature solution is independent of the element size; i.e., a chosen cutline was sampled uniformly from the models at different mesh sizes. The parameter σ was then computed according to equation 2:

𝜎 = √ 1 𝑁 ∑ (∆𝑇 𝑖 ) 2 𝑁 𝑖=1 (2)
where N is the sample number, ∆Ti is the temperature variation at depth xi=(i-1)D/N-1 along the chosen cutline where the mesh size of the model is refined. Here, D is the depth of the model. It demonstrated that a stable solution was obtained above 10 4 elements with a difference of about 1 °C or less (Figure 7).

Figure 7

Mesh convergence study for the a) Ambilobe and b) Ambanja models.

Boundary conditions

Groundwater flow conditions

Topography-driven groundwater flow affecting forced convection was reproduced by imposing lateral flow boundary conditions (López and Smith 1995;Smith and Chapman 1983;Tóth 2009;Woodbury and Smith 1985). Without any information about the water table level, the lateral boundaries of the model are assigned a constant hydraulic head equal to the topographic elevation (Figure 5). Each segment representing a river in the model valley is also assigned a constant hydraulic head value equal to its nodal elevation. An outflow boundary was assigned at the upper-end of the faults to allow fluids to exit at the upper-end of the faults. The base of the models was assumed impermeable.

Heat transfer conditions

A constant temperature was specified at the surface of the models and assumed to be a linear function of the altitude. Equation 3 establishes the relationship between the surface temperature (Ts) and altitude (z):

Ts= -0.0046𝑧 + 26.85 (3)
where -0.0046 and 26.85 °C are constants determined from average annual temperatures measured at four meteorological stations in Madagascar: Antsiranana, Mahajanga, Antananarivo and Tolagnaro (Table 4; Stalenberg et al. 2018). The lower boundary of the models, at a depth of 10 km, was assigned a second type condition of constant heat flow (Table 5), defined according to the Curie point depth (CPD; Tanaka et al. 1999) and the internal heat generation in the crust. Heat flow at the CPD (QCDP, mWm 2 ) was initially calculated using equation 4, which takes into account the surface (T0 =20 °C) and Curie (TCDP=580 °C) temperatures:

𝑄 𝑐𝑝𝑑 = 𝜆 𝑇 CPD-𝑇 0 𝑍 CPD -𝐴 𝑍 CPD 2 (4)
where A (10 -6 W m -3 ) is the average heat generation due to the radioactive decay of elements in the crust, which is the weighted average of the lower and upper crust heat generation. The internal heat generation in the upper and lower crusts was estimated at 1.5 × 10 -6 W m -3 and 0.72 × 10 -6 W m -3 , respectively (Table 2; Horton et al. 2016).

λ (W m -1 K -1 ) is the average thermal conductivity of the crust and ZCPD (m) is the CPD.

The range of the CPD was estimated at 10 km to 33 km, corresponding to the average Moho depth in northern Madagascar, using the centroid method of Tanaka et al. (1999) and according to the available world map of the CPD (Li et al. 2017).

The heat flow at 10 km (Q10, mWm -2 ) was then calculated with:

𝑄 10 = 𝑄 cpd + 𝐴 LC 𝐻 LC + 𝐴 UC (𝐻 UC -𝐻 model ) (5)
where HLC and HUC are the thickness of the lower and upper crusts, respectively, Hmodel is the thickness of the model, and ALC and AUC are the heat generation of the lower and upper crusts, respectively (Table 2).

Adiabatic heat transfer conditions were applied to the lateral boundaries.

Parameter sensitivity

The sensitivity of the model to the uncertain parameters was evaluated using the parametric sweep function implemented in COMSOL Multiphysics. The value of the parameters was set from the minimum to maximum. These parameters are basal heat flow, fault aperture and fault permeability (Tables 2 and5). Combining these three parameters with the different fault dip scenarios (Table 1) yielded 96 simulations for the Ambilobe 2D model and 81 for Ambanja.

2D subsurface temperature distribution and flow velocity

The different scenarios shown in Figures 8 and 10 are for base conditions using minimum, median and maximum heat flow values and average fault properties (Tables 1 and2). The simulation results confirm that the highest temperatures occur when heat flow is maximum, and that fluids and heat ascend along faults compared the fault-free base case (scenario 0).

Temperature profiles as a function of depth along fault lines (Figures 9 and 11 The magnitude of the Darcy velocity varies between 10 -9 and 10 -2 m s -1 along the faults (e.g. Figure 12). The simulated outflow velocity at the upper-ends of favorable faults (Table 6 and Figure 12) shows no difference or a 1-2 magnitude difference compared to the measured flow velocity at some hot springs, which varies between 10 -9 and 10 -7 m s -1 (Table 6; Gunnlaugsson et al. 1981). Simulated outflow velocity at upperends of fault (m s -1 )

6.47  10 -7 to 1.21  10 -6

7.07  10 -7 to 1.11  10 -6 9.63 10 -7 and 1.24  10 -6 to 2. 6 10 -7 2.14 10 -7 and 2.72  10 -7 to 1.24  10 -6

2.32  10 -8 to 1.19  10 -3 9.37  10 -5 to 9.74  10 -4

Figure 13

Variation of simulated flow velocity and temperature as a function of depth along the fault. Ambilobe 2D model, scenario 3, fault F2 and median heat flow value.

Influence of faults

Faults play a fundamental role in the movement of fluids and pore pressure distribution and therefore exert considerable influence over the temperature distribution. This is evident when comparing scenario 0 (without faults) to faulted scenarios (Figures 8 and10). The faults serve as conduits for fluid upflow, with Darcy velocity vectors converging upward along the faults toward the hot springs location. The temperature along faults can be up to 150 °C higher than adjacent areas (Figure 8 and 10). Heat transfer becomes conductive in the surrounding host rock, where bulk permeability is lower than the high-permeability zones because fluids drain towards the fault where convection occurs. Faults create permeability heterogeneities in the crust, which affect fluid flow and convective heat transfer, thereby influencing the location of hot springs [START_REF] Faulds | Characterizing structural controls of geothermal reservoirs in the Great Basin, USA, and Western Turkey: Developing successful exploration strategies in extended terranes[END_REF][START_REF] Faulds | Favorable tectonic and structural settings of geothermal systems in the Great Basin region, western USA: Proxies for discovering blind geothermal systems[END_REF], especially at their intersections (Belgrano et al. 2016;[START_REF] Person | Analytical and numerical models of hydrothermal fluid flow at fault intersections[END_REF]. This is consistent with the numerical modeling of deep fluid flow associated with crustal-scale normal faults in the south-central Canadian Cordillera [START_REF] Ferguson | What do aqueous geothermometers really tell us?[END_REF]). The latter reported that the permeability contrasts between host-rocks and faults acting as conduits for preferential fluid flow, have the greatest influence on the subsurface temperature distributions consistent with observed thermal springs. However, the presence of a fault does not necessarily imply the formation of a hot spring unless the favorable conditions are present that allow fluids and heat to ascend and reach the surface (e.g. López and Smith 1995;Thiébaud 2008) Simulation scenario 3 for the Ambilobe 2D model clearly illustrates favorable conditions for fluid ascension, with higher simulated temperatures near the upper-end of both faults (Figure 8; Table 1) and fault dips of < 90° facing the fluid flow. As shown in Figure 13, the fluid at the surface, reaches the highest Darcy velocity (10 -6 m s -1 ) which match the temperature of the hot spring (50 °C). A flow velocity of 10 -8 m s -1 is found at the depth having the reservoir temperature of 180 °C.

On the other hand, if the fault is facing away from the fluid flow direction, like F1 in scenarios 1 and 3 for the Ambanja graben model (Figure 10; Table 1), then the hydraulic gradient between both ends of the fault is smaller and the simulated temperature near the upper-end of the faults is lower. Similar results were reported by [START_REF] Moreno | Fault characterization and heat-transfer modeling to the Northwest of Nevado del Ruiz Volcano[END_REF], who simulated the impact of faults on the circulation of fluids originating from the basement underlying the Nevado del Ruiz volcanic complex of Colombia, where hot springs are present. However, scenarios 1 and 2 for the Ambilobe 2D model showed that fluids are able to rise to the surface along a fault with a dip facing away from the fluid flow direction (F2 and F1, respectively) if it is close to a river that acts as a discharge zone (Figure 8 and9). This is also seen in scenario 4. It suggests that both the fault dip and its distance from a discharge zone where the hydraulic head can be low are important characteristics in favor of fluid ascension.

Additionally, convective heat transfer along faults implies that ascending fluids are diluted as they mix with groundwater of shallower origin, to the degree that half of the ascending fluids will reach the upper 75 % of the fault length [START_REF] Ferguson | What do aqueous geothermometers really tell us?[END_REF]. This can explain the finding of significant mixing of geothermal fluids with shallower groundwater (Giggenbah diagram of Figure 5) and the absence of boiling hot springs in the Ambilobe and Ambanja geothermal fields due to cooling along the fault between the reservoir depth and the surface.

Parameter sensitivity

For each scenario, we investigated the variation in depth at which reservoir temperature is attained as a function of basal heat flow, fault aperture and permeability (Figures 15 and16). The results reveal that the depth of reservoir temperature decreases with increasing heat flow, fault aperture and fault permeability. For example, in scenario 1 (Figure 14a), with the minimum fault aperture (Df) value and maximum fault permeability (kf), the depth of reservoir temperature decreases from >2500 to 200 m along fault F2 when heat flow increases from the minimum to the maximum. For the same scenario, using maximum heat flow and mean fault permeability, the depth of reservoir temperature decreases from 1400 to 1100 m along fault F2 by increasing the value of fault aperture. In scenario 3 (Figure 14b), using minimum fault permeability and maximum fault aperture, the depth of reservoir temperature decreases from >2500 to 750 m along fault F2 as function of the heat flow value. For the same scenario, using maximum heat flow and mean fault permeability, the depth of reservoir temperature decreases from 800 to 400 m along fault F2 by increasing the fault aperture. The least favorable cases in scenarios 1 and 2 are when minimum heat flow fails to produce any points on the plots because the reservoir temperature is always deeper than 2500 m despite increases in fault aperture and permeability. The parameter sensitivity trend described for different fault dip scenarios in the Ambilobe model is the same for the Ambanja model (Figure 15). The least favorable cases correspond to scenarios 1 and 2 using the minimum heat flow value because the reservoir temperature is deeper than 2500 m despite increasing fault aperture and permeability.

We assessed the influence of basal heat flow, fault aperture and fault permeability given the uncertainty associated with those inputs. Varying any of these parameters results in significant differences in pressure, velocity and field temperature, enhancing or diminishing fluid pathways along the fault.

The basal heat flow used in our models provides a first estimation of the heat flow range for the geothermal region in northern Madagascar, which is comparable to the known heat flows of the proposed analogue geothermal systems for Ambilobe and Ambanja discussed below (Rajaobelison et al. 2020). In zones of extension and continental rifting, regional heat flow estimates are high, generally from 75 to 125 mWm -2 (Jaupart et al. (2015). The range of heat flow for the Great Basin (USA) is from 50 to ~120 mWm -2 , averaging roughly 90 mWm -2 [START_REF] Blackwell | Geothermal Map of North America[END_REF]. In Western Turkey, particularly in the Gediz graben, the heat flow varies from ~50 to 140 mWm -2 , with values exceeding ~100 mWm -2 (Yilmazer et al. 2010). In the present study, the basal heat flow affects the temperature distribution at depth and thus influences the geothermal gradient along permeable faults (Figure 9,8 and 10,[START_REF] Faulds | Favorable tectonic and structural settings of geothermal systems in the Great Basin region, western USA: Proxies for discovering blind geothermal systems[END_REF]. Moreover, a higher basal heat flow amplifies the convective regime when the fault permeability also increases (López and Smith 1995). In turn, the reservoir temperature is reached at shallower depths. However, the thermal spring distribution does not appear to be a function of heat flow. Indeed, high heat flow is associated with higher temperatures near the upper-ends of the fault, while hot spring temperatures can be found near the upper-ends of the fault in any of the favorable fault scenarios regardless of heat flow (Figures 2).

Fault aperture impacts heat transfer along the fault, which is amplified when fault permeability is high (kf > 10 -14 m 2 ), allowing for fluid ascension as shown in scenarios 1 and 2 for the Ambilobe 2D model of Figure 14a. Indeed, a small fault aperture (Df > 250 μm) combined with a moderate to high fault permeability (10 -10 -10 -7 m 2 ) allows hot fluids to be found at shallow depths, ˂ 2 km from the surface. Increasing fault permeability shifts the system from petrothermal to transitional, thereby making this parameter a critical factor in numerical models. By varying fault permeability in the model, the transmissivity to flow along faults becomes significant when fault permeability is higher than 10 -14 m 2 , attracting hot upwelling fluids at depth (Figure 14b and15b). An optimum combination of country-rock (host rock) and fault permeabilities is also needed to obtain subsurface temperature distributions consistent with observed thermal springs (Forster and Smith 1989). This means, the permeability difference between the country rock and the fault can increase or decrease springs temperatures [START_REF] Ferguson | What do aqueous geothermometers really tell us?[END_REF]Forster and Smith 1989;López and Smith 1995). In the present study, a permeability difference of 4 to 7 orders of magnitude between the country rock and the fault depresses the thermal field. The decrease is due to the downward flux of recharge water, while hightemperature discharge associated with upward flow occurs along the fault plane. In fact, the pressure that forces the water to the surface is mainly supplied by the difference in elevation between the recharge and discharge areas.

However, it is well known that fractures and faults can seal and re-open over time through multiple tectonic constraints (e.g. Ingebritsen and Manning 2010). The current seismicity (magnitude Mw > 5 in the study areas; (Bertil and Regnoult 1998) suggests that faults may maintain efficient permeability over time.

Discussion

Groundwater flow boundary conditions

The hydraulic head applied to the lateral boundaries of the model reproduces topographydriven groundwater flow from NNE to SSE in the Ambilobe model and from NW to SE in the Ambanja model. Surface recharge is neglected since the goal is to evaluate the impact of faults intercepting the regional flow path. At the river segment, the flow vectors are directed toward the surface and the river constitutes a discharge zone (Figures 8 and10).

The magnitude of ascending flow vectors increases from 10 -10 to 10 -2 m s -1 when faults intercept ascending flow vectors. The impact of the regional flow path was previously highlighted by Taillefer et al. (2018), who suggested that topography (i.e., an elevated water table) is the key factor controlling fluid pathways and the hot spring distribution in the Têt hydrothermal system of the Eastern Pyrénées in France. Similar conclusions can be made for Ambilobe and Ambanja, where the topography-driven groundwater flow affects forced convective heat transfer along the faults. The flow lines within the models show that the fluids originating from higher elevations reach considerable depth where they can be heated and redirected towards the surface through permeable faults (Figures 8 and 10). This also confirms that hydraulic heads and hydrological budgets control the regional groundwater flow systems in an extensional domain geothermal play, as mentioned by [START_REF] Moeck | Catalog of geothermal play types based on geologic controls[END_REF].

Model limitations

As a first step, a simplified model is necessary to understand the flow and heat transfer mechanisms influencing the temperature distribution in a geothermal system where crustal-scale faults control the fluid migration pathway. Fault and host rock permeabilities are assumed constant with depth, although they can decrease with increasing confining pressure (Rajaobelison et al. 2021). With this assumption, forced convection versus conduction is clearly favored. Indeed, permeability probably decreases with depth in the basement due to the contrasting mechanical properties of the minerals present in the rocks, such as mica and feldspar, and in faults due to sealing by the mineral precipitation (Stober and Bucher 2015).

The scarcity of structural information in the study area is an important limitation to the models. More structural field data, including dip angle and fracture density measurements, would better constrain the fault dip, aperture and permeability values in the numerical models. Moreover, statistical features summarized from field data are needed for the stochastic processes in discrete fracture models (e.g. Ivanova et al. 2014). It is currently difficult to anticipate fracture behavior at depth. Ggeophysical surveys (e.g. gravity and magnetotellurics) would help to image deep structures, build more accurate 2D models, and eventually develop 3D models around faults.

Even though these aspects limited our modeling approach, the present work successfully reproduces hot spring temperatures observed at surface, with consistent reservoir temperatures at reasonable depths.

Geothermal exploration

This study identified favorable fault characteristics that trigger fluid ascension and favor the formation of geothermal energy sources, knowledge that can assist geothermal exploration efforts to better identify exploration targets. More specifically, the target should be normal faults with ˂ 90° dips facing the regional fluid flow path driven by topography. Fault permeability must be higher than the host rock, corresponding to the critical range of permeability under which forced convection occurs, combined with a minimum aperture of 250 μm.

The simulation results also revealed that an average drilling depth of 1 to 2 km along favorable faults should reach the temperature interval of 150-200 °C (Figures 8 and10).

This finding can positively impact an exploration project since drilling cost increases nonlinearly with depth (e.g. Augustine et al. 2006 ).

Conclusions

This study presents, for the first time, 2D numerical models simulating groundwater flow and heat transfer according to a realistic topography and fault hypothesis in Madagascar.

The aim is to better understand heat transfer mechanisms involved with the geothermal systems of the Ambilobe normal fault zone and Ambanja graben in northern Madagascar.

Figure 1

 1 Figure 1Simplified geological maps showing a) the six main geodynamic

Figure 2

 2 Figure 2 (a) Location of the Ambilobe and Ambanja geothermal areas. (b) and (c)

Figure 3 :

 3 Figure 3: a) Field view of the Andavakoera fault line and its conjugated faults (F). b)

Figure 5 :

 5 Figure 5: Geometry and boundary conditions of the 2D numerical model based on the

  ) show which simulations best reproduce the near-surface temperature of the hot springs. The results also indicate the depth at which reservoir temperature is attained, assuming the chemical geothermometer provides a reliable reservoir temperature estimate. The most realistic scenarios are when hot spring temperatures are attained near the upper-end of the faults, and the reservoir temperature is at a reasonable depth, below approximately 2 km. The most realistic simulation cases for Ambilobe are maximum heat flow value and fluids rising along fault F2 in scenario 1, maximum heat flow value and fluids rising along fault F1 in scenario 2, and median or maximum heat flow values and fluids rising along fault F2 in scenarios 3 and 4 (Figures 8 and 9).

Figure 8

 8 Figure 8Simulated temperature distribution for the Ambilobe 2D model considering different fault dips and heat flux variations

Figure 9

 9 Figure 9Temperature profiles as a function of depth along fault lines for the Ambilobe 2D

Figure 10

 10 Figure 10Simulated temperature distribution for the Ambanja 2D model considering different fault dips and heat flux variations at

  Figure 12a) Ambilobe 2D model, scenario 3: flow velocity along fault F2, maximum

Figure 14

 14 Figure 14Depth variation of reservoir temperature along fault lines in the Ambilobe

Figure 15

 15 Figure 15Depth variation of reservoir temperature along fault lines in the Ambanja

  Numerical simulations are used to define the large-scale temperature field in such extensional-type geothermal systems. The resulting simulated temperature distribution appears to be characteristic of conductive heat transfer in a low-permeability rock matrix and forced convective heat transfer along faults, representing favorable conditions that facilitate fluid ascension. The most realistic scenarios are when the hot spring temperature is attained near the upper-end of the fault and the reservoir temperature is at a reasonable depth, shallower than ~2 km. The simulated flow velocity at the upper-ends of faults is comparable to those measured in surface hot springs.Our findings support previous interpretations of conceptual models in northernMadagascar, where conduction is the dominant heat transfer mechanism in lowpermeability rock matrix, except along faults that experience a strong hydraulic gradient and provide a path for fluids to rise with forced convective heat transfer as the dominant mechanism Numerical models used in this study can be a predictive tool to help define geothermal exploration strategies in similar normal fault zone and graben settings. Further hydrogeological assessment using stable isotope analyses on springs and groundwaters can help improve and support the simulated fluid circulation concept by defining infiltration altitudes and determining fluid origin. Deep borehole data, such as temperature logs, will also be required for the next steps to realistically constrain the temperature prediction and advance geothermal exploration in North Madagascar.

Figure 1

 1 Figure 1Simplified geological maps showing a) the six main geodynamic

Figure 2

 2 Figure 2 (a) Location of the Ambilobe and Ambanja geothermal areas. (b) and (c)

Figure 3 :

 3 Figure 3: a) Field view of the Andavakoera fault line and its conjugated faults (F). b)

Figure 5 :

 5 Figure 5: Geometry and boundary conditions of the 2D numerical model based on the

  ) show which simulations best reproduce the near-surface temperature of the hot springs. The results also indicate the depth at which reservoir temperature is attained, assuming the chemical geothermometer provides a reliable reservoir temperature estimate. The most realistic scenarios are when hot spring temperatures are attained near the upper-end of the faults, and the reservoir temperature is at a reasonable depth, below approximately 2 km. The most realistic simulation cases for Ambilobe are maximum heat flow value and fluids rising along fault F2 in scenario 1, maximum heat flow value and fluids rising along fault F1 in scenario 2, and median or maximum heat flow values and fluids rising along fault F2 in scenarios 3 and 4 (Figures 8 and 9).

Figure 8

 8 Figure 8Simulated temperature distribution for the Ambilobe 2D model considering different fault dips and heat flux variations

Figure 9

 9 Figure 9Temperature profiles as a function of depth along fault lines for the Ambilobe 2D

Figure 10

 10 Figure 10Simulated temperature distribution for the Ambanja 2D model considering different fault dips and heat flux variations at

Figure 11

 11 Figure 11Temperature profiles as a function of depth along fault lines for the Ambanja

  Figure 12a) Ambilobe 2D model, scenario 3: flow velocity along fault F2, maximum

Figure 14

 14 Figure 14Depth variation of reservoir temperature along fault lines in the Ambilobe

Figure 15

 15 Figure 15Depth variation of reservoir temperature along fault lines in the Ambanja

  Numerical simulations are used to define the large-scale temperature field in such extensional-type geothermal systems. The resulting simulated temperature distribution appears to be characteristic of conductive heat transfer in a low-permeability rock matrix and forced convective heat transfer along faults, representing favorable conditions that facilitate fluid ascension. The most realistic scenarios are when the hot spring temperature is attained near the upper-end of the fault and the reservoir temperature is at a reasonable depth, shallower than ~2 km. The simulated flow velocity at the upper-ends of faults is comparable to those measured in surface hot springs.Our findings support previous interpretations of conceptual models in northernMadagascar, where conduction is the dominant heat transfer mechanism in lowpermeability rock matrix, except along faults that experience a strong hydraulic gradient and provide a path for fluids to rise with forced convective heat transfer as the dominant mechanism Numerical models used in this study can be a predictive tool to help define geothermal exploration strategies in similar normal fault zone and graben settings. Further hydrogeological assessment using stable isotope analyses on springs and groundwaters can help improve and support the simulated fluid circulation concept by defining infiltration altitudes and determining fluid origin. Deep borehole data, such as temperature logs, will also be required for the next steps to realistically constrain the temperature prediction and advance geothermal exploration in North Madagascar.

Table 1

 1

	Fault scenarios
	Area	Faults	Dip	Scenario
			No fault	0
			70° NW; 85° SE	1
	Ambilobe	F1; F2	60° NW; 60° SE 80° SE; 75° SE	2 3
			80° NW; 85° NW	4
			No fault	0
			80° SW; 80° NE	1
	Ambanja	F1; F2	75° NE; 80° NE	2
			75° NW; 80° NW	3

=thermal conductivity, ρ= density, c= heat capacity, A= radiogenic heat production

	Table 2 Model input parameters
	Parameter				Value				Unit	Reference
	Area			Ambilobe			Ambanja
	Formation		Sediments		Basement		Basement
	Geological Units	Permian Sandstone	Eo Triassic shale	Triassic sandstone	Mica Shist Amphibolite	Ortho-gneiss	Syenite	Para-gneiss
	Mean ∅	0.07	0.08	0.07	0.11	0.05	0.001	0.02	Fractional
	Mean k	1.4  10 -13	1.0  10 -17	1.0  10 -16	1.0  10 -15	1.0  10 -16	1.0  10 -18	1.0 10 -15	m 2	Rajaobelison et al. (2021)
	Mean λ	3.6	2.5	3.5	2.6	3.4	2.6	2.7	W m -1 .K -1
	Mean ρ	2542	2382	2628	2735	2678	2700	2678	kg.m -3
	Mean c	904	875	919	879	846	875	830	J kg -1 K -1
	Mean A	1.14	0.72	1.25	0.58	0.38	1.35	0.38	10 -6 W m-3
	Thickness Model		4		6		10		km
	Depth of crust Upper Crust (UC) Lower Crust (LC)				33 20.5 12.5				km km km	Andriampenomanana et al. (2017)
	Mean A (UC) Mean A (LC)				1.5 0.72				10 -6 W m -3	Horton et al. (2016)
	A-Total Crust				1.2				10 -6 W m -3
	Fault Properties				F1; F2
	Aperture Df	Min. 250  10 -6		Mean 1.1  10 -3		Max. 2  10 -3		m	Current study
	Permeability kf		10 -14		10 -10		10 -7		m 2
	λ

Table 3

 3 Mesh statistics of Ambilobe and Ambanja 2D models

		Ambilobe	Ambanja	Unit
	Number of element	101 672	274 403
	Maximum element size	104	404	m
	Minimum element size	0.208	0.404	m
	Nodes	49 506	137 998

Table 4

 4 Data from selected meteorological stations in Madagascar

	Station	Antsiranana	Mahajanga	Antananarivo	Tolagnaro
	Altitude (m)	105	27	1268	22
	Average annual temperature (°C)	24.7	27	19.5	24

Table 5

 5 Parameters for the heat transfer boundary condition

	Parameters		Value		Unit
	Surface temperature		Ts=-0.0046z+26.85		o C
	Heat flow at base (Q10 )	Min 50	Mean 99	Max 148	mWm -2

Table 6

 6 Simulated flow velocity at favorable faults for the most realistic scenarios for Ambilobe and Ambanja models

	Ambilobe	Ambanja

  Besairie, H. Contribution à l'étude des sources minérales et des eaux souterraines de Madagascar. Travaux du Bureau Géologique, Service Géologique Antananarivo; 1959. Report No.: 92. Bisdom, K., Bertotti, G., & Nick, H. M. (2016). The impact of in-situ stress and outcropbased fracture geometry on hydraulic aperture and upscaled permeability in fractured reservoirs.Tectonophysics, 690, 63, doi:10.1016Tectonophysics, 690, 63, doi:10. /j.tecto.2016.04.006. .04.006.

	Blackwell D, D., & Richards, M. (2004). Geothermal Map of North America.
	Brogi, A., Alçiçek, M. C., Liotta, D., Capezzuoli, E., Zucchi, M., & Matera, P. F. (2021).
	Step-over fault zones controlling geothermal fluid-flow and travertine formation
	(Denizli		Basin,	Turkey).	Geothermics,	89,	101941,
	doi:https://doi.org/10.1016/j.geothermics.2020.101941.
	Chagnon-Lessard, N., Mathieu-Potvin, F., & Gosselin, L. (2020). Optimal design of
	geothermal power plants: A comparison of single-pressure and dual-pressure
	organic		Rankine	cycles.	Geothermics,	86,	101787,
			Environmental	Earth	Sciences,	70(8),	3603-3618,
	Faulds, J., Coolbaugh, M., Bouchot, V., Moeck, I., & Oguz, K. (2010). Characterizing
	structural controls of geothermal reservoirs in the Great Basin, USA, and Western Petroleum Geology, 78, 750-758, Turkey: Developing successful exploration strategies in extended terranes. In World Geothermal Congress Bali, 25-30 April 2010 pp. 11 doi:https://doi.org/10.1016/j.marpetgeo.2016.10.002.
	Moeck, I. S. (2014). Catalog of geothermal play types based on geologic controls. Faulds, J. E., & Hinz, N. H. Favorable tectonic and structural settings of geothermal Renewable and Sustainable Energy Reviews, 37, 867-882, systems in the Great Basin region, western USA: Proxies for discovering blind geothermal systems. In World Geothermal Congress, Melbourne, 19-25 2015 pp. doi:http://dx.doi.org/10.1016/j.rser.2014.05.032.
	6, Moreno, D., Lopez-Sanchez, J., Blessent, D., & Raymond, J. (2018). Fault characterization
	and heat-transfer modeling to the Northwest of Nevado del Ruiz Volcano. Journal
	of	South	American	Earth	Sciences,	88,	50-63,

doi:https://doi.org/10.1016/j.geothermics.2019.101787. Cherubini, Y., Cacace, M., Blöcher, G., & Scheck-Wenderoth, M. (2013). Impact of single inclined faults on the fluid flow and heat transport: results from 3-D finite element simulations. doi:10.1007/s12665-012-2212-z. Chorowicz, J., Bardintzeff, J. M., Rasamimanana, G., Chotin, P., Thouin, C., & Rudant, J. P. (1997). An approach using SAR ERS images to relate extension fractures to volcanic vents: examples from iceland and Madagascar. Tectonophysics, 271(3), 263-283, doi:http://dx.doi.org/10.1016/S0040-1951(96)00250-8. Corbel, S., Schilling, O., Horowitz, F. G., Reid, L. B., Sheldon, H. A., Timms, N. E., et al. Identification and Geothermal Influence of Faults in the Perth Metropolitan Area, Australia. In Workshop on Geothermal Reservoir Engineering, Stanford University, 30 January-1 February 2012 Vol. 37, pp. 8: SGP-TR-194 De Waele, B., Thomas, R. J., Macey, P. H., Horstwood, M. S. A., Tucker, R. D., Pitfield, P. E. J., et al. (2011). Provenance and tectonic significance of the Palaeoproterozoic metasedimentary successions of central and northern Madagascar. Precambrian Research, 189(1), 18-42, doi:https://doi.org/10.1016/j.precamres.2011.04.004. Ferguson, G., Grasby, S. E., & Hindle, S. R. (2009). What do aqueous geothermometers really tell us? Geofluids, 9(1), 39-48, doi:https://doi.org/10.1111/j.1468-8123.2008.00237.x. doi:https://doi.org/10.1016/j.jsames.2018.08.008.

  generate electricity using an optimized Organic Rankine Cycle (e.g.[START_REF] Chagnon-Lessard | Optimal design of geothermal power plants: A comparison of single-pressure and dual-pressure organic Rankine cycles[END_REF] Liu et al. 2017; Shi and Pan 2019).

	No exploration boreholes have been drilled in northern Madagascar, so the optimal depth
	for finding geothermal sources remains unknown, as do the geothermal gradient and heat
	flux density. Previous work classified geothermal play types and evaluated rock thermo-
	hydraulic properties from surface outcrops to better understand heat transfer mechanisms
	at depth (Rajaobelison et al. 2020; Rajaobelison et al. 2021). Crustal-scale faults have been
	Numerical modeling has been previously used to simulate the depth distribution of
	temperature to better understand fractured geothermal reservoirs (e.g. Bakhsh et al. 2016;
	Corbel et al. 2012; Hao et al. 2012; Holzbrecher et al. 2010). Faults that serve as conduits
	, ranges from 140 and increase permeability (e.g. Brogi et al. 2021; Ferrill et al. 2004; Liotta et al. 2021;
	to 200 °C, high enough to be used for scalding vanilla beans (~60 °C), heating greenhouses Zucchi 2020) are important targets for geothermal resource exploration, but if they act as

(> 40 °C; Lindal 1973)

, and generating electricity (> 80 °C;

Grasby et al. 2011; Tomarov and Shipkov 2017)

. Fluids from geothermal sources with temperatures below 120 °C can be used to identified as potential targets for geothermal exploration, but further work is needed to better understand their role in the subsurface heat transfer mechanisms that control the depth and extent of geothermal sources. Fault influence must be taken into account when exploring geothermal energy sources, especially the influence of the large-scale faults identified by surface mapping and shown on the geological maps of Madagascar by

Roig et al. (2012)

and

Tucker et al. (2014)

. Therefore, numerical modeling was used to simulate groundwater flow and heat transfer according to various fault scenarios affecting the geothermal energy sources associated with petrothermal systems in this extensional domain of Madagascar. barriers and decrease permeability (e.g.

Gibson 1998)

, they can compartmentalize the reservoirs, making it challenging to access and exploit the resource

(Loveless et al. 2014

).

Table 1

 1

	Fault scenarios
	Area	Faults	Dip	Scenario
			No fault	0
			70° NW; 85° SE	1
	Ambilobe	F1; F2	60° NW; 60° SE 80° SE; 75° SE	2 3
			80° NW; 85° NW	4
			No fault	0
			80° SW; 80° NE	1
	Ambanja	F1; F2	75° NE; 80° NE	2
			75° NW; 80° NW	3

=thermal conductivity, ρ= density, c= heat capacity, A= radiogenic heat production

	Table 2 Model input parameters
	Parameter				Value				Unit	Reference
	Area			Ambilobe			Ambanja
	Formation		Sediments		Basement		Basement
	Geological Units	Permian Sandstone	Eo Triassic shale	Triassic sandstone	Mica Shist Amphibolite	Ortho-gneiss	Syenite	Para-gneiss
	Mean ∅	0.07	0.08	0.07	0.11	0.05	0.001	0.02	Fractional
	Mean k	1.4  10 -13	1.0  10 -17	1.0  10 -16	1.0  10 -15	1.0  10 -16	1.0  10 -18	1.0 10 -15	m 2	Rajaobelison et al. (2021)
	Mean λ	3.6	2.5	3.5	2.6	3.4	2.6	2.7	W m -1 .K -1
	Mean ρ	2542	2382	2628	2735	2678	2700	2678	kg.m -3
	Mean c	904	875	919	879	846	875	830	J kg -1 K -1
	Mean A	1.14	0.72	1.25	0.58	0.38	1.35	0.38	10 -6 W m-3
	Thickness Model		4		6		10		km
	Depth of crust Upper Crust (UC) Lower Crust (LC)				33 20.5 12.5				km km km	Andriampenomanana et al. (2017)
	Mean A (UC) Mean A (LC)				1.5 0.72				10 -6 W m -3	Horton et al. (2016)
	A-Total Crust				1.2				10 -6 W m -3
	Fault Properties				F1; F2
	Aperture Df	Min. 250  10 -6		Mean 1.1  10 -3		Max. 2  10 -3		m	Current study
	Permeability kf		10 -14		10 -10		10 -7		m 2
	λ

Table 3

 3 Mesh statistics of Ambilobe and Ambanja 2D models

		Ambilobe	Ambanja	Unit
	Number of element	101 672	274 403
	Maximum element size	104	404	m
	Minimum element size	0.208	0.404	m
	Nodes	49 506	137 998

Table 4

 4 Data from selected meteorological stations in Madagascar

	Station	Antsiranana	Mahajanga	Antananarivo	Tolagnaro
	Altitude (m)	105	27	1268	22
	Average annual temperature (°C)	24.7	27	19.5	24

Table 5

 5 Parameters for the heat transfer boundary condition

	Parameters		Value		Unit
	Surface temperature		Ts=-0.0046z+26.85		o C
	Heat flow at base (Q10 )	Min 50	Mean 99	Max 148	mWm -2

Table 6

 6 Simulated flow velocity at favorable faults for the most realistic scenarios for Ambilobe and Ambanja models

	Ambilobe	Ambanja

Acknowledgments

The authors would like to thank the Canadian Francophonie Scholarship Program CFSP, which supported the first author with a scholarship. The UNESCO International geoscience program (IGCP) is also acknowledged for providing a framework that supported collaboration among the members of the IGCP-636 group involved in geothermal research and who participated as coauthors to this paper. Additional thanks go to Elie R. for discussions and invaluable help with modeling issues using COMSOL.

Acknowledgments

The authors would like to thank the Canadian Francophonie Scholarship Program CFSP, which supported the first author with a scholarship. The UNESCO International geoscience program (IGCP) is also acknowledged for providing a framework that supported collaboration among the members of the IGCP-636 group involved in geothermal research and who participated as coauthors to this paper. Additional thanks go to Elie R. for discussions and invaluable help with modeling issues using COMSOL.

Availability of data and materials

The relevant datasets analyzed in this study are all presented in the manuscript.

Availability of data and materials

The relevant datasets analyzed in this study are all presented in the manuscript.

Author contributions

MR: Data collection, investigation, conceptualization, methodology, Writing (original draft preparation), formal analysis. JR and MM: Supervision, resources, project leads. JR, CD and MR: Validation. JR, MM, CD and SL: Drafting and editing. All authors reviewed the manuscript and have read and agreed to the published version of the manuscript.

Competing interests

The authors declare that they have no competing interests.

Author contributions

MR: Data collection, investigation, conceptualization, methodology, Writing (original draft preparation), formal analysis. JR and MM: Supervision, resources, project leads. JR, CD and MR: Validation. JR, MM, CD and SL: Drafting and editing. All authors reviewed the manuscript and have read and agreed to the published version of the manuscript.

Competing interests

The authors declare that they have no competing interests.