Interactive comment on “ Simulation of seismic waves at the Earth crust ( brittle-ductile transition ) based on the Burgers model ”

General Comments: The paper entitled Simulation of seismic waves at the Earth crust based on the Burgers model is concerned with wave propagation in the crust transition based on a Burger model. The paper argues this Burger model is appropriate (physically and geologically) for such a brittle-ductile transition. The paper also gives the formulation of numerical simulation for 2D P-SV and SH and then verifies the accuracy against the known analytical solutions. The papers have theoretical and numerical examples to demonstrate wave propagation in such a Burger models. Overall, the paper is easy to follow and well written. There is however a number of points which need clarification before publish.


Introduction
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 (Meissner and Strehlau, 1982;Zappone, 1994;Simpson, 1999) and the availability of geothermal energy (Manzella et al., 1998).The BDT in the earth is generally viewed as a transition between two different constitutive behaviors, viscoelastic and plastic (Dragoni, 1990).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 (Brogi et al., 2003).
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 (Matsumoto and Hasegawa, 1996).The ductile medium mainly flows when subjected to deviatoric stress, while it does not show major flow under hydrostatic stress.The flow deformation is then mainly associated with the shear modulus of the medium.The deviatoric stress is determined by the octahedral stress, σ o , a scalar that is invariant under coordinate transformations and whose value determines the character of the flow.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  Carcione, 2007).σ, ǫ, µ and η represent stress, strain, shear modulus and viscosity, respectively, where η 1 describes seismic relaxation while η is related to plastic flow and processes such as dislocation creep.
when σ o exceeds the elastic octahedral-stress limit.Below this limit, there is gradual creep deformation when constant stress is applied.Then, if σ o is lower than the elastic limit, the material follows a viscoelastic stress-strain relation.If σ 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 Hueckel et al. (1994), based on a thermoplasticity theory, where the elastic domain is postulated as temperature dependent, shrinking with temperature.On the other hand, Arcay (2012) proposes a thermomechanical model based on a non-Newtonian viscous rheology and a pseudo-brittle rheology.
It is widely accepted that linear viscoelastic-plastic models are appropriate to describe the behavior of ductile and brittle media.Gangi (1981Gangi ( , 1983) ) used this type of model to fit data for synthetic and natural rock salt.The viscoelastic creep of salt has been described with a Burgers model by Carcione et al. (2006).Carcione and Poletto (2013) used the Burgers model to describe the BDT transition, including the presence of anisotropy and seismic attenuation.The Burgers model is shown in Fig. 1 (the Maxwell and Zener model are particular cases of this model).The simulations presented here consider both media, ductile and brittle, with the Zener model used to describe 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 modeling, inversion and interpretation (Zappone, 1994;Long and Silver, 2009;Juhlin and Lund, 2011).In particular, grid methods are required for simulating wave propagation in heterogeneous realistic models (e.g., Carcione et al., 2002;Seriani et al., 1992).In this work, we propose to simulate seismic wave propagation in heterogeneous media involving the brittleductile transition.The differential equations are solved in the time domain by using memory variables (Carcione, 2007).We assume isotropic media and plane strain conditions and obtain the differential equations of motion for 2-D 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 (e.g., Carcione, 2007).

The Burgers mechanical model
The constitutive equation, including both the viscoelastic and ductile behavior, can be written as a generalization of the 1-D stress-strain relation reported by Dragoni (1990) and Dragoni and Pondrelli (1991) to the 3-D anelastic case, replacing the Maxwell model by the Burgers model (Carcione et al., 2006;Carcione, 2007;Carcione and Poletto, 2013).
The Burgers model is a series connection of a dashpot and a Zener model as can be seen in Fig. 1.The usual expression in the time domain is the creep function (Carcione et al., 2006;Chauveau and Kaminski, 2008), where t is time and H (t) is the Heaviside function.The quantities τ σ and τ are seismic relaxation times, µ 0 is the relaxed shear modulus (see below) and η is the flow viscosity describing the ductile behavior related to shear deformations.The frequency-domain shear modulus µ can be obtained as µ = [F( χ )] −1 , where F denotes time Fourier transform and a dot above a variable denotes time derivative.It is formulated as where i = √ −1 and ω is the angular frequency.The relaxation times can be expressed as where τ 0 is a relaxation time such that ω 0 = 1/τ 0 is the center frequency of the relaxation peak and Q 0 is the minimum quality factor.The dependence of the quality factor as a function of frequency is given below in Eq. ( 24).The limit η → ∞ in Eq. ( 2) recovers the Zener kernel to describe the behavior of the brittle material, while τ σ → 0 and τ → 0 yield the Maxwell model used by Dragoni (1990) and Dragoni and Pondrelli (1991): (e.g., Carcione, 2007).For η → 0, µ → 0 and the medium becomes a fluid.Moreover, if ω → ∞, µ → µ 0 τ /τ σ and µ 0 is the relaxed (ω = 0) shear modulus of the Zener element (η = ∞).
The viscosity η can be expressed by the Arrhenius equation (e.g., Carcione et al., 2006;Montesi, 2007).It is related to the steady-state creep rate ˙ by where σ o is the octahedral stress.The creep rate can be expressed as (e.g., Gangi, 1983;Carcione et al., 2006;Carcione and Poletto, 2013), where A and n are constants, E is the activation energy, R = 8.3144 J mol −1 K −1 is the gas constant and T is the absolute temperature.The form of the empirical relation (Eq.6) is determined by performing experiments at different strain rates, temperatures and/or stresses (e.g., Gangi, 1983;Carter and Hansen, 1983).
In order to obtain the equations of motion to describe wave propagation it is convenient to consider the Burgers relaxation function (Carcione, 2007), where and In terms of the relaxation times and µ 0 , it is The complex shear modulus is It can be verified that Eqs. ( 2) and ( 11) coincide.

2-D 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 (e.g., Carcione, 2007), where λ is a Lamé constant, σ are stress components, v are particle-velocity components, ∂ 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 " * " denotes time convolution.
The convolutions have the form ψ * ∂ i v j and can be overcome by introducing memory variables.We obtain The stress-strain relation becomes zx .
On the other hand, the dynamical equations of motion are (Carcione, 2007), where ρ is the mass density and s i are source components.
The equations of motion are given by Eqs. ( 15), ( 16) and (19) in the unknown vector where M is a 13 × 13 matrix containing the material properties and spatial derivatives.
In view of the correspondence principle (e.g., Carcione, 2007), the complex and frequency-dependent P and S wave velocities are and the P and S wave quality factors are given by (e.g., Carcione, 2007).

Propagation of SH waves
The stress-strain relations describing shear motion in the (x, z) plane are (Carcione, 2007).
On the other hand, the dynamical equation of motion is where s is the source (Carcione, 2007).
Applying the same procedure as in the P -S case we obtain with The equations of motion are given by Eqs.(15)(16)(17)(18)(19) in the unknown vector σ = (v y , σ xy , σ zy , e (m) xy , e (m) zy ) and can be recast as Eq. ( 20) with matrix M of dimension 7 × 7.

Numerical solution
The formal solution to Eq. ( 20) with zero initial conditions is given by where exp(tM) 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 (Jain, 1984;Carcione, 2007).
The spatial derivatives are calculated with the Fourier pseudospectral method (Kosloff and Baysal, 1982;Carcione, 2007).This method consists 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 max dt/d min < 2.79, where c max is the maximum phase velocity, dt is the time step and d min is the minimum grid spacing.

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 where t p is the period of the wave (the distance between the side peaks is √ 6t p /π) and we take t s = 1.4t p .Its frequency spectrum is 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 −1 .Considering a Poisson medium, we obtain an S-wave velocity of c S = 3.464 km s −1 .We have λ + 2µ 0 = ρc 2 P and µ 0 = ρc 2 S .Assuming a density of ρ = 2700 kg m −3 , we have λ = µ 0 = 32.4GPa.The seismic quality factor is Q 0 = 40 and ω 0 = 2πf p .The numerical mesh has 231 × 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  b) and (d) to v z .At the source peak frequency, the P-and S-wave quality factors for η = 10 20 Pa s are 60 and 40, while those corresponding to η = 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 Fig. 3.In this case, we assume Q 0 = ∞, i.e., attenuation is solely due to the plastic viscosity.The quality factor for η = 2 × 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 (Albertsson et al., 2003).The temperature dependence is expressed by the Arrhenius steadystate power law, e.g., Eq. ( 15) in Carcione and Poletto (2013).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 (Violay et al., 2010(Violay et al., , 2012)).These parameters were determined by mechanical observations in laboratory stress-deformation experiments performed at different confining pressures with glass-free basaltic rock sam-ples (GBF) (Violay et al., 2010(Violay et al., , 2012)), in agreement with the results of Hacker and Christie (1992).In the absence of glass (silica) -which strongly influences the ductile behavior of the basaltic rock at lower temperatures -the glassyfree basalt presents rapid BDT variation at high temperatures.Figure 4 shows the temperature profile used in the calculation, with a steep gradient at depth, similar to that reported in Foulger (1995).
The purpose of the numerical experiment is to predict Pand S-wave propagation with dispersion and attenuation at high temperatures (Carcione and Poletto, 2013), and study the observability of the BDT by seismic reflection methods.This issue poses the problem of evaluating 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 Carcione and Poletto (2013), with an average lithostatic confining pressure of approximately 95 MPa Comparison between the analytical (solid line) and numerical (symbols) SH-wave solutions.The re normalised with respect to the amplitude of the higher viscosity.
21 for a Poisson medium with unrelaxed P wave velocity of 6 km s −1 , S wave velocity of 3.464 km s −1 and density ρ = 2600 kg m −3 .The transition in the GBF velocity functions is quite rapid at approximately 3.55 km in Fig. 6, corresponding to T = 1120 • C. Complete melting is obtained above 1300 • C. The sharp variation of the velocity in the transition zone (Fig. 7) is due to the combined effect of the rapid velocity variation in the BDT zone (Fig. 5) and to the steep temperature gradient (Fig. 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 fre-  quency 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-whiledrilling reverse-VSP experiment (Poletto and Miranda, 2004;Poletto et al., 2011).Figure 9 shows the synthetic commonshots, 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, Fig. 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 (cf.Fig. 7).Higher dispersion can be expected in glassy basalt (GB) (Violay et al., 2010(Violay et al., , 2012)).

Conclusions
The upper -cooler -part of the crust is brittle, while deeper zones present ductile behavior.In some cases, this brittleductile transition is a single seismic reflector with an associated reflection coefficient.The stress-strain relation and its physical implications have been analyzed 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 transi-  tion 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.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 behaviors.
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 examples demonstrate the observability of the brittleductile 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.

Fig. 2 .Figure 2 .
Fig. 2. 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 η = 10 9 Pa s.

Figure 3 .Fig. 4 .Figure 4 .
Figure 3.Comparison between the analytical (solid line) and numerical (symbols) SH-wave solutions.The fields are normalized with respect to the amplitude of the higher viscosity.

Fig. 6 .Figure 6 .
Fig. 6.Phase-velocity profiles as a function of depth at 10 Hz for the GBF, based on the temperatu profile of Figure 4.

Fig. 9 .Figure 9 .
Fig. 9. 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.

Fig. 10 .Figure 10 .
Fig. 10.Synthetic VSP in the melted zone due to a vertical source located at the surface (vertical c ponent).RP indicates the reflected wave and TP the transmitted P wave.

JFigure 11 .
Figure 11.An f −k plot of the VSP.The ellipses indicate the transmitted P wave and direct S wave (solid line) and the reflected wave (dashed line).In the upper part of the model, the melting is negligible.The velocity dispersion effects are moderate; however, they are not clearly evident because of wave superposition in the first ellipse.