## Abstract

No complete physically consistent model of earthquake rupture exists that can fully describe the rich hierarchy of scale dependencies and nonlinearities associated with earthquakes. We study mesh sensitivity in numerical models of earthquake rupture and demonstrate that this mesh sensitivity may provide hidden clues to the underlying physics generating the rich dynamics associated with earthquake rupture. We focus on unstable slip events that occur in earthquakes when rupture is associated with frictional weakening of the fault. Attempts to simulate these phenomena directly by introducing the relevant constitutive behaviour leads to mesh-dependent results, where the deformation localizes in one element, irrespective of size. Interestingly, earthquake models with oversized mesh elements that are ill-posed in the continuum limit display more complex and realistic physics. Until now, the mesh-dependency problem has been regarded as a red herring—but have we overlooked an important clue arising from the mesh sensitivity? We analyse spatial discretization errors introduced into models with oversized meshes to show how the governing equations may change because of these error terms and give rise to more interesting physics.

## 1. Introduction

One of the fundamental issues in seismology is the dynamics of rupture in the Earth’s crust. Why are some earthquakes small and others huge? What are the origins of earthquake complexity? Although the evidence for complexity is overwhelming, the dynamics of faulting is still poorly understood. Despite this complexity, several empirical scaling relations are universal, for example, the power-law frequency–magnitude distribution (Gutenberg–Richter law) of earthquakes. We demonstrate that it may be possible to isolate the nonlinear physical processes leading to the different length scales in earthquake-related phenomena. We show that mesh sensitivity in numerical models of earthquake rupture can unearth new, exciting physical phenomena and provide hidden clues in discovering a physical law capable of generating earthquake complexity.

Nonlinear behaviour is ubiquitous. Extracting universal behaviour from nonlinear systems and representing the underlying physical phenomena in theoretical and numerical models is a common research goal of all science. In contrast to fluid dynamics where a simple law, the Navier–Stokes equation, is generally used to represent the underlying physics, no such universal law or model underpins earthquake dynamics as yet (Sethna *et al.* 2001). While fluid motion is very complicated on the microscopic scale, the Navier–Stokes equation can describe most fluids on long time and size scales. Physics works because simple laws emerge on large scales. Can we derive such a law for earthquake dynamics?

Several mechanisms have been proposed to explain the wide range of scales of earthquake sizes. One explanation is that geometric and material heterogeneities in the Earth dominate the behaviour and are the source of spatio-temporal complexity at faults. Previous numerical simulations (Rice & Ben-Zion 1996; Zöller *et al.* 2005; Hillers *et al.* 2006) indicate that a range of earthquake event sizes can be produced by treating geometrical heterogeneities of fault zones as geometrical variations of frictional properties in continuum models. The range of size scales characterizing the fault heterogeneities can act as an effective tuning parameter of the fault dynamics. The implications of the structure and distribution of these heterogeneities for predicting the frequency–magnitude distribution of earthquakes remains an outstanding research question.

An alternative argument proposes that, because the power-law scaling of earthquakes implies regular behaviour over a huge range of sizes, this behaviour is likely to be independent of microscopic and macroscopic details (Sethna *et al.* 2001). Complexity is thought to arise primarily as a consequence of dynamical instabilities. Complexities evolve as the system evolves without imposing fixed heterogeneities or stochastic forcing. Several authors (Bak *et al.* 1987; Carlson *et al.* 1994; Langer *et al.* 1996) have concluded that realistic slip complexity including a wide range of event sizes is a generic outcome of inertial dynamics on a homogeneous fault governed by nonlinear friction.

Without question, material and geometrical heterogeneities play a role in producing the complexity displayed by earthquakes. The question is to what extent, and are other more general factors, such as dynamic instabilities, also involved in generating their scale-invariant complex behaviour? We present a third argument that helps unify these two perspectives. We begin by demonstrating that earthquake complexity arises in mesh-dependent discrete earthquake models (Rice 1993; Olsen-Kettle & Mühlhaus in press) with very nonlinear frictional laws. We then show that the source of slip complexity in discrete earthquake models may arise from nonlinear, scale-dependent terms present in these models (Olsen-Kettle & Mühlhaus in press). Scale-dependent terms may be present as geometrical heterogeneities such as asperities, giving further credence to the first argument mentioned above. Isolation and inclusion of the nonlinear, scale-dependent terms present in discrete earthquake models may lead to a universal law where complexity arises naturally (as in the second argument) in the equations and is not an artefact of the numerical mesh.

## 2. Earthquake-fault model

Models of earthquake faults with dynamically unstable frictional laws display strain-softening behaviour and strong mesh dependence, where the physics at short length scales can sensitively determine the behaviour at the largest scales. The mesh-dependence problem can be resolved in two ways. The first is simple—choose a mesh-element size that resolves the underlying physics. However, this becomes computationally prohibitive if the length scale is very small. The second method allows the governing equations to be regularized by implementing higher order continua theories. One can interpret a coarse mesh from the vantage point of a higher order continuum model. This work investigates if higher order perturbative terms present in models with oversized mesh elements produce the spatio-temporal slip complexity, as well as act to regularize the earthquake failure process.

In any numerical approach, there is the question of spatial and temporal resolution of the model space. A well-documented problem that arises in numerical simulations of both mechanical failure (Belytschko & Lasry 1989; de Borst *et al.* 2004) and earthquakes (Rice 1993; Rice & Ben-Zion 1996; Lapusta *et al.* 2000; Shaw & Rice 2000) is that they become extremely sensitive to spatial discretization owing to nonlinearities in the physics. Typically in earthquake studies, this dependency is removed by employing friction laws that do not show this mesh sensitivity and have well-defined continuum limits. An interesting question is whether or not a friction law should be formulated in such a way that it leads to a well-posed problem in elastic-wave analysis. Friction laws that are ill-posed in the continuum limit typically display richer spatio-temporal complexity (Rice 1993; Rice & Ben-Zion 1996; Lapusta *et al.* 2000; Olsen-Kettle & Mühlhaus in press). Until now, the mesh-dependency problem has been regarded as a red herring—but have we overlooked an important clue arising from the mesh sensitivity?

Classical continuum models suffer from excessive mesh dependence when strain-softening models are used in mechanical-failure studies. Continuous descriptions of material damage are prone to non-physical responses beyond a certain level of damage, and numerical simulations then become extremely sensitive to the spatial discretization. In contrast to earthquake studies, where the numerical results using coarser meshes are more realistic, in mechanical failure, this is not usually the case. Thus, mesh dependency in mechanical models has been more rigorously studied to remove the unrealistic effects (Mühlhaus & Vardoulakis 1987; Belytschko & Lasry 1989; Mühlhaus & Aifantis 1991; de Borst *et al.*1993, 2004; de Borst & Mühlhaus 2005; Pasternak & Mühlhaus 2005).

In treating mesh dependency in earthquake models, we draw on concepts from these studies and introduce a new tool for analysing discrete earthquake models. By analysing spatial discretization errors introduced into continuum-fault models, we are able to explain why the coarse meshes (discrete models) sometimes produce more interesting results than those obtained in the continuum limit. We also demonstrate how mesh dependency can be eliminated, and accurate results obtained by subtracting differential terms representing the discretization errors. This has important consequences in unifying discrete and continuum rupture models.

Originally, discrete models referred to discrete fault models using cellular automata or spring-loaded sliding blocks. Several authors (Bak *et al.* 1987; Carlson *et al.* 1994; Langer *et al.* 1996) have shown that these discrete fault models produce realistic slip complexity, including a wide range of event sizes. However, the application of these discrete fault models to describe rupture along a fault embedded in a surrounding elastic continuum has met with some objections. We introduce the concept of inherently discrete models first used by Rice (1993). In this paper discrete models are fault models formulated using a very strongly slip-weakening frictional law where the underlying elastodynamic equations cannot be represented in the continuum limit by the finite-element meshes employed here. In contrast, when we consider a smoother slip-weakening frictional law, it falls into the category of ‘continuum’ models, as they can be represented in the continuum limit by the finite-element mesh and thus do not exhibit mesh dependence.

## 3. Numerical results for earthquake model

Comparison of the mesh sensitivity of a discrete and a continuum fault model is illustrated in figure 1. Figure 1 shows the final slip profile along the fault at the end of an earthquake event for a discrete and a continuum model with varying mesh spacing: 20, 39 and 78 cm. Figure 1*a* clearly shows the mesh dependence of the discrete model where very distinct final slip profiles are generated using different mesh spacing.

The discrete and continuum models correspond to two different slip-weakening frictional laws governed by the coefficient of friction (*μ*_{f}) defined as (Olsen-Kettle & Mühlhaus in press)
3.1
and
3.2where *s* is the inelastic slip during each earthquake event, *D*_{b} is the characteristic slip distance, *μ*_{s} is the static coefficient of friction and *μ*_{d} is the dynamic coefficient of friction. Static re-strengthening occurs in our slip-weakening model when a point on the fault is no longer slipping inelastically (Olsen-Kettle *et al.* 2008). Both pulse-like and crack-like rupture can be produced since the fault can re-strengthen to its initial static value of friction after rupture. The use of nonlinear slip-weakening laws has already been advocated by others on the basis of purely seismic considerations (Abercrombie & Rice 2005; Rice *et al.* 2005; Wibberley & Shimamoto 2005; Chambon *et al.* 2006; Rice & Cocco 2006).

The elastic wave equation was solved, , subject to boundary conditions at the fault, *σ*⋅*n*^{1}=*F*_{n}*n*^{1}+*F*_{τ}*τ*^{1} on *Γ*^{1}, where ** u** is the displacement,

*ρ*is the density of the elastic medium and

*σ*is the stress tensor field given by

*σ*

_{ij}=

*λu*

_{k,k}

*δ*

_{ij}+

*μ*(

*u*

_{i,j}+

*u*

_{j,i}).

*F*

_{n}and

*F*

_{τ}are normal and tangential restoring forces at the fault interface given by the penalty method to prevent slip from occurring. The elastic wave equation was solved using the finite-element method (Olsen-Kettle

*et al.*2008).

The discrete model is considered discrete because it requires a mesh refinement less than 10 cm to resolve its cohesive zone in the continuum limit, whereas the continuum model only requires a mesh refinement less than 1 m to resolve its cohesive zone (Rice 1993; Day *et al.* 2005). A fractal distribution for *μ*_{s} along the fault was employed to investigate the effect of heterogeneities along the fault (Olsen-Kettle *et al.* 2008; Olsen-Kettle & Mühlhaus in press). The dynamic coefficient of friction was proportional to *μ*_{s}, where *μ*_{d}=*α*_{μ}*μ*_{s} and *α*_{μ} ranged from 0.6 to 0.75.

Figure 2*a* demonstrates that richly complex slip at the fault occurs over multiple earthquake events using the discrete model with a coarse mesh spacing of 78 cm. Whereas in figure 2*b*, only simple limit cycles of periodically repeating large earthquakes were produced using the continuum model and the same mesh spacing (Rice 1993). Figure 2*a* shows a wide range of earthquake sizes occurring for the discrete model in contrast to the continuum model. The slip in the discrete fault model is also more strongly influenced by the imposed spatial heterogeneities along the fault.

Analysis of the frequency–magnitude distribution of earthquake events shows that the earthquake sizes fit into three broad categories for the discrete model. The earthquake sizes relate to mesh size (when only one fault element ruptures), fault length (when the entire fault length ruptures) and mid-sized events (when only sections of the fault rupture), corresponding to the fractal distribution employed for *μ*_{s}. In contrast, the frequency–magnitude distribution for the continuum model shows a Gutenberg–Richter distribution of small events (when only one fault element ruptures) combined with enhanced statistics around a larger characteristic earthquake (when the entire fault length ruptures).

We have shown that mesh dependence in earthquake models gives rise to richly complex spatio-temporal slip over multiple earthquake cycles. To show why rich dynamics is only produced in the discrete models using coarse meshes, we demonstrate how spatial discretization errors can substantially influence the numerical solution of an analogous system: solving a nonlinear one-dimensional wave equation. In the next section, we investigate the origins of complexity in the numerical solution of a nonlinear equation.

## 4. Mesh dependence in nonlinear wave equations

We investigate numerical solutions of the following nonlinear one-dimensional wave equation:
4.1The quadratic nonlinearity in equation (4.1) is known to produce high gradients, shocks in fact (Tabor 1989), not unlike the steep gradients at the front of a propagating discontinuity within the earthquake fault zone. We show that numerical solutions of equation (4.1) also exhibit mesh dependence. We found that solutions using coarse meshes gave rise to soliton-like solutions and not shock solutions because of the presence of spatial discretization errors. The Fermi–Pasta–Ulam oscillator (Fermi *et al.* 1955) is another example of a well-known nonlinear problem that gives rise to soliton solutions.

We will plot the numerical solution of equation (4.1) using different mesh resolutions and methods to demonstrate the complexity that can arise using coarse meshes and low-order numerical methods. We can draw analogies between complexities arising through mesh sensitivity in this nonlinear wave equation and our earthquake model.

To solve equation (4.1), we introduce dimensionless variables *χ*,*η*,*τ*,*δ*, where *x*=*x*_{c}*η*/*δ*, and *u*=(*x*_{c}/*α*)*χ*, with
4.2We consider the effect of different length scales on the solution in the limit .

Comparisons of the numerical convergence of the solution of the nonlinear and linear (*α*=0 in equation (4.1)) wave equation are shown in figures 3 and 4, respectively. We consider periodic boundary conditions and solve for 0<*η*<*L* and *δ*=*L*/128. It is relatively simple to isolate the origin of mesh dependence in the nonlinear system in figure 3. Both the earthquake-fault model and the nonlinear one-dimensional model show mesh sensitivity. The mesh sensitivity arises in both cases because of the nonlinearities present in the governing equations that produce high gradients (shock fronts in the one-dimensional case and a propagating steep discontinuity in the earthquake model), making them extremely difficult to solve numerically and causing ill-posedness. We will show how the numerical solution using coarse meshes for the one-dimensional model in effect changes the governing equation because of the additional error terms that arise with these meshes. This means that a different partial differential equation is actually solved for the one-dimensional model using coarse meshes. In an analogous fashion, we propose that the numerical solution for fault rupture using coarse meshes may also change the underlying physics of the governing equations owing to spatial discretization errors present in the discrete models. In this instance, it would appear that the errors lead to more complex spatio-temporal slip at the fault.

Mesh dependence arises from the nonlinear term in the one-dimensional wave equation. If we look for an asymptotic solution to equation (4.2) of the form *χ*∼*ϕ*(*ζ*,*T*), where *ζ*=*η*−*δτ* and *T*=*δ*^{2}*τ*, we obtain a shock solution (Tabor 1989),
4.3However, we actually obtain soliton-like solutions for the coarse meshes in figure 3—not shock-wave solutions. We investigate why we obtain soliton-like solutions for mesh sizes with Δ_{η}≥*L*/256 by analysing the discretization error introduced in the central difference method we employ to approximate the spatial derivatives numerically. The central difference approximation for the spatial derivatives is only accurate to second order in the mesh spacing. This is typically used in finite-difference and finite-element codes (Olsen-Kettle *et al.* 2008).

We can analyse the error solving the nonlinear shock-wave equation using the central-difference method. To numerically solve equation (4.2), we discretize the spatial region 0<*η*<*L* with *n* mesh elements with spacing Δ_{η}=*L*/*n*. The discretized approximation to equation (4.2) is
4.4where
The reason soliton-like solutions arise from the numerical solution with coarse meshes becomes clear if we consider discretization errors (DEs) up to the order in equation (4.4). This term () and its associated length scale (Δ_{η}) introduces diffusion that acts to limit or smooth the steepness of the wave front. This process is similar to the localization limiters described in Belytschko & Lasry (1989); de Borst *et al.* (1993, 2004); Pasternak & Mühlhaus (2005); Mühlhaus & Vardoulakis (1987); Mühlhaus & Aifantis (1991); de Borst & Mühlhaus (2005). If we look for an asymptotic solution of the form *χ*∼*ϕ*(*ζ*,*T*), where *ζ*=*η*−*δτ*, *T*=*δ*^{5/2}*τ* and *U*=*ϕ*_{,ζ}, we obtain soliton-like solutions when this error term is included in the governing equation (Tabor 1989),
4.5Thus, more complex, soliton-like behaviour results from the numerical solution of the one-dimensional nonlinear wave equation using coarse meshes because the spatial DEs present lead to the solution of the soliton-wave equation (4.5) and not a shock-wave equation (4.3). Equation (4.5) contains an internal length scale (Δ_{η}) in contrast to equation (4.3).

This example clearly shows how different numerical solutions can be obtained when either coarse or fine meshes are employed and low-order finite difference approximations are applied to solve nonlinear partial differential equations. The numerical error introduced by the discretization of equation (4.2) are large enough to affect the solution when coarse meshes are used because these error terms become non-negligible and in effect change the underlying governing equations. Figure 3 clearly shows that soliton solutions arise for meshes where the spacing (Δ_{η}) becomes large enough that the DE influences the physics of the numerical solution. Likewise, for the earthquake model, we showed in figure 1*a* that very different slip profiles are produced when different mesh sizes are employed for the solution. Although this same phenomena of mesh sensitivity is occurring in both cases, the fault model and nonlinear one-dimensional model are completely different and not chosen to resemble each other physically. Both show how nonlinearities in the governing equations can make it extremely difficult to solve when low-order methods and coarse meshes are used. However, as we have seen in these two cases, numerical solutions with coarse meshes can sometimes produce more interesting physics because of the DEs, which can, in effect, change the underlying equations.

Figure 5 shows dramatically improved results if we solve the shock wave equation (4.2) with a higher order numerical method that removes the spatial DEs (up to ), 4.6This demonstrates the importance of spatial DEs when considering highly nonlinear systems and the inadequacy of the central-difference method in figure 3. However, figure 3 illustrates the complexity that can arise through spatial DEs. It is this unintentional complexity that we wish to exploit in earthquake models.

We can check that equation (4.5) is an accurate picture of the asymptotic behaviour of the numerical solution with coarse meshes by plotting the numerical solution of the nonlinear one-dimensional wave equation using progressively higher order methods in figure 6. 0CT refers to the numerical solution just using the central-difference method; 1CT refers to the numerical solution with one higher order corrective term removed ; 2CT refers to the numerical solution with two higher order corrective terms removed ; and 3CT refers to the numerical solution with three higher order corrective terms removed . Figure 6 gives us confidence that even the inclusion of just one higher order perturbative corrective term of order improves the numerical solution of equation (4.2) considerably and that our asymptotic analysis earlier for soliton solutions was reasonable.

Thus, we have shown that the numerical solution of nonlinear partial differential equations with coarse meshes and inadequate schemes can result in mesh dependence and inaccurate results. There are two solutions at hand. The first one requires more computational time and is simple to implement because the same low-order numerical schemes (e.g. central difference) can be used and mesh refinement is used to gain accuracy. The second method we demonstrated here is that higher order numerical schemes can give more accurate results with the same mesh resolution. However, as we demonstrated in both the earthquake-fault model and the nonlinear one-dimensional wave equation, numerical solutions with coarse meshes can produce richer physics owing to the nonlinearities arising from two sources, the first from the existing equations and secondly from spatial DEs present giving rise to nonlinear, mesh-dependent (or scale-dependent) terms.

## 5. Conclusions

Through this simple one-dimensional nonlinear model, we have shown that mesh sensitivity provides a window into discovering new physics underlying nonlinear problems. In both the earthquake model and our simplified one-dimensional model, we have shown that mesh dependence can unravel complex behaviour arising from the numerical solutions of highly nonlinear problems. The internal length scales introduced into the earthquake model through the mesh spacing could be physically interpreted as the dependence of the constitutive response on a characteristic length scale associated with the fault asperities. Cocco & Tinti (2008) proposed that the complex structure of real fault zones implies that earthquakes in such fault structures are dominated by scale-dependent processes. This could bridge the gap between existing theories on the source of earthquake complexity. Analysis of spatial DEs in earthquake models could help constrain a universal law capable of capturing earthquake complexity.

## Acknowledgements

This research is supported by the Australian Research Council Linkage project LP0562686 with the Queensland Department of Main Roads and University of Queensland and AuScope Ltd under the National Collaborative Research Infrastructure Strategy. Computations were performed on the Australian Earth Systems Simulator and the Australian Partnership for Advanced Computing national facility at the Australian National University through an award under the Merit Allocation Scheme. We are grateful for support from a University of Queensland Major Equipment and Infrastructure grant and L.O. is grateful for support from a University of Queensland New Staff grant.

## Footnotes

One contribution of 17 to a Theme Issue ‘Patterns in our planet: applications of multi-scale non-equilibrium thermodynamics to Earth-system science’.

- © 2010 The Royal Society