\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage[round]{natbib}
\let\cite\citep




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[ngerman,greek,english]{babel}










\begin{document}

\title{Effects of pre-existing~ structures on the seismicity of the Charlevoix
Seismic Zone~~~~~~~~~~~~}



\author[1]{ifadugba}%
\affil[1]{University of Memphis}%


\vspace{-1em}



  
  \date{May 07, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup








\sloppy


\section*{Abstract}

{\label{619939}}

\textbf{{[}Fill in.{]}}

\par\null

\section*{Introduction}

{\label{619939}}

The Charlevoix Seismic Zone (CSZ) is the most seismically active region
in eastern Canada. CSZ is within the Cambro-Ordovician St. Lawrence rift
zone in southeastern Quebec and overlapped by a Devonian impact crater
(Fig.~{\ref{749307}}). The impact crater has a 28
km-surface radius and is 15 km deep. The crater has a more damaged
interior zone with a surface radius of 18 km
(Fig.~{\ref{749307}}). The impact crater is
superimposed on three major rift faults trending N35\selectlanguage{ngerman}°E and dipping at
70° to the southeast (\citealp*{anglin1984};~\citealt{Rondot_1971}). The CSZ
is considered to pose a~high risk of seismic hazard due to its history
of generating large earthquakes
(\citealp*{Anglin1984};~\citealp*{Lamontagne1999}).~Earthquakes in the CSZ occur
within the volume bounded by the rift faults. Those of greater
magnitudes occur in the northeast from the crater
(Fig.~{\ref{749307}}). The earthquakes within and
beneath the crater form the background seismicity (Nuttli magnitude,
mN~\selectlanguage{english}[?] 3) of the CSZ~\citep{Baird2010}.

\par\null\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/vel-model-CSZ/Map-of-CSV-now3}
\caption{{Topography and seismicity of the Charlevoix Seismic Zone (CSZ) as well
as the location of the impact crater (outer cyan circle) and the more
damaged inner crater (inner cyan circle). Small circles are the
relocated epicenters~\protect\citep*{Powell2017} and their colors represent the
focal depths. The focal~mechanisms are for the earthquakes used
by~\protect\citet*{Mazzotti2010}~for the stress inversion. Large circles labeled as
NW and SE show orientations of SH\textsubscript{max}~from the stress
inversion of focal mechanism (SH\_S) and from borehole breakout
measurements (SH\_B) for the earthquake clusters in the northwest and
southeast of the crater center. Solid black lines mark the rift faults
known in the region: GNF, Gouffre Northwest Fault; SLF, Saint-Laurent
fault; CHF, Charlevoix Fault; and SSF, South Shore Fault
(\protect\citealp*{Rondot1971},~\protect\citealp*{Lamontagne1999}). The inset shows the location
of the CSZ in eastern Canada.~
{\label{749307}}%
}}
\end{center}
\end{figure}

\citet{Vlahovic_2003} conduct a three-dimensional (3D) P-wave velocity
study of the CSZ. The tomography models show several high-velocity
bodies at mid-crustal depths. A significant distribution of earthquakes
in the CSZ are~located around the high-velocity bodies~and~some larger
earthquakes along the northern edges of the high-velocity
bodies.~\citet{Vlahovic2003} suggest that the high-velocity bodies cause
stress amplification which can explain the local stress rotation. A
recent earthquake relocation and tomography study
of~\citet*{Powell_2017}~give a better constrain to the geometries of these
high-velocity bodies using a larger number of earthquakes in the CSZ. In
addition to the high-velocity bodies,~\citet*{Powell2017}~show~that the
relocated earthquakes occur along the rift faults and thus defined the
geometry of the rift faults. They also conclude that the dip of the
northernmost rift fault is 65\selectlanguage{ngerman}°SE while the other two faults dip 40°SE.
This is contrary to previous studies in the CSZ which believe that the
earthquakes occur within highly fractured volumes bounded by the rift
faults with a dip of 70°SE (\citealp{Baird2010};~\citealp*{Lamontagne1999}).~

\par\null

\citet*{Mazzotti2010} assess the stress state in the CSZ. They determine
the orientation of maximum principal horizontal stress
(SH\textsubscript{max}) from stress inversion of focal mechanisms in ten
seismic zones in central and eastern North America. The orientations of
SH\textsubscript{max} from focal mechanism stress inversion of the
entire earthquakes in the CSZ show a clockwise rotation of about 32°
relative to the SH\textsubscript{max} determined from borehole breakouts
(\citealp{Zoback1992};~\citealp*{Mazzotti2010}). Also, within the CSZ, there is
a local clockwise rotation of about 47° in the SH\textsubscript{max}
orientation between clusters of earthquakes in the SE compared to the
earthquakes in the NW part of the CSZ
(Fig.~{\ref{749307}},~\citealp*{Mazzotti2010}).~

\par\null

To explain the seismicity of the CSZ in terms of the interactions
between the impact structure and the rift faults,~\citet{Baird2010}
numerically model the partitioning and distribution of seismicity in the
CSZ. The geodynamic model is able to explain the distribution of the
earthquakes within the faulted volume bounded by the fault planes and
also suggests weak~rift faults with friction angle of 5°. The numerical
model is also able to explain the relatively higher number of
earthquakes in the NE part of the crater relative to the SW
part.~\citet{Baird2010}~provide an explanation to the rotation of the
observed stress orientation by plotting the principal stresses in the
seismogenic zone of CSZ and the slips on the rift faults on stereonets.
The former to represent the SH\textsubscript{max} orientation in the
entire~CSZ while the slips on the faults represent the
SH\textsubscript{max} orientation in the SW part of the CSZ. However,
the models could not explain the observed rotation of the stress
orientation using the principal stresses in the numerical model for the
NW and SE part of the CSZ.

\par\null

In this study, we present a new set of numerical models that include
more realistic relative spacings of the rift faults and correlate the
modeled stress distributions with the recently-relocated hypocenters.
The results of this research will be compared with the previous study on
the influences of rift faults on the regional seismicity
by~\citet{Baird2010}.~

~

\section*{MODEL SETUP}

{\label{698233}}

We use PyLith version 2.1.0
(\citealp{Aagaard_2013};~\citealp{Aagard2015};~\citealp{c2015}), an
open-source finite element code for modeling dynamic and quasi-static
tectonic deformation. PyLith is suitable to our models that should
incorporate three major rift faults and the damaged impact crater
because it can easily handle complicated fault geometries generated by
meshing software.~

\par\null

\subsection*{Model geometry}

{\label{900956}}

Our models are composed of three units, crust, crater, and rift faults,
reflecting the simplified geology of the region
(Fig.~{\ref{854521}}A).~ The crust is a 220 x 220 x 40
km box with the crater cut out. Embedded at the center of the crust, the
impact crater is a spherical bowl of 60 km-surface diameter with a
maximum depth of 15 km.

The three rift faults trend at N35°E across the crust except~for the
outermost 10 km-thick layer, which is required for numerical stability
at the fault edges. In other words, the rift faults are embedded in the
crust layer except~for the top surface edge that is exposed on the
surface. ~The three~faults correspond to the Gouffre River, St. Laurent,
and the Charlevoix faults (Fig.~{\ref{749307}}). The
locations of the rift faults relative to the crater are based on the
geologic maps by~\citet*{Lamontagne1999} and~\citet*{Rondot1971}. We create
two sets of models with different dip angles of the faults based on the
two known geometries of the rift faults. One set has a uniform dip of
70°SE for all three faults (\citealp*{Anglin1984};~\citealp{Baird2010})
while the other has 65°, 40° and 40° for the three planar faults, from
north to south, respectively (\citealp*{Powell2017}).~

\par\null

We conducted resolution tests in which the change in displacement and
stress solutions relative to a 1-km resolution model are monitored as a
function of element size. Displacement was almost constant for 2-km or
smaller element sizes. We also observe edge effects by comparing a 420 x
420 x 40 km model with the original 220 x 220 x 40 km model, finding
that the seismogenic parts around the crater show consistent
displacement and stress solutions for both domain sizes. Based on this
result, we choose the smaller domain for this study and set the element
size to be about 2 km.~

\par\null

\subsection*{Initial and boundary
conditions}

{\label{912602}}

Lithostatic stress is assumed initially. We include the gravity in order
to simulate the expected increase in differential stress
(\(\sigma_D\)=~\(\sigma_1-\sigma_3\), where~\(\sigma_{1}\)
and~~\(\sigma_{3}\)~are the greatest and least compressive principal
stresses) and decrease in the slip on the rift faults with depth.~Due to
the nonplanar geometry of the crater, computing initial lithostatic
equilibrium stresses is not trivial.~We first construct a model with the
density distribution from an adopted velocity model for the CSZ. The
initial stress is computed~\emph{a priori} assuming the static balance
in the vertical direction only. Any remaining deviations from the
perfect lithostatic equilibrium reflect in the displacement and stress
solutions. We use the updated stress distribution as the true initial
stress conditions for the lithostatic equilibrium in the model. We
verify that this approach guarantees numerically zero displacements at
zero simulation time as expected.~

\par\null

We apply velocity boundary conditions to sides perpendicular to x- and
y-axis (Fig.~{\ref{854521}}A) after rotating the model
domain such that the y-perpendicular sides are parallel to the direction
of the regional SH\textsubscript{max}~(\(\phi_{\max}\)) from the
borehole breakout measurements,~N55°E~\citep{Zoback1992}
(Fig.~{\ref{854521}}B). We apply a 1-m per year
compressive velocity in the x-axis direction. The magnitude of boundary
velocity does not represent any tectonic loading in the continental
interiors but has the sole purpose of increasing the~\(\sigma_D\)
up to~the Byerlee's law (\(\tau\) = 50 MPa +
0.6\(\sigma_n\),~\citealp*{Byerlee_1978}).
The~\(\sigma_3\)~corresponds to the lithostatic stress. If the
average density of crustal rock is 2500
kg/m\textsuperscript{3},~\(\sigma_3\) is 250 MPa at a depth of 10
km. Noting that few earthquakes occur at greater depths than 10 km
(\citealp{Baird2010};~\citealp*{Powell2017}), we increase the boundary
displacement until the \(\sigma_D\)reach 706 MPa at 10 km depth.
The bottom boundary is a free-slip surface and the top boundary is
traction-free.~

\par\null\par\null\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/model-geometry-now/model-geometry-now5-(1)}
\caption{{(A) Model domain with the crater (gray) and crust (red). The three rift
faults (black lines), dimensions, boundary conditions and~the total
amount of displacement are annotated. The outer 10 km-thick layer of
crust (green) is added to contain fault edges within the domain, a
requirement by PyLith. (B) The orientation of the model domain~(red
dotted box) and the associated coordinate axes (red arrows) relative to
the geographic reference frame. Cyan circles on the geographic map show
the inner and outer crater boundaries and block solid lines trace rift
faults (see Fig.~{\ref{749307}}).
{\label{854521}}%
}}
\end{center}
\end{figure}

\subsection*{Velocity models for material
parameters}

{\label{553867}}

We assume linear elasticity~for both crust and crater. To obtain elastic
moduli and density distributions, we consider three regional velocity
models and one recent tomography (\citealp*{Powell2017}). The regional
models are the one for the Saguenay region in Quebec
(\citealp{helmberger1990}), the 1D standardized halfspace velocity model of
eastern Canada (\citealp*{Lamontagne1999}), and the 1D velocity model
of~\citet*{Lamontagne1999} derived from seismic refraction studies
(Fig.~{\ref{398984}}).~\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/vel-model-CSZ4/layered-velocity-models}
\caption{{Velocity models used in this study. Red, blue and black lines are
velocity models for the Saguenay region (\protect\citealp{Somerville1990}),
a~standardized velocity of Eastern Canada (\protect\citealp*{Lamontagne1999}) and the
velocity model of CSZ (\protect\citealp*{Lamontagne1999}).
{\label{398984}}%
}}
\end{center}
\end{figure}

We also consider a non-uniform distribution of density and elastic
moduli inferred from the 3D Vp and Vs tomography by~\citet*{Powell2017},
of which Vp and Vs variations for the depth range from 12 to 14 km is
shown in Fig.~{\ref{740205}}. The tomography has a
resolution of 2 km. Since this tomography does not cover the entire
study area, we assume a halfspace velocity model (\citealp*{Lamontagne1999})
outside of the tomography's coverage. Smooth transition between the
tomography and the half-space model is achieved by linear
interpolation.~Densities (\(\rho\)~) are computed according to
the Gardner's principle (\citealp{Gardner_1974}):~\(\rho=0.31\ V_p^{0.25}\),
where~Vp is the P-wave velocity in m/s and density is given as
g/cm\textsuperscript{3}.~

\par\null

For simplicity, we assume that the crater has one half or one fourth of
the crust's bulk and shear moduli (e.g.,~\citealp{Baird2010}). Similarly,
the density of the crater is reduced by 10 \% from that of the velocity
models over the corresponding depth range. ~\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/TOMO-velocity-models/TOMO-velocity-models}
\caption{{Variation of Vp and Vs at the depths from 12 to 14 km from the 3D
seismic tomographic model by~\protect\citet*{Powell2017}.
{\label{740205}}%
}}
\end{center}
\end{figure}

\subsection*{Models with and without
faults}

{\label{504267}}

We treat the rift faults as cohesive-frictional planes, on which~slip
vectors are not prescribed but determined by the model. In order to
create relative motions across the rift fault planes, PyLith~decouples
the motion of the two sides of the faults by inserting cohesive cells on
the fault surfaces. We assign a coefficient of static friction
(\(\mu\)) and the cohesion (C) on the fault. We assume
that~\(\mu\) and C~are~uniform on the fault planes in each
model but model the consequences of \textbf{~}C changing from 0 to 30
MPa and~\(\mu\)~from 0.1 to 0.6.~

\par\null

The reference model has the model geometry described above but without
the rift faults. This is equivalent to a model with infinitely strong
faults and high cohesion. The initial stress state is close to the
initially assumed lithostatic stress after the~nonplanar correction due
to the crater. The maximum principal stress is parallel to the
applied~loading direction perpendicular to the x-axis of the model
geometry. The intermediate principal stress corresponds to the initial
stress in addition to the Poisson~effect in the y-direction due to the
applied loading while the~minimum principal stress corresponds to the
initial stress in the vertical direction. Hence, the maximum principal
stress is equivalent to SH\textsubscript{max} and their directions are
also equivalent. We determine the percentage change in differential
stress (\(\Delta\sigma_D'=100\left(\sigma_D-\sigma_{D,\text{ref}}\right) / \sigma_{D,\text{ref}}\)) and the change in the orientation of
SH\textsubscript{max}~(\(\Delta\phi_{\max}'=\phi_{\max}-\phi_{\text{max,ref}}\)).~The
SH\textsubscript{max}~orientation (\(\phi_{\max}\)) shows the
direction of the maximum principal~stress (\(\sigma_1\)) in the
model and it corresponds to the direction of stress borehole breakout
measurements. Since our model domain was rotated clockwise about the
center of the crater by 35\selectlanguage{ngerman}°, a zero-degree~\(\phi_{\max}\) is
parallel to the regional stress field. Positive values
of~\(\phi_{\max}\) the signify an anticlockwise rotation from the
direction of applied loading and vice versa.~\(\Delta\sigma_D'\) is a
proxy for earthquake generation potential of the CSZ.~

\par\null

The entire modeling procedure adopted in this study is summarized in
Fig.~{\ref{321193}}.

\par\null\par\null\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/vel-model-CSZ8/meth2}
\caption{{Flow chart summarizing the methodology of this study.
{\label{321193}}%
}}
\end{center}
\end{figure}

\subsection*{Numerical Models}

{\label{973116}}

We select and control model parameters such as the velocity models, dips
of the rift faults, crater-crust moduli factor, frictional coefficients
and cohesion.~

\par\null

We describe an overview of the effects of faults using a model with
three rift~faults, simplified to be planar, with a general strike of
N70\selectlanguage{ngerman}°E and dip 70°SE. We use the velocity model for~the Saguenay region
(\citealp{Somerville1990}), the coefficient of friction and cohesion of 0.3
and 0 MPa, respectively; and the crater-crust elastic moduli and density
ratios of 0.25 and 0.9, respectively (i.e., the~SD70R25 model in
Table~{\ref{table}}). We stop the simulations
when~~\(\sigma_D\)\emph{(=~}\(\sigma_1-\sigma_3\)\emph{)} is 850MPa
at 10~km~\citep*{Byerlee1978}. This is a new stopping criterion from the
200 MPa~used in~\citet{Baird2010}.~ We repeat the simulation with 200
MPa stopping criterion and the crater-crust elastic moduli of 0.1 for
comparison of our results with the results from~\citet{Baird2010}
(SD70R25M1 model in Table {\ref{table}}).

\par\null

We determine an optimal~\(\mu\) that explains the seismicity
of the CSZ (SD70R25 models in Table~{\ref{table}}). We
use the~geometry of the SD70R25 model
(Table~{\ref{table}}) and use the 1D velocity models
for the Saguenay region for all simulations. We assume
that~\(\mu\) is the same for the three rift faults and it is
uniform on the fault planes, including the part of the faults within the
crater region. We vary~\(\mu\) from 0.1 to 0.6 at a constant
cohesion of 0 MPa.~ We stop the simulations when~~\(\sigma_D\)is
850MPa at 10~km~\citep*{Byerlee1978}. As~\(\mu\) increases, the
rift faults are stronger, and hence the model approaches the reference
model. In other words, we expect the~\(\Delta\sigma_D'\) to approach zero
as~\(\mu\)~increases.~

\par\null

We then determine an optimal C that explains the seismicity of the CSZ
(SD70R25 models). We use the 1D velocity models for the Saguenay region
for all simulations. We use C values of 0, 3, 5, 10, 20 and 30 MPa at a
constant~\(\mu\) and crater-crust moduli ratio of 0.3 and
0.25, respectively. The C of 10, 20 and 30 MPa are to test the model
results at higher values of cohesion. We also assume that
the~\(\mu\) and C~are constant on the three rift faults.~

\par\null

Different dip angles for the rift faults could cause the stress
orientation partitions observed in the CSZ. There are two fault
geometries for the CSZ. One model as dip~angles of 70°SE for the three
rift faults (SD70R25 model) (\citealp{Baird2010},
Fig.~{\ref{305206}}A) and the second as dip angles of
65°SE, 40°SE and~40°SE for the three rift faults (from north~to south)
(SD65R25 model) (\citealp*{Powell2017}
,~Fig.~{\ref{305206}}B). We test this hypothesis by
using the two fault geometries to explain observed seismicity and stress
rotation. Based on the correlation of the model results with the
observed seismicity and stress rotations, we aim to suggest the most
appropriate fault geometries~for the~CSZ. We use the 1D velocity models
for the Saguenay region~\citep{Somerville1990} and a crater-crust
moduli~ratio of 0.25 for both models. We use the~\(\sigma_D\) of
850 MPa at 10 km depth as the stopping criterion, and the frictional
coefficient and cohesion of 0.3 and 0 MPa, respectively.

\par\null

After determining the optimal uniform~\(\mu\) and C for the
three rift faults, we determine the relative strength between the crater
and the surrounding~crust that explains the observed~seismicity. We use
the 1D velocity models for the Saguenay region for all simulations. We
use crater-crust moduli ratios of 0.25, 0.5 and 1.0 (SD70R25, SD70R50
and SD70R100 models, respectively). The moduli ratio of 1.0 signifies
that the crater has the same velocity model as the crust. Based on model
results described in previous sections, we use a~\(\mu\) and
C of 0.3 and 0 MPa, respectively.~

Previous numerical modelings in this study, as
used~in~\citet{Baird2010}, use the 1D velocity model of the
Sangeney~region for the CSZ. Availability of 1D and 3D velocity models
specific to the CSZ poses a need to test how the specific velocity
models affect the results we describe hitherto. We compare the results
of three velocity models. These velocities are the~1D velocity of
the~Sangeney model (SD70R25) (\citealp{Somerville1990}), the 1D velocity model
of the CSZ (LD70R25) (\citealp*{Lamontagne1999}) and the 3D Vp and Vs
tomography (TD70) (\citealp*{Powell2017} )
(Figs.~{\ref{398984}}
and~{\ref{740205}} ). The model setups have the same
fault dips (70°-70°-70°). The crater-crust moduli ratio,
~\(\mu\),~and C are 0.25, 0.3 and 0 MPa, respectively. We do
not apply the moduli ratio in the TOMO simulations.

\par\null

\emph{Following~}\citet{Baird2010}\emph{,~we use~}
\(\sigma_{D}\)\emph{~of 200 MPa at 10 km depth as the stopping
criterion.}

\par\null

The summary of all the models is presented in
table~{\ref{table}}.



\section*{}

{\label{621593}}\par\null\par\null\selectlanguage{english}
\begin{table}
\caption{{List of numerical models}}
\centering
\begin{tabular}{lcccccc}
\hline
Model & Velocity Model & Dip of Rift Faults & Moduli ratio &
    $\mu$ & C \\
% Run  & Time (min)  \\
\hline
    SNFR25  & Saguenay Region & No faults & 0.25 & 0.3 & 0 \\
    SD70R25  &  & $70^\circ-70^\circ-70^\circ$  & 0.25 & 0.3 & 0 \\
    SD70R25$^{a}$  &  & $70^\circ-70^\circ-70^\circ$  & 0.25 & 0.3 & 0 \\
    SD70R25M1 &  & $70^\circ-70^\circ-70^\circ$  & 0.25 & 0.1 & 0 \\
    SD65R25  &  & $65^\circ-40^\circ-40^\circ$  & 0.25 & 0.3 & 0\\
    SD70R50  &  & $70^\circ-70^\circ-70^\circ$  & 0.5 & 0.3 & 0 \\
    SD70R100   &  & $70^\circ-70^\circ-70^\circ$  & 1.0 & 0.3 & 0 \\
    LD70R25   & 1-D Model of L99 & $70^\circ-70^\circ-70^\circ$  & 0.25 & 0.3  & 0 \\
    TD70 & Tomography by PL17 & $70^\circ-70^\circ-70^\circ$ & - & 0.3 & 0 \\
    
\hline
\multicolumn{2}{l}{$^{a}$Model has }
\end{tabular}
\end{table}




\par\null

\section*{RESULTS}

{\label{973504}}

\subsection*{Reference Model}

{\label{527307}}

We observe that the effects of the weaker crater are dominant in our
reference model SNFR25 (see
Table~{\ref{table}}).~\(\sigma_D\)~is
smaller~\textbf{{[}by how much?{]}} in the crater region relative to the
surrounding crustal rock as shown by 5 and 10-km depth sections in
Fig.~{\ref{660170}}A.~\(\sigma_D\)is greater in
the direction perpendicular to the applied loading than in the parallel
direction. 15 and 20 km depth sections show that~\(\sigma_D\)is
greater at the central region below the crater
(Fig.~{\ref{660170}}A). The effects of the crater
disappears at greater depths (e.g., 25 km depth section in
Fig.~{\ref{660170}}A) and \(\sigma_D\) field
becomes more or less uniform.

\par\null

~\(\phi_{max}\), the orientation of SH\textsubscript{max}, shows a
four-lobe pattern of~alternating polarities around the crater
(Fig.~{\ref{660170}}B)~\textbf{{[}Do you think you can
plot some oriented line segments or something that visualize
\textbackslash{}phi\_\{max\}?{]}}, uniformly approaching~
0\(\degree\), the loading direction, away from the crater and
also with depth (Fig.~{\ref{660170}}B). In the SE
corner of the crater,~the~\(\phi_{\max}\) shows a maximum clockwise
(i.e., negative) rotation of about 18.3\selectlanguage{ngerman}° from the regional stress
direction (5 km panel in
Fig.~{\ref{660170}}B).~\(\phi_{\max}\) at 5 km and
10 km depths is close to~0\(\degree\) within the crater, being
subparallel to the regional stress orientation (\citealp*{Zoback1992}).~

\par\null\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/F1-diff-stress-model-results/F1-m3-c0-ref-model-now}
\caption{{(A) Differential stress (\(\sigma_D\)) and (B)
SH\textsubscript{max} orientation (\(\phi_{\text{max}}\)) in the reference
model (SNFR25) at different depths. The white dots on the panels in (A)
are hypocenters within 2km of each depth slice.
{\label{660170}}%
}}
\end{center}
\end{figure}

\subsection*{}

{\label{487798}}

\subsection*{Fault Models}

{\label{487798}}

\textbf{{[}I copied the following texts in italic from the previous
section because they seem to fit better here. They might be redundant
with the existing contents of this section, though.{]}}

\emph{In SD70R25M1 model,
the~}\(\sigma_D\)\emph{~and~}\(\phi_{\max}\)\emph{~are similar
to that of the SD70R25 model. However, there are subtle differences. For
example, we didn't observe the bifurcation
in~}\(\sigma_D\)\emph{below 15 km depth as observed in SD70R25.
The maximum rotation of the~}\(\phi_{\max}\)\selectlanguage{ngerman}\emph{increases to 20°
compared to the 18.3° in the SD70R25 model.~}

\emph{We also observe similar reference model results in the LD70R25 as
observed in SD70R25 model. The maximum value
of~}\(\phi_{\max}\)\emph{is 19.4}\(\degree\)\emph{.}

\par\null

\emph{Increasing the crater-crust moduli ratio to 0.5 (i.e. SD70R50
model), the~}\(\sigma_D\)\emph{has a lower magnitude compared to
previous models. For example, the range of~}\(\sigma_D\) \emph{in
the 5 and 10 km panels are -1100/-250 MPa and -870/-350 MPa for the~
LD70R25 and SD70R50 models, respectively. The~}\(\phi_{\max}\)
\emph{reduces to 10}\(\degree\) \emph{due to the increased
strength of the~crater region relative to the surrounding~crustal rock.}

\par\null

\emph{The~}\(\sigma_D\) \emph{and~}\(\phi_{\max}\) \emph{in
SD70R100 model show little or no lateral variations at all depths but
also increases with depth.~The value of~}\(\phi_{\max}\) \emph{is
near zero ranging from 0.08}\(\degree\)\emph{(25 km depth) to
1.3}\(\degree\)\emph{(5 km depth).~In other words, the
lateral~}\(\sigma_D\)\emph{and~}\(\phi_{\max}\)
\emph{variations~observed in previous reference models are due to the
weaker crater region relative to the surrounding crustal rock.}

\par\null

\emph{The patterns in the~}\(\sigma_D\)
\emph{and~}\(\phi_{\max}\) \emph{in TD70 model~are different from
all the layered velocity models. We observe a
lower~}\(\sigma_D\)\emph{in regions with lower stiffness and
immediate surroundings. The maximum value of~}\(\phi_{\max}\)
\emph{is 2}\(\degree\)\emph{which is significantly lower than
the~}\(\phi_{\max\ }\)\emph{in the layered models. We observe that
the~}\(\phi_{\max}\)\emph{is zero (i.e., parallel to the regional
stress orientation) and the value of~}\(\sigma_D\) \emph{is
uniform at regions with no lateral variation in the velocity model. This
observation corroborates~the idea that the lateral changes in
the~}\(\sigma_D\) \emph{and~the rotation in
the~}\(\phi_{\max}\)\emph{, in model without rift faults, are due to
lateral changes~in the velocity model.}

\par\null

\subsubsection*{1. SD70R25}

{\label{409297}}

In this model that is the same with the reference model except that it
includes the rift faults, smaller magnitudes of~\(\sigma_{D}\) in
the crater than in the surrounding crust are still prominent but~the
rift faults significantly modifies the~\(\sigma_{D}\) field in the
vicinity of the rift faults
(Fig.~{\ref{864176}}A).~\(\sigma_D\)at a
distance greater than 70 km from the center of the crater is similar to
that of the reference model~\textbf{{[}at each depth?{]}}. The regions
of increase in~\(\sigma_D\) correlate with the locations of the
earthquakes at each depth \textbf{{[}\textless{}- I don't see this in
Fig. 7A.{]}}.

\par\null

The field of~\(\phi_{\max}\)is similar to that of the reference
model in that it shows a four-lobe pattern of~alternating polarities
around the weaker crater and approaches the E-W direction, i.e.,
0\(\degree\), away from the crater and also with depth
(Fig.~{\ref{864176}}B). However, the rift faults
further create significant clockwise rotation in their vicinity at all
depths.~\(\phi_{\max}\)is about~\(-\)10° within the
crater at a depth of 5 km but is subparallel to the regional stress
orientation at 10 km depth.~\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/F1-fault-model-results/F1-m3-c0-fault-model}
\caption{{(A) Differential stress and (B) SH\textsubscript{max} orientation in the
fault model (SD70R25) at different depths. The white dots represent
earthquakes within 2km of each depth slice.
{\label{864176}}%
}}
\end{center}
\end{figure}

\subsubsection*{2. SD70R25M1}

{\label{111899}}

The~\(\sigma_D\) and~\(\phi_{\max}\) is similar to the results
of SD70R25 model even though the~\(\mu\) is now 0.1. However,
the~\(\sigma_D\) as a lower value compared to SD70R25 model at
corresponding depths. The~\(\phi_{\max}\) increase to a maximum of
25.6\(\degree\)compared to the 20.3\(\degree\)we observe
in SD70R25 model.\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/F1-m3-c0-fault-model/F1-m1-c0-fault-model}
\caption{{SD70R25M1
{\label{279029}}%
}}
\end{center}
\end{figure}

\par\null

\subsubsection*{3. SD65R25}

{\label{424236}}

In the~SD65R25 model, the~\(\sigma_D\)~and~\(\phi_{\max}\)
results are different from the results from SD70R25
(Fig.~{\ref{353153}}A). We observe some regions with an
increase in~\(\sigma_D\)but they do not correlate with the
locations of earthquakes below 15 km depth.

\par\null

The SD65R25 model shows a higher clockwise rotation within the crater at
the 10 km depth while the~\(\phi_{\max}\) in the SD70R25 model has
near-zero values (Fig.~{\ref{353153}}B).~The maximum
value of~\(\phi_{\max}\) increases to 26.5\(\degree\) compared
to the 20.3\(\degree\) in the~SD70R25 model. However, the two
models~show the four-lobe pattern due to their lower crater-crust moduli
ratios and the magnitude of the rotation decreases with depth in both
models.\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/F2-m3-c0-fault-model/F2-m3-c0-fault-model}
\caption{{SD65R25
{\label{353153}}%
}}
\end{center}
\end{figure}

\par\null

\subsubsection*{4. LD70R25}

{\label{850025}}

We also observe an increase in~\(\sigma_D\)~near the rift faults
in the LD70R25 model (Fig.~{\ref{383000}}). However,
there is a decrease in the~lateral extent of the~region with~a decrease
in~\(\sigma_D\) at 5km depth. The maximum value
of~\(\phi_{\max}\) is 21.7\(\degree\) similar to the SD70R25
model.\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/F13-LAM-m3-c0-fault-model/F13-LAM-m3-c0-fault-model}
\caption{{LD70R25
{\label{383000}}%
}}
\end{center}
\end{figure}

\subsubsection*{5. SD70R50}

{\label{788824}}

The value of~\(\sigma_D\) is lower than in the SD70R25 model
(Fig.~{\ref{739684}}). The maximum value
of~\(\phi_{\max}\)reduces to 13\(\degree\) due to the~increase
in the strength of the impact~crater relative to the SD70R25 model.\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/F1-r05-m3-c0-fault-model/F1-r05-m3-c0-fault-model}
\caption{{SD70R50
{\label{739684}}%
}}
\end{center}
\end{figure}

\subsubsection*{6. SD70R100}

{\label{837890}}

The~SD70R100 model shows the effect of the rift faults since
crater-crust moduli ratio is 1.0. In other words, there is no lateral
variation in the velocity model. The~\(\sigma_D\)
and~\(\phi_{\max}\) are different from models with the crater.
The~\(\sigma_D\)did not show a significant decrease in the center
of the model geometry (i.e. the `crater' region) but shows an increase
near the rift faults and with depth
(Fig.~{\ref{226986}}A).~The maximum value
of~\(\phi_{\max}\)reduces to 12.8\(\degree\) (5km depth in
Fig.~{\ref{226986}} B).~ The absence of the crater
(i.e. absence of the four-lobe pattern in \(\phi_{\max}\))
emphasizes~the change in the polarity of~\(\phi_{\max}\) with
depth.~\(\ \phi_{\max}\) is negative (clockwise) at 5 and 10 km depth,
change to a near-zero value at 15 km depth and then turns
positive~(anticlockwise) at 20 and 25 km depths. (negative)
(Fig.~{\ref{226986}} B).

\par\null\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/F1-r10-m3-c0-fault-model/F1-r10-m3-c0-fault-model}
\caption{{SD70R100
{\label{226986}}%
}}
\end{center}
\end{figure}

\par\null

\subsubsection*{7. TD70}

{\label{707293}}

As~in the corresponding reference model, the results
of~\(\sigma_D\) and~\(\phi_{\max}\) in TD70~model~are different
from that of the layered velocity models. We observe alternate increase
and decrease in the~\(\sigma_D\)~in regions adjacent to the rift
faults (Fig.~{\ref{236109}}A). The~\(\phi_{\max}\)
is similar to the results of SD70R100 and its maximum value is
12.4\(\degree\)(Fig.~{\ref{236109}}B).~

\par\null\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/F5-T-m3-c0-fault-model/F5-T-m3-c0-fault-model}
\caption{{TD70
{\label{236109}}%
}}
\end{center}
\end{figure}

\par\null

\subsubsection*{Overall effects of rift faults on stress
distribution}

{\label{912965}}

\subsubsection*{SD70R25}

{\label{457621}}

We compare the model SD70R25 (see Table~{\ref{table}})
with the fault-free reference model (described
in~{\ref{527307}} and~{\ref{409297}})
because the effects of faults on stress distributions are best
highlighted by the differences between them. We determine the~percentage
change in differential stress (\(\Delta\sigma_D'\)) and change in
orientation of the SH\textsubscript{max} (\(\Delta \phi_{\text{max}}^{\prime}\)) at depths
5, 10, 15, 20 and 25 km (Fig.~{\ref{353858}}).

Positive~\(\Delta\sigma_D'\), corresponding to the increased differential
stress in SD70R25 relative to the corresponding reference model,
spatially overlap with the seismicity in general and particularly well
at 10 km (Fig.~{\ref{353858}}A). However, some
earthquakes in the southwestern part of the crater~at a depth
shallower~than 15 km~fall in the region of~negative~\(\Delta\sigma_D'\).
Also, some earthquakes~at 10 km depth in the northeastern corner of the
crater overlap with negative~\(\Delta\sigma_D'\) regions.

\par\null

The~\(\Delta\phi_{\max}'\) isolates the effect of the rift faults
on~\(\phi_{\max}\) . It's value varies spatially and approaches zero
at distances above 70 km from the crater region as in previous
quantities (Fig.~{\ref{353858}}B). The rift
faults~predominantly rotate the original~~\(\phi_{\max}\)in the
clockwise direction.~It's absolute value decrease with depth and goes
negative, i.e., anticlockwise, at 20 and 25 km depths. Despite the
traverse of the rift faults through the crater~and its uniform~friction
coefficient and cohesion everywhere on the faults, the
~\(\Delta\phi_{\max}'\) shows that the rift faults have a lesser effect
on~~\(\phi_{\max}\) within the crater. This lower clockwise rotation
due to the rift faults within the crater could be due to the less
strength of rocks therein relative to the surrounding crustal rock. At
10 km depth within the crater, the~~\(\Delta\phi_{\max}'\)goes to zero and
shows a small anticlockwise rotation at some regions
(Fig.~{\ref{353858}}B).~

\par\null\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/F1-model-results-now/F1-m3-c0-ch-diff-shmax}
\caption{{Percentage change in differential stress (\(\Delta\sigma_{\max}'\)) and the
change in \(\phi_{\max}\) (\(\ \Delta\phi_{\max}'\)) of the SD70R25 model.
{\label{353858}}%
}}
\end{center}
\end{figure}

\subsection*{}

{\label{836526}}

\subsubsection*{}

\subsubsection*{SD70R25M1}

The~\(\Delta\sigma_D'\) and~\(\Delta\phi_{\max}'\)in the SD70R25M1 are
similar to the results of the~SD70R25 model,~but their maximum values
have increased to~45\% and 25.5\(\degree\).~\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/F1-m3-c0-ch-diff-shmax/F1-m1-c0-ch-diff-shmax}
\caption{{SD70R25M1
{\label{440320}}%
}}
\end{center}
\end{figure}

\par\null

\subsubsection*{SD65R25}

{\label{229865}}

The~\(\Delta\sigma_D'\)in SD65R25 model does not spatially correlate with
the observed seismicity as we observe in the SD70R25 model
(Fig.~{\ref{602131}}A). In this model, most of the
earthquakes in the southeastern cluster fell in regions with no
potential for generating earthquakes,~i.e., within
negative~\(\Delta\sigma_D'\ \). At 15 km depth, the southwestern cluster of
earthquakes consistently deviate from the locations of the rift faults
in the model.~

In the~\(\Delta\phi_{\max}'\), the~ SD65R25 model shows more stress rotation
than the SD70R25 model (Fig.~{\ref{602131}}B). The
regions of this stress rotations are~at or near the location of the rift
faults in both models although the different dip angles of the rift
faults create different spacing between them with depth. The
~\(\Delta\phi_{\max}'\) is skewed to the southeast in the direction
perpendicular to the strike of the rift faults. The magnitude of
the~\(\Delta\phi_{\max}'\) decreases away from the rift faults as observed
in~\(\phi_{\max}\) plots.~

\par\null\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/F2-m3-c0-ch-diff-shmax/F2-m3-c0-ch-diff-shmax}
\caption{{SD65R25
{\label{602131}}%
}}
\end{center}
\end{figure}

\subsubsection*{LD70R25}\label{ld70r25}

Despite the difference in the velocity models, the results of LD70R25 is
similar to that of the SD70R25. The maximum value of~\(\Delta\sigma'_D\)
and~\(\Delta\phi_{\max}'\) in both models are also comparable.

\par\null\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/F13-LAM-m3-c0-ch-diff-shmax/F13-LAM-m3-c0-ch-diff-shmax}
\caption{{LD70R25
{\label{837380}}%
}}
\end{center}
\end{figure}

\par\null

\subsubsection*{SD70R50}\label{sd70r50}

The region of negative~\(\Delta\sigma'_{\max}\)~in~SD70R50 is smaller compared
to the size in the~SD70R25 model and does not correlate with the
seismicity of the CSZ (Fig.~{\ref{997185}}A).~
The~\(\Delta\sigma_{\max}'\) decrease to 27\% due to the increased strength of
the crater. Also, the~\(\Delta\phi_{\max}'\) show a more clockwise rotation
within the crater region than that of the SD70R25
(Fig.~{\ref{997185}}B). Specifically, the non-zero
value of~\(\Delta\phi_{\max}'\)~is absent at 10 km depth.~ ~

\par\null\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/F1-r05-m3-c0-ch-diff-shmax/F1-r05-m3-c0-ch-diff-shmax}
\caption{{SD70R50
{\label{997185}}%
}}
\end{center}
\end{figure}

\subsubsection*{SD70R100}\label{sd70r100}

The~\(\Delta\sigma_{\max}'\) in the SD70R100 model predicts~earthquakes along
the rift faults which do not explain observed seismicity
(Fig.~{\ref{627344}}A). At 5 and 10 km depths, the
earthquakes fall at regions with a~negative value of~\(\Delta\sigma_{\max}'\)
. The~\(\Delta\phi_{\max}'\) does not correlate with observed stress
orientation (Fig.~{\ref{627344}}B). The value
of~\(\phi_{\max}\)decrease to a near-zero value at 15 km depth. The
values of~\(\Delta\phi_{\max}\) became negative in the fault region at 20
and 25 km.

\par\null\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/F1-r10-m3-c0-ch-diff-shmax/F1-r10-m3-c0-ch-diff-shmax}
\caption{{SD70R100
{\label{627344}}%
}}
\end{center}
\end{figure}

\subsubsection*{TD70}\label{td70}

The results in the TD70 is similar to the~results when the crater-crust
moduli ratio is 1.0 (Fig.~{\ref{886711}}).
The~\(\Delta\sigma_{\max}'\) also predicts~earthquakes along the rift faults
which do not explain observed seismicity. The~\(\Delta\phi_{\max}'\) with
the TD70 does not correlate with observed stress orientation.

\par\null\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/F5-T-m3-c0-ch-diff-shmax/F5-T-m3-c0-ch-diff-shmax}
\caption{{TD70
{\label{886711}}%
}}
\end{center}
\end{figure}

\subsection*{}

{\label{836526}}\par\null\par\null\par\null\par\null\par\null\par\null\par\null\par\null\par\null\par\null\par\null\par\null\par\null\par\null\par\null\par\null

\section*{Discussions}

\par\null

\subsubsection*{Effect of Dips of the Rift
Faults~}

{\label{344318}}

Different dip angles for the rift faults could also cause the stress
orientation partitions observed in the CSZ. There are two fault
geometries for the CSZ. One model as dip~angles of 70\selectlanguage{ngerman}°SE for the three
rift faults (F1) (\citealp{Baird2010},
Fig.~{\ref{305206}}A) and the second as dip angles of
65°SE, 40°SE and~40°SE for the three rift faults (from north~to south)
(F2) (\citealp*{Powell2017} ,~Fig.~{\ref{305206}}B). We
test this hypothesis by using the two fault geometries to explain
observed seismicity and stress rotation. Based on the correlation of the
model results with the observed seismicity and stress rotations, we aim
to suggest the most appropriate fault geometries~for the~CSZ. We use the
1D velocity models for the Saguenay region~\citep{Somerville1990} and a
crater-crust moduli~ratio of 0.25 for both models.
Following~\citet{Baird2010},~we use the~\selectlanguage{greek}σ\selectlanguage{english}\textsubscript{D} of 200 MPa
at 10 km depth as the stopping criterion, and the frictional coefficient
and cohesion of 0.1 and 0 MPa, respectively. We use depth slices at 10
and 15 km depth to highlight the effect of dip angles of the rift
faults. We use these~depth slices because of the higher number of
earthquakes at these depths (Fig.~{\ref{305206}} ).~

\par\null\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/effect-fault-dips-now/effect-fault-dips-now}
\caption{{This is a caption
{\label{305206}}%
}}
\end{center}
\end{figure}

The~\(\Delta\sigma_D'\)in F1 models correlates better with the observed
seismicity than the F2 models. However,~the F1 model could not explain
some~earthquakes at 10 km in the SW part of the crater
(Fig.~{\ref{305206}}A). In the F2 model, most of the
earthquakes in the SE cluster fell regions with no potential for
generating earthquakes,~i.e., within negative~\(\Delta\sigma_D'\ \).
Another observation in the F2 model is that, at 15 km depth, all the
earthquakes in the SW cluster consistently deviate from the locations of
the rift faults.~

\par\null

Both models have some similarities in the~\(\phi_{\max}\),~but the F1
model is more consistent with observed
stress~orientation~\citep*{Mazzotti2010}. They both~show the four-lobe
pattern due to their lower crater-crust moduli ratios. F2 model has a
relatively~higher stress rotation than in the F1 model, but the
magnitude of the rotation in both models decreases with depth in both
models. However, the F2 model shows a slightly higher clockwise rotation
within the crater at the 10 km depth while the~\(\phi_{\max}\) in the
F1 model has near-zero values, hence consistent with
observations~\citep*{Mazzotti2010}.

\par\null

Also~in the~\(\Delta\phi_{\max}'\), the~ F2 model shows more stress rotation
than the F1 model. The regions of this stress rotations are~at or near
the location of the rift faults in both models.~The magnitude of
the~\(\Delta\phi_{\max}'\) decreases away from the rift faults as observed
in~\(\phi_{\max}\) plots.

\par\null

\subsubsection*{Effect of friction coefficient (\(\mu\)) of
the rift
faults}

{\label{102288}}

We determine an optimal~\(\mu\) that explains the seismicity
of the CSZ. We use the~geometry of the F1 model
(Table~{\ref{table}}) and use the 1D velocity models
for the Saguenay region for all simulations. We assume
that~\(\mu\) is the same for the three rift faults and it is
uniform on the fault planes, including the part of the faults within the
crater region. We vary~\(\mu\) from 0.1 to 0.6 at a constant
cohesion of 0 MPa.~ We stop the simulations when~~\(\sigma_D\)is
850MPa at 10~km~\citep*{Byerlee1978}. This is a new stopping criterion
from the one we use hitherto. As~\(\mu\) increases, the rift
faults are stronger, and hence the model approaches the reference model.
In other words, we expect the~\(\Delta\sigma_D'\) to approach zero
as~\(\mu\)~increases. We create a depth slice
of~\(\Delta\sigma_D'\),~~\(\phi_{\max}\) and~\(\Delta\phi_{\max}'\) at 10
km depth for the different~\(\mu\).

The~\(\Delta\sigma_D'\) with a~\(\mu\) of 0.3 explains the
seismicity of the CSZ (Fig.~{\ref{889549}}A). As
we~expect, the values of~\(\Delta\sigma_D'\) decreases with an increase
in~\(\mu\) . Below a~\(\mu\) of 0.3, the region of
positive~\(\Delta\sigma_D'\) is small compared to the spatial extent of
the observed seismicity.~Both~\(\mu\) of 0.3 and 0.4 can
explain the seismicity of the CSZ. However, the model result
with~\(\mu\) of 0.4 predicts an occurrence~of earthquakes
throughout the entire length of the rift faults. This result is
different from the spatial extent of observed
seismicity.~\(\Delta\sigma_D'\)is numerically zero
when~\(\mu\)is 0.5 and 0.6, which signifies that the seismic
zone approaches the reference Byerlee's model.~

\par\null

The~\(\phi_{\max}\)supports a~\(\mu\)of 0.3 for the rift
faults (Fig.~{\ref{889549}}B). The~\(\phi_{\max}\)
results show a clockwise stress rotation in the entire seismic zone
at~\(\mu\) below 0.3. At~\(\mu\) above 0.3, the
model results approach the reference model. However,
when~\(\mu\)is 0.3, the~\(\phi_{\max}\) is subparallel
to the region stress orientation and clockwise rotation in the SE part
of the crater, hence, consistent with the observation.~

\par\null\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/vel-model-CSZ9/effect-of-friction-coef-new}
\caption{{Effect of friction coefficient on the~percentage change in differential
stress and SH\textsubscript{max} in the fault model relative to the
reference model. The texts on each plot show the friction coefficients.
Note the change in the scale of the plots.
{\label{889549}}%
}}
\end{center}
\end{figure}

\subsection*{}

{\label{940877}}

\subsubsection*{Effect of fault cohesion
(C)}

{\label{645413}}

We determine an optimal C that explains the seismicity of the CSZ. We
use the 1D velocity models for the Saguenay region for all simulations.
We use C values of 0, 3, 5, 10, 20 and 30 MPa at a
constant~\(\mu\) and crater-crust moduli ratio of 0.3 and
0.25, respectively. The C of 10, 20 and 30 MPa are to test the model
results at higher values of cohesion. We also assume that
the~\(\mu\) and C~are constant on the three rift faults. We
use depth slice at 10 km depth because of its~high frequency of
earthquakes (Fig.~{\ref{243126}}).

\par\null

A C of 0 MPa explains the seismicity of the CSZ
(Fig.~{\ref{243126}}A).
The~\(\Delta\sigma_D'\)\textsubscript{~}and~\(\phi_{\max}\) decrease
with an increase in fault cohesion but their spatial distribution could
explain observed seismicity. The~\(\Delta\sigma_D'\) is not affected by a
change in C between 0 and 5 MPa. Hence,~\(\mu\)~is an
essential fault parameter to explain the seismicity of the CSZ as it
affects the slope of the failure envelope.~However, the models could not
explain some earthquakes~in the SW corner of
the~crater.~\(\Delta\phi_{\max}'\) also decreases with an increase in C. As
C increases, the effect of the rift faults on stress rotations within
the crater decreases and concentrate in the fault region outside the
crater.

\par\null\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/effect-of-cohesion/effect-of-cohesion}
\caption{{Effect of cohesion on percentage change in differential stress,
SH\textsubscript{max} and the change in SH\textsubscript{max}.
{\label{243126}}%
}}
\end{center}
\end{figure}

\subsection*{~}

{\label{408620}}

\subsection*{Effects of crater-crust elastic
moduli~ratio}

{\label{484079}}

After determining the optimal uniform~\(\mu\) and C for the
three rift faults, we determine the relative strength between the crater
and the surrounding~crust that explains the observed~seismicity. We use
the 1D velocity models for the Saguenay region for all simulations. We
use crater-crust moduli ratios of 0.25, 0.5 and 1.0. The moduli ratio of
1.0 signifies that the crater has the same velocity model as the crust.
Based on model results described in previous sections, we use
a~\(\mu\) and C of 0.3 and 0 MPa, respectively. Since the
moduli ratio will have a~maximum effect within the crater, we create
depth slices at 5 and 10 km (Fig.~{\ref{435037}}). The
maximum depth of the crater is 15 km.~

\par\null\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/effect-moduli-ratio-now/F1-effect-of-moduli-ratio}
\caption{{Effect of crater-crust moduli ratio of the seismicity of the CSZ.
{\label{435037}}%
}}
\end{center}
\end{figure}

The moduli ratio of 0.25 shows the better correlation between
the~spatial extent of the observed seismicity and
the~\(\Delta\sigma_D'\). The value and spatial extent
of~\(\Delta\sigma_D'\)decreases in value and spatial extent as the moduli
ratio increases from 0.25 to 0.5.~ When the~crater is absent, i.e., a
crater-crust moduli ratio of 1.0, the value of~\(\Delta\sigma_D'\)is
generally positive and thus could not explain the observed seismicity.

\par\null

The~\(\phi_{\max}\) supports a moduli ratio of 0.25. At 10 km depth,
the~\(\phi_{\max}\) is near-zero within the crater when the moduli
ratio is 0.25 and deviates from the regional stress orientation as the
moduli ratio increases. At a moduli ratio of 1.0, the four-lobe pattern
disappears and the model gives a~net clockwise stress rotation. This
confirms that the four-lobed pattern is due to the lower crater-crust
moduli ratio. The~\(\phi_{\max}\) decreases to zero away from the CSZ
in all the models,~therefore, approach the regional stress orientation

\par\null\par\null

\subsection*{Effects of the velocity
models}

{\label{972731}}

Previous numerical modelings in this study, as
used~in~\citet{Baird2010} , use the 1D velocity model of the
Sangeney~region for the CSZ. Availability of 1D and 3D velocity models
specific to the CSZ poses a need to test how the specific velocity
models affect the results we describe hitherto. We compare the results
of three velocity models. These velocities are the~1D velocity of
the~Sangeney model (F1\_0.25) (\citealp{Somerville1990}), the 1D velocity
model of the CSZ (LAM\_0.25) (\citealp*{Lamontagne1999}) and the 3D Vp and Vs
tomography (TOMO) (\citealp*{Powell2017} )
(Figs.~{\ref{398984}}
and~{\ref{740205}} ). The model setups have the same
fault dips (70\selectlanguage{ngerman}°-70°-70°). The crater-crust moduli ratio,
~\(\mu\),~and C are 0.25, 0.3 and 0 MPa, respectively. We do
not apply the moduli ratio in the TOMO simulations. We create depth
slices at 10 and 15 km to investigate the effect of the velocity models
on the displacement and stress solutions within and under the crater
region. Many earthquakes occur at these two depths.

\par\null

The velocity model of~\citet*{Lamontagne1999}~is significantly different from
that of the Saguenay region. The velocity model of~\citet*{Lamontagne1999}
has higher values of P- and S-wave velocities,~and density than the
Saguenay~region at shallower depths up to 35km. At deeper depths, the
velocities of the Saguenay region are higher than that of LAM\_0.25. The
1D velocity model of LAM\_0.25 is similar to the halfspace
velocity~model for the Canadian Shield
(\citealp*{Lamontagne1999};~\citealp*{Powell2017}). The TOMO, on the other hand,
show some lateral variations in the P- and S-wave velocities. The TOMO
shows the geometry of a high-velocity body at mid-crustal depth.
However, the TOMO~ do not show a distinct decrease in velocities in the
crater region. Therefore, except for the lateral variation in
velocities, the TOMO is similar to a 1D velocity model with a
crater-crater moduli ratio of 1.0.

\par\null\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/effect-of-high-velocity-body1/effect-of-high-velocity-body}
\caption{{Effects of the velocity models.
{\label{798761}}%
}}
\end{center}
\end{figure}

\par\null

All the models show a decrease in~\(\Delta\sigma_D'\),
~\(\phi_{\max}\),~and~\(\Delta\phi_{\max}'\) with depth and with distance
away from the crater (Fig.~{\ref{798761}}).~

\par\null

The~\(\Delta\sigma_D'\)in the F1\_0.25~correlates with observed
seismicity at both 10 and 15 km depths
(Fig.~{\ref{798761}}A) whereas the
positive~\(\Delta\sigma_D'\) region extends farther at 15 km depth.~
Despite the difference in the velocity models used for the LAM\_0.25 and
F1\_0.25, the~\(\Delta\sigma_D'\)results are similar. This similarity in
the model results indicates that the 1D velocity from the refraction
study in the CSZ (\citealp*{Lamontagne1999}) is appropriate for stress modeling
in the CSZ.~The results in the TOMO is similar to the~results when the
crater-crust moduli ratio is 1.0
(Fig.~{\ref{435037}}C). TOMO predicts~earthquakes along
the rift faults which do not explain observed seismicity.~

\par\null

The~\(\phi_{\max}\) at 10 and 15 km depth using the 1D velocity
models correlate with the
observed~\(\phi_{\max}\)(Fig.~{\ref{798761}} B).
The~\(\phi_{\max}\) is parallel to the regional stress field within
the crater a 10 km depth. However, the~\(\phi_{\max}\) with the TOMO
do not correlate with observed stress orientation. However, near-zero
values of~\(\phi_{\max}\) at 15 km depth in the~TOMO could be due~to
the high-velocity mid-crustal bodies.

\par\null

We observe similar~\(\Delta\phi_{\max}'\)in the two 1D velocity models
(Fig.~{\ref{798761}} C). This observation further
supports that the 1D velocity models with a crust moduli ratio of 0.25
fits better with observed seismicity and~\(\phi_{\max}\). the stress
rotation in the TOMO models is due to the rift faults due to the
similarities between the~\(\Delta\phi_{\max}'\) and
the~\(\phi_{\max}\). ~

\par\null\par\null\par\null\par\null\par\null\par\null\par\null\par\null\par\null\par\null

\section*{}

{\label{436996}}

In other words, the majority of the Discussion and Conclusions sections
should be an interpretation of your results. You should:

\begin{itemize}
\tightlist
\item
  Discuss your conclusions in order of~\textbf{most to least important.}
\item
  \textbf{Compare~}your results with those from other studies: Are they
  consistent? If not, discuss possible reasons for the difference.
\item
  Mention any~\textbf{inconclusive results}~and explain them as best you
  can. You may suggest additional experiments needed to clarify your
  results.
\item
  Briefly describe the~\textbf{limitations~}of your study to show
  reviewers and readers that you have considered your experiment's
  weaknesses. Many researchers are hesitant to do this as they feel it
  highlights the weaknesses in their research to the editor and
  reviewer. However doing this actually makes a positive impression of
  your paper as it makes it clear that you have an in depth
  understanding of your topic and can think objectively of your
  research.
\item
  Discuss~\textbf{what your results may mean}~for researchers in the
  same field as you, researchers in other fields, and the general
  public. How could your findings be applied?
\item
  State how your results~\textbf{extend the findings}~of previous
  studies.
\item
  If your findings are preliminary, suggest~\textbf{future studies}~that
  need to be carried out.
\item
  At the end of your Discussion and Conclusions sections,~\textbf{state
  your main conclusions once again}.
\end{itemize}

\par\null\par\null

\subsection*{Discuss your conclusions in order of~most to least
important.}

{\label{796086}}

\subsubsection*{The coefficient of friction and cohesion of rift
faults}

{\label{218853}}

~The friction coefficient and cohesion of the rift faults in the CSZ are
0.3 and 0 MPa, respectively.~

\par\null

~ ~ ~ ~ ~ ~ Discuss the effects of friction coefficient and cohesion.

\par\null

\subsubsection*{Crater-crust moduli
ratio}

{\label{918530}}

A lower stiff crater, up to a quarter of the elastic moduli, relative to
the surrounding crust is required to explain the seismicity in the upper
10 km of the CSZ and the observed clockwise rotation of the SHmax in the
SE relative to the NW part of the crater. Below 10 km depth, the
reactivation of the rift faults explains the observed seismicity in the
CSZ.

\par\null

~ ~ ~ ~ ~ ~ Discuss the effects of moduli ratio

\par\null\par\null

\subsubsection*{Dips of the rift faults}

{\label{399239}}

Magnitude of SHmax are higher in the 65-40-40 fault dips than the
constant 70 fault dips by 6 deg but does not fit the seismicity. The
remaining SHmax could be due to change in the dip of the rift faults
with depth.

\par\null

~ ~ ~ ~ ~ ~ Discuss effects of dips

\par\null

\subsubsection*{}

{\label{135552}}

\subsubsection*{Velocity model of the
CSZ}

{\label{766326}}

1D velocity models of the CSZ with embedded crater and constant-dip rift
faults can explain the observed seismicity and more than 50\% of the
stress rotations (SHmax).

\par\null

~ ~ ~ ~ ~ ~ DISCUSS THE EFFECTS OF VELOCITY MODEL

~

\par\null

\subsubsection*{Direction of loading~}

{\label{162144}}

The observed clockwise rotation of the SHmax relative to the first order
direction of the tectonic stress is due to the relative angle between
the rift faults and the regional tectonic stress and not the damaged
zone.

\par\null\par\null\par\null\par\null\par\null\par\null\par\null

~

\par\null

\par\null\par\null\par\null\par\null

Our study shows that the rift faults are responsible for the stress
rotations observed in the CSZ (Fig.~{\ref{353858}}C).
The~\(\Delta\phi_{\max}'\)isolates the contribution of the rift faults to
the stress rotation observed in the~\(\phi_{\max}\).
(Fig.~{\ref{353858}}C).~\(\Delta\phi_{\max}'\) shows a
general clockwise stress rotation within the rift faults with a maximum
of 21.5\selectlanguage{ngerman}°. We observe a lower stress rotation~within the crater. The
\(\Delta\phi_{\max}'\) is near-zero within the crater at a depth of 10 km.

This study shows that a lower stiff crater, up to a quarter of the
elastic~moduli, relative to the surrounding crust is required to explain
the seismicity in the upper 10 km of the CSZ. Below 10 km depth, the
reactivation of the rift faults explains the observed seismicity in the
CSZ.

\par\null

Most of the earthquakes in the NW cluster are within the crater region
and occur at about 10 km depth (Fig.~{\ref{749307}}).
The subparallel~stress orientation in the crater region could explain
the observed stress inversion result from focal mechanisms in the NW
cluster of the CSZ. Also, the significant negative values of the
SH\textsubscript{max} in the SE part of the model (outside of the
crater) could explain the observed stress orientation in the SE cluster.

\par\null\par\null\par\null\par\null\par\null\par\null\par\null\par\null

\subsection*{Compare~your results with those from other studies: Are they
consistent? If not, discuss possible reasons for the
difference.}

{\label{536638}}

However, the distinct and lower strength damaged zone is not apparent in
the 3D tomography (\citealp*{Powell2017}).

\par\null

The observed clockwise rotation of the SH\textsubscript{max} in the SE
relative to the NW part of the crater is as a result of the relative
angle between the rift faults and the first order direction of the
tectonic stress (Fig.~{\ref{798761}}). A clockwise
rotation in the orientation of SH\textsubscript{max} from stress
inversion relative to SH\textsubscript{max} from borehole breakout is
also observed in the St Lawrence seismic zone~\citep*{Mazzotti2010}. St.
Lawrence seismic zone is also along the St. Lawrence River to the south
of the CSZ. This observation further suggests that the observed rotation
in the SH\textsubscript{max} may not due to the presence of the damaged
zone created by the crater. We propose that the rotation may be due to
the presence of the rift faults.

\par\null

This study shows that the friction coefficient and cohesion of the rift
faults are 0.3 and 0 MPa. This is contrary to the previous works in the
CSZ (\citealp*{Lamontagne1999};~\citealp{Baird2010}). For
example,~\citet{Baird2010} conclude that the rift faults are very weak
and a friction angle of 5°. This friction angle corresponds to
a~\(\mu\) of about 0.1. We reproduce the results
of~\citet{Baird2010}~using the~\(\mu\) and the same stopping
criterion of~\(\sigma_D\) of 200 MPa at 10 km depth. The model
results also support the very weak (\(\mu\) = 0.1) rift
faults. This shows that the stopping criterion plays a very important
role in the model results. In this study, the~\(\sigma_D\) is the
difference between the maximum and minimum principal stresses.~We derive
a~\(\sigma_D\)~of 706 MPa at 10 km depth without making any
assumptions in the intermediate principal stress. We use this as the
stopping criterion for our simulations.

\par\null\par\null\par\null\par\null\par\null\par\null\par\null\par\null\par\null\par\null

\subsection*{Mention any~inconclusive results~and explain them as best
you can. You may suggest additional experiments needed to clarify your
results.}

{\label{585578}}\par\null\par\null\par\null\par\null

\subsection*{Briefly describe the~limitations~of your study to show
reviewers and readers that you have considered your experiment's
weaknesses.}

{\label{468113}}\par\null\par\null\par\null\par\null\par\null

\subsection*{Discuss~what your results may mean~for researchers in the
same field as you, researchers in other fields, and the general public.
How could your findings be
applied?}

{\label{747341}}\par\null\par\null\par\null\par\null\par\null

\subsection*{State how your results~extend the findings~of previous
studies.}

{\label{692365}}\par\null\par\null\par\null\par\null\par\null

\subsection*{If your findings are preliminary, suggest~future
studies~that need to be carried
out.}

{\label{857425}}\par\null\par\null\par\null\par\null\par\null

\subsection*{}

{\label{932499}}\par\null\par\null\par\null\par\null\par\null\par\null\par\null\par\null\par\null\par\null\par\null\par\null\par\null\par\null\par\null\par\null\par\null\par\null\par\null\par\null





The distribution of earthquakes in the CSZ can be explained by the
change in differential stress in a~numerical model.







\par\null

~

\par\null\par\null\par\null\par\null\par\null

However, the less stiff crater is necessary to explain why the SHmax
orientation in the~NW part of the crater is parallel to the direction of
tectonic stress.

\par\null\par\null\par\null\par\null

2. We didn't see~ and/or the mid-crustal high-velocity body.~ Rotation
of the SHmax orientation relative to the regional stress orientation
observed in the New Madrid Seismic Zone has also been attributed to the
presence of high velocity denser mafic intrusions (Levandowski et al.,
2016).

\par\null

2b. Numerical models with layered velocity model and a crater-crust
moduli~ratio of 0.25 and 0.5 show the gap between the cluster of
earthquakes on fault A and B. This gap is also observed in the TOMO
models. However, there is a significant difference in the SHmax
orientation and~\selectlanguage{greek}Δ\selectlanguage{english}SHmax. The region with the mid-crustal high-velocity
body shows a near zero~\selectlanguage{greek}Δ\selectlanguage{english}SHmax and SHmax orientation. This means that the
mid-crustal high-velocity body acts to align thus rotate the SHmax
orientation anticlockwisely~with the direction of regional stress.~ ~

\par\null

3. Faults' friction coefficients:~This conclusion correlates with the
result of Baird et al. (2009).

\par\null

4. We can only explain part of the observed rotation. These earthquakes
may be explained with change in dip angle with depth.~A clockwise
rotation of about 20\selectlanguage{ngerman}° in the SHmax orientation. This is similar to the
32° observed by (Mazzotti and Townend, 2010).~Within the CSZ, we
observed a local clockwise rotation of about 10°. (Mazzotti and Townend,
2010) observed 46° clockwise rotation.~

5. We can not expalin~all the eartqhuake~distributions.~For a model with
crater-crust moduli ratio of 0.25, the~\selectlanguage{greek}Δσ\selectlanguage{english}D~correlates with the
seismicity of the CSZ. The model shows a positive~\selectlanguage{greek}Δσ\selectlanguage{english}D~on
the~northwestern part of the rift faults. The velocity model still
explains the seismicity and further explains the absence of earthquakes
on the northwestern side of the rift faults.

6. 70-70-70 dip angle of faults explains the observed seismicity
better.~The result patterns in these models explain the recorded
seismicity but the locations of the rift faults are inaccurate. The
SHmax orientation gives a higher value in the crater region compared to
the 70\selectlanguage{ngerman}°SE-dip models.

\par\null

7.~Therefore, the F2 models with lower dip angles (F2) are consistently
inadequate at explaining the seismicity and the SHmax orientation.

\par\null

8.~

\par\null

We will later change the stopping criterion to a \selectlanguage{greek}σ\selectlanguage{english}D~of 706 MPa at 10 km
depth to discuss the difference in the model results.~

~

The rift faults are responsible for the stress rotation on near and away
from the rift faults but not much effect within the crater.

Therefore, the SHmax orientation spatially overlaps the observed SHmax
from the stress inversion of focal mechanisms~\citep*{Mazzotti2010}.



\subsection*{Conclusions}

{\label{929830}}\par\null

\subsection*{At the end of your Discussion and Conclusions
sections,~state your main conclusions once
again.}

\par\null

Our work explains the change in the SH\textsubscript{max} orientation
inverted from earthquake focal mechanism compared to
SH\textsubscript{max} orientation determined from borehole breakout
measurements. We also consider and discuss the effects of the
high-velocity bodies imaged at mid-crustal depths~\citep*{Powell2017} on
the region's stress distribution.





\par\null

\subsection*{Acknowledgment}

{\label{459725}}

We thank the Computational Infrastructure for Geodynamics
(\href{http://geodynamics.org/}{http://geodynamics.org}) which is funded
by the National Science Foundation under awards EAR-0949446 and
EAR-1550901.






\section*{Results}
\subsection*{Reference Model}
We observe that the effects of the weaker crater are dominant in our reference model SNFR25 (see Table 1). $\sigma_D$ is smaller \textbf{[by how much?]} in the crater region relative to the surrounding crustal rock as shown by 5 and 10-km depth sections in Fig. 6A. $\sigma_D$ is greater in the direction perpendicular to the applied loading than in the parallel direction. 15 and 20 km depth sections show that $\sigma_D$ is greater at the central region below the crater (Fig. 6A). The effects of the crater disappears at greater depths (e.g., 25 km depth section in Fig. 6A) and $\sigma_D$ field becomes more or less uniform.

$\phi_{max}$, the orientation of SH$_{max}$, shows a four-lobe pattern of alternating polarities around the crater (Fig. 6B) \textbf{[Do you think you can plot some oriented line segments or something that visualize $\phi_{max}$?]}, uniformly approaching 0$^\circ$ , the loading direction, away from the crater and also with depth (Fig. 6B). In the SE corner of the crater, the $\phi_{\max}$ shows a maximum clockwise (i.e., negative) rotation of about 18.3$^\circ$ from the regional stress direction (5 km panel in Fig. 6B). $\phi_{\max}$ at 5 km and 10 km depths is close to 0$^\circ$ within the crater, being subparallel to the regional stress orientation (Zoback, 1992).\selectlanguage{english}
\begin{figure}[ht]
\centering
\includegraphics[width=30pc]{Figures/SD70R25/F1_m3_c0_ref_model_now.png}
\caption{{(A) Differential stress ($\sigma_D$) and (B) SH$_{max}$ orientation ($\phi_{\text{max}}$) in the reference model (SNFR25) at different depths. The white dots on the panels in (A) are hypocenters within 2km of each depth slice.}}
\label{figfour}
\end{figure}


\subsection*{Fault Models}
\textbf{[I copied the following texts in italic from the previous section because they seem to fit better here. They might be redundant with the existing contents of this section, though.]}

\selectlanguage{ngerman}\textit{In SD70R25M1 model, the $\sigma_D$  and $\phi_{\max}$ are similar to that of the SD70R25 model. However, there are subtle differences. For example, we didn't observe the bifurcation in $\sigma_D$  below 15 km depth as observed in SD70R25. The maximum rotation of the $\phi_{\max}$ increases to 20$^\circ$ compared to the 18.3° in the SD70R25 model.}

\textit{We also observe similar reference model results in the LD70R25 as observed in SD70R25 model. The maximum value of $\phi_{\max}$ is 19.4$^\circ$.}

\textit{Increasing the crater-crust moduli ratio to 0.5 (i.e. SD70R50 model), the $\sigma_D$ has a lower magnitude compared to previous models. For example, the range of $\sigma_D$ in the 5 and 10 km panels are -1100/-250 MPa and -870/-350 MPa for the  LD70R25 and SD70R50 models, respectively. The $\phi_{\max}$ reduces to 10$^\circ$ due to the increased strength of the crater region relative to the surrounding crustal rock.}

\textit{The $\sigma_D$ and $\phi_{\max}$ in SD70R100 model show little or no lateral variations at all depths but also increases with depth. The value of $\phi_{\max}$ is near zero ranging from 0.08$^\circ$ (25 km depth) to 1.3$^\circ$ (5 km depth). In other words, the lateral $\sigma_D$ and $\phi_{\max}$ variations observed in previous reference models are due to the weaker crater region relative to the surrounding crustal rock.}

\textit{The patterns in the $\sigma_D$ and $\phi_{\max}$ in TD70 model are different from all the layered velocity models. We observe a lower $\sigma_D$ in regions with lower stiffness and immediate surroundings. The maximum value of $\phi_{\max}$ is 2$^\circ$ which is significantly lower than the $\phi_{\max\ }$ in the layered models. We observe that the $\phi_{\max}$ is zero (i.e., parallel to the regional stress orientation) and the value of $\sigma_D$ is uniform at regions with no lateral variation in the velocity model. This observation corroborates the idea that the lateral changes in the $\sigma_D$ and the rotation in the $\phi_{\max}$, in model without rift faults, are due to lateral changes in the velocity model.}

\subsubsection*{SD70R25}
In this model that is the same with the reference model except that it includes the rift faults, smaller magnitudes of $\sigma_{D}$ in the crater than in the surrounding crust are still prominent but the rift faults significantly modifies the $\sigma_{D}$ field in the vicinity of the rift faults (Fig. 7A). $\sigma_D$ at a distance greater than 70 km from the center of the crater is similar to that of the reference model \textbf{[at each depth?]}. The regions of increase in $\sigma_D$ correlate with the locations of the earthquakes at each depth \textbf{[<- I don't see this in Fig. 7A.]}.

The field of $\phi_{\max}$ is similar to that of the reference model in that it shows a four-lobe pattern of alternating polarities around the weaker crater and approaches the E-W direction, i.e., 0$^\circ$, away from the crater and also with depth (Fig. 7B). However, the rift faults further create significant clockwise rotation in their vicinity at all depths. $\phi_{\max}$ is about -10° within the crater at a depth of 5 km but is subparallel to the regional stress orientation at 10 km depth.\selectlanguage{english}
\begin{figure}[ht]
\centering
\includegraphics[width=30pc]{Figures/SD70R25/F1_m3_c0_fault_model.png}
\caption{{(A) Differential stress and (B) SHmax orientation in the fault model (SD70R25) at different depths. The white dots represent earthquakes within 2km of each depth slice.}}
\label{figfour}
\end{figure}


\subsubsection*{SD70R25M1}
The $\sigma_D$ and $\phi_{\max}$ is similar to the results of SD70R25 model even though the $\mu$ is now 0.1. However, the $\sigma_D$ as a lower value compared to SD70R25 model at corresponding depths. The $\phi_{\max}$ increase to a maximum of 25.6$^\circ$ compared to the 20.3$^\circ$ we observe in SD70R25 model.\selectlanguage{english}
\begin{figure}[ht]
\centering
\includegraphics[width=30pc]{Figures/SD70R25M1/F1_m1_c0_fault_model.png}
\caption{{SD70R25M1}}
\label{figfour}
\end{figure}




\subsubsection*{SD65R25}
In the SD65R25 model, the $\sigma_D$ and $\phi_{\max}$ results are different from the results from SD70R25 (Fig. 9A). We observe some regions with an increase in $\sigma_D$ but they do not correlate with the locations of earthquakes below 15 km depth.

The SD65R25 model shows a higher clockwise rotation within the crater at the 10 km depth while the $\phi_{\max}$ in the SD70R25 model has near-zero values (Fig. 9B). The maximum value of $\phi_{\max}$ increases to 26.5$^\circ$ compared to the 20.3$^\circ$ in the SD70R25 model. However, the two models show the four-lobe pattern due to their lower crater-crust moduli ratios and the magnitude of the rotation decreases with depth in both models.\selectlanguage{english}
\begin{figure}[ht]
\centering
\includegraphics[width=30pc]{Figures/SD65R25/F2_m3_c0_fault_model.png}
\caption{{SD70R25M1}}
\label{figfour}
\end{figure}


\subsubsection*{LD70R25}
We also observe an increase in $\sigma_D$ near the rift faults in the LD70R25 model (Fig. 10). However, there is a decrease in the lateral extent of the region with a decrease in $\sigma_D$ at 5km depth. The maximum value of $\phi_{\max}$ is 21.7$^\circ$ similar to the SD70R25 model.\selectlanguage{english}
\begin{figure}[ht]
\centering
\includegraphics[width=30pc]{Figures/LD70R25/F13_LAM_m3_c0_fault_model.png}
\caption{{LD70R25}}
\label{figfour}
\end{figure}



\subsubsection*{SD70R50}
The value of $\sigma_D$ is lower than in the SD70R25 model (Fig. 11). The maximum value of $\phi_{\max}$ reduces to 13$^\circ$ due to the increase in the strength of the impact crater relative to the SD70R25 model.\selectlanguage{english}
\begin{figure}[ht]
\centering
\includegraphics[width=30pc]{Figures/SD70R50/F1_r05_m3_c0_fault_model.png}
\caption{{SD70R50}}
\label{figfour}
\end{figure}



\subsubsection*{SD70R100}
The SD70R100 model shows the effect of the rift faults since crater-crust moduli ratio is 1.0. In other words, there is no lateral variation in the velocity model. The $\sigma_D$ and $\phi_{\max}$ are different from models with the crater. The $\sigma_D$ did not show a significant decrease in the center of the model geometry (i.e. the 'crater' region) but shows an increase near the rift faults and with depth (Fig. 12A). The maximum value of $\phi_{\max}$ reduces to 12.8$^\circ$ (5km depth in Fig. 12 B). The absence of the crater (i.e. absence of the four-lobe pattern in $\phi_{\max}$) emphasizes the change in the polarity of $\phi_{\max}$ with depth. $\phi_{\max}$ is negative (clockwise) at 5 and 10 km depth, change to a near-zero value at 15 km depth and then turns positive (anticlockwise) at 20 and 25 km depths. (negative) (Fig. 12 B).\selectlanguage{english}
\begin{figure}[ht]
\centering
\includegraphics[width=30pc]{Figures/SD70R100/F1_r10_m3_c0_fault_model.png}
\caption{{SD70R100}}
\label{figfour}
\end{figure}



\subsubsection*{TD70}
As in the corresponding reference model, the results of $\sigma_D$ and $\phi_{\max}$in TD70 model are different from that of the layered velocity models. We observe alternate increase and decrease in the $\sigma_D$in regions adjacent to the rift faults (Fig. 13A). The $\phi_{\max}$ is similar to the results of SD70R100 and its maximum value is 12.4$^\circ$ (Fig. 13B).\selectlanguage{english}
\begin{figure}[ht]
\centering
\includegraphics[width=30pc]{Figures/TD70/F5_T_m3_c0_fault_model.png}
\caption{{TD70}}
\label{figfour}
\end{figure}



\subsection*{Overall effects of rift faults on stress distribution}

\subsubsection*{SD70R25}
We compare the model SD70R25 (see Table 1) with the fault-free reference model (described in 4.1 and 4.2.1) because the effects of faults on stress distributions are best highlighted by the differences between them. We determine the percentage change in differential stress ($\Delta\sigma_D'$) and change in orientation of the SH$_{max}$ ($\Delta \phi_{\text{max}}^{\prime}$) at depths 5, 10, 15, 20 and 25 km (Fig. 14).

Positive $\Delta\sigma_D'$, corresponding to the increased differential stress in SD70R25 relative to the corresponding reference model, spatially overlap with the seismicity in general and particularly well at 10 km (Fig. 14A). However, some earthquakes in the southwestern part of the crater at a depth shallower than 15 km fall in the region of negative $\Delta\sigma_D'$. Also, some earthquakes at 10 km depth in the northeastern corner of the crater overlap with negative $\Delta\sigma_D'$ regions.

The $\Delta\phi_{\max}'$ isolates the effect of the rift faults on $\phi_{\max}$. It's value varies spatially and approaches zero at distances above 70 km from the crater region as in previous quantities (Fig. 14B). The rift faults predominantly rotate the original $\phi_{\max}$ in the clockwise direction. It's absolute value decrease with depth and goes negative, i.e., anticlockwise, at 20 and 25 km depths. Despite the traverse of the rift faults through the crater and its uniform friction coefficient and cohesion everywhere on the faults, the $\Delta\phi_{\max}'$ shows that the rift faults have a lesser effect on $\phi_{\max}$ within the crater. This lower clockwise rotation due to the rift faults within the crater could be due to the less strength of rocks therein relative to the surrounding crustal rock. At 10 km depth within the crater, the $\Delta\phi_{\max}'$ goes to zero and shows a small anticlockwise rotation at some regions (Fig. 14B).\selectlanguage{english}
\begin{figure}[ht]
\centering
\includegraphics[width=30pc]{Figures/SD70R25/F1_m3_c0_fault_model.png}
\caption{{Percentage change in differential stress ($\Delta\sigma_{\max}'$) and the change in $\phi_{\max}$ ($\Delta\phi_{\max}'$) of the SD70R25 model.}}
\label{figfour}
\end{figure}


\subsubsection*{SD70R25M1}
The $\Delta\sigma_D'$ and $\Delta\phi_{\max}'$ in the SD70R25M1 are similar to the results of the SD70R25 model, but their maximum values have increased to 45$\%$ and 25.5$^\circ$.\selectlanguage{english}
\begin{figure}[ht]
\centering
\includegraphics[width=30pc]{Figures/SD70R25M1/F1_m1_c0_fault_model.png}
\caption{{SD70R25M1}}
\label{figfour}
\end{figure}


\subsubsection*{SD65R25}
The $\Delta\sigma_D'$ in SD65R25 model does not spatially correlate with the observed seismicity as we observe in the SD70R25 model (Fig. 16A). In this model, most of the earthquakes in the southeastern cluster fell in regions with no potential for generating earthquakes, i.e., within negative $\Delta\sigma_D'$. At 15 km depth, the southwestern cluster of earthquakes consistently deviate from the locations of the rift faults in the model.

In the $\Delta\phi_{\max}'$, the SD65R25 model shows more stress rotation than the SD70R25 model (Fig. 16B). The regions of this stress rotations are at or near the location of the rift faults in both models although the different dip angles of the rift faults  create different  spacing between them with depth. The $\Delta\phi_{\max}'$ is skewed to the southeast in the direction perpendicular to the strike of the rift faults. The magnitude of the $\Delta\phi_{\max}'$ decreases away from the rift faults as observed in $\phi_{\max}$ plots.\selectlanguage{english}
\begin{figure}[ht]
\centering
\includegraphics[width=30pc]{Figures/SD65R25/F2_m3_c0_ch_diff_shmax.png}
\caption{{SD65R25}}
\label{figfour}
\end{figure}


\subsubsection*{LD70R25}
Despite the difference in the velocity models, the results of LD70R25 is similar to that of the SD70R25. The maximum value of $\Delta\sigma'_D$ and $\Delta\phi_{\max}'$ in both models are also comparable.\selectlanguage{english}
\begin{figure}[ht]
\centering
\includegraphics[width=30pc]{Figures/LD70R25/F13_LAM_m3_c0_ch_diff_shmax.png}
\caption{{LD70R25}}
\label{figfour}
\end{figure}



\subsubsection*{SD70R50}
The region of negative $\Delta\sigma'_{\max}$ in SD70R50 is smaller compared to the size in the SD70R25 model and does not correlate with the seismicity of the CSZ (Fig. 18A). The $\Delta\sigma_{\max}'$ decrease to 27$\%$ due to the increased strength of the crater. Also, the $\Delta\phi_{\max}'$ show a more clockwise rotation within the crater region than that of the SD70R25 (Fig. 18B). Specifically, the non-zero value of $\Delta\phi_{\max}'$ is absent at 10 km depth.\selectlanguage{english}
\begin{figure}[ht]
\centering
\includegraphics[width=30pc]{Figures/SD70R50/F1_r05_m3_c0_ch_diff_shmax.png}
\caption{{SD70R50}}
\label{figfour}
\end{figure}



\subsubsection*{SD70R100}
The $\Delta\sigma_{\max}'$ in the SD70R100 model predicts earthquakes along the rift faults which do not explain observed seismicity (Fig. 19A). At 5 and 10 km depths, the earthquakes fall at regions with a negative value of $\Delta\sigma_{\max}'$. The $\Delta\phi_{\max}'$ does not correlate with observed stress orientation (Fig. 19B). The value of $\phi_{\max}$ decrease to a near-zero value at 15 km depth. The values of $\Delta\phi_{\max}$ became negative in the fault region at 20 and 25 km.\selectlanguage{english}
\begin{figure}[ht]
\centering
\includegraphics[width=30pc]{Figures/SD70R100/F1_r10_m3_c0_ch_diff_shmax.png}
\caption{{SD70R100}}
\label{figfour}
\end{figure}



\subsubsection*{TD70}
The results in the TD70 is similar to the results when the crater-crust moduli ratio is 1.0 (Fig. 20). The $\Delta\sigma_{\max}'$ also predicts earthquakes along the rift faults which do not explain observed seismicity. The $\Delta\phi_{\max}'$ with the TD70 does not correlate with observed stress orientation.\selectlanguage{english}
\begin{figure}[ht]
\centering
\includegraphics[width=30pc]{Figures/TD70/F5_T_m3_c0_ch_diff_shmax.png}
\caption{{TD70}}
\label{figfour}
\end{figure}












\par\null\par\null\par\null\par\null\par\null\par\null

\subsection*{Model Stopping Criterion~}\label{model-stopping-criterion}

We stop the simulation when the differential stress (\(\sigma_D=\sigma_1-\sigma_3\))
profile conforms with the Byerlee's law for brittle
rocks~\citep*{Byerlee1978}. For a given failure envelope
(\(\tau=C_o\ +\ \mu\sigma_n\)),~\(\sigma_D\) can be determined
using~\(\sigma_D=\frac{\left(2\sin[\tan^{-1}(\mu)](C_o+\mu\sigma_3)\right)}{\mu(1-\sin[\tan^{-1}(\mu)])}\), where~\(\mu\) is the friction
coefficient, Co~is the cohesion, and~\(\sigma_3\) is the minimum
principal stress. In our model,~\(\sigma_3\) corresponds to the
lithostatic stress. The~\(\sigma_3\) is 250 MPa at a depth of 10
km if the average density of crustal rock is 2500 kg/m3. For a Byerlee's
law (\(\tau\) = 50 MPa +
0.6\(\sigma_n\),~\citet*{Byerlee1978}), the~\(\sigma_D\)~is
706 MPa at 10 km depth without making any assumptions in the
intermediate principal stress. We will use this as the stopping
criterion for our simulations.

\par\null\par\null\par\null

To differentiate the model names, we adopt ``model\_moduli ratio'' e.g.
F1\_0.25 is a F1 model with a crater-crust moduli ratio of 0.25.

\par\null\par\null\par\null\par\null

~

\selectlanguage{english}
\FloatBarrier
\nocite{*}

\bibliographystyle{plainnat}
\bibliography{bibliography/converted_to_latex.bib%
}

\end{document}

