%% Copernicus Publications Manuscript Preparation Template for LaTeX Submissions
%% ---------------------------------
%% This template should be used for the following class files: copernicus.cls, copernicus2.cls, copernicus_discussions.cls
%% The class files, the Copernicus LaTeX Manual with detailed explanations regarding the comments
%% and some style files are bundled in the Copernicus Latex Package which can be downloaded from the different journal webpages.
%% For further assistance please contact the Publication Production Office (production@copernicus.org).
%% http://publications.copernicus.org


%% Differing commands regarding the specific class files are highlighted.


%% copernicus.cls
%\documentclass[se]{copernicus}

%% copernicus2.cls
%%\documentclass[journal abbreviation]{copernicus2}

%% copernicus_discussions.cls
%\documentclass[sed, hvmath, online]{copernicus_discussions}        
\documentclass[sed, hvmath]{copernicus_discussions}  %% ? tolto 'online', messo come esempio nel manuale Aronne

%% added Aronne
\usepackage{xcolor}
\usepackage{epstopdf}
\bibliographystyle{copernicus_discussions}
%% commands from original paper:
\newcommand{\ii}{{\rm i}}
\newcommand{\myfrac}[2]{\displaystyle \frac{#1}{#2}}
\newcommand{\BM}[1]{\mbox{\boldmath $ #1 $}}
\newcommand{\fo}{\rm}
%% end commands from original paper

\begin{document}

\linenumbers

\title{Simulation of seismic waves at the Earth crust (brittle-ductile transition) based on the Burgers model \\
   \large{Short title: Wave simulation based on the Burgers model}}


\author[]{J.~M.~Carcione}
\author[]{F.~Poletto}
\author[]{B.~Farina}
\author[]{A.~Craglietto}

\affil[]{Istituto Nazionale di Oceanografia e di Geofisica Sperimentale (OGS), 
Borgo Grotta Gigante 42c, 34010 Sgonico, Trieste, Italy.}

%% The [] brackets identify the author to the corresponding affiliation, 1, 2, 3, etc. should be inserted.

%\runningtitle{Simulation of seismic waves at the Earth crust (brittle-ductile transition) based on the Burgers model
%Short title: Wave simulation based on the Burgers model}
\runningtitle{Wave simulation based on the Burgers model}

\runningauthor{J.~M.~Carcione et al.}

\correspondence{Jos\'e M. Carcione\\ (jcarcione@inogs.it)}


\received{}
\pubdiscuss{} %% only important for two-stage journals
\revised{}
\accepted{}
\published{}

%% These dates will be inserted by the Publication Production Office during the typesetting process.

\firstpage{1}

\maketitle  %% Please note that for the copernicus2.cls this command needs to be inserted after \abstract{TEXT}


\begin{abstract}
The Earth crust presents two dissimilar rheological behaviours depending on the in-situ stress-temperature conditions. 
The upper,  cooler, part is brittle while deeper zones are ductile. Seismic waves may reveal the presence of 
the transition but a proper characterization is required. We first obtain a stress-strain relation including the 
effects of shear seismic attenuation and ductility due to shear deformations and plastic flow.
The anelastic behaviour is based on the Burgers mechanical model to describe the effects of seismic attenuation 
and steady-state creep flow. The shear Lam\'e constant of the brittle and ductile media depends on the in-situ stress 
and temperature through the shear viscosity, which is obtained by the Arrhenius equation and the octahedral stress 
criterion. 
The P- and S-wave velocities decrease as depth and temperature increase due to the geothermal gradient, an effect which 
is more pronounced for shear waves. 
We then obtain the P-S and SH equations of motion recast in the velocity-stress formulation, including memory variables 
to avoid the computation of time convolutions. The equations correspond to isotropic anelastic and inhomogeneous media 
and are solved by a direct grid method based on the Runge-Kutta time stepping technique and the Fourier pseudospectral 
method.
The algorithm is tested with success against known analytical solutions for different shear viscosities. A realistic 
example illustrates the computation of surface and reverse-VSP synthetic seismograms  in the presence of an abrupt 
brittle-ductile transition. 
%% keywors not requested 
% \keywords{brittle-ductile transition \and Burgers model \and wave simulation \and attenuation \and Fourier method.}
\end{abstract}

%% only used for copernicus2.cls
%%\abstract{
%% TEXT
%% \keywords{TEXT}}

\introduction  %% \introduction[modified heading if necessary]
The seismic characterization of the brittle-ductile transition (BDT) is essential in earthquake seismology and 
geothermal studies, since it plays an important role in determining the nature and nucleation depth of earthquakes 
\citep{meissner1982,zappone1994,simpson1999} and the availability of geothermal energy \citep{manzella1998}.
The BDT in the Earth is generally viewed as a transition between
two different constitutive behaviours, viscoelastic and plastic \citep{dragoni1990}.
There is evidence that the K horizon in the upper crust of Central Italy corresponds to a shear plane separating 
the brittle crust from the ductile crust \citep{brogi2003}.

The viscosity of the crust is a fundamental factor in defining the properties of the BDT interface.   
The contrast in properties at the transition is mainly due to the dissimilar shear rigidity with much lower values in 
the ductile medium \citep{matsumoto1996}.
The ductile medium mainly flows when subjected to deviatoric stress, while does not show major flow under 
hydrostatic stress. 
Then, the flow deformation is mainly associated with the shear modulus of the medium.
The deviatoric stress is determined by the octahedral stress, $\sigma_o$,
a scalar that is invariant under coordinate transformations
and whose value determines the character of the flow.
{\fo 
When the stress vector associated with the normal to the octahedral plane
is generated, its components in the principal directions are the eigenstresses
(or principal stresses). Alternatively, it has two components -- one normal to
the plane (which has a magnitude
equal to the mean stress) and one tangential to
the plane which has a magnitude equal to the
octahedral stress (and the latter is
proportional to the magnitude of the deviatoric stress).
The rock starts to yield when $\sigma_o$ exceeds the elastic
octahedral-stress limit. Below this limit, there is
gradual creep deformation when constant stress is applied. Then, if
$\sigma_o$ is lower than the elastic limit, the material
follows a viscoelastic stress-strain relation. If $\sigma_o$ exceeds 
this limit, steady-state
flow and failure occurs (Carcione and Poletto, 2013).
}

The flow viscosity is a function of temperature and pressure, determined by the geothermal gradient and the lithostatic 
stress, respectively.  
An alternative constitutive equation is proposed by \citet{hueckel1994}, based on a thermo-plasticity theory, 
where the elastic domain is postulated as temperature dependent, shrinking with temperature. 
On the other hand, \citet{arcay2012} proposes a thermomechanical model based on a non-Newtonian viscous rheology and
a pseudo-brittle rheology.

It is widely accepted that linear visco-elastic-plastic models are appropriate to describe the
behavior of ductile media. \citet{gangi1981,gangi1983} used this type of models to fit data for synthetic and 
natural rocksalt.
The viscoelastic creep of salt has been described with a Burgers model by \citet{carcione2006}.
\citet{carcione2013} used the Burgers model to describe the BDT transition, including the presence of anisotropy 
and seismic attenuation.
The Burgers model is shown in Figure 1 (the Maxwell and Zener model are particular cases of this model).
Here, the Zener model is used to model the viscoelastic motion with no plastic flow, obtained as the limit of
infinite plastic viscosity. Seismic wave losses are solely due to shear deformations. 

Computational geophysics is essential to study the structure of the Earth on the basis of seismic forward modelling, 
inversion and interpretation \citep{zappone1994,long2009,juhlin2011}.
In particular, grid methods are required to simulate wave propagation in heterogeneous realistic models 
\citep[e.g.,][]{carcione2002,seriani1992}.
In this work, we propose to simulate seismic wave propagation in heterogeneous media involving the brittle-ductile 
transition. The differential equations are solved in the time domain by using memory variables \citep{carcione2007}.
We assume isotropic media and plane strain conditions and obtain the differential equations of motion 
for 2D P-S and SH waves. The equations are recast in the velocity-stress formulation, requiring eight and four memory 
variables in the first and second cases when using one shear relaxation mechanism.
The equations are solved by a direct grid method based on the Runge-Kutta and the Fourier methods, corresponding to the 
time and spatial discretizations \citep[e.g.,][]{carcione2007}.

\section{The Burgers mechanical model}

The constitutive equation including both the viscoelastic and ductile behaviour, can be written as a generalization of 
the 1D stress-strain relation reported by \citet{dragoni1990} and \citet{dragoni1991} to the 3D anelastic case, 
replacing the Maxwell model by the Burgers model \citep{carcione2006,carcione2007,carcione2013}.

The Burgers model is a series connection of a dashpot and a Zener model as can be seen in Figure 1. The usual expression 
in the time domain is the creep function 
\begin{equation} \label{6}
\chi = \left( \frac{t}{\eta} + \frac{1}{\mu_0} \left[ 1 - \left(  1 - \frac{\tau_\sigma}{\tau_\epsilon}  \right) \exp (- 
t / \tau_{\epsilon} )  \right] \right) H (t)  
\end{equation}
\citep{carcione2006,chauveau2008}, where $t$ is time and $H(t)$ is the Heaviside function. 
The quantities $\tau_{\sigma}$ and $\tau_{\epsilon}$ are seismic relaxation times, 
$\mu_0$ is the relaxed shear modulus (see below) and $\eta$ is the flow viscosity describing the ductile behaviour 
related to shear deformations.
The frequency-domain shear modulus $\mu$ can be obtained as 
$\mu = \left[ {\cal F} ( \dot \chi ) \right]^{-1}$, where ${\cal F}$ denotes time Fourier transform and a dot above a 
variable denotes time derivative. It gives 
\begin{equation} \label{1}
\mu = \frac{\mu_0 ( 1 + \ii \omega \tau_{\epsilon})} 
{ 1 + \ii \omega \tau_{\sigma} - \myfrac{\ii \mu_0}{ \omega \eta} ( 1 + \ii \omega \tau_{\epsilon}) } , 
\end{equation}
where $\ii = \sqrt{-1}$ and $\omega$ is the angular frequency. 
The relaxation times can be expressed as 
\begin{equation} \label{2}
\tau_{\epsilon} = \frac{\tau_0}{Q_0}
\left( \sqrt{ {Q_0}^2 +1} +1 \right) ,
\ \ \
\tau_{\sigma} =
\tau_{\epsilon} - \frac{2 \tau_0}{Q_0} ,
\end{equation}
where $\tau_0$ is a relaxation time such that $\omega_0 = 1/\tau_0$ is the center frequency
of the relaxation peak and $Q_0$ is the minimum quality factor. 
{\fo 
The dependence of the qualiy factor as a function of frequency is given below in equation (\ref{n5}).
}

The limit $\eta \rightarrow \infty$ in equation (\ref{1}) recovers the Zener kernel to describe the behaviour of the 
brittle material, while $\tau_{\sigma} \rightarrow 0$ and $\tau_{\epsilon} \rightarrow 0$ yield the Maxwell model used 
by \citet{dragoni1990} and \citet{dragoni1991}:
\begin{equation} \label{3}
\mu = \mu_0 \left( 1 - \frac{ \ii \mu_0}{\omega \eta} \right)^{-1}
\end{equation}
\citep[e.g.,][]{carcione2007}. For $\eta \rightarrow 0$, $\mu \rightarrow 0$ and the medium becomes a fluid. 
Moreover, if $\omega \rightarrow \infty$, $\mu \rightarrow \mu_0 \tau_\epsilon / \tau_\sigma$ and 
$\mu_0$ is the relaxed ($\omega$ = 0) shear modulus of the Zener element  ($\eta = \infty$). 

The viscosity $\eta$ can be expressed by the Arrhenius equation \citep[e.g.,][]{carcione2006,montesi2007}.
It is related to the steady-state creep rate $\dot \epsilon$ by
\begin{equation} \label{4}
\eta = \frac{\sigma_o}{2 \dot \epsilon} ,
\end{equation}
where $\sigma_o$ is the octahedral stress. 
The creep rate can be expressed as
\begin{equation} \label{5}
\dot \epsilon =  A  \sigma_o^{n} \exp ( - E / R T )
\end{equation}
\citep[e.g.,][]{gangi1983,carcione2006,carcione2013}, where $A$ and $n$ 
are constants, $E$ is the activation energy, $R$ = 8.3144 J/mol/$^{\rm o}$K  is the gas constant and $T$
is the absolute temperature. The
form of the empirical relation (\ref{5}) is
determined by performing experiments at different strain rates,
temperatures and/or stresses \citep[e.g.,][]{gangi1983,carter1983}.

In order to obtain the equations of motion to describe wave propagation it is convenient to consider the Burgers 
relaxation function 
\begin{equation} \label{19}
\psi (t ) = [ A_1 \exp ( - t / \tau_1 ) - A_2 \exp (-t / \tau_2 )] H (t) 
\end{equation}
\citep{carcione2007}, where
\begin{equation} \label{20}
\tau_{1,2} = - \frac{1}{\omega_{1,2}} \ \ \ {\rm and} \ \ \
A_{1,2} = \myfrac{\mu_1 \mu_2 + \omega_{1,2} \eta_1 \mu_2}{\eta_1 ( \omega_1 - \omega_2)} .
\end{equation}
and 
\begin{equation} \label{21}
(2 \eta \eta_1) \omega_{1,2} = - b \pm \sqrt { b^2 - 4 \mu_1 \mu_2 \eta \eta_1},  
\ \ \ \
b = ( \mu_1 + \mu_2 ) \eta + \mu_2 \eta_1 .
\end{equation}
In terms of the relaxation times and $\mu_0$, it is
\begin{equation} \label{22}
\mu_1 = \myfrac{\mu_0 \tau_\epsilon}{ \tau_\epsilon - \tau_\sigma} , \ \ \
\mu_2 = \mu_0 \myfrac{\tau_\epsilon}{\tau_\sigma} , \ \ \
\eta_1 = \mu_1 \tau_\epsilon .
\end{equation}
The complex shear modulus is 
\begin{equation} \label{23}
\mu = {\cal F} (\dot \psi ) = \ii \omega \left(  \frac{A_1 \tau_1}{1 + \ii \omega \tau_1} - \frac{A_2 \tau_2}{1 + \ii 
\omega \tau_2}\right) .
\end{equation}
It can be verified that equations (\ref{1}) and (\ref{23}) coincide. 

\section{2D Propagation of P-S waves}

Let us consider plane-strain conditions and propagation in the (x,z)-plane. 
The simplest P-S stress-strain relation, with shear loss and flow, are
\begin{equation} \label{24}
\begin{array}{l}
\dot \sigma_{xx} = \lambda (\partial_x v_x + \partial_z v_z ) + 2 \dot \psi \ast \partial_x v_x , \\ \\
\dot \sigma_{zz} = \lambda (\partial_x v_x + \partial_z v_z ) + 2 \dot \psi \ast \partial_z v_z , \\ \\
\dot \sigma_{xz} = \dot \psi \ast (\partial_x v_z + \partial_z v_x )   
\end{array}
\end{equation}
\citep[e.g.,][]{carcione2007}, where $\lambda$ is a Lam\'e constant, 
$\sigma$ are stress components, $v$ are particle-velocity components,
$\partial_i$ indicates a spatial derivative with respect to the variable $x_i$, $i$ = 1,2,3 ($x_1 = x$, $x_2 = y$ and $x_3$ = z), 
and ``$\ast$" denotes time convolution.

The convolutions have the form $\dot \psi \ast \partial_i v_j$ and can be overcome by introducing memory variables. 
We obtain, 
\begin{equation} \label{25}
\dot \psi \ast \partial_i v_j = A_1 ( \partial_i v_j + e_{ij}^{(1)} ) - A_2 ( \partial_i v_j + e_{ij}^{(2)} ) , \ \  i,j 
=1,3
\end{equation}
where 
\begin{equation} \label{26}
e_{ij}^{(m)}= g_m H \ast  \partial_i v_j, \ \ \  g_m = - \frac{1}{\tau_m} \exp ( - t / \tau_m ) , \ \ m = 1,2 ,
\end{equation}
which satisfies 
\begin{equation} \label{27}
\dot e_{ij}^{(m)}= -\frac{1}{\tau_m} ( \partial_i v_j + e_{ij}^{(m)} ) .
\end{equation}
The stress-strain relation becomes
\begin{equation} \label{28}
\begin{array}{l}
\dot \sigma_{xx} = ( \lambda + 2 A_1 -2 A_2) \partial_x v_x + \lambda \partial_z v_z  + 2 (A_1 e_{xx}^{(1)} - A_2 
e_{xx}^{(2)} ) , \\ \\
\dot \sigma_{zz} = ( \lambda + 2 A_1 -2 A_2) \partial_z v_z + \lambda \partial_x v_x  + 2 (A_1 e_{zz}^{(1)} - A_2 
e_{zz}^{(2)} ) , \\ \\
\dot \sigma_{xz} = (A_1 - A_2) (\partial_x v_z + \partial_z v_x ) + A_1 e_{xz}^{(1)} - A_2 e_{xz}^{(2)} + A_1 
e_{zx}^{(1)} - A_2 e_{zx}^{(2)} . 
\end{array}
\end{equation}

On the other hand, the dynamical equations of motion are
\begin{equation} \label{29}
\begin{array}{l}
\dot v_x =  \myfrac{1}{\rho} ( \partial_ x \sigma_{xx}  + \partial_ z \sigma_{xz} )  + s_x , \\ \\
\dot v_z =  \myfrac{1}{\rho} ( \partial_ x \sigma_{xz}  + \partial_ z \sigma_{zz} )  + s_z 
\end{array}
\end{equation}
\citep{carcione2007}, where $\rho$ is the mass density and $s_i$ are source components.

The equations of motion are given by equations (\ref{27}), (\ref{28}) and (\ref{29})  in the unknown vector 
${\bf v} = (v_x , v_z , \sigma_{xx}, \sigma_{zz}, \sigma_{xz},  e_{ij}^{(m)} )^\top$. In matrix notation 
 \begin{equation} \label{30}
\dot {\bf v} = {\bf M} \cdot {\bf v} +{\bf s}  ,
\end{equation}
where ${\bf M}$ is a 13 $\times$ 13 matrix containing the material properties and spatial derivatives. 

In view of the correspondence principle \citep[e.g.,][]{carcione2007}, the complex and frequency-dependent 
P- and S-wave velocities are
\begin{equation} \label{n2}
v_P (\omega ) = \sqrt{ \frac{\lambda + 2 \mu (\omega )}{\rho}} , \ \ \ {\rm and} \ \ \
v_S (\omega) = \sqrt{ \frac{\mu (\omega )}{\rho}} , 
\end{equation}
respectively. 

For homogeneous waves in isotropic media, the phase velocity and attenuation 
factors are given by 
\begin{equation} \label{n3}
c_{P(S)} = \left[ {\rm Re} \left( \frac{1}{v_{P(S)}} \right) \right]^{-1}  
\end{equation}
and 
\begin{equation} \label{n4}
\alpha_{P(S)}  = - \omega {\rm Im} \left( \frac{1}{v_{P(S)}} \right) ,  
\end{equation}
and the P- and S-wave quality factors are given by 
\begin{equation} \label{n5}
Q_{P(S)} =  \frac{{\rm Re}(v_{P(S)}^2)}{{\rm Im}(v_{P(S)}^2)}   
\end{equation}
\citep[e.g.,][]{carcione2007}.

\section{Propagation of SH waves}

The stress-strain relations describing shear motion in the (1,3)-plane are
\begin{equation} \label{7}
\begin{array}{l}
\dot \sigma_{xy} = \dot \psi \ast \partial_x v_y  , \\  \\
\dot \sigma_{zy} = \dot \psi \ast \partial_z v_y 
\end{array}
\end{equation}
\citep{carcione2007}.

On the other hand, the dynamical equation of motion is 
\begin{equation} \label{15}
\dot v_y =  \frac{1}{\rho} ( \partial_ x \sigma_{xy}  + \partial_ z \sigma_{zy} )  + s ,
\end{equation}
where $s$ is the source \citep{carcione2007}.

Applying the same procedure as in the P-S case we obtain
\begin{equation} \label{16}
\begin{array}{l}
\dot \sigma_{xy} = (A_1 - A_2 )  \partial_x v_y  + A_1 e_{xy}^{(1)} - A_2 e_{xy}^{(2)}  , \\ \\
\dot \sigma_{zy} = (A_1 - A_2 )  \partial_z v_y  + A_1 e_{zy}^{(1)} - A_2 e_{zy}^{(2)}  ,  
\end{array}
\end{equation}
with 
\begin{equation} \label{17}
\begin{array}{l}
\dot e_{xy}^{(m)}= -\myfrac{1}{\tau_m} ( \partial_x v_y + e_{xy}^{(m)} ) , \\ \\
\dot e_{zy}^{(m)}= -\myfrac{1}{\tau_m} ( \partial_z v_y + e_{zy}^{(m)} ) , \ \ \ m=1,2 .
\end{array}
\end{equation}
The equations of motion are given by equations (\ref{15}), (\ref{16}) and (\ref{17})  in the unknown vector 
${\BM \sigma} = (v_y , \sigma_{xy}, \sigma_{zy},  e_{xy}^{(m)}, e_{zy}^{(m)} )^\top$ and can be recast 
as equation (\ref{30}) with matrix ${\bf M}$ of dimension 7 $\times$ 7. 
                                                                      
\section{Numerical solution}

The formal solution  to equation (\ref{30}) with zero initial conditions is given by 
\begin{equation} \label{n1}
{\bf v} (t) = \int_0^t \exp ( \tau {\bf M} ) \cdot {\bf s} (t - \tau) d \tau , 
\end{equation}
where $\exp ( t {\bf M})$ is the evolution operator of the system. The numerical solution is based on a Taylor expansion 
of this operator up to the fourth order, and the Runge-Kutta algorithm is used \citep{Jain1984,carcione2007}. 

The spatial derivatives
are calculated with the Fourier pseudospectral method \citep[]{Kosloff1982,carcione2007}. 
This method consist of a spatial discretization and 
calculation of spatial derivatives using the fast Fourier transform. 
It is a collocation technique in which a continuous
function is approximated by a truncated series of trigonometric functions,  
wherein the spectral (expansion) coefficients are chosen such that
the approximate solution coincides with the exact solution at the discrete set
of sampling or collocation points.
The collocation points are defined by equidistant sampling points.
Since the expansion functions are periodic, the Fourier method
is appropriate for problems with periodic 
boundary conditions.  
The method is infinitely accurate up to the 
maximum wavenumber of the mesh, that corresponds to a spatial wavelength of 
two grid points. 

The stability condition of the Runge-Kutta method is 
$c_{\rm max} \ dt / d_{\rm min} < 2.79$, where $c_{\rm max}$ is the maximum phase velocity, 
$dt$ is the time step and $d_{\rm min}$ is the minimum grid spacing.

\section{Examples}

First, we test the numerical code against an analytical solution for P-S waves in homogeneous media (Appendix A). 
To compute the transient responses, we use a Ricker wavelet of the form: 
\begin{equation} \label{Rick}
w (t) = \left ( a - \frac{1}{2} \right) \exp ( - a) , \ \ \
a = \left[ \frac{ \pi (t - t_s) }{t_p} \right]^2 , 
\end{equation}
where $t_p$ is the period of the wave (the distance between the side peaks is $\sqrt{6} t_p/ \pi$) and 
we take $t_s = 1.4 t_p$. Its frequency spectrum is 
\begin{equation} \label{Rickw}
W ( \omega ) = \left( \frac{t_p}{\sqrt{\pi}} \right) \bar a \exp ( - {\bar a} - \ii \omega t_s )  , \ \ \
\bar a = \left( \frac{\omega}{\omega_p} \right)^2 , \ \ \ \omega_p = \frac{2 \pi}{t_p} . 
\end{equation}
The peak frequency is $f_p = 1/t_p$. 

The rock is described by an unrelaxed P-wave velocity of 
$c_P$ = 6 km/s. Considering a Poisson medium, we obtain an S-wave velocity of $c_S$ = 3.464 km/s.
We have $\lambda + 2 \mu_0 = \rho c_P^2$ and $\mu_0 = \rho c_S^2$.
Assuming a density of $\rho$ = 2700 kg/m$^3$, we have 
$\lambda = \mu_0$ = 32.4 GPa. The seismic quality factor is $Q_0$ = 40 and $\omega_0 = 2 \pi f_p$. 
The numerical mesh has 231 $\times$ 231 grid points and a grid spacing of $dx$ = $dz$ = 30 m. The source is a vertical 
force with $f_p$ = 10 Hz and the receiver is located at  $x$ = $z$ = 1.2 km from the source. 
The solution is computed using a time step $dt$ = 1 ms.  Figure 2 shows the comparison between the numerical and 
analytical PS-wave solutions for $\eta =10^{20}$ Pa s (a-b) and  $\eta =10^9$ Pa s (c-d), where (a) and (c) correspond 
to $v_x$ and (b) and (d) to $v_z$.
At the source peak frequency, the P- and S-wave quality factors for $\eta = 10^{20}$ Pa s are 60 and 40, while 
those corresponding to $\eta = 10^9$ Pa s  are 2.9 and 1.8, respectively, i.e., very strong attenuation. 
The results for the SH wave are displayed in Figure 3. In this case, we assume $Q_0 = \infty$, i.e., attenuation is 
solely due to the plastic viscosity. The quality factor for $\eta = 2 \times 10^{10}$ Pa s is 39. 

Next, we present examples showing the capabilities of the full-waveform simulation algorithm
for the seismic characterization of crustal rocks at very high temperatures, 
as those encountered at depths where 
hydrothermal fluids are present at supercritical conditions \citep{albert2003}. 
The temperature dependence is
expressed by the Arrhenius steady-state power law, e.g., Eq. 15 in \citet{carcione2013}.
We use realistic Arrhenius constants $A=10^{30}$ MPa$^{-n}$s$^{-1}$, $n=3.5$ and
activation energy $E=990$ kJ/mol
for the rheology of the Icelandic crust \citep{violay2010,violay2012}. 
These parameters were determined by mechanical observations in laboratory stress-deformation experiments
performed at different confining pressures with
glass-free basaltic rock samples (GBF) \citep{violay2010,violay2012},
in agreement with the results of \citet{hacker1992}.
In the absence of glass (silica)---that strongly influences the ductile behaviour of the basaltic rock at lower 
temperatures---the
glassy-free basalts present rapid BDT variation at high temperature.
Figure 4 shows the temperature profile used in the calculation, with a steep gradient at depth,
similar to that reported in \citet{foulger1995}.

The purpose of the numerical experiment is to predict P-  and S-wave propagation 
with dispersion and attenuation at high 
temperatures \citep{carcione2013}, and study the observability of the BDT
by seismic reflection methods.
This issue poses the problem to evaluate suitable transition zones and gradients 
for the crustal rock properties
to get reflections in the frequency range typical of seismic exploration. 
In the following examples, 
we assume propagation in a uniform, isotropic medium, under lithostatic pressure in the absence of tectonic stress.
Figures 5 and 6 show the phase velocities calculated at a frequency of 10 Hz
using the approach of \citet{carcione2013}, 
with an average lithostatic confining pressure of approximately 95 MPa
for a Poisson medium with unrelaxed P-wave velocity of 6 km/s, S-wave velocity of 3.464 km/s and density $\rho = 2600$ 
kg/m$^3$. 
The transition in the
GBF velocity functions is quite rapid at approximately 3.55 km in Figure 6,
corresponding to $T$ = 1120 $^\circ$C. Complete melting is obtained above 1300 $^\circ$C. 
The sharp variation of the velocity in the transition zone (Figure 7) is due to the combined effect of the rapid
velocity variation in the BDT zone (Figure 5) and to the steep temperature gradient (Figure 4).
Figure 8 shows the dispersion of the P-wave velocity calculated at frequencies of 3, 10 and 30 Hz.

The full-waveform synthetic data are calculated using a vertical source with a Ricker wavelet with maximum frequency of 
50 Hz.
The source is assumed to be located at the free surface and the grid spacing of the numerical mesh is 10 m. 
The signals are recorded by a horizontal line of
receivers at the surface (shot gather) and by a vertical array of receivers. The latter experiment simulates a 
seismic-while-drilling reverse-VSP experiment \citep{poletto2004,poletto2011}. 
Figure 9 shows the synthetic common-shots, where (a) corresponds to the vertical component and (b) to 
the horizontal component. 
In both shot gathers, we can see the reflections of the transition interface for P waves (RP) (a) and S waves (RS) (b).
We remark the fact that the only change in the model is the temperature profile. 
Figure 10 shows the VSP recorded at zero offset, from 0.2 km depth to the bottom of the model.
We can observe the direct transmitted arrival (TP) in the deeper melted zone below 3.5 km,
and the reflection of  the P wave (RP) at approximately 3.5 km depth.
Finally, Figure 11 shows the $f$-$k$ plot of the VSP signal calculated every 10 m in depth. 
This plot confirms that
the dispersion in this example is moderate due to the thinness of the GBF transient zone (cfr. Figure 7).
Higher dispersion can be expected in glassy basalts (GB) \citep{violay2010,violay2012}.

\conclusions
The upper -- cooler -- part of the crust is brittle, while deeper zones present ductile behaviour. In
some cases, this brittle-ductile transition is a single seismic
reflector with an associated reflection coefficient. The stress-strain relation and its physical implications have been 
analysed in a previous work. Here, we have developed a full-waveform algorithm to simulate temperature-dependent 
propagation of
seismic waves in geothermal and magmatic crustal rheologies in general and the brittle-ductile transition in particular.
This abrupt transition is believed to be
the lower limit of seismicity and may be an indication of
geothermal activity, since its reflectivity may reveal the
presence of partial melting and/or overpressured fluids.

The method uses the Burgers viscoelastic model and the Arrhenius equation to calculate the flow viscosity. 
{\fo Existing viscoelastic codes, based on the Maxwell, Kelvin-Voigt and Zener models, cannot be used, because they fail to model both the brittle and ductile behaviours.
}
The time convolutions appearing in the stress-strain relations are circumvented by introducing memory variables and the numerical 
algorithm is based on the Fourier pseudospectral method to compute the spatial derivatives. The modeling technique, 
developed for P-SV and SH waves, is successfully tested against known analytical solutions. 

The example demonstrate the observability of the brittle-ductile transition using surface-seismic and VSP methods under 
appropriate conditions. Similar simulations, using this forward modeling technique, can be performed, including 
estimations of P/S velocity relations, dispersion and attenuation related to temperature profiles, temperature-gradient 
variations, pressure conditions and many aspects of rock physics in hot geothermal reservoirs with hydrothermal fluids 
and inhomogeneous media. 

%\section{HEADING}
%TEXT

%\subsection{HEADING}
%TEXT

%\subsubsection{HEADING}
%TEXT

%\conclusions  %% \conclusions[modified heading if necessary]
%TEXT

\appendix
\section{\\ \\ \hspace*{-7mm} Analytical solution in 2D homogeneous media}        %% Appendix A
The P-S Green's function corresponding to the wave field generated by an impulsive vertical 
force of strength $F_0$ is given by
\begin{equation} \label{Eas1w}
G_x ( x, z, \omega , v_P , v_S ) =  \left( \frac{F_0}{2 \pi \rho} \right) \frac{x z}{r^2} 
[ F_1 (x, z, \omega, v_P , v_S ) + F_3 (r, \omega ,  v_P , v_S) ] ,
\end{equation}
\begin{equation} \label{Eas3w}
G_z ( x, z, \omega , v_P , v_S ) =  \left( \frac{F_0}{2 \pi \rho} \right) \frac{1}{r^2} 
[ z^2 F_1 (x, z, \omega, v_P , v_S ) - x^2 F_3 (r, \omega ,  v_P , v_S) ] ,
\end{equation}
where $r = \sqrt{x^2 + z^2}$, 
\begin{equation} \label{G1w}
F_1 ( r ,\omega , v_P , v_S ) =  \frac{ \pi \omega}{2} \left[ \frac{1}{v_P^2} H_0^{(2)} \left( 
\frac{\omega r}{v_P} \right) + \frac{1}{\omega r v_S} H_1^{(2)} \left( \frac{\omega r}{v_S} \right) - 
\frac{1}{\omega r v_P} H_1^{(2)} \left( \frac{\omega r}{v_P} \right) \right] , 
\end{equation}
\begin{equation} \label{G3w}
F_3 ( r ,\omega , v_P , v_S ) =  - \frac{\pi \omega}{2} \left[ \frac{1}{v_S^2} H_0^{(2)} \left( 
\frac{\omega r}{ v_S} \right) - \frac{1}{\omega r v_S} H_1^{(2)} \left( \frac{\omega r}{v_S} \right) + 
\frac{1}{\omega r v_P} H_1^{(2)} \left( \frac{\omega r}{v_P} \right) \right] , 
\end{equation}
and $H_0^{(2)}$ and $H_1^{(2)}$ are
the zero- and first-order Hankel functions of the second kind \citep{eason1956,carcione2007}. 

The 2D viscoelastic particle velocities can then be expressed as 
\begin{equation} \label{Gr2d1}
v_i ( x, z, \omega ) =  \left\{ 
\begin{array}{l}
W ( \omega) G_i ( x, z, \omega , v_P , v_S ) , \ \ \ \ \omega \ge 0 , \\ 
W^\ast ( \omega) G_i^\ast ( x, z,  - \omega , v_P , v_S ) , \ \ \ \ \omega < 0 ,  \ \ \ i = 1,3  
\end{array} 
\right.
\end{equation}
(hermiticity), where the superscript ``$\ast$" denotes complex conjugate and $G_x$ and $G_z$ are assumed to be 
zero at $\omega = 0$. A numerical inversion to the time domain by a discrete Fourier 
transform yields the desired time-domain solution.

On the other hand, the SH Green's function is
\begin{equation} \label{SH1}
G ( x, z, \omega , v_S ) =  \pi \omega H_0^{(2)} \left( 
\frac{\omega r}{v_S} \right)  
\end{equation}
\citep[][Section 6.4]{carcione2007}, and 
\begin{equation} \label{SH2}
v_y ( x, z, \omega ) =  \left\{ 
\begin{array}{l}
W ( \omega) G ( x, z, \omega , v_S ) , \ \ \ \ \omega \ge 0 , \\ 
W^\ast ( \omega) G^\ast ( x, z, - \omega  , v_S ) , \ \ \ \ \omega < 0  .  
\end{array} 
\right.
\end{equation}

%\subsection                               %% Appendix A1, A2, etc.

\begin{acknowledgements}
We thank Josep de la Puente for useful comments.
\end{acknowledgements}


\newpage

\begin{thebibliography}{}

\bibitem[Albertsson et al.(2003)]{albert2003}
Albertsson, A.J.,  Bjarnason, O., Gunnarsson, T., Balizus, C. and 
Ingason, K.: Part III: Fluid Handling and Evaluation,
In: Iceland Deep Drilling Project, Feasibility Report OS-2003-007, 33 pp., 2003.

\bibitem[Arcay(2012)]{arcay2012}
Arcay, D.: 
Dynamics of interplate domain in subduction zones: influence of
rheological parameters and subducting plate age,
Solid Earth, 3, 467--488, 2012.

\bibitem[Brogi et al.(2003)]{brogi2003}
Brogi, A., Lazzarotto, A., Liotta, D. and Ranalli, G.:
Extensional shear zones as imaged by reflection seismic lines: the
Larderello geothermal field (central Italy), Tectonophysics, 363, 127--139, 2003.

\bibitem[Carcione(2007)]{carcione2007}
Carcione, J.M.: Wave fields in real media: Wave
propagation in anisotropic, anelastic, porous and electromagnetic media. 
Handbook of Geophysical Exploration, vol. 38, Elsevier (2nd edition, revised and extended), 2007.
 
\bibitem[Carcione et al.(2006)]{carcione2006}
Carcione, J.M., Helle, H.B. and Gangi, A.F.: Theory of borehole stability when
drilling through salt formations, Geophysics, 71, F31--F47, 2006.

\bibitem[Carcione et al.(2002)]{carcione2002}
Carcione, J.M., Herman, G. and ten Kroode, F.P.E.:
Seismic modeling,
Geophysics, 67, 1304--1325, 2002.

\bibitem[Carcione and Poletto(2013)]{carcione2013}
Carcione, J.M. and Poletto, F.: Seismic rheological model and reflection coefficients
of the brittle-ductile transition, Pure and Applied Geophysics, 
DOI 10.1007/s00024-013-0643-4, 2013.

\bibitem[Carter and Hansen(1983)]{carter1983}
Carter, N.L. and Hansen, F.D.:
Creep of rocksalt, Tectonophysics, 92, 275--333, 1983.

\bibitem[Chauveau and Kaminski(2008)]{chauveau2008}
Chauveau, B. and Kaminski, E.:
Porous compaction in transient creep regime and
implications for melt, petroleum, and CO2 circulation, J. Geophys. Res, 113, B09406, doi:10.1029/2007JB005088, 2008.

\bibitem[Dragoni(1990)]{dragoni1990}
Dragoni, M.: Stress relaxation at the lower dislocation edge of great shallow earthquakes,
Tectonophysics, 179, 113--119, 1990.

\bibitem[Dragoni and Pondrelli(1991)]{dragoni1991}
Dragoni, M. and Pondrelli, S.:
Depth of the brittle-ductile transition in a transcurrent boundary zone,
Pure and Applied Geophysics, 135, 447--461, 1991.

\bibitem[Eason et al.(1956)]{eason1956}
Eason, G., Fulton, J. and Sneddon, I.N.:
The generation of waves in an infinite elastic solid by variable body
forces, Phil. Trans. Roy. Soc. London, Ser. A, 248, 575--607, 1956.

\bibitem[Foulger(1995)]{foulger1995}
Foulger, G.R.:
The Hengill geothermal area, Iceland -- variation of temperature gradients deduced from the maximum 
depth of seismogenesis, J. Volc. Geotherm. Res., 65, 119--133, 1995.

\bibitem[Gangi(1981)]{gangi1981}
Gangi, A.F.:
A constitutive equation for one-dimensional transient and steady-state flow of
solids. Mechanical Behavior of Crustal Rocks,
Geophysical Monograph AGU, 24, 275--285, 1981.

\bibitem[Gangi(1983)]{gangi1983}
Gangi, A.F.:
Transient and steady-state deformation of synthetic rocksalt,
Tectonophysics, 91, 137--156, 1983.

\bibitem[Hacker and Christie(1992)]{hacker1992}
Hacker, B.R. and Christie, J.M.:
Experimental deformation of glassy basalt,
Tectonophysics, 200, 79--96, 1992.

\bibitem[Hueckel et al.(1994)]{hueckel1994}
Hueckel, T., Peano, A. and Pellegrini, R.:
A constitutive law for thermo-plastic behaviour of rocks: an analogy with clays, 
Surv. Geophys., 15, 643--671, 1994.

\bibitem[Jain(1984)]{Jain1984}Jain, M. K.: 
Numerical solutions of differential equations, Wiley Eastern Ltd., 1984.  

\bibitem[Juhlin and Lund(2011)]{juhlin2011}
Juhlin, C. and Lund, B.:
Reflection seismic studies over the end-glacial Burtr\"ask fault,
Skellefte\aa, Sweden, 
Solid Earth, 2, 9--16, 2011. 

\bibitem[Kosloff and Baysal(1982)]{Kosloff1982}
Kosloff, D., and Baysal, E.: Forward modeling by the Fourier method: Geophysics, 47, 1402--1412, 1982. 

\bibitem[Long and Silver(2009)]{long2009}
Long, M.D. and Silver, P.G.:
Shear Wave Splitting and Mantle Anisotropy: Measurements, Interpretations, and New 
Directions,
Surv. Geophys., 30, 407--461, 2009.

\bibitem[Manzella et al.(1998)]{manzella1998}
Manzella, A., Ruggieri, G., Gianelli, G. and Puxeddu, M.:
Plutonic-geothermal systems of southern Tuscany: a review of the
crustal models, Mem. Soc. Geol. It., 52, 283--294, 1998.

\bibitem[Matsumoto and Hasegawa(1996)]{matsumoto1996}
Matsumoto, S. and Hasegawa, A.:
Distinct S-wave reflector in the midcrust beneath
Nikko-Shirane volcano in the northeastern Japan arc,
J. Geophys. Res., 101(B2), 3067--3083, 1996.

\bibitem[Meissner and Strehlau(1982)]{meissner1982}
Meissner, R. and Strehlau, J.:
Limits of stresses in continental crusts and their relation to
the depth-frequency distribution of shallow earthquakes, Tectonics, 1, 73--89, 1982.

\bibitem[Montesi(2007)]{montesi2007}
Montesi, L.G.J.:
A constitutive model for layer development in shear zones near the
brittle-ductile transition, Geophys. Res. Lett., 34, 
L08307, doi:10.1029/ 2007GL029250, 2007.

\bibitem[Poletto et al.(2011)]{poletto2011}
Poletto, F., Corubolo, P., Schleifer, A., Farina, B., Pollard, J. and Grozdanich, B.: 
Seismic while drilling for geophysical exploration in geothermal wells, GRC Transactions, 35, 1737--1741, 2011.

\bibitem[Poletto and Miranda(2004)]{poletto2004}
Poletto, F. and Miranda, F.:
Seismic while drilling. Fundamentals of drill-bit seismic for
exploration, Pergamon, Amsterdam, 2004.

\bibitem[Seriani et al.(1992)]{seriani1992}
Seriani, G., Priolo, E., Carcione, J.M. and Padovani, E.:
High-order spectral element method for elastic wave 
modeling: High-order spectral element method for elastic wave modeling, 
62nd Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1285--1288, 1992.
                      
\bibitem[Simpson(1999)]{simpson1999}
Simpson, F.: Stress and seismicity in the lower crust: a challenge
to simple ductility and implications for electrical conductivity
mechanisms, Surv. Geophys., 20, 201--227, 1999.

\bibitem[Violay et al.(2010)]{violay2010}
Violay, M., Gibert, B., Mainprice, D., Evans, B., Pezard, P.A., Flovenz, O.G.
and Asmundsson, R.: The brittle-ductile transition in experimentally deformed
basalt under oceanic crust conditions: Evidences for presence of permeable reservoirs at supercritical
temperatures and pressures in the Icelandic crust, Proceeding World Geothermal Congress, 
Bali, Indonesia, 25-29 April, 2010. 

\bibitem[Violay et al.(2012)]{violay2012}
Violay, M., Gibert, B., Mainprice, D., Evans, B., Dautria, J.M., Azais, P.
and Pezard, P.A.: An experimental study of the brittle-ductile transition of
basalt at oceanic crust pressure and temperature conditions,
Geophys. Res., 117, 1--23, 2012.

\bibitem[Zappone(1994)]{zappone1994}
Zappone, A.:
Calculated and observed seismic behaviour of exposed upper mantle,
Surv. Geophys.,  5, 629--642, 1994.

\end{thebibliography}


%% Literature citations
%% command                        & example result
%% \citet{jones90}|               & Jones et al.\ (1990)
%% \citep{jones90}|               & (Jones et al., 1990)
%% \citep{jones90,jones93}|       & (Jones et al., 1990, 1993)
%% \citep[p.~32]{jones90}|        & (Jones et al., 1990, p.~32)
%% \citep[e.g.,][]{jones90}|      & (e.g., Jones et al., 1990)
%% \citep[e.g.,][p.~32]{jones90}| & (e.g., Jones et al., 1990, p.~32)
%% \citeauthor{jones90}|          & Jones et al.
%% \citeyear{jones90}|            & 1990






%% FIGURES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% ONE-COLUMN FIGURES

%f
%\begin{figure}[t]
%\vspace*{2mm}
%\begin{center}
%\includegraphics[width=8.3cm]{fm1.eps}
%\end{center}
%\caption{TEXT}
%\end{figure}



%% TWO-COLUMN FIGURES

%f
%\begin{figure*}[t]
%\vspace*{2mm}
%\begin{center}
%\includegraphics[width=12cm]{fm2.eps}
%\end{center}
%\caption{TEXT}
%\end{figure*}

\newpage

\begin{figure}
\vspace{1cm}
\includegraphics[width=9cm,height=6cm]{fm1.eps}
\vspace{1cm}
\caption{
Mechanical representation of the Burgers viscoelastic model for shear deformations \citep[e.g.,][]{carcione2007}.
$\sigma$, $\epsilon$,  $\mu$ and $\eta$ represent stress, strain, shear modulus and viscosity, respectively, 
where $\eta_1$ describes seismic relaxation while $\eta$ is related to plastic flow and processes such as  
dislocation creep.
}
\end{figure}

\begin{figure}
\vspace{1cm}
\includegraphics[width=12cm,height=9cm]{fm3.eps}
\vspace{1cm}
\caption{
Comparison between the analytical (solid line) and numerical (symbols) PS-wave solutions. The fields are normalised.
The amplitude in (c) and (d) are much lower than in (a) and (b) due to attenuation arising from the plastic viscosity. 
The S-wave has disappeared for $\eta = 10^9$ Pa s.}
\end{figure}

\begin{figure}
\vspace{1cm}
\includegraphics[]{fm4.eps}
\vspace{1cm}
\caption{
Comparison between the analytical (solid line) and numerical (symbols) SH-wave solutions. The fields are normalised with 
respect to the amplitude of the higher viscosity.}
\end{figure}

\begin{figure}
\begin{center}
\vspace{1cm}
\includegraphics[width=10cm, height=10cm]{fm5.eps}
\vspace{1cm}
\caption{Temperature profile with a sharp increase at about
3.5 km depth.
\label{efig1}}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\vspace{1cm}
\includegraphics[width=11cm,height=8cm]{fm6.eps}
\vspace{1cm}
\caption{Phase-velocity profiles as a function of temperature at 10 Hz.
\label{efig2}}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\vspace{1cm}
\includegraphics[width=11cm,height=8cm]{fm7.eps}
\vspace{1cm}
\caption{Phase-velocity profiles as a function of depth at 10 Hz
for the GBF, based on the temperature profile of Figure 4.
\label{efig3}}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\vspace{1cm}
\includegraphics[width=11cm,height=8cm]{fm8.eps}
\vspace{1cm}
\caption{Vertical gradient of the P and S-wave velocities in the transition zone at 10 Hz.
\label{efig4}}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\vspace{1cm}
\includegraphics[width=11cm,height=8cm]{fm9.eps}
\vspace{1cm}
\caption{Phase-velocity profiles for P waves at 3, 10 and 30 Hz.
\label{efig5}}
\end{center}
\end{figure}
\newpage

\begin{figure}
\begin{center}
\vspace{1cm}
\hspace{-1cm}
\includegraphics[width=14cm,height=9cm]{fm10.eps}
\vspace{1cm}
\caption{Synthetic shots calculated from a 2D model with a vertical source at the surface, where 
a) shows the vertical component and 
b) the  horizontal component.
In both panels we can see the reflections of the P-waves (RP) and
S-waves (RS) from the transition zone. 
\label{efig6}}
\end{center}
\end{figure}

\newpage

\begin{figure}
\begin{center}
\vspace{1cm}
\includegraphics[width=10cm,height=9cm]{fm11.eps}
\vspace{1cm}
\caption{Synthetic VSP in the melted zone due to a vertical source located at the surface (vertical component). RP 
indicates the reflected wave and TP the transmitted P wave.
\label{efig7}}
\end{center}
\end{figure}

\newpage

\begin{figure}
\begin{center}
\vspace{1cm}
\includegraphics[width=10cm,height=9cm]{fm12.eps}
\vspace{1cm}
\caption{An $f$-$k$ plot of the VSP shows that the dispersion effects are moderate. 
\label{efig8}}
\end{center}
\end{figure}




%% TABLES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% ONE-COLUMN TABLE

%t
%\begin{table}[t]
%\caption{hfjhfhjf kjgt}
%\vskip4mm
%\centering
%\begin{tabular}{column = lcr}
%\tophline

%\middlehline

%\bottomhline
%\end{tabular}
%\end{table}


%% TWO-COLUMN TABLE

%t
%\begin{table*}[t]
%\caption{TEXT}
%\vskip4mm
%\centering
%\begin{tabular}{column = lcr}
%\tophline

%\middlehline

%\bottomhline
%\end{tabular}
%\end{table*}


%% The different columns must be seperated with a & command and should
%% end with \\ to identify the column brake.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% If figures and tables must be numbered 1a, 1b, etc. the following command
%% should be inserted before the begin{} command.

\addtocounter{figure}{-1}\renewcommand{\thefigure}{\arabic{figure}a}

\end{document}