Visco elastic theory and assumptions pdf

We use controlled meshing to gradually change the size of the elements in the areas of interest. A coarser element discretization size is assigned to the deep parts of both mantles from about 50 km Figure 5b. Mesh size is chosen with the following considerations: 1 ill-conditioned problems can be avoided when constructing the relation between fault slip and surface deformation; 2 fault creep and variation of material properties will be better resolved close to the oceanic slab; 3 highly distorted elements, which would result in nonconverging solutions and numerical errors, can be avoided close to the trench; and 4 computational time and cost can be saved without affecting resolution accuracy by using coarse elements in deep mantle.

We use this convergence velocity as the reference for simulating back slip rate of fully locked plate interface and calculating obtained locking degree.

The published GPS observations in the study area are mainly based on survey-mode GPS data and therefore do not include the vertical interseismic displacements.

Most of the GPS data are located near the coast in the fore arc, with some data distributed sparsely in the back arc. The use of these data allows us to directly compare our viscoelastic locking prediction with previously published elastic models. The blue curve with solid squares represents the elastic models. The red curve with solid circles represents the viscoelastic models. Therefore, the viscoelastic effects in the interseismic deformation are very likely already being observed in the modern geodetic data, especially for the back-arc region.

The residuals between the GPS velocities and the predicted velocities from the elastic and viscoelastic models 55 and 45 km locking termination depths, respectively are plotted in Figures 8a and 8b.

Elastic and Viscoelastic Inversion of Locking Degree As viscoelastic effects are contained in the geodetic data, it is necessary to make a viscoelastic inversion based on GPS data. In order to reasonably decrease computational cost, we group nearby fault nodes as nonoverlapping patches with a size of about 20 km2 [Masterlark, ; Masterlark and Hughes, ].

An example of three patches in map view is shown in Figure S7. In this way, we achieve an accurate FEM resolution with dense nodes along the fault but a faster calculation of the inversion with larger fault patches.

During the inversion, no back slip constraints are imposed near the updip limit i. The residuals close to the coast point to the trench indicating elastic model overestimation of the deformation there, while the residuals in the back arc point landward showing that the elastic model underestimates the deformation there.

For the viscoelastic model the residuals in the back arc are considerably less than for the elastic model Figure 9c. For the viscoelastic model, the back slip distribution from the inversion does not need deeper back slip on the plate interface and the downdip limit of the locked zone is restricted to a maximum of 50 km depth Figure 9c.

Hence, the back slip pattern of the viscoelastic model is more patchy in both dip and strike direction. The red vectors are published GPS data used in our case study. The blue vectors are predicted deformation from the inverted back slip distribution. The white vectors with black outlines are the residuals from the GPS observations and predictions of inverted slip. In all the four panels, the gray solid contours are the isodepths of the subduction interface the values are given in kilometer.

Discussion 7. The model produces a time-dependent interseismic deformation; eventually, reaching a steady deformation rate after the relaxation time of the viscoelastic materials has elapsed.

In order to determine the optimal viscosity value for the continental mantle, we calculate the average WRMS values of physical locking depths i. Hence, we use this optimal value in the following modeling sections 7.

Purely elastic models cannot produce a long-wavelength deformation signal large enough to be observed in the back arc, restricting the interseismic compression to mainly the fore arc. This has a clear impact when geodetic data are used to invert the depth of the locking zone. Thus, elastic models incorrectly need a deeper downdip limit of the locked zone to reproduce the observed deformation [see Wang et al.

Indeed, the assumed back-arc shortening and sliver motion can be corrected for with a joint modeling of locking and microplate motion [e. While it is very likely that some signal in the back-arc may be due to the long-term geological shortening, as observed in the geological record [e.

Following this relaxation time, the fault locking related deformation becomes dominant. However, the decay to a steady deformation is reached at different times depending the distance to the trench. One implicit assumption for performing a linear viscoelastic inversion in our study is the quasi- time-independent behavior of surface velocity in the late stage of the interseismic period.

Without considering an initial earthquake, the Maxwell material exhibits constant viscoelastic response under constant interseismic loading. Hence, the velocity of surface displacement from an ideal viscoelastic seismic cycle model would remain constant after reaching the relaxation of the interseismic stress.

In a numerical model the effect of initiating back slip on the fault will be diminished after about 20 times the Maxwell time about years in a relaxed simulated system [Hu et al. Therefore, the stabilized velocities on the model surface in the late stages of simulation of for our model are likely due to the constant interseismic loading in a relaxed system Figures S3 and S4.

The viscoelastic response of this earthquake diminishes to near zero 60 to 80 years after the earthquake Figures S5 and S6. After years of simulation time, the differences of the velocities between the earthquake and nonearthquake models are only a few millimeters per year, a value much less than uncertainties in the GPS velocity vectors.

This time that we calculated for previous earthquake effects becoming negligible is consistent with previous numerical studies [e. Thus, by not starting with an initial great earthquake and simulating a total of years of constant interseismic loading, the considered surface velocities in the last stage of simulation can capture well the quasi-time-independent behavior of steady state viscoelastic interseismic deformation in a relaxed system. One simple assumption that could be made is that the late interseismic deformation that we observe is repeated over many cycles.

To achieve this assumption numerically, we can spin-up the model [e. Hence, we do not deem it necessary to make an alternative model which considers the viscoelastic deformation of the previous earthquake.

Another physical inconsistency related to elastic models is the need for very deep locking [e. The modeling approach that we adopt for the case study of the Peru-North Chile subduction margin does not consider the motion of possible microplates sliver motion and back-arc shortening. The residuals of the elastic inversion show a pattern similar to the elastic forward model, i.

Without back-arc shortening and sliver motion corrections, the elastic model needs a larger average back slip magnitude with little variation in lock- ing degree in the along-strike direction Figure 9b when compared to the vis- coelastic inversion results Figure 9d. The residuals of the viscoelastic model Figure 9d are much less than the residuals of the elastic model Figure 9b , and the along margin Figure Comparison of our optimal viscoelastic locking map with the slip segmentation of highly locked patches distributions of the Iquique earthquake derived from our own FEM for the viscoelastic model is increased, inversion , Tocopilla earthquake [Schurr et al.

The gray solid contours are the isodepths of the subduction interface the values are given in kilometer. We use a checkerboard synthetic back slip model to evaluate the model resolution and to demonstrate which features of the locking map can currently be resolved by the GPS data Figure S The elastic and viscoelastic models have similar recovery patterns. Both models give the best resolution under the coastline with resolution rapidly decreasing toward the trench and decreasing more gradually in the downdip direction Figures S10a and S10c.

Our favored locking model is chosen to be the viscoelastic model and is shown in Figure 10 along with the slip distributions of the most recently observed large earthquakes in this margin.

Correlation of Historical Earthquake Slip Distributions With Locking Patterns From Viscoelastic Models The spatial and temporal span of the GPS data is an important factor to consider in interpreting the relationship of the derived apparent locking degree and the slip distributions of historical ruptures. The older set of GPS observations of Kendrick et al. Therefore, the locking state around Iquique earthquake rupture zone was determined by GPS observations from the late stage of the interseismic period, whereas the locking state around Tocopilla and Antofagasta earthquake rupture zones was obtained by GPS observations that are potentially affected by postseismic signals.

The slip distribution shown is obtained from inversion of data corresponding to the main shock and Mw 7. For the Iquique-Pisagua earthquake, there is a very good spatial correlation between a highly locked region and the slip extents of the main shock and largest aftershock ruptures.

Moreover, the patch of locking extends in depth and southward in agreement with the aftershock propagation direction after the Iquique main event [Schurr et al.

The Tocopilla rupture zone is thought to have been highly coupled up to a depth of 50 km before the event [Chlieh et al. This moderate size Mw 7. Hence, the locking degree near the rupture area of the earthquake appears to be low, while the degree of locking updip of this which has not ruptured yet is higher. Therefore, the locking pattern around Mejillones Peninsula at the spatial limits of the Antofagasta and Tocopilla events is not likely to be contaminated by any postseismic relaxation signal and represents an interface in the interseismic state.

This pattern may be related to a long-term creeping barrier [e. By means of synthetic modeling we reveal the pitfalls of inverting viscoelastic interseismic deformation with elastic models: the most notable pitfall being that using a purely elastic model to invert horizontal GPS velocities for locking degree results in an overestimation of the true locking depth.

Therefore, previous locking estimations using a purely elastic model may be overestimating the contribution to the GPS velocities caused by microplate motions. The locking pattern of the viscoelastic model has a better spatial correlation with the slip distributions of the Mw 8.

Knepley, and C. Solid Earth, , —, doi Klotz, and C. Shaoyang Li is a , — Kendrick, R. Smalley, B. Brooks, R. Allmendinger, and B. Isacks , On the strength of interplate coupling and the rate of back arc research school GeoSim. Moreno and Jonathan Bedford have Brooks, B. The authors would Brooks, B. Kogan, G. Steblov, G. Hilley, V. Levin, and E. Apel , Interseismic coupling and asperity distribution along the open source PyLith software for their Kamchatka subduction zone, J.

Chlieh, M. De Chabalier, J. Ruegg, R. Armijo, R. Dmowska, J. Campos, and K. Feigl , Crustal deformation and fault slip during the []. We would like to thank the doi Pardo , Reappraisal of great historical earthquakes in the northern Chile and southern Peru seismic gaps, Nat.

Pardo, L. Dorbath, C. Dorbath, H. Haessler, L. Rivera, A. Cisternas, and L. Ponce , Determination of seismogenic matplotlib module [Hunter, ]. Contreras-Reyes, E. Du, Y. Aydin, and P. Segall , Comparison of various inversion techniques as applied to the determination of a geophysical deformation model for the Borah Peak earthquake, Bull. Freed, A. Micksch, and S. Haberland, C. Rietbrock, D. Lange, K. Bataille, and T. Dahm , Structure of the seismogenic zone of the south-central Chilean margin revealed by local earthquake traveltime tomography, J.

Harris, R. Hayes, G. Wald, and R. Johnson , Slab1. Herman, W. Barnhart, K. Furlong, S. Riquelme, H. Benz, E. Bergman, S. Barrientos, P. Earle, and S.

Samsonov , Continuing megathrust earthquake potential in Chile after the Iquique earthquake, Nature, , — Hetland, E. Hager , Interseismic strain accumulation: Spin-up, cycle invariance, and irregular rupture sequences, Geochem. Kley, E. Klosko, S. Stein, T. Dixon, and E. Norabuena , Consistency of geologic and geodetic displacements during Andean orogenesis, Geophys. Hu, Y. Wang, J. He, J. Klotz, and G. Burgmann, J. Freymueller, P. Banerjee, and K. Wang , Contributions of poroelastic rebound and a weak volcanic arc to the postseismic deformation of the Tohoku earthquake, Earth Planets Space, 66 1 , 1—10, doi Hunter, J.

Hyndman, R. Ivins, E. Kanda, R. Simons , An elastic plate model for interseismic deformation in subduction zones, J. Kendrick, E. Bevis, R. Smalley, and B. Kirby, S. Kronenberg , Rheology of the lithosphere: Selected topics, Rev. Li, S. Moreno, M. Rosenau, D. Melnick, and O.

Masterlark, T. Hughes , Next generation of deformation models for the M9 Sumatra-Andaman earthquake, Geophys. McCaffrey, R. Zones, — Socquet, C. Vigny, D. Carrizo, S. Peyrat, A. Delorme, E. Maureira, M. Valderas-Bermejo, and I. Mitrovica, J. Bolte, J. Klotz, and D. Melnick , Impact of megathrust geometry on inversion of coseismic slip from geodetic data: Application to the Chile earthquake, Geophys. Rosenau, and O. Our PIC finite element method is based closely on the standard finite element method, and is a direct development of the material point method of Sulsky et al.

The standard mesh is used to discretize the domain into elements, and the shape functions interpolate node points in the mesh in the usual fashion.

The problem is formulated in a weak form to give an integral equation, and the shape function expansion produces a discrete matrix equation. For the discretized problem, these integrals occur over sub-domains elements and are calculated by summation over a finite number of sample points within each element.

For example, in order to integrate Eq. Constraints on the values of wp come from the need to integrate polynomials of a minimum degree related to the degree of the shape function interpolation, and the order of the underlying differential equation e. Hughes, These Lagrangian points carry the history variables including the director orientation which are therefore directly available for the element integrals without the need to interpolate from nodal points.

Moresi et al. Numerical simulations We present an example of a simulation of folding of a layer of anisotropic visco-elastic material sandwiched between two isotropic layers of equal viscosity on each side Fig.

To accommodate the shortening of the system, the outermost isotropic layers are compressible. In benchmarking, this sandwich of incompressible and compressible embedding material was found to give good agreement with analytic results assuming an infinite domain Moresi et al. The first two ex- amples are purely viscous. For case 1 Fig. We observe that the growth rate is higher for smaller wavenumber in the infinitesimal deformation limit, and this persists with finite amplitude de- formation.

For case 2 Fig. However, we also observe that low wavenumber modes are excited in the finite deformation limit irrespective of the perturbation wave- number. The perturbation causes a noticeable secondary variation in the interface deflection.

Initial geometry for the folding experiment. Layer 1 is compressible, viscous gM background material, layer 2 is identical to layer 1 but incompressible see text for an explanation , layer 3 is the test sample: visco-elastic l; lS ; g; gS with a director orientation n. The anisotropic layer contains small perturbations to the otherwise horizontal internal layering.

Case 1: evolution of folding in anisotropic viscous layer. Isotropic embedding material has viscosity 1, layer has shear viscosity 1, normal viscosity Case 2: evolution of folding in anisotropic viscous layer. Isotropic embedding material has viscosity 1, layer has shear viscosity 10, normal viscosity Case 3: evolution of folding in anisotropic visco-elastic layer. Case 4: evolution of folding in anisotropic visco-elastic layer. Couple stresses We conclude the main body of this paper with an excursion into the potential influence of couple stresses on the mechanical behaviour of layered materials.

We restrict ourselves to viscous materials, mainly for algebraic convenience; the extension to more general constitutive relationships is reasonably straight- forward however.

Couple stresses are significant in situations where the gradient of ni changes strongly over a short dis- tance limiting case: disclination. In such cases we have to take the variations of the normal stresses across the layer cross sections into consideration e. The couple stress theories see e. Generalised continuum models——so called gradient plasticity models——motivated by dislocation structures in single polycrystalline metals were proposed by Aifantis , In these models the stress tensor is symmetric and non-local hardening effects e.

In the present case the couple stress enhancement leads naturally to the superposition of visco-elastic bending stiffness on our standard continuum model 6.

In connection with layered materials the internal length scale introduced by the couple stresses is proportional to the layer thickness ranging from microns to kilometers in geological applications and to the differences between the viscosities and shear moduli governing pure and simple shear respectively see e. In layered materials the explanation why the stress tensor is non-symmetric in couple stress materials is straight forward: in a continuum description the stresses rep- resent average values over multiples of the layer thickness.

In bending the shear stress obtained by aver- aging normal to the layering is different in general from the shear stress parallel to the layering. The latter may even be zero——for instance, in the case of a stack of perfectly smooth playing cards a standard continuum model would break down in this case. Within the framework of a couple stress theory one considers the variation of the normal stress across the layer thickness in much the same way as in the standard engineering beam and plate theories , introduces statically equivalent couple stresses balancing the difference between the shear stresses Fig.

The couple stress tensor l moment per unit area is conjugate in the rate of energy to the rate of curvature j. In the context of layered materials, a natural choice for the rate of curvature reads: Fig.

Evolution of folding in viscous layer. Isotropic embedding material has viscosity 1, layer has, normal viscosity In this example the internal length was chosen relatively large hence the Cosserat terms act on the larger wavelength, increase the larger length scale. The effect of couple stresses in connection with higher values of the shear viscosity and in particular in the context of visco-elasticity requires further investigation.

There are a number of complications associated with the application of couple stress theories in finite element analyses. For instance rotations cannot independently varied on surfaces where the normal component of the displacements or velocities are prescribed. In this case the velocity gradients on the surface have to be decomposed into normal and surface parallel components.

The surface parallel part will——after application of the surface divergence theorem——produce a contribution to the stress traction see e. Another difficulty arises from the fact that the volume integrals in 16 contain second order velocity gradients so that the shape functions in a finite element model must be continuously differentiable across element boundaries C1 continuity.

Standard finite element programs mostly support C0 continuity only. In 19 P is the so-called penalty parameter. If the energy supply is bounded then we expect that the original constraint 18 will be satisfied in the limit as P!

In the relaxed form of the governing equations we have recovered the equations of the full unconstrained Cosserat continuum where the rotations ui are independent degrees of freedom. In our case this is true as long as P is finite. The independence of the rotations means that the non-independence of surface gradients and the C1 continuity problem have disappeared. In finite element calculations the velocities and rotations are approximated independently and since the highest order derivative of both velocities and rotations are of the first order in the power balance both may be approximated by using the same type of shape function, which needs to be C0 continuous only.

For illustration of the model we consider the simple shearing of an infinite, viscous layer x1 ; 0 6 x2 6 h. First we consider the convergence of the penalty scheme for increasing values of the penalty parameter P. The analytic- and the one-point numerical solution coincide to the first three digits, however the finite element solution based on full integration diverges for increasing values of the penalty parameter: the velocity on top tends to zero.

Finite element model: eight by twelve four noded quad- rilaterals; sixteen particles per element. Periodic boundaries on the sides, i. The latter produces to a positive definite contribution to the argument of the penalty stiffness depending quadratically on x2. In Figs. Conclusions We have presented a simple formulation for the simulation of large, visco-elastic deformations in layered systems. The influence of the bending stiffness of the individual layers is considered within the framework of a couple stress theory.

The combination of the basic model with a large deformation, PIC finite element method allows the simulation of a diverse range of crustal deformation problems. By way of examples we have given a realistic treatment of folding and simple shear processes, which includes the mechanical in- fluence of fine-scale.

The model is relatively simple in its present form but still gives a useful insight into the physical processes involved in certain types of folding processes involving simple shear. One of the most interesting results occurs for purely viscous, layered simulations where low-wavenumber folding is induced even for very low viscosity contrasts between embedded and embedding media. In the past, the very large viscosity contrasts required to produce Biot-type folding in purely viscous media have led people to discount the possibility that viscous buckling occurs at all in geology.

If couple stresses are considered the deformation of the layer is no longer homogeneous; we observe, depending on the viscosity ratio the formation of a slowly deforming core, rapid shearing concentrates around the layer surfaces.

Corresponding results are shown in Figs. Appendix A. Derivation of the anisotropy tensor K In preparation for the derivation of constitutive relations for our layered material we define some auxiliary relationships. From the definition of K see A. Couple stresses We model the layered material as a superposition of thin plates of thickness t.

Basic constitutive relationships are see e. In this case the bending stiffness in C. Appendix D. On the microstructural origin of certain Inelastic models. Aifantis, E. The physics of plastic deformation. Plasticity 3, — Biot, M. Bending of an infinite beam on an elastic foundation. ASME J. A 59, 1—7. Folding instability of a layered visco-elastic medium under compression. A , — Theory of internal buckling of a confined multi layered structure.

The Mechanics of Incremental Deformations. Theory of similar folding of the first and second kind. Caporali, A. Buckling of the lithosphere in western Himalaya: constraints from gravity and topography data. Chapple, W. Fold shape and rheology: the folding of an isolated viscous-plastic layer.

DeGennes, P. The Physics of Liquid Crystals, Second ed. Clarendon Press, Oxford. Fleck, N. Strain gradient plasticity. Fletcher, R. Wavelength selection in the folding of a single layer with power law rheology.

Coupling of diffusional mass transport and deformation in tight rock. Hobbs, B. The influence of chemical migration upon fold evolution in multi-layered materials, vol. In: Krug, H. Hunt, G. Folding processes and solitary waves in structural geology. Johnson, A. Folding of Viscous Layers. Columbia University Press, New York. Leroy, Y. Stability of a frictional cohesive layer on a viscous substratum: Variational formulation and asymptotic solution.