A nonlinear and time-dependent visco-elasto-plastic rheology model for studying shock-physics phenomena

We present a simple and efficient implementation of a viscous creep rheology based on diffusion creep, dislocation creep and the Peierls mechanism in conjunction with an elasto-plastic rheology model into a shock-physics code, the iSALE open-source impact code. Our approach is based on the calculation of an effective viscosity which is then used as a reference viscosity for any underlying viscoelastic (or even visco-elasto-plastic) model. Here we use a Maxwell-model which best describes stress relaxation and is therefore likely most important for the formation of large meteorite impact basins. While common viscoelastic behavior during mantle convection or other slow geodynamic or geological processes is mostly controlled by diffusion and dislocation creep, we showed that the Peierls mechanism dominates at the large strain rates that typically occur during meteorite impacts. Thus, the resulting visco-elasto-plastic rheology allows implementation of a more realistic mantle behavior in computer simulations, especially for those dealing with large meteorite impacts. The approach shown here opens the way to more faithful simulations of large impact basin formation, especially in elucidating the physics behind the formation of the external fault rings characteristic of large lunar basins.


Introduction
Viscoelastic behavior of material has been studied intensively over the last decades. Studies range from engineering and industrial applications (e.g. Darabi et al., 2011;Larsen et al., 2009;Morian et al., 2014;Hong, 2011) as well as biological (e.g. Grooman et al., 2012;Rouse Jr., 2004), medical (Gennissson et al., 2010;Streitberger et al., 2011;Teran et al., 2010) and geological objectives (Lange et al., 2014;Remus et al., 2011;Wang et al., 2012;Nield et al., 2014;Theofanous, 2011). The concept of viscoelasticity was first proposed by Maxwell in 1867 (Maxwell, 1867). Later on, progress in the understanding of crystal structures of geological materials resulted in improved mathematical and physical models [Karato, 2010], such as the diffusion creep model [Nabarro, 1948] or a power--law dislocation creep model [Weertman, 1955]. These and other improvements to the description of viscoelastic behavior contributed to a better understanding of (geo)dynamical problems, such as plate tectonics, mantle convection, mountain ridge formation or dynamic earthquake ruptures. Over the last decades, extensive efforts have been made to

The iSALE hydrocode
The hydrocode iSALE (Elbeshausen et al., 2009;Elbeshausen and Wünnemann, 2011;Wünnemann et al., 2006;Collins et al., 2004; http://www.isale--code.de) is a multi--material shock physics code (historically called "hydrocode"). It solves the Navier--Stokes equations in a compressible fashion on a Cartesian staggered mesh by using finite differences and finite volumes in two and three dimensions. The solver follows the scheme of Hirt et al., 1974 which allows the calculation of flows at nearly arbitrary speeds. The simulations can be performed in an Eulerian or Lagrangian fashion as well as by using an Arbitrary Lagrangian Eulerian (ALE) technique. For calculations of meteorite impacts, however, the large strains would greatly distort a Lagrangian mesh, making the Eulerian approach the most reasonable one. The code additionally consists of a number of different equations of state, such as ANEOS [Thompson and Lauson, 1972] or Tillotson [Tillotson, 1962], allowing computations with a realistic thermodynamic behavior for many geomaterials. In addition different constitutive models are required to account for material damage and failure [Collins et al., 2004;Ivanov et al., 1997], porosity [Wünnemann et al., 2006], dilatancy [Collins, 2015], acoustic fluidization [Wünnemann and Ivanov, 2003], thermal softening [Ohnaka, 1995], or low--density weakening. The underlying constitutive model, however, is either elasto--plastic or purely viscous (simple linear and time--independent viscosity). A visco--elastic or even visco--elasto--plastic rheology has not been considered yet. The development and implementation of this extension is presented in the subsequent sections.

The Peierls mechanism -description
The rheology of rocks depends on a large number of constitutive and environmental factors, including mineralogy, fluid content and chemistry, mineral grain size (see e.g. Riedel and Karato, 1997), melt fraction, temperature, pressure, and differential stress conditions [Bürgmann and Dresen, 2008]. The same is true for the viscous behavior of e.g. the Earth's interior: While the upper crust is usually assumed to be in a frictional equilibrium with active faults that might limit the strength [Bürgmann and Dresen, 2008], a pressure dependent increase of frictional strength with depth is counteracted by thermally activated creep processes which reduce the viscous strength with increasing temperature and depth [Goetze & Evans 1979;Goetze, 1978;Hirt and Kohlstedt, 2003;McBirney and Murase, 1984;Melosh, 1980]. A simple linear viscous rheology is not sufficient for a realistic calculation of mantle behavior. We therefore use a more sophisticated model that has been previously presented in several papers (e.g. Kameyama et al., 1999;Kawazoe et al., 2009;Karato and Wu, 1993). This model consists of three different regimes depending on stress and temperature: Grain--size dependent diffusion creep, a power-law dislocation creep (see also e.g. Faul et al., 2011), and the exponential Peierls mechanism for larger stresses. While deformation caused by classic geodynamical mantle convection is dominated by diffusion and dislocation creep, the Peierls mechanism becomes essential for deformation at larger stresses and strain rates, such as during meteorite impact and the subsequent collapse of the initial deep crater.
In our approach the overall strain rate due to inelastic deformation ε v consists of a sum of viscous strain rates due to the three creep mechanisms: Diffusion creep ( ! ε 1 ), Dislocation (or power--law) creep ( ! ε n ), and Peierls creep ( ! ε P ).

Diffusion creep accounts for deformation of material by diffusion of vacancies in the crystal lattice. It is
therefore a plastic deformation, which gives rise to ductile, not brittle, behavior. It is generally a linear function of stress. The strain rate due to diffusion creep depends not only on temperature T and deviatoric stress σ (this is usually taken to equal the second invariant of the stress tensor, as the appropriate generalization of shear stress) but also on grain size a: where a 0 is the reference grain size (here: a 0 =1 mm), R is the gas constant, η 0 is a reference viscosity, and E 1 is an activation enthalpy. m is a material--dependent parameter describing the grain size--dependency of diffusion creep, which ranges between 2 for volume diffusion (Nabarro--Herring creep) and 3 for grain boundary diffusion (Coble creep). Diffusion creep is mainly important at lower stresses, is highly temperature dependent and results in very small strain rates.

Dislocation creep (power--law creep)
describes material deformation as the result of dislocations moving through the crystal lattice (see also Carrez and Cordier, 2010). It is best described by a power--law function of the deviatoric stress, in which the strain rate depends on temperature and stress but is independent of the grain size: Here E n is the activation enthalpy of non--Newtonian creep, σ c is the critical stress, and n is the exponent for the power--law dependence of this creep that usually ranges between 3 and 5. This mechanism usually predominates at intermediate stresses and higher temperatures.
The Peierls mechanism is a special type of dislocation creep that is thermally activated and dominates only at high deviatoric stress. The creep law for Peierls mechanism is given by where A P and q are material--dependent parameters, E P is the activation enthalpy for the Peierls mechanism, and σ P is the Peierls stress. The Peierls creep dominates for stresses above about 500 MPa, although this depends upon the material: Mantle materials such as olivine are much more susceptible to Peierls creep than typical crustal minerals such as plagioclase.
The parameters of the equations shown above are strongly material dependent. Since our study focuses on the mantle rheology we chose values that best reflect dry olivine, as published in Kameyama et al., 1999. See Table 1 below for the parameters. Expanding equations (1)--(4), the overall viscous strain rate is given by The deformation map in Figure 1 shows the relationship between strain rate, stress, and temperature as given by equation (5). We can clearly see that the diffusion creep only accounts for strain rates which are (usually) much smaller than those expected for large--scale planetary impact events, except for much slower viscous relaxation long after the crater forms. Thus, it is likely that diffusion creep does not play an important role for this study, where typical strain rates range from 1×10 --3 to 10 s --1 . For performance reasons equation (2) might be neglected if applied for planetary scale impact simulations. In the framework of this study, however, we included this mechanism for completeness. Instead we ignored a possible pressure dependence in the Peierls regime (as e.g. suggested in Kawazoe et al., 2009), because it is much smaller than the temperature and strain rate dependence for the simulations intended within the framework of our study. Indeed, if long--term crater relaxation is the goal of a simulation, then iSALE is a poor choice, because the inclusion of inertial accelerations in iSALE requires very small time steps that would be impractical for deformation occurring over millions of years.

Coupling the Peierls mechanism with an elasto--plastic model
The time dependence of the response of a viscoelastic system is similar to that of an electrical circuit. From a mathematical perspective both systems can be described by an identical set of ordinary differential equations. A convenient approach to this problem is the so--called "spring--dashpot" models, where the elastic rheology is represented by an elastic spring, described by σ =k ⋅ε Here σ and ε represent the spring force and displacement, and the spring constant k plays the role of the shear modulus µ. The viscous rheology can be visualized by a Newtonian dashpot where η is the viscosity. These two basic elements can be combined in different ways to simulate more complex rheologic behavior.
Following this approach, different rheologic models are in general use: The Maxwell--configuration connects one spring and one dashpot in series, so that the strain is the sum of each element while the stress is the same in both. The Kelvin--Voigt puts a viscous and an elastic element in parallel so that the strain is equal in both elements while the total stress is the sum of the individual stresses. A Standard Linear Solid combines a Maxwell element and a spring in parallel. To describe a more complex rheology such as that of polymers, mathematical solutions have been presented that rely upon arbitrary combinations of multiple elements [Wiechert, 1893;Tschoegl, 1989;Brinson and Brinson, 2008].
We used a Maxwell--configuration in conjunction with a plastic (Bingham) flow model to connect the viscous creep rheology as described above to an elasto--plastic model to allow visco--elastic or even visco-elasto--plastic material behavior. When using elasticity, viscosity, and plasticity in series, every element in the entire system is experiences the same stress while the strain rates are additive: The subscripts T, v, e, and p refer to the total, viscous, elastic, and plastic contributions, respectively. The constitutive equation for viscoelasticity in a Maxwell--configuration is directly inherited from equations (6), (7), and (8): The advantage of utilizing a Maxwell--model is that is permits relaxation of stresses, while allowing arbitrarily large strains, consistent with observations of the behavior of geologic materials. For example, under the condition of a strain that is suddenly applied and then held fixed, the stress relaxes by where M T η µ = is the Maxwell decay time.

Implementation of visco--elasto--plasticity in a hydrocode
The constitutive equation for viscous creep, equation (5), shows how to calculate the strain rate from a given stress. However, the hydrocode we utilize for our study, iSALE, works in the opposite way: it updates the stresses from a given strain rate. This creates two complications that must be solved: 1. The constitutive equation (5) must be inverted to compute the stress from a given strain rate. Because this equation is a transcendental equation it must be solved numerically, which is very time--consuming. 2. Because we calculate the stress contributions from the total strain rates and we assume that each element acts in series, we need to determine the correct partitioning of the total strain rate ( ! ε T ) into an elastic strain rate ( ! ε e ), a viscous strain rate ( ! ε v ) and a plastic strain rate ( ! ε p ).
In the series configuration the same stress acts on each element, which responds with its own strain rate, according to the size of each term in equation (5). Please note: when we use deviatoric stress and strain, these equations hold separately for each component of the stress and strain.
In the following, we present an approach to deal with these issues. This approach to the computation of visco--elasto--plastic behavior requires four different steps which are described in this section:

Derive an effective viscosity
For viscous material behavior, we consider diffusion creep, dislocation creep, and Peierls creep. The strain rate due to each mechanism is additive, which results in the main equation (5). Note that most experimenters use pure shear, so they measure 2 ij ε ⋅ and not ij ε . Thus, our basic equation is Considering the constitutive equation for the viscous element we can easily compute eff η at any working stress σ: Here we use the second invariant of the old stress tensor as a reference stress.

Calculate the effective stress
At the beginning of each time step we know the previous stress 0 ( ) ij t σ . In our approach viscosity only contributes to the stress of "intact" material, i.e. material whose stress does not exceed the yield envelope ( II σ <Y). In this case, no plastic strain occurs ( ! ε p = 0 ). For the current (i.e. new) time step, we therefore know Using the constitutive equations (6) and (7) we obtain Please note that here we use the effective viscosity eff η η = . The differential equation for σ is given by where M T η µ = is the Maxwell decay time. This results in the general solution given by to first order in t Δ . Now we can use equation (17) to calculate the effective stress from the effective viscosity.

Partitioning viscous and elastic strain rates
By using the derived effective stress, we can now compute the elastic and viscous strain easily. The elastic strain rate is given by and the viscous strain rate is obtained by

Plasticity -rheology above the yield envelope
Finally, we compare the new stress (here: the second invariant II σ ) with the yield envelope. If II σ < Y, no plastic strain occurs and the total strain rate is partitioned into elastic and viscous strain rates as given in the equations (18) and (20) above. If the stress exceeds the yield limit ( II Y σ ≥ ), then stresses relax to the yield envelope and, once the new stress is known, we can use again the constitutive equations to compute the new (reduced) elastic and viscous strain rates ! ε e and ! ε v at a stress on the yield envelope: The remaining part contributes to the plastic strain rate

Results and validation
To test and explore our rheology model we performed numerical simulations of a 1 m sized block sheared with a velocity of 1 m/s and recorded the resulting stresses and effective viscosities. While this test becomes quite unrealistic after strains exceed a few tens of percent, it provides a stringent test of the stability of our numerical method. The block is composed of dunite. To calculate the thermodynamic behavior we used the Analytical Equation of State (ANEOS; see Thompson and Lauson, 1972) with tabularized data for dunite. For the viscous calculations we used the same parameters as those listed in Table  1. We assumed a pre--heated block at an initial temperature of T=1700 K to focus on the regime where the Peierls creep becomes important (see also Figure 1). This temperature also better reflects the thermal conditions in the (Earth) mantle.
The elastic stress contribution e σ is calculated according to equation (17) by where the elastic shear modulus µ is computed from the speed of shear waves and density and, thus, on the thermodynamic behavior of the material. At the conditions described above, it is on the order of 130 GPa.
The plastic behavior is simulated by limiting stresses to the yield envelope and reducing strain rates properly (see section 4.4). The iSALE hydrocode contains different methods for calculating the yield envelope of a given material. For these tests we utilized the so--called rock--model (see Collins et al., 2004). In this model, the yield envelope comprises different stress paths for intact and damaged material. While the yield envelope for completely damaged material is calculated by using a simple Mohr--Coulomb approach, intact material is represented by a Lundborg relationship. Between these two states a linear interpolation is used to retrieve an adequate yield limit. The parameters used for this plasticity model are listed in Table 2. 1000. J kg --1 K --1 Shear strength at initial condition (intact material) 10 MPa Coefficient of internal friction (intact material) 1.1 Shear limit (intact material) 2.5 GPa Shear strength at initial condition (damaged material) 10 KPa Coefficient of internal friction (damaged material) 0.8 Shear limit (damaged material) 2. GPa We calculated the shearing of the block for different material rheologies: (a) purely elastic, (b) elasto-plastic, (c) visco--elastic and (d) visco--elasto--plastic. The results are shown in Figure 2ff. Stresses and strain rates are always represented by their XY--component. Figure  2 shows the time evolution of stress for all four different rheologies. While a purely elastic rheology results in a significant and continuous stress increase, stresses in the elasto--plastic regime are significantly lower. In the latter case the resulting stresses are always above the yield envelope. Thus, when considering elasto--plastic behavior, plastic work and plastic strain occurs at any time during deformation.
A comparison between the elastic and visco--elastic behavior reveals an identical stress evolution in the very early stage. Later on, however, stresses begin to relax in the viscoelastic case. The rate of the stress decay is defined by the Maxwell time M T η µ = which describes the time required to decrease the stress to 1/e of the initial stress. In laboratory experiments this measure is often used to obtain insights into the viscous behavior of a given material. In our approach, however, the Maxwell time is itself time-dependent, because the effective viscosity η varies with time, as shown in Figure  3. In this experiment, the Maxwell time quickly declines from 35 s in the beginning to less than 30 ms after 0.1 s from the beginning of the shearing process. In the same time the effective viscosity decreases over three orders of magnitude from > 10 12 Pa--s to 1 GPa--s.

Figure 2 Evolution of stress during shearing for elastic, elasto--plastic, visco--elastic and visco--elasto--plastic material behavior.
The rollover of the elastic stresses after about 0.1 s is due to the extreme distortion of the initial 1 x 1 x 1 m cube of dunite (the strain is 10% at 0.1 s). At 1 s the total strain is 100% and continues to increase thereafter as the block is distorted from a cube into a long parallelepiped. Our algorithm nevertheless faithfully follows the stress evolution in this distorted block. Note that in our rheological model the stress is the same in each element of the system, so the stress plotted acts equally on the elastic, viscous and plastic elements.
Please note that the magnitude of these contributions is strongly dependent on the size of the time increment ( t Δ ), as equations (17) and (28) illustrate. Also note the negative sign in the viscous term of this equation indicating that the viscous contribution counteracts the elastic stresses. As is also evident in Figure 3, the rheology is initially dominated by elastic (or elasto--plastic) behavior. The viscous stresses are significantly lower than the elastic stresses, but increase rapidly. Thus, very shortly after the onset of shearing, viscous stresses become more prominent than the elastic contributions. The stress therefore decays. This happens roughly when the Maxwell time and the current time step are the same (at the intersection between the Maxwell decay time and the decay limit in Figure 3). The dominance of the viscous contributions results in a quick decay of the total stresses. While the elastic stresses increase slowly, the viscous stresses decrease after reaching their maximum. Thus, the decrease of stress is also damped by time, in agreement with equation (10).
When considering a full visco--elasto--plastic rheology, the material initially behaves as purely elasto-plastic (see Figure 2) because there is not time for viscous strains to accumulate. The resulting stresses in the visco--elastic and visco--elasto--plastic scenarios are exactly the same until the block is strongly deformed and the Bingham yield stress is exceeded. Despite stress relaxation in the viscous element (and the corresponding accumulation of viscous strain), the resulting stresses are above the yield envelope during large parts of the deformation. Stresses are, thus, reduced (see section 4.4) and part of both the viscous and elastic strain is converted into plastic strain. At some later time, however, the viscous flow relaxes the stress below the yield envelope. From this moment on, the resulting stress veers away from the yield envelope and the material behaves as purely visco--elastic.
All three regimes considered by us to calculate a realistic viscous behavior - Diffusion creep, dislocation creep and the Peierls mechanism -depend on temperature. Each of these regimes is activated or dominant at different temperature conditions, as Figure 1 illustrates. Thus, we repeated the experiments described above and systematically varied the initial temperature of the block material. Figure  5 shows the stress evolution for both purely elastic (dashed lines) and visco--elastic (solid lines) material behavior and for different temperatures. It is evident that the resulting stresses decrease with increasing temperatures. If only the viscous behavior is considered, however, the stresses are much lower and even the temperature dependence on stress is much more pronounced. Furthermore, also the separation of the visco--elastic path from the purely elastic path occurs much earlier for higher temperatures, as becoming visible in Figure 5 (right), which shows the stress evolution for the very early stage. This observation is also demonstrated in Figure  6, showing the effective viscosity as a function of time. The initial viscosity is much higher for lower temperatures. At room temperature, the resulting initial viscosity 0 ( ) t η ≈ is even above 10 60 Pa--s and the material behaves therefore more solid--like until viscosity starts to decrease after some time due to the stress dependence of the viscosity. The onset of a significant and nonlinear decrease of the viscosity depends strongly on temperature and occurs much earlier for higher temperatures. At a temperature of 2000 K the initial viscosity drops down to 10 GPa--s and viscous effects are therefore more prominent right from the beginning of the deformation. Despite the different temperature conditions all viscosities asymptotically approach the same value at a later stage (here ≈10 GPa s). The resulting stresses for T=2000 K are a bit lower which can be explained by the onset of changes in the thermodynamic behavior (beginning phase transitions) at conditions close to the melting point of dunite.

Summary and conclusion
We have presented a simple and efficient way of implementing a viscous creep rheology based on Diffusion creep, Dislocation creep and the Peierls mechanism in conjunction with an elasto--plastic rheology model into a shock--physics code. The three regimes are combined continuously, i.e. there is no threshold initiating a sudden change from one regime to the next. While the viscoelasticity of the mantle as utilized for mantle convection or other slow geodynamic or geological processes is mostly controlled by diffusion and dislocation creep, we showed that the Peierls mechanism becomes the dominant mechanism at the large stresses and strain rates occurring during meteorite impacts.
Our approach to implementing the model into the shock--physics--code iSALE is based on the calculation of an "effective viscosity" which is then used as a reference viscosity for any underlying viscoelastic (or even visco--elasto--plastic) model. Here we use a Maxwell--configuration, where elastic, viscous and plastic strain rates sum up to the total strain and the stress contributions at the viscous, elastic and plastic elements are identical). The Maxwell--model best describes stress relaxation (see equation (10)) in natural materials, and so is likely to be important for the formation of large meteorite impact basins, for which this code was intended. However, this model does not consider the full range of observed creep phenomena accurately. Under constant stress conditions, strain in this model increases linearly with time. Most natural materials, however, exhibit a more complex relationship and at lower temperatures strain rates usually decrease with time, in a regime known as primary creep. The Kelvin--Voigt approach, where the viscous and elastic elements are arranged in parallel (and, thus, each element encounters the same strain rate, but the resulting stresses sum up to the total stress), simulates creep at small strains much more accurately, but it is less accurate for predicting large strain behavior. In the last decades more complex rheologies have been developed and proposed, such as Standard--Linear--Solid (where a Maxwell--model is combined with an additional elastic element in parallel) or a Generalized Maxwell Model or Maxwell--Wiechert model (arbitrary number of Maxwell--elements in parallel, see Wiechert, 1893;Tschoegl, 1989;Brinson and Brinson, 2008). Our approach of calculating an effective viscosity can be easily applied to these more complex configurations, which we intend to do in future work.