\makeatletter
\@ifundefined{HCode}
{\documentclass[screen, CRGEOS, Unicode, Thematic]{cedram}
%\usepackage[sort&compress,square,comma,authoryear]{natbib}
\usepackage[square,comma,authoryear]{natbib}
\newenvironment{Table}{\begin{table}}{\end{table}}
\renewcommand\bibsection{\section*{\refname}}
\newenvironment{noXML}{}{}
\def\thead{\noalign{\relax}\hline}
\def\endthead{\noalign{\relax}\hline}
\def\tbody{\noalign{\relax}\hline}
\def\botline{\\\hline} 
\def\tabnote#1{\vskip4pt\parbox{.97\linewidth}{#1}}
\def\tsup#1{$^{{#1}}$}
\def\tsub#1{$_{{#1}}$}
\RequirePackage{etoolbox}
\usepackage{upgreek}
\def\jobid{crgeos20221173}
%\graphicspath{{/tmp/\jobid_figs/web/}}
\graphicspath{{./figures/}}
\newcounter{runlevel}
\usepackage{multirow} 
\def\nrow#1{\@tempcnta #1\relax%
\advance\@tempcnta by 1\relax%
\xdef\lenrow{\the\@tempcnta}}
\def\morerows#1#2{\nrow{#1}\multirow{\lenrow}{*}{#2}}
\def\xmorerows#1#2{{#2}}
\let\MakeYrStrItalic\relax
\def\refinput#1{}
\def\back#1{}
\def\myRef{\ref}
\csdef{Seqnsplit}{\\}
\DOI{10.5802/crgeos.213}
\datereceived{2022-12-01}
\daterevised{2023-03-27}
\dateaccepted{2023-04-06}
\ItHasTeXPublished
\def\ndash{\text{--}}
\makeatletter
\g@addto@macro{\UrlBreaks}{\UrlOrds}
\gappto{\UrlBreaks}{\UrlOrds}
\makeatother
}
{
\PassOptionsToPackage{authoryear}{natbib}
\documentclass[crgeos]{article}
\usepackage[T1]{fontenc} 
\def\CDRdoi{10.5802/crgeos.213}
\makeatletter
\def\CDRsupplementaryTwotypes#1#2{}
 \def\selectlanguage#1{} 
 \def\xmorerows#1#2{\morerows{#1}{#2}}
}
\makeatother

\dateposted{2023-05-15}
\begin{document}

\begin{noXML}

%\TopicEF{Geodynamics of continents and oceans - A tribute to Jean
%Aubouin}{G\'eodynamique des continents et des Oc\'eans - Hommage \`a
%Jean Aubouin}

\shortrunauthors

\title{Heat flow, thermal anomalies, tectonic regimes and
high-temperature geothermal systems in fault zones}

\author{\firstname{Laurent} \lastname{Guillou-Frottier}\CDRorcid{0000-0002-2698-1447}\IsCorresp}
\address{BRGM, Georesources Division, 45060 Orl\'{e}ans, France}
\address{ISTO, UMR 7327, Universit\'{e} d'Orl\'{e}ans, CNRS, BRGM, 45071, 
Orl\'{e}ans, France}
\email[L. Guillou-Frottier]{l.guillou-frottier@brgm.fr}

\author{\firstname{Ga\'{e}tan} \lastname{Milesi}\CDRorcid{0000-0001-8658-9837}}
\address{GeoRessources, CNRS, 
Universit\'{e} de Lorraine, Labcom CREGU, 
54506 Vand\oe{}uvre-l\`{e}s-Nancy, France}
\email[G. Milesi]{gaetan.milesi@univ-lorraine.fr}

\author{\firstname{Vincent} \lastname{Roche}\CDRorcid{0000-0003-1482-3115}}
\address{ISTEP, UMR 7193, Sorbonne Universit\'{e}, 
CNRS-INSU, 75005 Paris, France}
\address{Laboratoire des Fluides Complexes et leurs R\'{e}servoirs -- IPRA, 
E2S-UPPA, TotalEnergies, CNRS, 
Universit\'{e} de Pau et des Pays de l'Adour, 
UMR 5150, 64012 Pau, France}
\email[V. Roche]{vincent.roche@sorbonne-universite.fr}

\author{\firstname{Hugo} \lastname{Duwiquet}\CDRorcid{0000-0003-0278-3005}}
\address{IFP Energies nouvelles, 1 et 4 avenue de Bois-Pr\'{e}au, 
92852 Rueil-Malmaison, France}
\email[H. Duwiquet]{hugo.duwiquet@ifpen.fr}

\author{\firstname{Audrey} \lastname{Taillefer}\CDRorcid{0000-0001-8946-6875}}
\address{Universit\'{e} Grenoble Alpes,  Universit\'{e} Savoie Mont
Blanc, CNRS, IRD,  Universit\'{e} Gustave Eiffel, ISTerre, 38000
Grenoble, France}
\email[A. Taillefer]{audrey.taillefer@gmail.com}

\keywords{\kwd{Geothermal energy}
\kwd{Heat flow}
\kwd{Permeability}
\kwd{Fault zones}
\kwd{Fault geometry}}

\begin{abstract} 
The potential of high-temperature
$({>}150~\text{\textdegree}\mathrm{C})$ geothermal systems in crustal
fault zones (fault cores and hundreds of meters wide networks of
interconnected fractures in the damage zone) is underestimated. Based
on numerical models, we show that topography-driven,
poroelasticity-driven as well as buoyancy-driven forces play a
significant role in the establishment of shallow  (1--4~km) thermal
anomalies in fault zones. We investigate the role of permeability,
topography, fault dip, tectonic regime and fault geometry on the
amplitude of thermal anomalies. Favorable conditions include:  (i)~a
damage zone thickness  $>$ 100~m,  (ii)~a minimum cumulative
displacement of 100--150~m and  (iii)~fault zone lengths of at least
one kilometer. Based on these parameters, we propose new potential
targets for the geothermal exploration of fault zones in Western
Europe.
\end{abstract}

\maketitle

\clearpage

%\vspace*{3.3pc plus 1pt minus 1pt}

\twocolumngrid

\end{noXML}

\section{Introduction}\label{sec1}
In a global energy transition context, the potential of geothermal
power, i.e. electricity generation from geothermal energy, might be
underestimated. The installed capacity of geothermal power plants has
increased by 27\% between 2015 and 2020  \citep{Huttrer2021}. The
world's total geothermal power reached 15.9 GWe (electric gigawatt) in
May 2020 and the regular increase since 2010 would imply a total world
installed capacity of around 20~GWe (the equivalent of the average
power of  20 nuclear reactors) in 2025  \citep{Huttrer2021}. However,
this prediction is based on present-day geological models of geothermal
reservoirs, which are generally assumed to be located in specific
tectonic and geodynamic contexts  (e.g. extensional tectonics or
volcanic areas). The three main countries where geothermal exploration
has recently increased are Indonesia, Kenya, and Turkey. In these and
other countries where geothermal energy is exploited, surface
manifestations such as geysers, hot springs, or fumaroles are observed
and represent direct evidence of an underlying geothermal reservoir
(e.g., Larderello, Italy; Wairakei, New Zealand; The Geysers,
California; Olkaria, Kenya; Germencik, Turkey). However, in this study,
we focus on specific geological systems that may host geothermal
reservoirs without showing surface evidence. Permeable structures such
as crustal fault zones (fault cores and hundreds of meters wide
networks of interconnected fractures in the damage zone) are
particularly discussed. These structures may represent unconventional
geothermal resources for power production  \citep{Duwiquet2022}, since
reservoir permeability would be sufficiently favorable to avoid
stimulation  [see  \citealp{Jolieetal2021} for a detailed review of
conventional and unconventional geothermal  resources].

In the following, a ``high-temperature'' geothermal system refers to a
geothermal reservoir where the temperature is greater than 
150~\textdegree C. Various definitions are used in the literature 
\citep[e.g.,][]{Jolieetal2021} and our choice corresponds to a
threshold value of temperature where exploration of such systems in
France has to be officially declared. Actually,  \citet{Lundetal2022}
emphasized that electric power can be generated from geothermal systems
where the temperature is as low as  90~\textdegree C, but our study
will be focused on higher-temperature geothermal systems with better
energy performance.

Several high-temperature geothermal systems at exploitable depths 
(1--4~km) do not exhibit surface manifestations (e.g., Landau in the
Upper Rhine graben). Actually, ``blind'' or ``hidden'' geothermal
systems were discovered accidentally during oil, gas, or mining
activities  [e.g., the geothermal reservoirs in the Imperial  Valley of
southern California, or those in the Basin and Range province; 
\citealp{FauldsHinz2015,Gargetal2010,Dobson2016}]. Likewise, the
Soultz-sous-For\^{e}ts thermal anomaly (Upper Rhine graben, France) was
first evidenced by temperature measurements performed in petroleum
boreholes  \citep{HaasHoffmann1929}. Recently, a blind geothermal
system in western Nevada has been discovered through a ``geothermal
play fairway analysis'', where gravity, magnetics, resistivity, and
geothermometry data were combined before deciding to drill
temperature-gradient boreholes  \citep{Craigetal2021}. Interestingly,
their predicted geothermal power reaches 16~MWe. From this example and
many others, one can infer that the presence of a magmatic activity is
not a prerequisite for the formation of shallow thermal anomalies. To
summarize, exploitable thermal anomalies at shallow depths are probably
much more present than previously\unskip\break 
assumed.

Promising areas where geothermal systems could be exploited can also be
inferred from the presence of anomalously high subsurface temperature
gradients. However, direct measurements of underground temperatures
depend on the availability of sufficiently deep (hundreds to thousands
of meters) and open boreholes. When such temperature measurements are
available, downward extrapolation to greater depths may however be
biased by several physical and geological processes, as detailed later.
When heat flow or underground temperature data are not available,
geothermal potential becomes difficult to assess, except if
multidisciplinary studies are conducted  \citep{Craigetal2021}. 

However, it is well known that temperature is not the only fundamental
parameter to make an exploitable geothermal reservoir. Permeability and
flow rate are essential. In a sense, the role of permeability was
neglected in the last decades, and artificial techniques were developed
to enhance the efficiency of fluid circulation (the ``Enhanced
Geothermal Systems~(EGS)'' concept). 
Yet, instead of looking for anomalously hot areas where
permeability could be enhanced, it might be worth searching for
anomalously permeable zones, such as crustal fault zones 
\citep{Bellangeretal2016,Duwiquetetal2019,Duwiquetetal2021} where deep
hot fluids would naturally rise to shallow crustal levels. This
concept, which prioritizes permeable zones instead of hot zones, may
strongly increase the number of favorable geological systems where
blind geothermal reservoirs would be present. We therefore focus our
study on hydraulically active parts of fault zones, which are
recognized as preferential pathways for fluid flow.

The main objective of this study consists in opening different ways of
thinking and exploring geothermal energy. More precisely, this study is
aimed at (i)~describing the various physical processes at the origin of
the formation of geothermal reservoirs; (ii)~suggesting new
perspectives for geothermal exploration through different views of
high-temperature geothermal systems; (iii)~estimating possible targets
in fault zones for the development of geothermal energy in Europe in
the next decades. After describing several mechanisms at the origin of
thermal anomalies in the shallow crust  (Section~\ref{sec2}), we
present tectonic and geodynamic processes that can affect the thermal
regime of the crust at different scales (in space and time), from
mantle to crustal units  (Section~\ref{sec3}). Then a literature review
focuses on crustal permeability  (Section~\ref{sec4}) and highlights
the role of damage zone geometry on the heat transport mechanisms and
on the conditions for the onset of buoyancy-driven convection 
(Section~\ref{sec5}). We discuss the role of tectonic regimes in the
establishment of thermal anomalies in permeable fault zones. The effect
of local topography, the roles of fault dip and fault thickness, and
the role of tectonic regimes are then investigated by new numerical
results based on previously published numerical procedures 
\citep{Tailleferetal2017,Tailleferetal2018,Guillou-Frottieretal2020,Duwiquetetal2022}. 
Finally, we describe physical processes occurring at more local scales,
in particular the interplay between hydrothermal fluid flow
(buoyancy-driven and/or topography-driven) and local stress regime
(inducing a poroelasticity-driven flow), which may result in
potentially exploitable geothermal reservoirs. Whatever the cause
leading to fluid flow (topography, buoyancy, poroelasticity) our study
is mainly focused on the establishment and the amplitudes of thermal
anomalies at exploitable depths.

\section{Heat flow, temperature anomalies, and geothermal
energy}\label{sec2}

\subsection{Heat flow components}\label{sec2.1}
The thermal regime of the crust is mainly controlled by heat
conduction. By definition, the measured surface heat flow includes the
heat flow coming from the mantle and the crustal component due to
radiogenic heat production  \citep{JaupartMareschal2011}. If mantle
heat flow and thermal properties of crustal rocks (thermal
conductivity, thermal diffusivity, and heat production rates) were well
constrained, then crustal temperatures would be predicted with small
uncertainties. However, this is not the case, and uncertainties in
mantle heat flow values or in heat production rates may under- or
over-estimate crustal temperatures by several tens of\ \textdegree C at
a depth of a few kilometers. 

Mantle heat flow is generally deduced from surface heat flow values and
models of crustal composition. For an average  40~km--thick stable
continental crust, the average surface and mantle heat flow values
amount respectively to 67 and  $30~\mathrm{mW}{\cdot}\mathrm{m}^{-2}$,
with an average crustal heat production providing
$37~\mathrm{mW}{\cdot}\mathrm{m}^{-2}$  [\citealp{Lucazeau2019};  see 
\citealp{Jaupartetal2015} for more detailed statistics on continental
heat  flow]. Below cratonic areas, mantle heat flow is much lower,
around  $10\ndash15~\mathrm{mW}{\cdot}\mathrm{m}^{-2}$ [see Table~7 in 
\citealp{Jaupartetal2015}]. Outside cratonic areas, the mantle heat
flow may exceed $50~\mathrm{mW}{\cdot}\mathrm{m}^{-2}$ 
\citep[e.g.][]{Lucazeauetal1984}.

When radiogenic heat production of crustal rocks is not anomalously
high, the crustal component of surface heat flow is typically around 30
to  $40~\mathrm{mW}{\cdot}\mathrm{m}^{-2}$ [see average crustal
compositions detailed in  \citealp{Jaupartetal2016}]. However, when
surface heat flow is around 
$80\ndash100~\mathrm{mW}{\cdot}\mathrm{m}^{-2}$ or even higher, the
crustal component is necessarily larger due to higher heat production
rates, and the mantle heat flow component probably exceeds 
$30~\mathrm{mW}{\cdot}\mathrm{m}^{-2}$. Even when several geophysical
datasets are taken into account to infer the most probable crustal
models, the uncertainty on mantle heat flow value can reach 15 to 20\% 
\citep{Jaupartetal2015}. In addition, when the temperature dependence
of thermal properties  \citep{ClauserHuenges1995,Whittingtonetal2009}
is not considered, the inferred crustal temperatures may be strongly
underestimated  \citep{Braun2009}. 

Apart from these uncertainties in the context of a thermally conductive
regime, the occurrence of fluid flow in permeable rocks (advective or
convective regimes) may strongly disturb the isotherms that would have
been inferred from a purely conductive regime. For example, the
uprising of crustal fluids through permeable fault zones may induce the
formation of kilometer-sized positive thermal anomalies exceeding
several tens  of\ \textdegree C 
\citep[e.g.][]{Wanneretal2019,Guillou-Frottieretal2020}. In the
following, a positive/negative thermal anomaly corresponds to
temperature values above/below an average conductive temperature--depth
profile.

The knowledge of crustal temperatures is necessary to understand
tectonic regimes and processes. Rock rheology is strongly
temperature-dependent, and the deformation style may evolve from
brittle to ductile for a temperature increase  \citep{Burov2011}. In
addition to spatial variations of crustal temperatures, rheological
contrasts and their tectonic consequences add a temporal variation over
a long timescale. Processes such as sedimentation, erosion, faulting
events, exhumation and topography building, affect crustal temperatures
over several tens of Myrs, thus involving transient evolution of
crustal temperatures that could even lead to partial melting and
magmatism  \citep{JaupartMareschal2011,Chenetal2018}. Because heat
diffusion is a long-term heat transfer mechanism, such large-scale
processes may still affect the present-day temperatures of the upper
crust, leading to the formation of thermal anomalies. 

Although the typical time scale of geothermal systems is much shorter
(tens to hundreds of thousands of years), it is fundamental to consider
past tectonic processes as a potential source of present-day thermal
anomalies. Indeed, the shape of the isotherms in the upper crust 
$({<}$10~km) is particularly controlled by topography 
\citep{Braun2002,Foekenetal2007}. Therefore, rapid relief evolution due
to fault activity can affect the isotherm distribution 
\citep{EhlersFarley2003}. In this case, isotherm disequilibrium is
slower than deformation and thermal anomalies can be generated by past
and no-longer active\unskip\break 
tectonics. 

\subsection{Heat flow and crustal temperatures}\label{sec2.2}
Interestingly, a positive thermal anomaly in the shallow crust is not
necessarily correlated with a high surface heat flow, and a high
surface heat flow is not\unskip\break 
always associated with an underlying positive
thermal anomaly. Yet, geothermal exploration is generally focused
either on areas where surface manifestations are present, or in areas
showing abnormally high surface heat flow. 

Terrestrial surface heat flow represents the energy loss at the surface
of the Earth, through a unit surface
$(\mathrm{W}{\cdot}\mathrm{m}^{-2})$. Expressed as the product of the
vertical temperature gradient by the thermal conductivity, it is thus
not directly correlated to temperature:
{\begin{equation}\label{eq1}
\vec {q}=-\lambda\cdot \vec {\nabla }T
\end{equation}}\unskip\noindent
where $\vec {q}$  ($\mathrm{W}{\cdot}\mathrm{m}^{-2}$) is the surface
heat flow vector, from hot to cold (upward), while the vertical
temperature gradient $\vec {\nabla
}T\ (\mathrm{K}{\cdot}\mathrm{m}^{-1})$ is directed downward, and where
$\lambda\ (\mathrm{W}{\cdot}\mathrm{m}^{-1}{\cdot}\mathrm{K}^{-1})$ is
thermal conductivity. An elevated surface heat flow value can thus be
explained by (i)~a high thermal conductivity with a normal temperature
gradient, (ii)~a high-temperature gradient with a normal thermal
conductivity, or (iii)~a high thermal conductivity and a
high-temperature gradient. An example of a high thermal conductivity
where the temperature gradient is not anomalous is shown in 
Figure~\ref{fig1}, left case (green subvertical body). This example
corresponds to the Thompson Nickel belt, Manitoba 
\citep{Guillou-Frottieretal1996}, where subvertical and narrow bodies
of quartzites (elevated thermal conductivity) do not disturb the
temperature gradient but result in an anomalously high heat flow (label
1). In this case, the abnormal surface heat flow value is not
associated with a thermal anomaly at depth (see the undisturbed white
contours in the colored image in  Figure~\ref{fig1}).

\begin{figure}
\includegraphics{fig01}
\caption{\label{fig1}Case of heat refraction due to lithological
contrasts. From top to bottom: geometries, boundary conditions, and
physical properties; temperature field with isotherms in white contours
separated by  20~\textdegree C; surface heat flow; horizontal
temperature profile at a depth of 4~km. A conductive and narrow body
(in green on the top image) does not disturb the isotherms (white
contours) but results in an elevated surface heat flow (label 1). On
the opposite, the large and insulating body (in pink on the top image)
results in an important uplift of the isotherms, but a small heat flow
anomaly (label 2). However, a temperature of almost  150~\textdegree C
is present at a depth of 4~km (label 3).}
\end{figure}

On the contrary, when large thermally insulating bodies are present
(e.g., a porous sedimentary layer or a thick sequence of ignimbrites,
pink horizontal body on top of  Figure~\ref{fig1}), the isotherms are
uplifted and a positive thermal anomaly develops  (Figure~\ref{fig1},
label 3). However, the product of the low thermal conductivity by an
elevated temperature gradient may result in a normal surface heat flow 
(Figure~\ref{fig1}, label 2). In this case, the underlying thermal
anomaly  (150~\textdegree C at a depth of 4~km, label 3) cannot be
inferred by looking at surface heat flow data only.

Furthermore, thermal conductivity is also influenced by various factors
such as temperature, the dominant mineral phase, porosity and
saturating fluid, and anisotropy  \citep{ClauserHuenges1995}. All these
factors may significantly affect the bulk thermal conductivity of the
rocks. For example, anisotropy of thermal conductivity should be
considered in sedimentary layers, especially in shales and clay-rich
rocks  \citep[e.g.,][]{Davisetal2007}. 

One consequence of these simple heat refraction effects is that
geothermal exploration may be underdeveloped if thermal features are
considered through heat flow data only. To estimate correctly the
underground temperatures, any analysis of heat flow data should be
accompanied by an analysis of thermal conductivity data.

\subsection{Thermal features at Soultz-sous-For\^{e}ts (upper Rhine
graben, France)}\label{sec2.3}
The geothermal site of Soultz-sous-For\^{e}ts is a striking example of
what is addressed in this study. Over the  20~km surrounding the area,
high surface heat flow values range between 121 and 
$183~\mathrm{mW}{\cdot}\mathrm{m}^{-2}$  \citep{Lucazeau2019}, probably
due to a large-scale anomalously high mantle heat flow, associated with
the formation of the European Cenozoic Rift System 
\citep{Dezesetal2004}. However, this high surface heat flow is not
accompanied by hot springs or other surface manifestations of
underlying geothermal\unskip\break 
activity. 

The presence of an anomalous surface temperature gradient in the area
was established at the beginning of the XXth century when tens of
boreholes were drilled in the Pechelbronn oil field, a few kilometers
east of Soultz-sous-For\^{e}ts. With the available bottom-hole
temperature measurements,  \citet{HaasHoffmann1929} established a
temperature map at a depth of 400~m and discovered the thermal anomaly
of Soultz-sous-For\^{e}ts (around  55~\textdegree C at a depth of
400~m). Before deep (up to five kilometers) boreholes were drilled at
Soultz-sous-For\^{e}ts, the available subsurface temperature gradient
measured at depths of a few hundred meters was exceptionally high
$(110~\text{\textdegree}\mathrm{C}{\cdot}\mathrm{km}^{-1})$. Thus, the
existence of a high-temperature geothermal reservoir at a shallow depth
(estimated around  170~\textdegree C at a depth of 1.5~km, black star
in  Figure~\ref{fig2}a, or between 140 and  180~\textdegree C at a
depth of 1900~m, black rectangle in  Figure~\ref{fig2}a) was then
suggested  \citep{GerardKappelmeyer1987}. In the late  80's and in the
90's, when the first shallow boreholes at Soultz were deepened to
greater depths, it turned out that a temperature of only 
140~\textdegree C was reached at a depth of 2~km, and only 
200~\textdegree C at a depth of 5~km  [\citealp{Genteretal2010};
Figure~\ref{fig2}a]. The regularly decreasing temperature gradient from
1~km to 3.5~km depth represents the vertical extent of a nearly
isothermal reservoir, where hydrothermal convection dominates and
homogenizes temperatures 
\citep{ClauserVillinger1990,Guillou-Frottieretal2013}. Similar
temperature--depth profiles (an elevated surface temperature gradient
and a nearly isothermal reservoir over several kilometers) were indeed
measured in many other high-temperature geothermal systems (see
examples of such temperature profiles in\unskip\break  
Figure~\ref{fig2}b). 

\begin{figure*}
\includegraphics{fig02}
\caption{\label{fig2}Temperature profiles in geothermal systems: 
(a)~measured temperature profile (grey) in the GPK2 borehole at
Soults-sous-For\^{e}ts (Alsace) and modifications of the temperature
gradient (dashed lines), after  \citet{Genteretal2010}. The black star
and the black rectangle indicate the predictions in 1987 (see text).
(b)~General shape of temperature profiles in different geothermal
systems: Fenton Hill  \citep{Harrisonetal1986}, Svartsengi and
Nesjavellir  \citep{Steingrimsson2013}, Germencik 
\citet{InancTureyenetal2014}, Bouillante 
\citep{Guillou-Frottier2003}, Tiwi  \citep{Mooreetal2000}, Salton
Sea  \citep{Sassetal1986}, Kakkonda  \citep{Doietal1998}.}
\end{figure*}

With the example of Soultz-sous-For\^{e}ts in mind, where no hot
springs or fumaroles are visible, one may postulate that the discovery
of new high-temperature geothermal systems is still plausible, as soon
as hydrothermal convection can occur in the shallow crust. The
occurrence of hydrothermal convection requires that the host rock must
be sufficiently permeable. Since the density of hydrothermal fluids is
temperature-dependent, whatever the fluid composition, hot and light
fluid at a depth of several kilometers would naturally rise to shallow
levels, where less permeable rocks will enclose the geothermal
reservoir. 

Beside the occurrence of buoyancy-driven convection  [or ``free''
convection, as detailed in  \citealp{Alt-Eppingetal2021}], non-Rayleigh
convection and topography-driven (``forced'') convection may play a
significant role in fluid flow and thermal regimes. In this study,
although we put forward the role of damage zone geometry 
(Section~\ref{sec4.3}) on the occurrence of buoyancy-driven convection
(see  Section~\ref{sec5.2}), the other driving forces are also
considered (topography-driven and poro-elasticity-driven forces).

\section{Large-scale (10--1000~km) thermal anomalies: geodynamical
control}\label{sec3}

\subsection{Large-scale heat flow anomalies}\label{sec3.1}

\begin{figure*}
\includegraphics{fig03}
\vspace*{3pt}
\caption{\label{fig3}Three examples of interpolated maps of surface
heat flow data (top), from  \citet{ArtemievaMooney2001}, and synthetic
west-east cross-sections of geodynamic contexts (bottom). Long
wavelength heat flow anomalies (from red to purple) are visible in
southeast Australia (a), from the French Massif Central to the Rhine
graben (b), and in western North America (c).}
\vspace*{5pt}
\end{figure*}

Although lithological contrasts may result in small-scale discontinuous
heat flow variations  (Figure~\ref{fig1}), long wavelength heat flow
anomalies can be evidenced by interpolation maps such as the one
published by  \citet{ArtemievaMooney2001} (Figure~\ref{fig3}).
Anomalous heat flow values  [higher than the average continental heat
flow of  $67~\mathrm{mW}{\cdot}\mathrm{m}^{-2}$, 
\citealp{Lucazeau2019}] can be delineated over several hundreds of
kilometers. In southeast Australia  (Figure~\ref{fig3}a), a 
${\sim}$250~km large zone shows high heat flow values
$(92\pm10~\mathrm{mW}{\cdot}\mathrm{m}^{-2})$, which are interpreted as
resulting from high heat production of basement rocks  [an average
value of  $6~\upmu\mathrm{W}{\cdot}\mathrm{m}^{-3}$ is obtained,
instead of less than  $3~\upmu\mathrm{W}{\cdot}\mathrm{m}^{-3}$ outside
the area;  \citealp{Neumannetal2000}].

\begin{figure*}
\includegraphics{fig04}
\caption{\label{fig4}Effect of mantle heat flow on temperature--depth
profiles for a three-layer crust enriched in radioactive elements (see
Section~\ref{sec3.2}). The crustal geotherms on the left  (a) are
zoomed at shallow depth in (b), where the red arrow indicates the
temperature domain for high-temperature geothermal systems;
(c)~transient effect of an instantaneous increase of  $Q_m$ (from  20
to $50~\mathrm{mW}{\cdot}\mathrm{m}^{-2}$ during 10~Myrs).}
\vspace*{4pt}
\end{figure*}

In the French Massif Central  (Figure~\ref{fig3}b), a broad zone of
heat flow values exceeding  $100~\mathrm{mW}{\cdot}\mathrm{m}^{-2}$ can
be delineated over a distance of  ${\sim}$300~km 
\citep{LucazeauVasseur1989}.  \citet{Lucazeauetal1984} showed that
crustal heat production alone cannot explain this long wavelength
anomaly. They suggest that an additional mantle heat flow component,
reaching  $25\ndash30~\mathrm{mW}{\cdot}\mathrm{m}^{-2}$ (resulting in
a maximum mantle heat flow around 
$60~\mathrm{mW}{\cdot}\mathrm{m}^{-2}$) and associated with the recent
Cenozoic activity, corresponds to the signature of a mantle diapir.
According to recent thermo-mechanical models of plume--lithosphere
interactions, small-scale mantle upwellings (or ``baby-plumes'') might
have played a significant role in continental rifting 
\citep{Koptevetal2021}. 

The third example (Figure~\ref{fig3}c) illustrates the long wavelength
heat flow anomaly along the western coast of North America, below which
the Farallon oceanic plate had subducted.  \citet{ErkanBlackwell2008}
detailed the different subduction-related geodynamical events that
could explain the large-scale positive and negative heat flow
anomalies. In particular, the geothermal systems of the Basin and Range
Province, Nevada, would be associated with subduction-induced
extensional tectonics and with surface heat flow higher than 
$80~\mathrm{mW}{\cdot}\mathrm{m}^{-2}$ \citep{WisianBlackwell2004}.
Slab tears and their thermal consequences might also have played a role
in establishing this large-scale heat flow anomaly, promoting the
melting of the ascending asthenosphere, particularly in the northeast
portion of Baja California region  \citep{Ferrarietal2018}. At smaller
scales, the high surface heat flow values in western Anatolia, Turkey,
could be related to underlying subduction-related mantle flow (slab
rollback and slab tear processes), inducing high mantle heat flow 
\citep[e.g.,][]{Rocheetal2018,Rocheetal2019}. Similar processes may be
at the origin of the extensional tectonics in Tuscany, Italy, where the
Larderello-Travale geothermal field is characterized by an high
regional surface heat flow exceeding 
$100~\mathrm{mW}{\cdot}\mathrm{m}^{-2}$ over more than 100~km 
\citep{Boninietal2014}.

\subsection{Effect of mantle heat flow on crustal
temperatures}\label{sec3.2}
The previous examples suggest that deep thermal anomalies that occurred
in the past have diffused towards the surface, leaving a thermal
imprint that can be deciphered today.  
Figures~\ref{fig4}a
and~b illustrate the effect of mantle heat flow $(Q_m)$ on
steady-state crustal temperatures for a three-layer crust, where heat
production rates are twice as large as the ones described for an
averaged crustal composition. A 35~km thick crust is made of a 13~km
thick upper crust (thermal 
conductivity of 
$3.0~\mathrm{W}{\cdot}\mathrm{m}^{-1}{\cdot}\mathrm{K}^{-1}$,  heat
production of $3.2~\upmu\mathrm{W}{\cdot}\mathrm{m}^{-3})$, an
intermediate 11~km thick crust
($2.5~\mathrm{W}{\cdot}\mathrm{m}^{-1}{\cdot}\mathrm{K}^{-1}$ and 
$1.5~\upmu\mathrm{W}{\cdot}\mathrm{m}^{-3}$) and a 10~km thick lower
crust  ($2.0~\mathrm{W}{\cdot}\mathrm{m}^{-1}{\cdot}\mathrm{K}^{-1}$
and  $0.4~\upmu\mathrm{W}{\cdot}\mathrm{m}^{-3}$). For this model, the
crustal\unskip\break 
component reaches $62.1~\mathrm{mW}{\cdot}\mathrm{m}^{-2}$. A
mantle heat flow of  $20~\mathrm{mW}{\cdot}\mathrm{m}^{-2}$, leading to
a surface heat flow $(Q_{s})$ of 
$82.1~\mathrm{mW}{\cdot}\mathrm{m}^{-2}$, is not sufficient to get 
150~\textdegree C in the first 6~km of the crust. However, a  $Q_{m}$
value of  $50~\mathrm{mW}{\cdot}\mathrm{m}^{-2}$ 
$(Q_{s}=112.1~\mathrm{mW}{\cdot}\mathrm{m}^{-2})$ allows such high
temperatures, but only below a depth of 4~km.  Figure~\ref{fig4}c
illustrates the transient effect of a thermal pulse coming from the
mantle. Considering the steady-state temperature profile of 
Figure~\ref{fig4}a with a mantle heat flow value of 
$20~\mathrm{mW}{\cdot}\mathrm{m}^{-2}$, an instantaneous thermal pulse
of $30~\mathrm{mW}{\cdot}\mathrm{m}^{-2}$  (providing a mantle heat
flow of  $50~\mathrm{mW}{\cdot}\mathrm{m}^{-2}$) is applied during
10~Myrs. Mantle heat flow is then reduced instantaneously to its
initial value of $20~\mathrm{mW}{\cdot}\mathrm{m}^{-2}$. It results
that the induced temperature change at a depth of 5~km does not exceed 
24~\textdegree C. For a thermal pulse of 
$50~\mathrm{mW}{\cdot}\mathrm{m}^{-2}$, temperatures exceeding 
150~\textdegree C are found below a depth of 4100~m.

In summary, whatever the chosen crustal models, in the absence of
shallow magmatic intrusions, and for reasonable values of heat
production rates (from   $0.5~\upmu\mathrm{W}{\cdot}\mathrm{m}^{-3}$ in
the lower crust to  $4.0~\upmu\mathrm{W}{\cdot}\mathrm{m}^{-3}$ in the
upper crust), temperatures above  150~\textdegree C are difficult to be
reached at shallow depths  (1--4~km). If heat transfer occurs only by
conduction, then abnormally high heat production rates are required, as
for some Australian granites. Otherwise, high temperatures 
$({>}150~\text{\textdegree}\mathrm{C})$ at shallow depths can only be
reached if fluid circulation is taken into account.

\section{Small-scale (0.1--10~km) thermal anomalies: structural
control}\label{sec4}
This section is dedicated to abnormally permeable zones of the crust,
particularly fault zones within which fluid circulation induces thermal
anomalies at shallow depths  (Section~\ref{sec5}). Here, we first focus
on crustal permeability and transient effects on permeability as
detailed in  \citet{Taillefer2017}, and then on fault zone architecture
and damage zone properties.

\subsection{Crustal permeability}\label{sec4.1}
Permeability is a key parameter that controls the occurrence of fluid
flow (whether driven by topography or buoyancy) and thus the location
and the amplitude of thermal anomalies 
\citep[e.g.,][]{Tailleferetal2018,Guillou-Frottieretal2020,Alt-Eppingetal2021}.
Crustal permeability is considered to be strongly depth-dependent 
\citep{Andersonetal1985,IngebritsenManning1999,Shmonovetal2003,SaarManga2004,StoberBucher2007,EarnestBoutt2014,Ranjrametal2015},
and recent studies showed that the permeability at depth is the main
parameter that controls the longevity of a geothermal system 
\citep[e.g.][]{IngebritsenGleeson2017}. Exponential first-order
relationships were established between permeability and depth 
\citep{Andersonetal1985,Shmonovetal2003,SaarManga2004}, and the few
hydraulic data on crystalline rocks  \citep{StoberBucher2007} seem to
confirm this dependence (red curve in  Figure~\ref{fig5}). However,
based on large data sets, recent studies showed that permeability
values in the shallow crust (${<}$2.5~km) are mainly controlled by
lithology  [\citealp{Neuzil2019},  Figure~\ref{fig5}] and tectonic
setting  \citep{Ranjrametal2015}. Unfortunately, permeability evolution
in the targeted depth window for geothermal energy (1--4~km), is still
not well constrained due to the lack of data 
\citep{StoberBucher2007,Ranjrametal2015}. 

\begin{figure*}
\includegraphics{fig05}
\vspace*{5pt}
\caption{\label{fig5}Range of permeability  (m$^{2}$) vs. depth (km)
for different lithologies in the brittle continental crust \citep[modified
after][]{Neuzil2019}. Green and blue are site
($10^{-2}$--$10^{2}~\mathrm{km}^{2}$) and regional 
($10^{3}$--$10^{6}~\mathrm{km}^{2}$) clay and shale data, other colors
are sedimentary/crystalline rocks and metamorphic/crystalline rocks 
\citep{TownendZoback2000}, except brown for basalt and basaltic
andesite  \citep{SaarManga2004}. The red curve represents hydraulic
data from  \citet{StoberBucher2007}.}
\vspace*{10pt}
\end{figure*}

The high sub-surface permeability  (Figure~\ref{fig5}) is related to
supergene alteration processes and the fractured decompression zone 
\citep{Dewandeletal2006}. Interestingly, rocks affected by weathering
can preserve these high permeability characteristics when subsequently
buried  \citep{Lachassagneetal2011}. This is the case at
Soultz-sous-For\^{e}ts, where the permeability of the altered and
fissured horizon at the top of the granitic basement is preserved at
depth  \citep{Genteretal2010}.

Fractures can increase the permeability of protoliths by several orders
of magnitude  \citep{SonneyVuataz2009,Coxetal2015}. Because of their
interconnection, the set of fractures is also called a ``fracture
network'', whose permeability is dependent on several factors such as
the fracture morphology, the spatial organization of the fracture
network (density and connectivity), the fillings, and the ongoing
tectonic activity that keeps the fractures open  \citep{Crider2015}.
Fractures can be closed under confining pressure, or cementation, if
they are not held open by tectonic stresses or fluid pressure 
\citep{Laubachetal2004}. Larger stress magnitudes in tectonically
active regions may produce, at depth, fracture apertures larger than
expected  \citep{EarnestBoutt2014}. Fracturing associated with the
fault zones will have a significant impact on creating or 
maintaining\unskip\break
permeability, particularly at the depths targeted for deep geothermal
energy (1--4~km).

\subsection{Transient effects on permeability variations}\label{sec4.2}
The fractures, the fluid pressure, and the fluid composition play a key
role in the permeability evolution of the geothermal system 
\citep{Renardetal2000}. Numerous studies on hydrothermal systems show
that the permeability varies in time as well as space 
\citep[e.g.,][]{Coxetal2015}. It is therefore relevant to consider the
concept of ``dynamic permeability''  \citep{IngebritsenGleeson2017} and
detail the transient effects that can impact a geothermal system. In
particular, mineralization processes, or fault activity and consequent
fluid flow may decrease or increase permeability, respectively, as
detailed below.

The whole process of fracture sealing by mineralization plays an
important role in the lifespan of a geothermal system, and can
significantly reduce its permeability 
\citep{Renardetal2000,Griffithsetal2016}. Rock/fluid interaction will
be responsible for the nature of mineralization, and mineralization
can completely block fluid circulation  \citep{StanekGeraud2019}.
\citet{Renardetal2000} show that at the scale of a fracture,
pressure/solution mechanisms are faster (a few years) than fracture
cementation (from a hundred years to several million years), allowing
fluid circulation to be maintained over time. In fault zones that
localize fluid circulation, the presence of phyllosilicates and clays
is frequently reported  \citep[e.g.,][]{NishimotoYoshida2010}, and can
change the permeability properties 
\citep[e.g.,][]{Yieldingetal1997,Leclereetal2012}. 

Fault activity can increase permeability. Major tectonic setting
changes can lead to migration or modification of the main parameters of
the geothermal system: spring location, flow rates and temperatures 
\citep{Husenetal2004,IngebritsenManning2010,Coxetal2015,Tailleferetal2021}. 
\citet{Fauldsetal2010} showed that, in the Basin and Range (western
USA), geothermal systems are more likely associated with active faults.
However, the macro-seismicity is not a critical parameter for the
preservation of a permeable fracture network, if a favorably oriented
stress regime can maintain the permeability required for fluid flow 
\citep{Bartonetal1995,TownendZoback2000}. In this case, geothermal
systems reach a state of equilibrium allowing constant hydrothermal
activity  \citep{Tailleferetal2018}. 

The fluid flow itself influences permeability with fluid overpressures
in low-permeability or sealed zones \citep{Sibson1992}. When these
overpressures are located along a fault, they can initiate slip on the
fault and associated earthquakes  \citep{Cox2016}. Such a fault-valve
behavior  \citep{Sibson1992} allows post-failure discharge of fluids.
Importantly, in most cases, this leads to the upward migration of
fluids along the fault. In addition, the activity of the fault related
to these overpressures allows the re-opening of closed or sealed
fractures, and the generation of new fractures 
\citep{Schoenballetal2020}.

Finally, climatic factors can also influence fluid flow in a geothermal
system. Water from numerous geothermal systems has a meteoric origin
\citep[e.g.,][]{GrasbyHutcheon2001,SonneyVuataz2009,Diamondetal2018},
and therefore depends on climatic conditions. Fluid flow and
temperature of the water can be correlated with rainfall
\citep{Ratouisetal2017}. It is important to note that rainfall
correlation is not systematic for a long geothermal cell 
\citep{Tailleferetal2018}. Recent studies by 
\citet{Alt-Eppingetal2021,Alt-Eppingetal2022} demonstrated that flow
patterns can change during periods of glaciation. For example,
non-Rayleigh convection during glacial periods would evolve towards
topography-driven flow during interglacial periods. Glaciation can also
have several implications:  (1)~the loading effect of glaciers causes
local variations in stress  \citep{Neuzil2012};  (2)~the melting of ice
and the associated infiltration of cold fluids leads to a cooling of
the massif  \citep{Marechaletal1999};  (3)~the thickness of ice on the
relief can reduce the infiltration of meteoric water \break
\citep{Volpietal2017}. 

\subsection{Fault zone architecture and damage zone
thickness}\label{sec4.3}
The fault zones constitute one of the most permeable zones in the
Earth's crust (see  Sections~\ref{sec4.1} and~\ref{sec4.2}). It is
therefore relevant to detail the main characteristics and the
fracturation volume associated with these structural objects. Two main
zones are distinguished inside a deformed and altered rock volume in a
fault zone:  (i)~one, or several Core Zone(s) (CZ) characterized by the
presence of fault rocks frequently associated with clay or thin matrix 
$({<}0.1~\mathrm{mm})$ and well-developed cementation. These zones are
mostly recognized as a barrier for fluid circulation when the fault is
not active  \citep[e.g.][]{Sibson2000,Fossen2016}.  (ii)~two, or
several Damage Zone(s) (DZ), characterized by an increase in fracture
density compared to the protolith 
\citep[e.g.][]{Sibson2000,Faulkneretal2003}. The damage due to the
fault activity increases the permeability of the protolith by several
orders of magnitude, thus allowing fluid circulation 
\citep{MitchellFaulkner2012}. The damage zone is composed of an inner
part, highly fractured, commonly associated with intense fluid
circulation and an outer part, less fractured, but still more permeable
than the protolith  \citep[e.g.,][]{Choietal2016}.

The degree of maturity of the fault exerts a strong control on the
thickness of the DZ that localizes the maximum permeability. 
Figure~\ref{fig6} shows that the thickness of the DZ is correlated to
the fault displacement, a proportional relationship exists for
displacement until  ${\sim}150~\mathrm{m}$, regardless of the kinematic
of the fault and the nature of the protolith 
\citep{Benseetal2013,Mayolleetal2021}. Above 150~m, DZ thickness is
limited to a few hundred meters 
\citep[e.g.,][]{MitchellFaulkner2009,Benseetal2013,Mayolleetal2019}.
Characterization of the DZ thickness for faults with displacement above
hundred meters shows more uncertainty due to the presence of secondary
faults and the difficulty to find a continuous outcrop for fracture
quantification  \citep{Choietal2016}. Several models are proposed for
the bending of the DZ thickness observed in  Figure~\ref{fig6}: the
thickness of the brittle crust  \citep{AmpueroMao2017}, segmentation of
the fault  \citep{Schlischeetal1996}; accumulation of displacement and
fault damage thickness by incorporating damages from the former fault
segments  \citep[][and references
therein]{Mayolleetal2019,Mayolleetal2021}. Along pluri-kilometric
faults, the geometry of the fault itself has a strong control on the
permeability. Indeed, fault relays, fault terminations and fault
intersections are more favorable zones to channel fluid circulation due
to an increase in fracture density 
\citep{CurewitzKarson1997,FauldsHinz2015}. 

\begin{figure}
{\vspace*{4pt}}
\includegraphics{fig06}
{\vspace*{4pt}}
\caption{\label{fig6}Total fault zone thickness as a function of
displacement, modified from  \citet{Mayolleetal2019}. The damage zone
thickness is defined as the total adjacent damaged rock volume on both
sides of the fault. The bending of the linear scaling law trend is
highlighted in light brown.}
{\vspace*{4pt}}
\end{figure}

\section{Fluid circulation in fault geothermal reservoirs}\label{sec5}
Fluid flow in a permeable medium can be due to different driving forces
\citep{Alt-Eppingetal2021}. In our study, the permeable medium consists
of a finite width fault zone. Thus, the classical Rayleigh number
expression does not apply (see  Section~\ref{sec5.2}). Although we
mainly focus our study on buoyancy-driven convection, topography-driven
and poro-elasticity-driven forces will be also investigated 
(Sections~\ref{sec5.3} and~\ref{sec5.4}). The effects of fault
thickness, the impact of tectonic regimes or the role of the
poro-elasticity force are discussed but not compared to each other.
This  section is thus dedicated to illustrate independently different
causes leading to the establishment of shallow thermal anomalies.

\subsection{Fluid flow in fault zones and related thermal
anomalies}\label{sec5.1}
Apart from geological evidence of fluid circulation within fault zones
(see  Section~\ref{sec4}), theoretical and numerical studies brought
quantitative aspects related to the required conditions for thermal
convection to occur  \citep{ForsterSmith1989,MalkovskyMagri2016}. When
buoyancy-driven flow occurs within a fault zone, temperature anomalies
develop and may result in large volumes of high 
$({>}150~\text{\textdegree}\mathrm{C})$ temperatures at shallow depths
(1--4~km). For example,  \citet{Wanneretal2019} suggested that large
ellipsoidal thermal plumes are generated by fault-hosted orogenic
geothermal systems, such as the Grimsel Pass case, Swiss Alps, where
hot upwellings would provide hundreds of PJ  $(10^{15}~\mathrm{J})$ of
anomalous heat per km depth.  \citet{Guillou-Frottieretal2020}
illustrated similar convective patterns where temperature anomalies may
exceed  70~\textdegree C at a depth of 500~m. These studies demonstrate
that temporal and spatial variations of permeability play a key role in
the morphology and amplitude of temperature anomalies. They also show
that temperature anomalies are larger and higher for vertical fault
zones  \citep[see also][]{Duwiquetetal2019} and for wide damage zones.
Because damage zone thickness may exceed several hundreds of meters, it
is worth studying through simple theoretical considerations how fault
zone thickness influences thermal convection features.

\subsection{Role of fault zone thickness}\label{sec5.2}
\citet{MalkovskyMagri2016} studied the conditions for which thermal
convection occurs within a fault zone of thickness  $2\delta$ and
height  $H$, embedded into impermeable rocks. They extended the
constant viscosity case developed by  \citet{MalkovskyPek1997} by
considering a temperature-dependence of the fluid viscosity (as for the
case of hydrothermal fluids), characterized by a dimensionless
parameter,  $\gamma$. Their linear stability analysis provides a new
expression for the critical Rayleigh number $Ra^{\mathrm{crit}}$, where
the fault aspect ratio  $\Delta=\delta/H$ appears:
{\begin{eqnarray}
Ra^{\mathrm{crit}}(\Delta ,\gamma )&=&(1+0.493\gamma +0.12\gamma ^{2})
\nonumber\\
&&\times\,\left[\left(\frac{8.19}{\Delta }\right)^{5/4}+
(4\uppi^{2})^{5/4}\right]^{4/5}\label{eq2}
\end{eqnarray}}\unskip\noindent
It can be seen that for an infinitely wide medium and for a constant
viscosity  $(\gamma= 0)$, the $4\uppi^{2}$ value  \citep{Lapwood1948}
is recovered. For an average temperature gradient  $(\gamma\sim2)$,
the critical Rayleigh number can be calculated for different $\Delta$
values  (Figure~\ref{fig7}a). Equation~(\ref{eq2}) indicates that for
wide fault zones (i.e. large $\Delta$ or  $\delta$ values), 
$Ra^{\mathrm{crit}}$ will be small. In other words, all other
parameters being equal, thermal convection will occur more easily
within large fault zones. This is illustrated in  Figure~\ref{fig7}b,
where $\Delta$ varies from 0.005 to 0.05: the critical Rayleigh number
decreases from 4059 to 455. If the computed Rayleigh number equals
1000, then convection will not develop in the thin part of the fault
zone. On the opposite, within the wide part, convection will be
well-developed and isotherms can be strongly raised up (red curve 
in\unskip\break 
Figure~\ref{fig7}b).

\begin{figure}
\includegraphics{fig07}
\vspace*{-4pt}
\caption{\label{fig7}(a)~Role of the fault thickness  $(2\delta)$ (or
half-thickness over height ratio, $\Delta$) on the critical Rayleigh
number, after the theoretical law of  \citet{MalkovskyMagri2016};
(b)~sketch of a fault zone whose thickness increases from 20~m to
200~m: the critical Rayleigh number decreases by a factor 9, and
isotherms can be strongly distorted within the wide part.}
\vspace*{-8pt}
\end{figure}

When the Rayleigh number expression is considered [see
Equation~(12a) in  \citealp{MalkovskyMagri2016}], the critical
permeability above which convection occurs can be estimated for
different $\Delta$ ratios. The fault zone height is here fixed at 5~km,
and a temperature difference of  200~\textdegree C is assumed, while
other physical properties are well-known. For a fault thickness of
100~m  ($\delta=50$~m; $\Delta=0.01$), the critical Rayleigh number
equals 2050  (Figure~\ref{fig7}a) and the critical permeability above
which thermal convection occurs equals 
$2.6\times10^{-14}~\mathrm{m}^{2}$. For a fault thickness of 500 and
1000~m, the critical permeability decreases to 
$5.9\times10^{-15}~\mathrm{m}^{2}$, and to 
$1.4\times10^{-15}~\mathrm{m}^{2}$,\unskip\break 
respectively. 

However, these estimates depend on several hypotheses, such as the
chosen expression of the Rayleigh number  [see  Equations~(12a)
and~(12c) in  \citealp{MalkovskyMagri2016}]. Another way to
understand qualitatively the role of fault thickness on the temperature
field is to run simple 2D numerical experiments where realistic
physical parameters are chosen and where the thickness of the damage
zone varies. This 2D simplification is valid as soon as fault length
can be considered pluri-kilometric. The numerical procedure and set-up
are detailed in the Appendix~\ref{appA}. Thermal convection occurs
within the fault damage zone, and the maximum temperature at different
depths is recorded.  Figure~\ref{fig8}a illustrates this maximum
temperature at a depth of 3~km in a fault zone permeability--fault
zone width diagram. For example, for a fault width of 800~m where
permeability equals  $5\times10^{-15}~\mathrm{m}^{2}$,  120~\textdegree
C can be reached at a depth of 3~km. To get  150~\textdegree C for the
same fault width, permeability must be larger than 
$3\times10^{-14}~\mathrm{m}^{2}$.  Figure~\ref{fig8}b illustrates the
depths for which 100,  120 and 140~\textdegree C are observed, for four
different fault widths. For large fault thicknesses, temperatures above
100~\textdegree C are observed at depths lower than 2~km. On the right
of light orange areas (conductive regimes), maximum temperatures are
controlled by the convective regimes, where downwellings or upwellings
can decrease or increase the maximum recoded\break temperature.

\begin{figure}
\includegraphics{fig08}
\vspace*{-3pt}
\caption{\fontsize{9.1}{10.92}\selectfont\label{fig8}Results 
from 2D numerical models of thermal
convection within a fault-damage zone (see details in the text and
Appendix~\ref{appA}).  (a)~Maximum temperature observed at a depth of
3~km. For a 600~m wide damage zone, the temperature reaches 
150~\textdegree C for a permeability of
$6\times10^{-14}~\mathrm{m}^{2}$. If the fault-damage zone is twice as
wide (1200~m),  150~\textdegree C is reached for a permeability of 
$7.5\times10^{-15}~\mathrm{m}^{2}$. Curves for  120~\textdegree C, 
180~\textdegree C, and  200~\textdegree C are also shown.  (b)~Depths
of maximum temperatures observed (100, 120 and  140~\textdegree C) for
different fault widths as a function of fault zone permeability. Areas
in light orange correspond to conductive regimes.}
\end{figure}

\subsection{Impact of tectonic regimes}\label{sec5.3}
The three basic tectonic regimes, that is, extensional (normal),
compressional (reverse), and strike-slip, result from different stress
fields and express themselves in the field by typical geometric and
geomorphologic features (e.g., fault dip, form of the scarp
relief). Some of these differences are here tested in a 2D parametric
study considering a fault zone and its related topography in only one
of the fault compartments  (Figure~\ref{fig9}a). Mechanical stresses
are not explicitly taken into account in these models, where only
thermal and hydraulic processes are coupled, but their effects on fault
geometry and geomorphology are reproduced.

\begin{figure*}
\vspace*{-3pt}
\includegraphics{fig09}
\vspace*{-5pt}	
\caption{\fontsize{9.1}{10.92}\selectfont\label{fig9}Parametric 
study on the geometric features of a
scarp relief related to a fault offset. Colors associated with typical
fault dips (green $=$ reverse  30\textdegree, blue  $=$ strike-slip
90\textdegree, normal  $=$ 60\textdegree) are common to the three
figures.  (a)~Sketch of the numerical model, the tested parameters are
in red. (b)~Influence of the fault dip on the temperature distribution
(color chart), along a profile at a depth of 4000~m below the flat
right compartment (dashed black line). Thick and thin lines on the
sketches correspond to the isotherms and flow lines (colored according
to the temperature), respectively. Dashed colored lines correspond to
the position of the core zone at the surface.  (c)~Comparison of the
maximum temperatures reached at a depth of 4000~m for the three tested
parameters.}
\end{figure*}

\subsubsection{Numerical model}\label{sec5.3.1}
The numerical model is a 2D-block consisting of two compartments of the
basement with a depth-decreasing homogeneous permeability. A fault zone
composed of one 10~m-thick core zone with a constant permeability 
$K_c$, surrounded by two damage zones with a constant permeability $K_d$,
separates the two compartments. The choice of a constant permeability
$K_d$ may be questionable since it may depend on the tectonic regime.
However, based on the 
\citeauthor{Scibek2020}
database
\citeyearpar{Scibek2020}, we estimate an
average permeability of fault zones around  $10^{-14}~\mathrm{m}^{2}$,
whatever the tectonic regime. We thus fixed this value to allow fluid
circulation within damage zones for all experiments. The thickness of
the total fault damage zone (DZT, two half damage zones) is tested for
the values DZT $=$ 50, 300, and 1000~m. The fault dip is tested through
three typical dip values according to the main tectonic regimes:
reverse  (30\textdegree), strike-sip (90\textdegree), and normal
(60\textdegree). The left basement compartment (i.e., the hanging wall
of the reverse fault or the footwall of the normal fault) represents a
scarp relief associated with the fault offset which elevation (Scarp
Relief Elevation, SRE) is tested for the values  SRE $=$ 200, 1000, and
3000~m. The model is examined along a horizontal temperature profile at
4000~m below the surface without topography (i.e. the right
compartment). A mixed thermal boundary condition is applied to the
topographic surface. This mixed thermal condition relates the surface
heat flow with surface temperature through a heat transfer coefficient 
\citep[see details in][and in the 
Appendix~\ref{appA}]{Tailleferetal2018}. This allows variations in
surface temperatures, as illustrated by, e.g., the presence of hot
springs  \citep[see also][]{Magrietal2015}. The variables and
parameters used in the models are available in the Appendix~\ref{appA}.

\subsubsection{Numerical results}\label{sec5.3.2}

\paragraph{Influence of the fault dip}
The numerical results of  Figure~\ref{fig9}b represent the temperature,
the main isotherms, and flow lines colored according to the fluid
temperature inside the model for the three fault dips (FD)
configurations. For greater readability, only the intermediate
combination is presented in\unskip\break 
Figure~\ref{fig9}b, i.e. SRE $=$ 1000~m,
DZT $=$ 300~m, but all the parametric combinations were tested, and
show similar patterns. The flow lines in  Figure~\ref{fig9}b show that
cold fluids infiltrate at the model surface and then go down to depths
of several kilometers under the effect of the pressure gradient. They
warm during their descent under the effect of the geothermal gradient,
and reach a temperature of  150~\textdegree C around a depth of 4000~m.
These warm fluids are caught by the permeable damage zones where they
can easily circulate and reach the surface. In any fault dip
configuration, both the sketches and the temperature profiles at a
depth of 4000~m show that the deformed isotherms at the fault vicinity
create a positive thermal anomaly  (Figure~\ref{fig9}b). However, the
amplitude and the form of this anomaly depend on the fault dip. The
thermal anomaly associated with the reverse fault 
$(\mathrm{FD}=30\text{\textdegree})$ has the widest lateral amplitude
and is off-center with respect to the core zone. Conversely, the
strike-slip  $(\mathrm{FD}=90\text{\textdegree})$ and the normal 
$(\mathrm{FD}=60\text{\textdegree})$ fault anomalies are centered on
the fault core zones. The temperature anomaly related to the normal
fault is asymmetric, with the temperature in the left (relief)
compartment being more elevated than in the right (flat) compartment.
The fault dip, i.e., the tectonic regime, could then strongly influence
the temperature distribution at the fault zone\break vicinity.

\paragraph{Relative importance of the tested parameters on the maximum
temperature}
The maximum temperature reached along the horizontal temperature
profile at a depth of 4000~m gives an estimate of the relative
importance of the three tested geometrical parameters: the total damage
zone thickness DZT, the scarp relief elevation SRE, and the fault dip
FD  (Figure~\ref{fig9}c). The greatest influence on the maximum
temperatures is related to the scarp relief elevation (with a maximum
difference of  34~\textdegree C between SRE $=$ 200~m and SRE $=$
3000~m): the higher the scarp elevation relief, the higher the
temperature. Fault dip and damage zone thickness have a similar
influence on the maximum temperature, with a maximum difference of 
5~\textdegree C between normal and reverse fault (SRE  $=$ 200~m, DZT
$=$ 1000~m), and between a damage zone of 50 and 1000~m (SRE  $=$
3000~m, strike-slipe fault), respectively. It seems that there are no
clear rules about which fault dip or damage zone thickness is more
adapted for reaching the highest temperatures at a depth of 4000~m.
However, additional results (see Appendix~\ref{appA}) show that at
shallow depths (2000~m), fault dip has a stronger influence on the
maximum temperatures  (8~\textdegree C of difference), with reverse
faults recording the maximum temperatures for all the parameter
combinations. At these shallower depths, the damage zone thickness has
a poor influence on the maximum temperature. Nevertheless, the impacts
of the fault dip and damage zone thickness on the maximum temperature
are weak compared to those due to the topography.

\subsection{Role of the poroelasticity-driven force}\label{sec5.4}
Tectonic deformation influences fluid flow in different geological
contexts  \citep{Bethke1985,RowlandSibson2004,Eldursietal2021}. Within
a basement fault zone, such as the one described above (i.e., an
anomalously permeable zone), a single stress  $(S_{H\mathrm{max}})$
impacts the general convective dynamics and affects the thermal field
\citep{Duwiquetetal2021}. The distribution of positive and negative
temperature anomalies in these permeable areas depends, among others,
on tectonic regimes, whose poroelasticity-driven force has been
identified  \citep{Duwiquetetal2022}. These effects were highlighted
both by quantifying a very small permeability variation (less than
1\%), but also by recording fluid pressure variations related to the
tectonic regimes implemented.

\begin{figure*}
\includegraphics{fig10}
\caption{\label{fig10}3D numerical model of hydrothermal convection
within a 300~m wide fault zone embedded in impermeable rocks. A fixed
temperature and a fixed pressure are imposed at the upper surface. The
no-flow and insulating conditions apply to lateral boundaries. A fixed
heat flux and a no-flow condition are imposed at the bottom boundary.
Results correspond to a steady-state regime.  (a)~Boundary conditions
and result for the experiment without applying tectonic stresses; 
(b)~illustration of stress applications as boundary conditions for the
compressional tectonic regime (with 
$S_{H\mathrm{max}}\geq S_{h\mathrm{min}}\geq S_V$);
(c,d): thermal (red color for temperature greater than  150~\textdegree
C) and velocity fields (colored arrows along streamlines);  (e)~fluid
pressure variations for each experiment.}
\end{figure*}

Based on poroelasticity assumptions, we propose to investigate how the
poroelasticity-driven force, caused by tectonic regimes, can have an
impact on the fluid flow for a 0.3~km fault thickness (here the fault
core is not taken into acount). This fault thickness is typical of
natural fractured systems such as the T\^{e}t fault 
\citep{Tailleferetal2018,Milesietal2020} or the Seferihisar-Bal\c{c}ova
fault-hosted geothermal system  \citep{Magrietal2012}. Considering a
simplified 3D geometry and realistic physical properties 
\citep[see][for details]{Duwiquetetal2022} the approach consists first
in generating a numerical model without tectonic regime application
(the benchmark experiment,  Figure~\ref{fig10}a). Thereafter, we
applied two different tectonic regimes  (Figure~\ref{fig10}b,
compression and extension,  Figure~\ref{fig10}c,d). A
comparison between the benchmark experiment and the other simulations
will be discussed, especially in terms of fluid pressure variation 
(Figure~\ref{fig10}e). This comparison allows 
us to\vadjust{\pagebreak} 
understand the
different dynamic\unskip\break 
convective patterns 
obtained in the different
numerical experiments, and how the poroelasticity-driven force could
explain these differences. Note that this experiment is similar to the
ones discussed in  \citet{Duwiquetetal2022} except that fault thickness
is here reduced to 0.3~km.

To investigate the influence of tectonic regimes on fluid flow, we used
the Andersonian assumption, which is regularly used in geomechanical
reservoir studies  \citep{Zobacketal2003}. The principal stresses are
expressed with vertical $(S_{v})$, maximum horizontal
$(S_{H\mathrm{max}})$, and minimum horizontal
$(S_{h\mathrm{min}})$\unskip\break
components. Note that the relative stress magnitudes determine the
modeled tectonic regime: 

\begin{itemize}
\item Compressional (reverse/thrust faulting), with 
$S_{H\mathrm{max}}\geq S_{h\mathrm{min}}\geq S_{V}$,
\item Extensional (normal faulting) with 
$S_V\geq S_{H\mathrm{max}}\geq S_{h\mathrm{min}}$.
\end{itemize}

As the model is aligned with the principal stresses, pure normal
stresses are applied on the lateral boundaries. For the compressional
tectonic regime, the fault is perpendicular to the $S_{H\mathrm{max}}$
stress while for the extensional tectonic regime, the fault is
perpendicular to  $S_{h\mathrm{min}}$. In these models no shear is
applied on the vertical boundaries. The depth-dependence of stresses is
extracted from drilling measurements in the French Massif Central 
\citep{CornetBurlet1992} and is illustrated in  Figure~\ref{fig10}b.
This implementation follows the method detailed in
\citet{Duwiquetetal2022}.

In the benchmark experiment (i.e. without tectonic regime application),
the fluid shows an upward movement at the center of the fault zone 
(Figure~\ref{fig10}a) and two downward movements on either side of the
fault zone. Fluid flow velocities range from 2 to 
$17\times10^{-9}~\mathrm{m}{\cdot}\mathrm{s}^{-1}$. The velocity is
highest in the downward movements. The upward movements induce a rise
of the isotherms. For example, the 150~\textdegree C isotherm target
for geothermal energy is found at 1.2~km depth  (Figure~\ref{fig10}a),
while in a pure conductive regime with no fluid flow, the 
150~\textdegree C isotherm was observed at 5~km depth. Between the
basement and the fault zone, no fluid pressure variation exists 
(Figure~\ref{fig10}e). However, within the fault, there is a slight
increase in fluid pressure, which is related to upward fluid movement.
In light of these elements, buoyancy-driven forces alone act on fluid
convection. After tectonic regime application  (Figure~\ref{fig10}b),
the results are different  (Figure~\ref{fig10}c,d).

In the compressional tectonic regime (Figure~\ref{fig10}c), the fluid
pattern is characterized by two downward and two upward movements. The
fluid velocities vary from 1 to 
$13\times10^{-9}~\mathrm{m}{\cdot}\mathrm{s}^{-1}$,  and the fastest
velocities are found in the downward movements. The upward movements
raise the  150~\textdegree C isotherm at a depth of 1.4~km. This
convective dynamics differs from the benchmark experiment, 
(Figure~\ref{fig10}a), as well as from the extensional tectonic regime.
For the extensional tectonic regime  (Figure~\ref{fig10}d), the fluid
patterns follow a downward and two upward movements. Fluid flow
velocities range from 1 to 
$20\times10^{-9}~\mathrm{m}{\cdot}\mathrm{s}^{-1}$, and as in the
previous cases, the downward movements focus the fastest fluid flows.
The temperature rise follows the upward movement and is focused on
either side of the fault, and the  150~\textdegree C isotherm reaches a
depth of 1.7~km.

The three numerical experiments show three different convective
dynamics and three spatial isotherms distributions. Figure~\ref{fig10}e
shows the fluid pressure variations for all experiments. As previously
described, for the benchmark experiment, only a small variation is
noticed within the fault zone. After stress application, for the
compressional tectonic regime, a lateral fluid pressure variation is
present  (Figure~\ref{fig10}e), and this trend is found above the
values of the benchmark experiment. The difference in fluid pressure
between the basement and the fault is  1.45~MPa. For the extensional
system, a variation between the basement and the fault is also observed
and the trend of this variation is below the benchmark experiment. The
difference in fluid pressure between the basement and the fault is 
0.95~MPa. In the benchmark experiment, the only force impacting the
fluid flow is the buoyancy force, whereas the implementation of a
tectonic regime changes the way the fluids circulate. Note that the
pressure difference generated between the basement and the fault guides
the fluids from the high-pressure zones to the low-pressure zones. This
force has a similar impact as the effects of topography on fluid flow 
\citep{ForsterSmith1989,LopezSmith1995}. This effect is facilitated by
higher permeability values in the fault than in the basement. 

\section{Discussion and perspectives}\label{sec6}

\subsection{From large-scale mantle heat flow anomalies to fault-scale
thermal anomalies}\label{sec6.1}
Surface manifestations or areas of anomalously high surface heat flow
have long guided geothermal exploration. It seems that favorable areas
are located in regions where surface heat flow exceeds 
$80~\mathrm{mW}{\cdot}\mathrm{m}^{-2}$  [as in western North America, 
\citealp{WisianBlackwell2004}], a value probably correlated with
anomalously high mantle heat flow. Yet, the anomalous mantle heat flow
is not necessarily transferred homogeneously at the surface, and the
crustal heterogeneities may disturb and redistribute anomalous
temperatures in the shallow crust. In particular, permeable zones of
the crust contain crustal fluids that can circulate, by
pressure-driven, poroelasticity-driven, or buoyancy-driven forces, and
transfer deep hot fluids to shallow depths. 

In Western Europe, the anomalous large-scale surface heat flow 
(Figure~\ref{fig3}b) is probably at the origin of the European Cenozoic
Rifts System  [ECRIS,  \citealp{Dezesetal2004}] but associated
geothermal systems are sparsely distributed, probably because
permeability along associated fault-damage zones is heterogeneous. It
may be worth investigating all other ECRIS-related grabens, as it was
partly done in the past 
\citep{Guillou-Frottieretal2010,Calcagnoetal2014,Freymarketal2017}.
Hence, geodynamic settings and the history of mantle events in the last
30~Myrs have to be taken into account in exploration strategies, in
particular for the blind geothermal systems. 

The most permeable zones correspond to crustal fault zones, where the
thickness of the damaged zones may reach hundreds of meters. Indeed,
our results on the role of fault thickness on the critical Rayleigh
number  (Figure~\ref{fig7}) imply that wide crustal fault zones may
host geothermal systems at shallow depths. As shown in the
Appendix~\ref{appA} and in Figure~\ref{fig8}a, temperatures above 
150~\textdegree C may be present at a depth of less than 3~km and over
a width of several hundreds of meters  (Figure~\ref{figA1}b). Wide
fault zones and their surroundings may thus correspond to the most
promising geological targets hosting exploitable geothermal reservoirs
(given that the flow rate is sufficiently high), without using
extensive hydraulic fracturing. However, uncertainties in fault rock
permeability must be considered since several physical and chemical
processes may control the efficiency of fluid circulation (see 
Section~\ref{sec4}).

\subsection{Fluid flow processes within and around fault
zones}\label{sec6.2}
The numerical models performed here   (Figures~\ref{fig8}--\ref{fig10})
consider permeability as a static and spatially variable parameter
based on a conceptual model of the ``uni core'' fault zone, deduced
from field observations on the Punchbowl Fault zone (USA) 
\citep{ChesterLogan1986}. For larger systems, such as the Carboneras
Fault Zone (Spain)  \citep{Faulkneretal2010} or the Pontgibaud Fault
Zone (France)  \citep{Duwiquet2022},  the spatial variation of
permeability corresponds to the conceptual model of the ``multi-core''
fault zone  \citep{Faulkneretal2003}, where permeability can be
represented as spatially heterogeneous  \citep{Duwiquetetal2019}. Here,
our new numerical results shed light on other important factors that
may also control the geothermal potential. In particular, fault dip,
topography, and tectonic regime affect the underlying thermal anomaly
features  (Figures~\ref{fig9} and~\ref{fig10}).

\subsubsection{Effects of topography and tectonic
regimes}\label{sec6.2.1} 
The effects of topography and tectonic regimes on fluid flow are
similar. According to Darcy's law, the behavior of fluid flow is
related to the pore fluid pressure gradient; pore fluid always flows
from high to low pressure. Pressure gradients are generated from
topographic highs to topographic lows, and induce meteoric fluid flow
from peaks to valleys  \citep{ForsterSmith1989,Tailleferetal2018}. The
poroelasticity-driven force is related to the tectonic regime itself 
\citep{Duwiquetetal2022}, and, as for the topography-driven flow,
fluids are guided from high pressure (the basement) to low pressure
zones (the fault)  (Figure~\ref{fig10}e). In summary, when topography
and tectonic regimes are considered, fluid flow is not only driven by
buoyancy forces, but also by pressure-driven and poroelasticity-driven
forces. Thermal consequences of the presence of these three different
forces may be subtle. 

For example, there may be a competition between the infiltration of
cold fluids in the subsurface damage zone and the upwelling of warm
fluids from deep levels. Downward circulation of superficial cold
fluids is facilitated when the damage zone is particularly permeable,
or when topography is important, explaining why large damage zones
could be a disadvantage in some configurations of  Figure~\ref{fig9}c.
For example, the case of a strike-slip fault zone where SRE $=$ 3000~m
shows a decreasing temperature when fault zone thickness increases,
contrary to what theoretical considerations (with no topography)
suggest  (Figure~\ref{fig7}). This infiltration effect also explains
the asymmetry of the temperature anomaly observed in Figure~\ref{fig9}
as a function of the position of the fault-related relief. This should
be considered in three dimensions because lateral fluid flow also
occurs along the fault  \citep[see][]{Tailleferetal2018}. However, the
infiltration effect caused by the unfavorable conjunction of
permeability and topography decreases with depth because the buoyancy
effect on warm fluids takes over. We note that 
\citet{ForsterSmith1989} and  \citet{LopezSmith1995} already emphasized
these complex interactions between topography-driven fluid flow in
relatively permeable host rocks and buoyancy-driven fluid flow in 
fault\unskip\break zones.

\subsubsection{Effect of fault dip}\label{sec6.2.2}
We find that the vertical dip (90\textdegree) concentrates a higher
temperature anomaly than the  30\textdegree\ and 60\textdegree\ dips 
(Figure~\ref{fig9}b). These results were also found for other 
model\vadjust{\pagebreak}
configurations where the topography\unskip\break 
was not taken into account 
\citep{Duwiquetetal2019,Guillou-Frottieretal2020}. According to 
Figure~\ref{fig9}b, temperatures above  160~\textdegree C would be
present at a depth of 4~km for the three dips, and temperatures above 
150~\textdegree C would extend laterally (outside the fault zone) over
several kilometers. In  Figure~\ref{fig10}, such temperatures are
reached within the fault zone between 1.4 and 1.7~km, depending on the
tectonic regime. These shallower depths can be due to the absence of
topography, to the vigorous convection not present in 
Figure~\ref{fig9}, or to both phenomena. Indeed, the results from 
Figure~\ref{fig8} (2D models without topography),  Figure~\ref{fig9}
(2D models with topography and dip variation) and  Figure~\ref{fig10}
(3D models without topography) are difficult to compare, and
permeability distributions within the host rocks are not the same.
However, physical processes controlling fluid flow can be discussed.

\subsection{High-temperature geothermal systems in fault
zones}\label{sec6.3}
Although complex interactions between different fluid flow processes
may exist at exploitable depths (1--4~km), we have outlined some key
features of fault zones enabling us to consider them as potential hosts
of high-temperature geothermal systems. The width of the damage zone
must be sufficiently large to allow fluid convection to occur easily.
If the fault zone is surrounded by an important relief, then
large-scale cold downwellings could prevent hot temperatures to reach
shallow depths. Subvertical fault zones are also more favorable for
allowing important thermal anomalies at shallow depths.

\begin{figure}
\includegraphics{fig11}
\caption{\label{fig11}Maximum displacement--length diagram for faults
using data from several localities and fault settings 
\citep{Fossen2016}. The length of fault associated with displacement in
the range of  100--150~m associated with a thickness of damage zone of
100~m or more (see Figure~\ref{fig6}) is highlighted in light grey.
Areas in light orange correspond to the faults where heat conduction is
dominant.}
\end{figure}

Consequently, and following the main ideas of this study, we can define
potential targets for geothermal exploration using the main geometrical
parameters of fault zones. We focus here our analysis on thick fault
zones where free convection is possible. We do not include geothermal
plays where forced convection may lead to significant thermal anomalies
[e.g., Grimsel Pass, Swiss Alps,  \citealp{Wanneretal2019}]. The
conjunction of free and forced convection within permeable fault zones
in mountainous regions could also enhance the number of geothermal
targets. The difficulty is that permeability data are even less
collected in mountainous regions than in flat ones. This potential
should be examinated in a future study. 

\begin{figure*}
\includegraphics{fig12}
\caption{\label{fig12}Example of a normal fault with a damage zone
thickness of 600~m, where permeability ranges between $10^{-14}$ and 
$10^{-13}~\mathrm{m}^{2}$. The other geometrical parameters necessary
to generate convection cells and high temperatures at shallow depths
are indicated  \citep[modified after][]{Fossen2016}.}
\end{figure*}

Damage zone width is essentially controlled by the displacement on the
fault whatever
the fault kinematics  [\citealp{Benseetal2013},
Section~\ref{sec4.3}]. We must consider fault displacement greater or
equal to  100--150~m to reach a damage zone thickness of a hundred
meters or more. This critical thickness is necessary to reasonably
consider the generation of convection cells  (Figure~\ref{fig6}). A
power law defines the relationship between the maximum displacement on
the fault and the fault length  \citep{SchultzFossen2002}. This law
shows that a fault displacement of  100--150~m is associated with
faults lengths between 1.5 and 10~km  (Figure~\ref{fig11}). The
estimation of these parameters is based on a mean permeability in the
range of  $10^{-14}$--$10^{-13}~\mathrm{m}^{2}$, therefore damage zone
permeabilty changes can influence the geometrical parameters necessary
to generate fluid convection in a fault zone. The main geometrical
parameters for faults able to host geothermal systems are:  (i)~a
damage zone thickness $\geq$ 100~m, (ii)~a minimum\unskip\break 
displacement of
100--150~m and (iii)~a kilometric scale fault. These parameters are
illustrated in  Figure~\ref{fig12}, where a damage zone thickness of
600~m is chosen.

The database of  \citet{Scibek2020} contains numerous information on
the permeability of fault zones around the world, associated with fault
zone parameters available: protolith, kinematic, length, displacement,
damage zone thickness. Despite the important number of referenced sites
(511 sites), it remains incomplete and spatially distributed in a
heterogeneous way. We have sorted the database according to the  3 main
fault geometrical criteria to select the sites that can be potential
targets for geothermal exploration. The results highlight that one
third of the sites (173 sites) have favorable fault geometrical
settings for geothermal exploration. 

Table~\ref{tab1} shows a selection of a few fault zones, whose length
is larger than 1~km, whose displacement is larger than 100~m, and where
maximum fault thickness is larger than 150~m. Data belonging to the
``Geothermal reservoir'' category illustrates some well-known
geothermal systems providing several tens of MW (Dixie Valley in
Nevada; Kakkonda, Japan; Bal\c{c}ova, Turkey). In the lower half of 
Table~\ref{tab1}, other fault zone categories show cases where
permeability is probably sufficiently large because the fault width
exceeds several hundreds of meters (Gryphon, Saskatchewan; Eklutna,
Alaska). Some systems where permeability is unknown are also indicated
since their fault widths approach or even exceed 1~km (Sandwich,
Illinois; Carboneras, Spain).

\begin{table*}
\caption{\label{tab1}Selection of a few fault zones from the database
of \citet{Scibek2020}, where fault length is greater than 1000~m, fault
displacement greater than 100~m, and fault thickness greater than
150~m\vspace*{-2pt}}
\fontsize{9}{10}\selectfont\tabcolsep=3pt
\begin{tabular}{cccccccc}  
\thead
\raisebox{1pt}{\morerows{1}{\parbox[t]{1.3cm}{\centering Site number}}} & 
\xmorerows{1}{Short name} & \xmorerows{1}{Country/State} & 
\xmorerows{1}{Category} & 
\raisebox{1pt}{\morerows{1}{\parbox[t]{1.8cm}{\centering Max. thickness (m)}}} & 
\multicolumn{2}{c}{Permeability (m$^{2}$)} &
\raisebox{1pt}{\morerows{1}{\parbox[t]{2.5cm}{\centering Power production (MWe)}}}
\\\cline{6-7}
{} & {} & {} & {} & {} & {DZ} & {Bulk} & {} \\ 
\endthead
167 & Chingshui & Taiwan & Geothermal reservoir & 260 & / & 2.5 $\times$ 10$^{-14}$ & 4.2 \\ 
180 & Bal\c{c}ova & Turkey & Geothermal reservoir & 1000 & / & 4 $\times$ 10$^{-13}$ & 72 \\ 
139,1 & Ogiri & Japan & Geothermal reservoir & 150 & / & 6 $\times$ 10$^{-13}$ & 30 \\ 
154 & Los Azufres & Mexico & Geothermal reservoir & 1000 & / & 4 $\times$ 10$^{-13}$ & 161.5 \\ 
176 & St Gallen & Switzerland & Geothermal reservoir & 1000 & / & 5 $\times$ 10$^{-14}$ & Abandoned \\ 
186 & Dixie & Nevada & Geothermal reservoir & 1000 & / & 10$^{-12}$ & 70.9 \\ 
141,1 & Kakkonda & Japan & Geothermal reservoir & 1000 & / & 2 $\times$ 10$^{-14}$ & 50 \\ 
328 & Gryphon & Saskatchewan & Engineering & ${>}$500 & 10$^{-15}$ & / &  \\ 
370 & Eklutna & Alaska & Engineering & 1000 & / & 3 $\times$ 10$^{-12}$ &  \\ 
265,1 & Sandwich & Illinois & Groundwater supply & 800 & / & / &  \\ 
247 & Gloucester & Ontario & Groundwater supply & 300 & / & 2 $\times$ 10$^{-9}$ &  \\ 
33,1 & Neodani & Japan & Active fault/processes & 200 & / & 2.3 $\times$ 10$^{-13}$ &  \\ 
51 & Lansj\"{a}rv & Sweden & Active fault/processes & 150 & / & 10$^{-13}$ &  \\ 
40,1 & Alpine & New Zealand & Active fault/processes & 100~s & / & 6 $\times$ 10$^{-14}$ &  \\ 
47 & Carboneras & Spain & Active fault/processes & 2000 & / & / &  \\ 
55 & San Andreas & California & Active fault/processes & 500 & / & 2 $\times$ 10$^{-14}$ &
\botline
\end{tabular}
\tabnote{The site number comes from the database. DZ $=$ Damage zone.
``Bulk'' refers to data that cannot be relocated in fault core or
damaged zones. The first 7 lines correspond to geothermal systems for
which power production is indicated.\vspace*{-4pt}}
\end{table*}

\subsection{Potential geothermal systems in European fault
zones}\label{sec6.4} 
At the scale of Europe, fault geometrical parameters applied to the
databse of  \citet{Scibek2020} highlight 40 fault zones with a
favorable geometry for hosting geothermal systems. We also added three
crustal fault zones that have been previously mentioned in this study
(i.e., the Pontgibaud and T\^{e}t fault systems in France, and
the fault zone of Larderello in Italy). Fault kinematics are
represented in  Figure~\ref{fig13}a. The maps in  Figure~\ref{fig13}
show the spatial distribution of geothermal power plants and the
selected fault zones with their kinematics  (Figure~\ref{fig13}a) and
with their associated permeabilities  (Figure~\ref{fig13}b). 

\begin{figure}
\includegraphics{fig13}
\caption{\label{fig13}(a)~Geothermal fields installed and exploited in
Europe (star symbols)  (\url{https://www.thinkgeoenergy.com/map/}).
Selected fault zones from the database of  \citet{Scibek2020} with
geometrical criteria (length greater than 1000~m, displacements or
damage zone thickness greater than 100~m) favorable for geothermal
exploration are also indicated. Note the different types of faults. 
(b)~Map showing the main faults permeabilities from the database of 
\citet{Scibek2020} and geothermal systems. Permeabilities for the
Larderello geothermal system and Pontgibaud system are from 
\citet{BertaniCappetti1995} and  \citet{Duwiquetetal2021}
respectively.}
\end{figure}

Some countries do not contain any new data (Germany, Poland) but this
is only due to the lack of geometrical data in the database for these
sites. For example, the Unterhaching fault zone  [Germany, site number
101.1  in  \citealp{Scibek2020}] would have a high bulk permeability of
$10^{-13}~\mathrm{m}^{2}$, but there is no information on the fault
length, displacement or the damage zone thickness. Consequently, these
maps must be considered as a first step to define new potential targets
for geothermal exploration.

The maps in Figure~\ref{fig13} highlight that only a few fault zones
currently host exploited geothermal reservoirs, but some of them
control the location of several power plants. This is the case for the
Menderes Massif (Turkey) where at least 45 power plants are located
close to the low-angle normal fault systems with a total capacity of
1450~MWe. The identified new sites tend to show some similarities even
though they are located in different tectonic contexts. Further, these
maps show that all tectonic contexts may be interesting, in particular
the extensional regime represented by normal faults. However, this may
also be due to a bias in the incomplete database. Our selection, based
on purely geometrical criteria, shows that fault permeabilities range
from  $10^{-15}$ to  $10^{-12}~\mathrm{m}^{2}$, a favorable range for
high-temperature geothermal systems. Based on these observations and
criteria, we suggest that the geothermal potential in Europe is still
underestimated.

\subsection{Perspectives for geothermal exploration in fault
zones}\label{sec6.5} 
Beside being heterogeneous, the  \citet{Scibek2020} database lacks
numerous permeability data of European fault zones. Consequently, 
Figure~\ref{fig13} should not be considered as a predictivity map since
several additional features must be taken into account before being
usable to predict the presence of an exploitable geothermal system.
Indeed, information on fault zone geometry and permeability data should
be combined with the regional tectonic regime, the local stress field
and the surrounding topography. In addition, it is important to recall
that the permeability database of fault zones in Europe is not
exhaustive, (since in constant evolution). We however believe that our
study can help geothermal exploration towards a more efficient
targeting, where permeability enhancement techniques would not be
necessary.

We focused our study on thermal anomalies greater than  150~\textdegree
C at a depth lower than 4~km. It must be noted that the recent study by
\citet{Lundetal2022} indicates that ``\textit{electric power from
geothermal energy is now being produced from resources with
temperatures as low as  90~\textdegree C, using the organic Rankine
cycle process in binary power units}''. Consequently, our selection
must be considered as a lower limit since temperatures around 
100~\textdegree C at a depth of 1--2~km should be easily found in wide
fault zones  (Figure~\ref{fig8}b), thus increasing the number of fault
zones that could be considered as promising targets for geothermal
exploration.

\section*{Conflicts of interest}
Authors have no conflict of interest to declare.

\section*{Acknowledgements}
This research has been supported by the GERESFAULT project (grant no.
ANR-19-CE05-0043-01), GM acknowledges funding from the program TelluS
of the Institut National des Sciences de l'Univers, CNRS. Discussions
with Patrick Ledru and Vincent Bouchot helped to improve the ideas
developed in this study. We thank reviewers Domenico Montanari and
Christoph Wanner for their positive and constructive reviews. We also
thank Olivier Fabbri for additional suggestions that help to improve
the final version of the manuscript.

\CDRGrant[ANR]{ANR-19-CE05-0043-01}

\back{}

\appendix{}

\gdef\appendixlabel{Appendix}

\setcounter{equation}{0}
\renewcommand{\theequation}{\thesection.\arabic{equation}}

\setcounter{table}{0}
\renewcommand\thetable{A\arabic{table}}

\setcounter{figure}{0}
\def\thefigure{\Alph{section}\arabic{figure}}

\section{Numerical modeling: model set-up, equations and physical
parameters}\label{appA}
All numerical models of this study have been computed with the Comsol
Multiphysics{\texttrademark} software. Several different benchmark
experiments have been performed, and previously published results were
reproduced by our numerical procedure 
\citep{Eldursietal2009,Garibaldietal2010,Tailleferetal2017,Guillou-Frottieretal2020,Duwiquetetal2022}.
Figures~\ref{fig8},~\ref{fig9} and~\ref{fig10} have different model
setups and their peculiarities are described below.

\begin{figure*}
\includegraphics{fig14}
\vspace*{4pt}
\caption{\label{figA1}Details on the numerical modeling leading to the
results illustrated in Figure~\ref{fig8}. For these models, a fixed
temperature condition is imposed at the surface. (a)~Model set-up,
parameters, and boundary conditions; (b) temperature field (isotherms in
black contours, separated by 25~\textdegree C) for a fault width of 800 m
and a fault permeability of 
$5 \times 10^{-14}$~m$^{2}$, 3700 yr
after the beginning of the experiment; (c) horizontal temperature
profile at a depth of 3~km (profile A--B in b), at the time when
temperature perturbation is maximum.}
\vspace*{6pt}
\end{figure*}

\subsection{Models of Figure~\myRef{fig8}}\label{appA1}
The series of numerical models reported in Figure~\ref{fig8} have been
performed in 2D, where a permeable fault zone of variable width and
variable permeability is embedded into a low permeable host rock. Model
set-up, geometry and boundary conditions are illustrated in Figure~\ref{figA1}.
Physical properties and other details are also indicated. The heat
equation, Darcy law and mass conservation are coupled through the
velocity field and the temperature-dependence of fluid properties. An
identical numerical procedure as that described in
\citet{Guillou-Frottieretal2020} was used.

\subsection{Models of Figure~\myRef{fig9}}\label{appA2}
The models of Figure~\ref{fig9}, partly dedicated to the role of
topography, account for the altitude-dependence of the surface pressure
(see expression $p_0$ and $p$ in Table~\ref{tabA1}). A mixed thermal
boundary condition is imposed at the surface (details below), allowing
to get cold temperatures (5--10~\textdegree C) at high altitudes, hot
temperatures above the fault zone (thermal springs up to 80~\textdegree
C) and 10~\textdegree C elsewhere  \citep[see details
in][]{Tailleferetal2018}. Fluid and rock properties, tested parameters,
initial conditions, and boundary conditions are detailed in
Table~\ref{tabA1}. Additional results are shown in Figure~\ref{figA2}.

\begin{figure*}
\includegraphics{fig15}
\caption{\label{figA2}Additional results with horizontal temperature
profiles at a depth of 2000~m.}
\end{figure*}

\begin{table*}
\caption{\label{tabA1}Parameters and variables used in numerical models
of Figure~\ref{fig9}}
\begin{tabular}{lclc}  
\tbody
\multicolumn{4}{c}{Fluid variables and parameters} \\ 
Fluid dynamic viscosity  & ${\mu}$ & 
$2.414 \times 10^{-5}
\exp(570/((T-273.15)+133))$, $T$ in K & Pa${\cdot}$s \\ 
Fluid density & ${\rho}$f & \parbox[t]{7cm}{1036.5 $-$ 0.14167 ${\times}$ 
($T-273.15$)~$-$~0.0022381 ${\times}$ ($T-273.15$) ${\times}$ ($T-273.15$), 
$T$ in K} & kg/m$^{3}$\vspace*{2pt} \\ 
Thermal capacity & Cpf & 4180 & J/(kg${\cdot}$K) \\ 
Thermal conductivity & ${\lambda}$f & 0.60 & W/(m${\cdot}$K) \\ 
Thermal expansion coefficient & ${\alpha}$ & 10$^{-4}$ & 1/K\vspace*{7pt} \\ 
\multicolumn{4}{c}{Rock variables and parameters} \\ 
Mass density & ${\rho}$s & 2650 & kg/m$^{3}$ \\ 
Porosity & ${\phi}$ & 0.10 &  \\ 
Thermal capacity & Cps & 1000 & K/(kg${\cdot}$K) \\ 
Thermal conductivity & ${\lambda}$s & 2.50 & W/(m${\cdot}$K) \\ 
Thermal expansion coefficient & ${\beta}$ & 
$7\times10^{-6}$ & 1/K \\ 
Basement permeablity & Kb &  10$^{-16}\exp^{(y-\mathrm{SRE})/2500}$  $y$ in m, $-7000 < y < \mathrm{SRE}$ & m$^{2}$ \\ 
Core zone permeability  & $K_c$ & 10$^{-20}$ & m$^{2}$ \\ 
Damage zone permeability  & $K_d$ &  
$2\times10^{-14}$ & m$^{2}$\vspace*{7pt} \\ 
\multicolumn{4}{c}{Geometric features (tested parameters)} \\ 
Scarp relief elevation & SRE & 200, 1000, 3000 & m \\ 
Fault dip  & FD & 30, 90, 60 & deg \\ 
\parbox[t]{4.5cm}{Total Damage zone Thickness ($2 \times \frac{1}{2}$ DZ)} & DZT & 50, 300, 1000 & m\vspace*{7pt} \\ 
\multicolumn{4}{c}{Model variables and parameters} \\ 
Ambiant temperature & Text & 20 & \textdegree C \\ 
Heat flow at the model base & Q$_0$ & 0.07 & W/m$^{2}$ \\ 
\parbox[t]{4.5cm}{Mixed thermal condition at the model surface} & qs & 0.05 ${\times}$ (Text-$T$),
$T$ in~\textdegree C & W/m$^{2}$\vspace*{2pt} \\ 
Initial pressure field & 
$p$ & \parbox[t]{7.3cm}{10$^{5}$ ${\times}$ (($1-0.006$~$y$/288.15)$^{5.255}$) $+$ ($1000\times 9.8\times (-y)$), $y$ in m} & Pa \vspace*{2pt}\\ 
Surface pressure & $p_{0}$ & 10$^{5}$ ${\times}$ (($1-0.006$~$y$/288.15)$^{5.255}$), $y$ in m & Pa \\ 
Initial temperature field & 
$T$ & Text $+$ 0.03 ${\times}$ ($-y$), $y$ in m & \textdegree C \\ 
Thermal insulation  & \multicolumn{3}{c}{Lateral boundaries} \\ 
No flow conditon & \multicolumn{3}{c}{Lateral $+$ model base boundaries}
\botline
\end{tabular}
\end{table*}

{Details on the mixed thermal boundary condition fixed at the
topographic surface of the numerical model, issued from
\citet{Tailleferetal2018}.}

A theoretical constant heat flow 
$q_0 = 0.07~\mathrm{W}{\cdot}\mathrm{m}^{-2}$ is imposed at
the base of the model. The lateral limits of the model are thermally
insulated. In order to allow for temperature variations at the surface
(cold areas at high altitudes but also warm zones where hot springs
emerge), a mixed thermal boundary condition like the one used in 
\citet{Magrietal2015} and  \citet{Tailleferetal2018} is applied to the
topographic surface:
{\begin{equation}
Q = h (T_0-T)
\end{equation}}\unskip\noindent
where $Q$ (W${\cdot}$m$^{-2}$) is the surface heat flow, $h$
(W${\cdot}$m$^{-2}{\cdot}$K$^{-1}$) a coefficient of heat transfer and $T_0$
(\textdegree C) the reference temperature at the surface (chosen to be
20~\textdegree C). In order to impose a realistic thermal boundary
condition, a sensitivity analysis of the heat transfer coefficient
$h$ was performed. The condition applied to the topographic
surface needs to satisfy (1)~a realistic temperature at the surface
depending on the temperature (5--20~\textdegree C) and (2)~potential
convective flows allowing hot fluid emergence. The $h$ value was
therefore calibrated so as to combine realistic temperatures for
(i) the fluid infiltrated at the model surface and (ii) the hot fluids
that rise up to the model surface. To obtain consistent temperatures at
the surface (between 5 and 20~\textdegree C) while emergence
temperatures may reach 70~\textdegree C (which are observed values for
hot springs in mountain ranges), the optimum $h$ value appears
to range between 0.017 and 0.025~W${\cdot}$m$^{-2}{\cdot}$K$^{-1}$. 

Using this value, the few areas at high altitude (above 2000 m)
exhibit low temperatures (5--7~\textdegree C) while all other areas
show surface temperature of 10 ${\pm}$ 2~\textdegree C. The initial
conditions for the model are a hydrostatic pressure and a purely
conductive thermal state, from which a steady state is computed.

\subsection{Models of Figure~\myRef{fig10}}\label{appA3}
The numerical procedure for the models of Figure~\ref{fig10}, as well
as the benchmark for the 3D experiment, have been described in 
\citet{Duwiquetetal2022}.

\def\ensuremath#1{#1}

\bibliographystyle{crgeos}
\bibliography{crgeos20221173}
\refinput{crgeos20221173-reference.tex}

\end{document}
