<figure><img src="image_1.png" /></figure>

HAL
open science

# 3D finite-element mesh for modeling large-scale surface deformation induced by subduction megathrust earthquakes: Application to Chile

Hugo Boulze, Jean-Didier Garaud, Emilie Klein, Luce Fleitout, Christophe

Vigny, Vincent Chiaruttini

► To cite this version:

Hugo Boulze, Jean-Didier Garaud, Emilie Klein, Luce Fleitout, Christophe Vigny, et al.. 3D finite-element mesh for modeling large-scale surface deformation induced by subduction megathrust earthquakes: Application to Chile. Seismica, 2025, 4 (2), 10.26443/seismica.v4i2.1544. hal-05378163

HAL Id: hal-05378163

https://hal.science/hal-05378163v1

Submitted on 25 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A 3D finite-element mesh for modeling large-scale surface deformation induced by subduction megathrust earthquakes: Application to Chile

Hugo Boulze ID¹*, Jean-Didier Garaud ID², Emilie Klein ID¹, Luce Fleitout ID¹, Christophe Vigny¹, Vincent Chiaruttini ID²

¹Laboratoire de géologie - CNRS UMR 8538, École normale supérieure - PSL University, Paris, France, ²DMAS, ONERA, Université Paris-Saclay, 92320, Châtillon, France

**Author contributions:** Conceptualization: H.B., E.K, L.F., C.V.. Methodology: H.B, J-D.G.. Software: H.B., J-D.G., V.C.. Validation: H.B., J-D.G.. Formal Analysis: H.B.. Investigation: H.B.. Writing - Original draft: H.B., J-D.G., E.K, L.F, C.V.. Project administration: E.K., L.F., C.V..

**Abstract** Megaearthquakes (Mw > 8) cause continental-scale, long-lasting surface deformation, mainly due to viscoelastic relaxation of the asthenosphere. To investigate the links between this deformation and the slip history along subduction interfaces—including earthquakes, postseismic slip, and interseismic coupling—large 3D spherical finite-element meshes are required. This technical report introduces the various steps to build Chile_Mesh_v1.0, a customizable mesh for the Chilean subduction zone, designed as a robust platform for testing various viscoelastic rheologies. It spans ~8500 km in longitude, ~7300 km in latitude, encompassing the entire South American plate, and from the surface to 2900 km depth. Special care was taken to reproduce the complex slab geometry, especially in flat-slab regions such as the Pampean and Peruvian segments, following the Slab2 model. We show that accurately modeling both coseismic and postseismic deformation over large scales requires realistic meshed domains, extending down to the Core-Mantle boundary and thousands of kilometers from the trench. In some cases, depth-reduced meshes can be used to model viscoelastic postseismic deformation, but they fail to simultaneously capture coseismic deformation accurately. We hope this open-access mesh proves valuable for researchers studying subduction dynamics in Chile and supports the development of similar models for other regions.

**Production Editor:**

Yen Joe Tan

**Handling Editor:**

Wenbin Xu

**Copy & Layout Editor:**

Kirsty Bayliss

**Signed reviewer(s):**

Francisco

Ortega-Culaciati

Haipeng Luo

**Received:** February 19, 2025

**Accepted:**

July 26, 2025

**Published:**

September 5, 2025

## 1 Introduction

Megathrust earthquakes suddenly induce large stress variations in Earth's mantle. In the shallow part of the sub-lithospheric mantle, called the asthenosphere, these stresses are relaxed through viscous creep. This asthenospheric viscoelastic relaxation induces deformation of the overriding plate, which is both large-scale (thousands of kilometers from the subduction trench) and long-lasting (for at least decades for large events). Finite-element numerical methods have been widely used to explore the parameters that control this viscoelastic relaxation of the asthenosphere during the postseismic phase of the seismic cycle in Chile, Japan, South-East Asia and elsewhere (e.g. Freed and Lin, 2001; Khazaradze et al., 2002; Hu et al., 2004; Freed et al., 2006; Wang et al., 2007; Pollitz et al., 2008; Garaud et al., 2009; Suito and Freymueller, 2009; Hu and Wang, 2012; Moreno et al., 2011, 2012; Hu and Wang, 2012; Wang et al., 2012; Trubienko et al., 2013; Hu et al., 2014; Sun et al., 2014; Trubienko et al., 2014; Li et al., 2015; Bedford et al., 2016; Klein et al., 2016; Freed et al., 2017; Klein et al., 2017; Li et al., 2017; Melnick et al., 2017; Suito, 2017; Sun et al., 2018; Agata et al., 2019; Peña et al., 2019, 2020; Boulze et al., 2022; Nield et al., 2022; Hor-

mazábal et al., 2023; Liu et al., 2023; Lovery et al., 2025; Marsman et al., 2025). Depending on the context of the study, meshes differ in dimensionality (2D or 3D), extent of the meshed domain (lateral and depth), or modeled geophysical zones. The vast majority of the discretized domains extend less than 2000 km from the trench. Most of them investigate the case of an asthenosphere with uniform viscosity from the base of the lithosphere down to the bottom of the mesh situated at depths of 140, 340, 400, or 500 km, where either fixed, free-slip, or free boundary conditions are applied. One notable advantage of employing finite-element models, in contrast to analytical approaches, is their ability to accommodate the presence of strong lateral variations in viscosities, such as the presence of the subducting slab, low-viscosity zones (e.g., wedges or channels) or lateral viscosity gradients within the asthenosphere (Pollitz, 2003). These lateral variations are crucial elements in understanding postseismic and interseismic deformation (e.g. Trubienko et al., 2013, 2014; Klein et al., 2016; Li et al., 2017).

In the study of postseismic deformation following a megaearthquake, a key challenge is to distinguish between different sources of deformation: postseismic slip at seismogenic depths (afterslip), deeper slip along the plate interface (sometimes modeled with low-

<footnote>*Corresponding author: boulze@geologie.ens.fr</footnote>
viscosity channels, Klein et al., 2016), and viscoelastic relaxation in the asthenosphere. These three sources produce deformation that decays differently with distance from the earthquake, making it difficult to separate their contributions if only near-field data are considered. This highlights the need for models of postseismic deformation that are able to accurately reproduce both the elastic response to slip on the interface and the viscoelastic response in the asthenosphere at large distances from the source.

Immediately after an earthquake, postseismic deformation is concentrated near the rupture zone, but over time it spreads across a much wider region. Accurately tracking this deformation at larger spatial scales requires computational meshes that cover broad lateral extents and possibly significant depths. Moreover, modeling the far-field deformation from ancient megathrust earthquakes, such as the 1960 Valdivia event, similarly requires the capability to compute deformation over large distances with precision.

Beyond the transient velocity changes observed in the years following an individual earthquake, understanding deformation throughout the entire seismic cycle is also crucial. Such deformation provides valuable constraints on plate interface coupling, the history of past earthquakes, and consequently, seismic hazard. Computing deformation over the full seismic cycle requires accounting for the viscoelastic response of the system to slip history on the interface, corrected for convergence velocity (backslip), over periods long enough that the interface slip history no longer induces significant deformation (Trubienko et al., 2013).

These considerations motivated us to develop a mesh that discretizes a large domain in all three directions: depth, along-strike, and trench-perpendicular.

Developing a finite-element model of the entire Chilean subduction zone is particularly relevant because it is one of the most seismically active regions in the world. Specifically, during the last decade, three megathrust earthquakes have broken different segments of the subduction zone: the Maule (2010, M<sub>w</sub>8.8), Iquique (2014, M<sub>w</sub>8.1), and Illapel (2015, M<sub>w</sub>8.3) earthquakes. Moreover, the largest registered earthquake, Valdivia (1960, M<sub>w</sub>9.5), broke close to 1000 km of the subduction interface in Patagonia (Cifuentes, 1989; Barrientos and Ward, 1990). Contrary to other subduction zones (e.g. Japan, Indonesia), the presence of a continent adjacent to the subduction trench allows continuous geodetic monitoring of deformation (Klein et al., 2022). To enable the study of viscoelastic deformation at the scale of the South American continent, we have developed a new 3D spherical finite-element mesh, called *Chile_Mesh_v1.0*. It describes the entire Chilean subduction interface, capturing local variations in slab geometry (e.g. Peruvian and Pampean flatslabs) as described by the Slab2.0 model Hayes et al. (2018). The mesh covers the entirety of the South American continent, extending approximately 8500 km in longitude, 7300 km in latitude, and reaching from the Earth's surface to the Core-Mantle Boundary at a depth of 2900 km. The subduction interface itself is meshed continuously from Lima (~80.20° W, ~10.20° S) to the

Chilean Triple Junction (~76.30° W, ~47.40° S). This large-scale mesh enables the study of viscoelastic deformation associated with seismic cycles along the entire Chilean subduction zone. The mesh has been developed and tested using modules from Zset/Zebulon (http://zset-software.com/, last access: 15 June 2025) and MMG (https://www.mmgtools.org/, last access: 15 June 2025).

This paper describes which geophysical zones are integrated into the mesh (Sect. 2) and the methodology we developed to mesh the entire Chilean subduction (Sect. 3). We illustrate the use of this mesh on the simple case of coseismic and postseismic motions triggered by a coseismic slip located somewhere along the subduction interface (Sect. 4). Then, we perform different tests on the extent of the bounding box and their associated boundary conditions (Sect. 5). Finally, we discuss how the mesh can be modified by users (Sect. 6) and future improvements of *Chile_Mesh_v1.0* (Sect. 7).

## 2 Mesh setup: Definition of geophysical zones

In this section, we detail the parameter values (e.g. depth, thickness) selected for the geophysical zones included in the mesh.

**Peru-Chile trench.** To be consistent with the slab topography described by Slab2.0 (Hayes et al., 2018), the longitude-latitude path of the trench is defined as the shallowest points of Slab2.0 between the North of Lima ~(80.20° W, 10.20° S) and the Chilean Triple Junction point ~(76.32° W, 47.37° S). The trench depth is set to 6 km depth in the mesh (see Fig. S1).

**Bathymetry between the trench and the coast.** The distance between the trench and the coast varies along the subduction zone (Fig. 1) but is on average about 100 km (Fig. S1). Therefore, in the mesh, the bathymetry rises from 6 km depth at the trench to sea level at the coast over a distance of 100 km. However, in some places, the actual distance is larger than 100 km and the interpretation of the model predictions in these areas must be made with caution (e.g. ~200 km in the Arica Bend (~71.5° S, ~19° S)).

**Oceanic and continental lithospheres.** Since rheology is temperature-dependent, we use the thermal definition of the lithosphere. In the mesh, the oceanic lithosphere is represented as a layer of constant thickness (50 km Sodoudi et al., 2011) that extends from a depth of 4 km (representing the average depth of the oceanic seafloor, as shown in Fig. S1) down to a depth of 54 km. We also assume a constant thickness for the continental lithosphere, extending from the surface of the Earth to 70 km depth. Indeed, the continental lithosphere in front of the Chilean subduction appears similar to a mature continental lithosphere (Celli et al., 2020) for which a thickness of the order of 70 km seems appropriate.
<figure><img src="image_3.png" /><figcaption>**Figure 1**

a) Portion of the Andean subduction modeled in our 3D mesh. Contour lines (every 20 km) depict the depth of the modeled slab from the trench down to 400 km depth. North of 40° S, the surface of the slab is Slab2.0 (Hayes et al., 2018), South of 40° S, the surface is extrapolated (shaded area). CTJ stands for the Chile Triple Junction point between the Nazca, South American, and Antarctic plates. Extensions of the earthquakes mentioned in the paper are represented with white dashed lines. Inset globe: Tectonic context. N = Nazca plate; SA = South American plate. The arrow indicates the average convergence direction between the Nazca and South American plates in the studied region. The red frame delimits the mesh area (bounding box of the mesh). The blue frame delimits the area represented in the main map, where the subduction interface is meshed. b) Sketch in a 3D flattened frame showing the different geophysical zones modeled for a normal-dipping slab part of the Chilean subduction. Drawing not to scale. c) Top view of *Chile_Mesh_v1.0*. d) Lateral view of a mesh clipped along the profile PP' shown in c, highlighting the variation of slab geometry. The black line depicts the coast of South America. For color codes, refer to b, e) View of the slab meshed in *Chile_Mesh_v1.0*</figcaption></figure>

**Accretionary prism.** The Chilean subduction zone is an erosive continental margin North of Valparaíso (~71.3° W, ~33° S), resulting in minimal sediment accumulation in the trench. South of Valparaíso, the subduction zone is accretionary, with several thousand meters of sediment accumulated since the last deglaciation (e.g Glodny et al., 2006). In the mesh, we introduce a prism zone all along the subduction, which extends 25 km from the trench towards the continent (e.g. von Huene et al., 1997; Contreras-Reyes et al., 2010; Maksymowicz et al., 2018). This prism can be 'activated' in the accretionary region by applying a sedimentary rheology, or 'deactivated' in the erosive region by assigning an elastic crust rheology.

**Slab.** The surface topography of the slab in the mesh

is based on the Slab2.0 model (Hayes et al., 2018). We chose to use this widely adopted model because it provides a continuous slab geometry over a broad region. However, like any model, Slab2.0 has limitations and uncertainties, particularly in southern Chile, where the geometry is poorly constrained due to a lack of seismicity. In Slab2.0, the slab topography is not fully defined south of 40° S due to a lack of seismicity to image it. Still, there is a subduction slab in this region because this southern part of Chile experienced the M<sub>w</sub> 9.5 Valdivia earthquake, which ruptured ~1000 km of the subduction interface in 1960 (Cifuentes, 1989; Barrientos and Ward, 1990). This earthquake must be included in the seismic cycle models to study the deformation of the South American continent. Therefore, a slab must be modeled in this region. For this purpose, we extrapolate the surface of Slab2.0 southwards to the Chilean Triple
Junction (CTJ, near 48° S) (shaded area in Fig.1a). The extrapolation consists in taking the slab profile from Slab2.0 at 40°S — one of the last latitudes where the slab geometry is fully defined from the trench down to 400 km depth — and then copying this profile along the trench path south of 40°S. Because the slab is a continuation of the oceanic lithosphere, we assign it a thickness of 50 km.

Emerging new slab models, such as Contreras-Reyes et al. (2025), show differences of up to 10 km with Slab2.0, which can affect offshore and near-trench predictions. However, a sensitivity test (Fig. S10) shows that a uniform 10 km depth offset of the slab surface between the trench and 40 km depth has a negligible impact on predicted deformation. For coseismic displacements, differences are maximal at 160 km for both horizontal and vertical. Considering a *M*<sub>w</sub>9 earthquake, it represents respectively 29 and 27 mm, which is negligible, given the multi-metric displacements generated there by such an earthquake. They fall below 1 mm for horizontal and vertical further than 1000 km and 600 km (respectively) from the trench, which is also negligible, given the pluri-centimetric displacements at that distance from the source. For the postseismic velocities 5 yrs after the earthquake, the difference reaches a maximum of 3.5 mm/yr at around 130 km, but it drops below 1 mm/yr beyond 300 km from the trench. For the vertical component, the difference remains below 1 mm/yr across the entire profile. Again, it is negligible given the pluri-millimeter velocities along the profile.

**Asthenosphere.** The viscosity stratification of the upper mantle plays a crucial role in controlling the pattern of postseismic deformation. However, both the thickness of the low-viscosity asthenosphere beneath the lithosphere and the very existence of a viscosity gradient with depth in the upper mantle remain subjects of debate.

Postseismic deformation is often modeled with an asthenosphere about 100 km thick (e.g. Pollitz et al., 2006; Broerse et al., 2015), corresponding to the thickness of the low-quality factor Q (implying low short-term viscosity) and low-velocity seismic zone (e.g. Resovsky et al., 2005; Dalton et al., 2008). Trubienko et al. (2014) demonstrated that the shape of the ratio of postseismic to coseismic deformation (the *post/co* ratio) as a function of distance from the trench provides valuable information about asthenospheric thickness: the thicker the asthenosphere, the farther from the trench the maximum of the *post/co* ratio is located. They also show that observations following recent major earthquakes are most consistent with an asthenosphere extending to depths of approximately 200–300 km. Additional arguments for a strong increase of viscosity with depth in the upper mantle come from comparisons between long-term viscosity inferred beneath cratons by glacial isostatic adjustment studies (e.g. Argus et al., 2021) and the viscosity required for small-scale convection to match the observed heat flow at the base of the plates (Dumoulin et al., 1999). However, many finite-element models of postseismic deformation assume a constant viscosity

extending to depths of 400–500 km, typically corresponding to the bottom boundary of the mesh in these studies (e.g. Hu et al., 2004; Moreno et al., 2011; Li et al., 2017; Peña et al., 2019).

To address these issues, in the mesh, the asthenosphere is divided into two layers: an upper layer and a lower layer, providing the possibility to modify the effective thickness of the asthenosphere. The upper layer has different thicknesses because the oceanic and continental lithospheres do not have the same thickness. Thus, it extends from 54 to 200 km depth on the ocean domain and from 70 to 200 km depth on the continental domain. The lower layer is identical in thickness and depth range for both the oceanic and continental domains, extending from 200 to 270 km.

**Channel.** The presence of a low-viscosity channel located immediately above the slab was suggested by Klein et al. (2016) after the Maule earthquake. Such a channel has a significant influence on vertical surface deformation. Therefore, we define a ~15 km thick channel above the slab. To explore the potential depth this channel could reach, we divide it into four segments: CH1 (extending from approximately 35 km to 70 km depth), CH2 (from 70 km to 200 km depth), CH3 (from 200 km to 270 km depth), and CH4 (from 270 km to the bottom of the slab). Each channel segment can be activated or deactivated by changing its assigned viscosity.

**Deeper mantle.** Below 270 km depth, the mantle is divided into two layers: the deep upper mantle and the lower mantle. We define the deep upper-mantle layer between 270 km and 670 km depth. The 670 km depth corresponds to the discontinuity between the upper mantle and the lower mantle (Ringwood and Iribune, 1988; Deuss et al., 2006). The lower mantle is defined as the layer extending from 670 km to 2900 km (Core-Mantle boundary).

**Extent of the meshed domain.** To cover the South American continent, the mesh bounding box extends from 105° W to 20° W in longitude (~8500 km) and from 13° N to 60° S in latitude (~7300 km). The bounding box ranges from 0 to 2900 km depth. The significance of a large mesh domain is addressed in Sect. 5.

## 3 Meshing methodology

In this section, we describe the step-by-step process involved in constructing the Chilean subduction mesh (Fig. 2).

**Creation of the subduction interface (Fig. 2a).** We first build a slice of the subduction in a 2D flat frame perpendicular to the trench direction, including some geophysical zones: the slab (with a constant dipping angle at this stage), the channel, and the accretionary prism. The 2D slice is then meshed and referenced in a 2D spherical reference frame where the center of the Earth is defined as the center of a sphere with a radius of 6371 km (mean of the Earth's radius).

**Extrusion into a 3D spherical frame (Fig. 2b).** The 3D spherical subduction interface is created by extruding the 2D slice along the trench path. As in the previous
<figure><img src="image_4.png" /><figcaption>**Figure 2** Sketch of the meshing pipeline. **a-2D** sketch of the subduction interface (2D, flat frame). **b-3D** extrusion of (a) along the path of the Peru-Chile trench (3D, spherical frame). **c-Subduction** interface after the projection of the slab on Slab2.0 (3D, spherical frame). **d-e-f-Insertion** of the subduction interface into the mesh bounding box (3D, flat frame). **g-Insertion** of the 4 km, 54 km, 70 km, 200 km, 270 km and 670 km planes (3D, flat frame). **h-Insertion** of vertical planes to delimit continental and oceanic domains. **i-Reassignment:** Creation of the final geophysical zones. **j-Remeshing** process and final conversion to 3D spherical frame.</figcaption></figure>

step, the dip of the slab remains constant throughout this step.

**Projection of the slab on Slab2.0 (Fig. 2c).** Using an iterative process, we deform and mesh the geophysical zones to align the slab surface with Slab2.0. The iterative approach is crucial in our methodology. Previous attempts to achieve deformation in a single step resulted in poor-quality elements (e.g., high aspect ratios), particularly in regions with rapid changes in slab curvature

(e.g. in flat-slab zones where the slab dip is not constant with depth). At the end of the process, we verify the slab and channel thicknesses. While the slab's average thickness (~48 km, Fig. S2) is well-preserved with minimal dispersion, the channel thickness is increased on average by 36% compared to the initial 2D specifications (11 km before the projection, 15 km after) and exhibits greater variability (Sect. S2). This point is discussed in more detail in Sect. 7.
**Insertion of the subduction interface into the bounding box (Fig. 2d-e-f).** We first create a bounding box in a 3D flat frame (Fig. 2d). Then, we triangulate the surface of the deformed subduction interface from the previous step (first converted into a 3D flat frame) to connect it with the top corners of the bounding box (Fig. 2e). At the end of this step, all geophysical zones are remeshed (Fig. 2f).

**Insertion of the horizontal and vertical planes (Fig. 2g-h).** Using the Zcracks module, we create layers in the bounding box by inserting horizontal planes at specific depths: 4, 54, 70, 200, 270 and 670 km (Fig. 2g). Note that non-plane boundaries can also be inserted at this stage, for example, for creating thinner (below the volcanic chain) or thicker (cratons) lithospheres. We separate the oceanic and the continental parts of the lithosphere and of the asthenosphere by introducing two vertical planes between the extremities of the subduction interface to the corners of the bounding box (Fig. 2h).

**Reassignment and Remeshing. (Fig. 2i-j)** Because the mesh has been sliced with horizontal and vertical planes, we need to reassemble the different geophysical zones in the mesh, using different techniques (e.g. renaming, connectivity between zones). We then remesh all the geophysical zones using MMG software. The mean size *s* of an element edge is defined with distance *d* from the fault plane: refined elements near the trench (in Chile-Peru coast), where stress and strain gradients are more important, coarser further away from the trench, to limit computation costs. The mesh's sDomap function is therefore divided into three parts: near-field elements (*d* < ~ 130 km, *s* ~ 25 km); mid-field elements (~ 650 km > *d* > ~ 130 km, *s* ~ 60 km); far-field elements (*d* > ~ 650 km, *s* ~ 250 km). We then create the subduction fault plane, inserting a crack in the mesh by duplicating the nodes between 6 km and 95 km depths. Finally, the mesh is referenced in a 3D spherical frame (Fig. 2j).

At the end of the meshing process, our preferred mesh has 137 330 nodes (i.e. 411 990 degrees of freedom) and 775 129 linear tetrahedral elements (Fig. 1c-d-e). Calculations require ~ 9 GB of RAM using Zset's direct linear solver.

# 4 Example of co- and post-seismic modeling

In this section, we test *Chile_Mesh_v1.0* through a numerical simulation of a hypothetical earthquake. We first describe the boundary conditions and how the fictive earthquake is modeled. Then, using specific rheological parameters, we analyze the resulting coseismic and postseismic deformation.

## 4.1 Boundary conditions

The mesh includes six sides: four lateral sides (North, South, East, and West) and two base sides, with one representing the Earth's surface and the other the Core-Mantle boundary (hereafter referred to as the bottom of the mesh). In this simulation, all four lateral sides

are fixed. The bottom of the computational domain has a free-slip condition (differences with using fixed or free conditions are negligible). A Robin-Winkler condition is introduced to simulate vertical gravity for the Earth's surface. This condition is a surface force equal to $-\rho g U_r$, where $\rho$ is the material density, $g$ the gravity acceleration, and $U_r$ the induced displacement oriented towards the center of the Earth. We set $\rho = 3 \text{ g/cm}^3$ (density of mantle rocks) for the following numerical simulations. This technique of 'simplified gravity' neglects several physical terms (e.g. self-gravitation, density variations induced by interface distortion) present in the elastic gravitational free oscillations (Backus, 1967) or Glacial Isostatic adjustment theories (e.g. Wu and Peltier, 1982; Cathles, 2015). However, it has been shown through comparison based on spectral methods that this is a good approximation (Pollitz, 1997). The impact of this 'simplified gravity' is shown in supplementary materials (Fig. S4). We also assume a purely continuous model: the interface between two geophysical zones of different rheologies does not have a specific behavior (e.g. no friction laws).

## 4.2 Earthquake simulation

In the following, earthquake slip model refers to a static slip distribution along the fault-plane. The earthquake is imposed through specific multi-point constraints (MPC) by suddenly forcing a relative motion along the fault plane between the subducting and the overriding plates. For that purpose, we use the duplicate node technique (Melosh and Raefsky, 1981): along the dislocation, the nodes are duplicated into two groups that move relative to each other. When studying real earthquakes, the earthquake slip model is obtained through the inversion of different surface observations of deformation generated by the earthquake (so-called coseismic displacement): e.g. GNSS, InSAR, seismological data (accelerograms and seismograms), and tsunami heights. For a large earthquake, its associated coseismic slip model is usually complex with heterogeneous slip. While it is crucial to consider slip variability to understand the deformation pattern in the near-field, its impact diminishes significantly with distance, playing a minor role in the mid-field and an even smaller role in the far-field. For numerical testing, we simulate a hypothetical Mw9 earthquake characterized by a uniform slip, where all nodes slip by the same amount and in the same direction across a 400 km rupture, spanning latitudes from 34° S to 38° S and depths from 6 km to 70 km (Fig. 3a and Fig. S3). Each node slips by 9 m parallel to the 12° E direction, which is the convergence direction between the Nazca and South American plates at this latitude (Vigny et al., 2009).

## 4.3 Calculating the coseismic and postseismic displacements

We perform a numerical simulation by applying the uniform coseismic slip previously defined on the fault plane and computing the induced coseismic and postseismic deformation at the surface of the Earth. For
<figure><img src="image_5.png" /><figcaption>**Figure 3** **a-** Coseismic displacement field produced by a fictive M<sub>w</sub>9 earthquake. **b-** Postseismic velocity field 5 yrs after the fictive M<sub>w</sub>9 earthquake. For **a** & **b**, the color scales depict the vertical motion and the arrows the horizontal motion. Note that scales are different between coseismic displacements and postseismic velocities. The vertical color scales are saturated to make far-field displacements visible. **c-d-e-** Synthetic time-series of displacements for three points located at different distances from the trench (yellow dots in **a**. and **b.**) aligned along a profile at 36° S (orange dashed line). NF: near-field point (~200 km from the trench); MF: mid-field point (~500 km from the trench); FF: far-field point (~1500 km from the trench).</figcaption></figure>

numerical testing, we use simple Maxwell rheologies for the various geophysical zones considered to be viscoelastic (Table S1). In this test, the asthenosphere (3 × 10<sup>18</sup> Pa.s) extends from the base of the continental and oceanic lithospheres down to 270 km depth (Fig. 1b). We set the viscosity in the channels equal to that of the asthenosphere.

We run the simulation over 100 yrs, computing the state of the model (i.e. stresses at integration points, displacements at nodes) following an evolving and adaptive time step that varies based on the time elapsed since the earthquake. In the early postseismic phase, time steps are small (in the order of days to weeks) to better capture the non-linear transient deformation, gradually increasing as time goes on (up to 10 yrs). The finite-element results (displacements) are interpolated through time at specific points. For real earthquakes, these points typically correspond to GPS station locations. For synthetic events, as in this study, they can be selected on grids or along profiles. Here, we interpolate displacements on a 0.5° regular grid across South America and along a 0.1° profile at 36° S (the average latitude of the uniform coseismic slip distribution). For both the profile and grid, we choose to keep only points located inland.

**Coseismic displacements (Fig. 3a).** On the continent, coseismic horizontal displacements exhibit a rather simple pattern: they all point seaward (to the West) and simply decrease in amplitude with distance to the trench. Vertical displacements follow a more complex pattern: overall, they also decrease with distance to the trench but alternate between uplift along the coast, subsidence in the Andes mountain range, and uplift again further away inland.

**Postseismic displacements (Fig. 3b).** Facing the rupture zone, horizontal postseismic velocities (calculated 5 years after the earthquake) display a pattern similar to that of the coseismic displacements: velocities are oriented seaward (Fig. 3b). However, contrary to the coseismic displacements, horizontal postseismic velocities do not simply decrease with distance to the trench: they first increase with distance before decreasing further inland. At the North and South of the rupture zone, horizontal postseismic velocities describe rotations: clockwise in the North and counterclockwise in the South. The vertical postseismic velocities decrease with distance to the trench, alternating between regions of subsidence and uplift. The shape of these regions and the precise position of the uplift/subsidence transition depend on the chosen rheological parameters. For
example, Klein et al. (2016) showed that a subduction channel with a ~ $10^{17}$ Pa.s viscosity extending from 55 to 135 km depth is needed to fit the uplift in the Andes after the Maule earthquake.

# 5 Impact of the mesh extent

Finite-element meshes used to simulate earthquake-related deformation typically cover a limited volume of the Earth's crust and mantle around the seismically active region of interest. To avoid numerical artifacts, it is critical that the meshed domain be large enough that the simulation results remain unchanged even if a larger domain were used. In this technical report, we do not address the validity of specific physical parameters, such as asthenospheric thickness and viscosity or the choice of rheological models (e.g., Maxwell, Burgers). Our focus is on determining the optimal extension of the finite-element mesh required to ensure numerical accuracy. A key question we explore is whether the mesh extent can be significantly reduced for certain types of applications without compromising the reliability of the results. Given that most existing studies use mesh extents much smaller than the one considered here, our analysis also serves to assess the spatiotemporal limits of validity for simulations based on these reduced meshes.

In this section, we investigate the impact of the extent of the meshed domain on the modeled surface deformation. For all numerical experiments, we compare the results obtained using the Earth realistic mesh *Chile_Mesh_v1.0* with those from simulations employing reduced mesh extents, using identical rheological parameters within the overlapping domain. Due to the large extent of its meshed domain, *Chile_Mesh_v1.0* is assumed to produce results free from numerical artifacts associated with the insufficient extent of the meshed domain.

To perform these experiments, we employ the synthetic $M_w$9 coseismic slip distribution defined in Sect. 4.2. For each test, we analyze coseismic displacements due to elastic deformation and postseismic velocities due to viscoelastic deformation in the asthenosphere. Both horizontal and vertical components are considered; the horizontal component refers to the norm of the horizontal displacement vector. All tests in this study are conducted using Maxwell rheology, although Burgers rheology is often preferred for modeling the postseismic phase (Klein et al., 2016; Llovery et al., 2025). While the choice of rheology may influence the time-dependent evolution of velocity, it does not significantly affect the model's sensitivity to boundary conditions.

We focus on postseismic velocities at two key times: 5 yrs and 95 yrs after the earthquake to limit the number of figures. Examining the effects of slip on the interface long after an earthquake is crucial, as a megathrust earthquake that occurred several decades ago can still influence deformation over a wide area. For example, deformation from the 1960 Valdivia earthquake remains detectable today in Argentina (Klotz et al., 2001).

Furthermore, the long-term response to slip on the interface is a major factor in calculating interseismic velocity, as this integrates the asthenosphere's response to both backslip (Savage, 1983) and past earthquakes over extended periods (Trubienko et al., 2013). It is therefore essential to test the mesh to ensure that deformation induced by slip events that occurred a long time ago is accurately modeled.

## 5.1 Depth extension of the mesh

We investigate the impact of the mesh depth. In these tests, we consider two types of depth dependence for the viscosity. First, as mentioned in the introduction, most studies assume a uniform asthenosphere beneath the lithosphere, extending to the base of the mesh situated at depths varying between 140 km and 500 km. To assess the impact of boundary conditions on this commonly used configuration, we test models with uniform viscosity down to the base of the mesh, choosing two depths for the bottom boundary: 270 km (**model A**) and 670 km (**model B**). Second, we examine the case where viscosity increases with depth. Specifically, we consider a mesh extending to 670 km depth, but with a viscosity of $3 \times 10^{18}$ Pa.s above 270 km and $3 \times 10^{20}$ Pa.s below (**model C**). We compare model A and C with the full-depth model (**model R1**, down to 2900 km) featuring an asthenosphere down to 270 km, presented in Sect. 4. We compare model B with a full-depth model (reference **model R2**) featuring an asthenosphere down to 670 km depth. Here is a resume of the 5 tested models (see also on top of Fig. 4):

* **Model A:** *Chile_Mesh_v1.0* truncated at 270 km depth with an asthenosphere down to 270 km depth

* **Model B:** *Chile_Mesh_v1.0* truncated at 670 km depth with an asthenosphere down to 670 km depth

* **Model C:** *Chile_Mesh_v1.0* truncated at 670 km depth with an asthenosphere down to 270 km depth

* **Reference model R1:** full-depth mesh (i.e. down to 2900 km), with an asthenosphere down to 270 km (same as in Sect. 2)

* **Reference model R2:** full-depth mesh, with an asthenosphere down to 670 km

Detailed rheologies for all models are provided in Tables S1 and S2. For models A, B, and C, we test three different boundary conditions (hereafter referred to as bc) for the bottom of the mesh (expressed in spherical coordinates below):

1. fixed: $U_r = U_\theta = U_\varphi = 0$, displacement vector is zero in all directions

2. free: $\sigma_{rr} = \tau_{r\theta} = \tau_{r\varphi} = 0$, the components of the stress tensor normal and tangential to the bottom boundary vanish

3. free-slip: $U_r = \tau_{r\theta} = \tau_{r\varphi} = 0$, no radial displacement, and the tangential components of the stress tensor vanish
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4** Panel A. Schematics of the different tested models. Panel B. Horizontal (a.) and vertical (b.) predicted coseismic displacement for model A (red lines) and model R1 (gray line). Horizontal (c.) and vertical (d.) predicted coseismic displacements for model B and C (same displacements, blue lines) and model R1 and R2 (same displacements, gray line). Coseismic displacements are equal, whatever the boundary conditions, between models R1 and R2, and between models B and C, since the viscous asthenosphere is not playing a role in the determination of the coseismic offset, purely elastic. Boundary conditions at the bottom of models A, B, and C vary: i) fixed (solid line), ii) free (dotted line), and iii) free-slip (dashed line). Coseismic displacements are calculated along an inland profile at 36° S (dashed orange line in Fig. 3a).</figcaption></figure>

**Effect on coseismic displacements.** At first view, the boundary condition at the base of the mesh seems to have limited impact on the horizontal and vertical coseismic displacements (Fig. 4). The free bc induces the largest differences between A, B, C, and the refer-

ence models: e.g., overestimation for model A both in horizontal and vertical components. At large distances to the trench (> 1000 km), the relative differences for the three bc become significant (Fig. S5) for the three meshes: differences with the reference models reach
<figure><img src="image_8.png" /><figcaption>**Figure 5** Horizontal (a, b) and vertical (c, d) predicted postseismic velocities for model A (red lines) and R1 (gray line). Panels on the left (a, c) show velocity at +5 yrs after the earthquake, and panels on the right (b, d) at +95 yrs. Boundary conditions at the bottom of model A vary: i) fixed (solid line), ii) free (dotted line), and iii) free-slip (dashed line). Postseismic velocities are computed along an inland profile at 36°S (dashed orange line in Fig. 3).</figcaption></figure>

several cm, which is about the same order of magnitude as the displacements (Note that models B and C have the same coseismic displacements, since the elastic moduli of all the geophysical layers are the same and the viscous behavior of the different configurations does not play any role in the prediction of the coseismic displacement).

**Effect on postseismic velocities.** For model A (Fig. 5), the fixed bc is the closest to the reference model, exhibiting non-negligible discrepancies. At 5 yrs, differences reach up to 10 mm/yr in both components: this configuration also does not fully capture the viscoelastic response of the system. At 5 yrs, the *free* and *free-slip* bc yield similar horizontal velocity patterns, but their decay beyond 800-1000 km from the trench is significantly slower than in the reference model R1. For the vertical component, discrepancies are more pronounced; for example, both *free* and *free-slip* bc predict subsidence around 400 km from the trench, which is not visible in the reference model R1. At 95 yrs, these two boundary conditions clearly underestimate horizontal velocities between 400 and 1000 km (detailed differences in Fig. S6).

For model B (Fig. 6), the fixed bc features minor differences with reference R2, unlike the free or free-slip bc: also with a relatively thick asthenosphere, the free-slip

or free bc predict postseismic deformations very different from model R2 (detailed differences in Fig. S7), both in terms of amplitude and of pattern. This is also true at short distances to the trench concerning the vertical or the horizontal velocity at 95 yrs. For a fixed bc the postseismic deformations are similar.

For model C (Fig. 7), differences compared to the reference model R1 are minor. At 5 yrs, some differences appear mainly for the free bc for both horizontal (up to 2.5 mm/yr at 700 km) and vertical (up to 2 mm/yr at 1000 km) components. For all other bc and times, velocity differences remain below 1 mm/yr (detailed differences in Fig. S8).

**Effect on post/co ratios.** For all models, we calculate the *post/co* ratio (cumulative postseismic displacement over 5 yrs, normalized by the coseismic displacement) as a function of the distance to the trench (Fig. 8) (Trubienko et al., 2013, 2014; Boulze et al., 2022). As mentioned earlier, the observed *post/co* ratios describe for all subduction earthquakes a similar bell shape with a maximum between 600-800 km from the trench (Trubienko et al., 2014; Boulze et al., 2022). The modeled slope of *post/co* curve at large distances to the trench is characteristic of the thickness of the asthenosphere (Trubienko et al., 2014). For a thick asthenosphere, the slope stays positive, as can be seen for the reference
<figure><img src="image_9.png" /><figcaption>**Figure 6** Same caption as Fig. 5 but for model B (light blue lines) and model R2 (gray line).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7** Same caption as Fig. 5 but for model C (blue lines) and model R1 (gray line).</figcaption></figure>
model R2 of (Fig. 8b). All models A, B, and C show divergence from their respective reference models, regardless of the applied bottom boundary condition. What distinguishes them is the distance from the trench at which these deviations become significant. Model A (Fig. 8a), which has the smallest depth extent, diverges the earliest, with noticeable differences appearing as close as 300 km from the trench for all boundary conditions. In Model B (Fig. 8b), divergence occurs between 300 km and 500 km, depending on the boundary condition: around 300 km for the free-slip case, and between 400–500 km for the fixed and free conditions. Model C (Fig. 8c) maintains good agreement with the reference up to greater distances. Divergence begins at approximately 500 km for the free boundary, 600 km for the fixed condition, and as far as 900 km for the free-slip case.

**Conclusion of the test.** Reducing the mesh depth has a significant impact on the predicted deformation. All models that exclude the lower mantle, and in the case of Model A, also part of the deeper upper mantle, introduce substantial differences, regardless of the bottom boundary condition applied. These discrepancies arise because shallower models fail to capture the long-wavelength response of the deeper mantle, which plays a critical role in reproducing the coseismic signal at large distances from the trench. To accurately model both coseismic displacements and their postseismic evolution, at least the upper portion of the lower mantle must be included in the computational domain.

## 5.2 Lateral extension of the mesh

In this section, we evaluate the impact of lateral reduction of the mesh bounding box in both continental and oceanic domains. Because there are no pre-existing geophysical zones whose exclusion would simplify this analysis, we introduce fictive zones within the mesh using Zset's mesher module. To simulate the lateral reduction, we deactivate these zones in the displacement calculations by assigning them properties of infinite elasticity (PREM bulk and shear moduli from Table S3 increased by a factor of $10^4$) and viscosity ($3 \times 10^{30}$ Pa.s). We consider two reduced models, which will be both compared to the reference model (full-lateral extent, R1, Sect. 5.1):

* **Model D: Chile_Mesh_v1.0** reduction of the continental domain beyond 1550 km from the trench (instead of ~5000 km), with a fictive zone between 57° W and 20° W (Fig. 9, inset 2, red zone).

* **Model E: Chile_Mesh_v1.0** reduction of the oceanic domain beyond approximately ~1000 km from the trench (instead of ~3000 km), with a fictive zone between 105° W and 85° W (Fig. 9, inset 3, blue zone).

In both cases, these fictive zones extend from the surface to 2900 km depth. The boundary conditions and rheologies of activated zones remain as described in Sect. 4.

**Effect on coseismic displacements.** Although coseismic displacements appear similar (Fig. 9a-b), model D

<figure><img src="image_11.png" /><figcaption>**Figure 8** a) Eastward post/co ratios at 5 yrs predicted for model R1 (gray line) and model A (red lines), with three different bottom boundary conditions: (i) fixed (solid), (ii) free (dotted), and (iii) free-slip (dashed). b) Same as a), but for model R2 (gray) and model B (light blue lines). c) Same as a., but for model R1 (gray) and model C (blue lines) with three different bottom boundary conditions</figcaption></figure>

shows differences of 10 mm at 100 km from the trench (negligible given the multi-metric scale of displacements), increasing to 49 mm at 1500 km, which becomes significant where displacements are only a few centimeters. For the vertical component, model D shows differences of only a few millimeters in the near field, which remain negligible relative to the metric deformation. However, at 1400 km, the difference exceeds 6 mm
(more than 500% relative to the reference model, which
induces just 1 mm of displacement). In model E, differ-
ences are negligible at both 100 km (17 mm on a multi-
metric displacement) and 1400 km, where they remain
below 5 mm, insignificant relative to multi-centimetric
offsets. Vertical differences are below 1 mm through-
out. See Fig. S9-a-b for detailed differences.

**Effect on postseismic velocities.** At 5 years (Fig.9c-d), model D predicts significantly slower horizontal velocities compared to the reference model R1, with differences exceeding 8 mm/yr at 1400 km from the trench, over 50% slower than R1. In the vertical component, however, model D remains within 1.5 mm/yr of R1. At 95 yrs, model D continues to show substantially slower horizontal velocities, with a difference of about 3 mm/yr at 1400 km, corresponding to nearly a 100% relative difference. In contrast, model E exhibits negligible deviations from R1 in both horizontal and vertical components at both 5 and 95 yrs (Fig.9e-f), with differences below 1 mm/yr. See Fig. S9-c-d-e-f for detailed differences.

**Effect on post/co ratios.** Like in Sect. 5.1, we compute the *post/co* ratio predicted by models D, E, and R1 (Fig. 10) at 5 yrs after the earthquake. The *post/co* ratio of model E slightly diverges from model R1, showing the limited impact on the continental deformation of reducing the mesh on the oceanic side. On the contrary, the *post/co* ratio of model D diverges significantly from the reference beyond 700-800 km to the trench.

**Conclusions of the test.** The model truncated at ~1550 km from the trench in the continental part is unsuitable for analyzing deformation beyond ~ 700-800 km. The model truncated on the oceanic part has a limited impact on predicted deformation. Ensuring a sufficiently large model bounding box is essential to avoid perturbations in estimating both coseismic and postseismic deformation.

# 6 Adaptability of the mesh

The mesh *Chile_Mesh_v1.0* is designed to study large-scale and long-term deformation throughout the seismic cycle accounting for far-field deformation induced by slip history along the various segments of the Chile subduction. However, it remains highly customizable for user-specific needs. The various geophysical zones defined in the mesh enable the testing of multiple rheological combinations. Furthermore, in *Chile_Mesh_v1.0*, we chose to assign constant thickness and depth parameters to each geophysical zone along the entire subduction. Without repeating the insertion of horizontal planes described in Sect. 3, these parameters can be first-order modified by manually adjusting the mesh nodes.

Although a synthetic earthquake was used for this
study, the mesh is also suitable for computing Green's
functions (e.g. dip and strike components) for each
node on the fault plane. This functionality enables the
inversion of coseismic slip distribution using geodetic
data, such as GNSS and InSAR. In addition, depending
on the finite-element software used, the mesh can be
modified or used in different ways, providing flexibility

for various applications and modeling needs.

**Using Zset/Zebulon.** Users can easily adjust rheological parameters by modifying the associated *.mat files for each geophysical zone. Using Zset's mesher module, new geophysical zones can be defined based on latitude, longitude, and depth conditions, even within existing elements (an example is provided in Sect. 5.2), and so without reconstructing the entire mesh. In regions with coarser elements, finite element precision can be improved, if needed, by using a meshing tool such as MMG.

**Using other finite-element software.** The mesh is fully compatible with other finite-element solvers/meshers than Zset/Zebulon, including Gmsh (https://gmsh.info/, last access: 15 June 2025), Abaqus (https://www.3ds.com/products/simulia/abaqus, last access: 15 June 2025), or PyLith (https://github.com/geodynamics/pylith, last access: 15 June 2025). Users can define their boundary conditions (e.g. gravity) and rheological parameters or create (or modify) geophysical zones in their preferred software.

# 7 Discussions and conclusion

## 7.1 Impact of the extent of the meshed domain on coseismic and postseismic deformation

As shown in Sect. 5, the size of the Earth's domain en-
compassed by the mesh has a non-negligible impact
on both postseismic and coseismic deformation. Ac-
cordingly, the large dimensions of the *Chile_Mesh_v1.0*
model are essential for accurately modeling deforma-
tion across large spatial scales and long timescales, par-
ticularly to capture the full complexity of slip history
along the Chilean subduction interface.

Many studies have employed 2D or 3D models with domains significantly smaller than that of *Chile_Mesh_v1.0* to investigate the seismic cycle (including coseismic, postseismic, and interseismic deformation) in regions such as Chile, Peru, Japan, Indonesia, and elsewhere (e.g. Freed and Lin, 2001; Khazaradze et al., 2002; Hu et al., 2004; Freed et al., 2006; Wang et al., 2007; Pollitz et al., 2008; Suito and Freymueller, 2009; Hu and Wang, 2012; Moreno et al., 2011, 2012; Wang et al., 2012; Hu et al., 2014; Sun et al., 2014; Li et al., 2015; Bedford et al., 2016; Freed et al., 2017; Li et al., 2017; Melnick et al., 2017; Suito, 2017; Sun et al., 2018; Agata et al., 2019; Peña et al., 2019, 2020; Hormazábal et al., 2023; Liu et al., 2023; Loveroy et al., 2025; Marsman et al., 2025). The numerical validity of these studies is left to the reader's judgment, based on both the sensitivity tests presented in Sect. 5 and the specific characteristics of each study: e.g., such as the distance from the trench of the observational data used, the spatial and temporal extent of the computational domain, and the applied boundary conditions. We provide below general recommendations to guide future modeling efforts.

For the lateral extent of the mesh, the farther the
<figure><img src="image_12.png" /><figcaption>**Figure 9** Horizontal (a) and vertical (b) coseismic displacements for model D (red dashed line), model E (blue line), and model R1 (gray line). Horizontal (c, e) and vertical (d, f) postseismic velocities for models D, E, and R1. Same color code. Panels (c–d) show velocities at 5 yrs, and (e–f) at 95 yrs after the earthquake. The colored area in the insets for models D and E indicates the deactivated regions of the mesh: Only the grey area is used in the computations. Both coseismic and postseismic deformations are evaluated along an inland profile at 36°S (dashed orange line in Fig. 3).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 10** Eastward *post/co* ratio at 5 yrs on the East component as a function of distance to the trench for model D (red dashed line) and E (blue solid line) compared to model R1 (gray solid line). The colored area in the insets for models D and E indicates the deactivated regions of the mesh: Only the grey area is used in the computations. Both coseismic and postseismic deformations are evaluated along an inland profile at 36°S (dashed orange line in Fig. 3 of the main paper).</figcaption></figure>

observation points are from the model boundary in the continental part (> 1000 km), the more accurate the modeled deformation becomes. The position of the lateral boundary on the oceanic side has limited influence on deformation within the continent, but it does affect deformation on the oceanic plate. A wide mesh in the oceanic direction also enables the inclusion of GNSS data from islands located on the subducting plate (Nazca, in our case) at various distances from the

trench, providing valuable constraints on seismic cycle processes.

Recommendations regarding the depth extent of the mesh:

**Coseismic deformation.** To obtain reliable estimates of coseismic deformation, the depth of the meshed domain must significantly exceed the horizontal distance between the trench and the farthest observation sites (e.g., GNSS stations or InSAR pixels). For instance, if geodetic data extend up to 1000 km from the trench, the meshed domain should reach depths greater than 1000 km (ideally up to or into the lower mantle). Deep meshes are therefore critical for accurately capturing both far-field coseismic deformation and the afterslip component of postseismic deformation.

**Postseismic deformation due to viscoelastic relaxation.** Accurate modeling of viscoelastic postseismic deformation can be achieved using one of the following configurations:

(i) a full-depth mesh that extends into the lower mantle, as in *Chile_Mesh_v1.0*,

(ii) a moderately deep mesh that includes a high-viscosity layer beneath the asthenosphere down to the mesh base (e.g., Model C in Sect. 5.1),

(iii) a shallower mesh that ends at the base of the asthenosphere, provided a fixed boundary condition is applied at its bottom (e.g., Models A and B with fixed bc in Sect. 5.1).

In contrast, shallow meshes truncated at the asthenosphere boundary (without any higher viscous layer below) with free-slip or free boundary conditions at the
bottom (models A and B in Sect. 5.1) fail to constrain the viscoelastic response correctly and should be avoided.

## 7.2 Modeling the seismic cycle using Chile_Mesh_v1.0

We have shown in this paper that *Chile_Mesh_v1.0* can be used to model coseismic and postseismic deformation. While interseismic deformation is not explicitly addressed here, the mesh can also be used for such modeling. Interseismic modeling remains a topic of debate, out of the scope of this paper, particularly regarding whether elastic or viscoelastic rheologies should be used. In either case, computing the interseismic deformation essentially involves calculating the response to a unit slip along the fault plane, and then using backslip theory (Savage, 1983) to determine the interseismic deformation.

Estimating the interseismic signal using viscoelastic rheologies requires computing the postseismic deformation over time-scales exceeding several decades after the earthquake (Trubienko et al., 2013): This requires a mesh large enough to avoid perturbation in the determination of the postseismic deformation over long time-scales (plots for 95 yrs in Sect. 5). Moreover, in the case of Chile, studying interseismic deformation in the Maule region requires accounting for the huge postseismic effects of the 1960 Valdivia earthquake even 60 years later (Moreno et al., 2011).

In the case of the Chilean subduction zone, *Chile_Mesh_v1.0* offers a robust framework to ensure that technical modeling choices, such as the extent of the meshed domain and boundary conditions, do not bias the interpretation of long-term deformation processes.

## 7.3 Future versions

The current version is named *Chile_Mesh_v1.0* to reflect its potential for future improvements. Subsequent iterations of the mesh will incorporate user feedback, minor local enhancements, and new constraints on geophysical zones: e.g., improved slab topography models. As mentioned in Sect. 6, the surface of the slab can locally be modified by moving nodes of the mesh. Another improvement could be to allow regional modifications of the thickness of some geophysical zones, such as the continental/oceanic lithospheres and asthenospheres, the slab, or the channel. Currently, with the existing element size, the channel is only one element thick, leading to thickness variability all along the subduction interface (see Sect. S2). A finer element size could be employed to ensure a more uniform thickness across the subduction zone, particularly during the projection process described in Sect. 3. However, this approach cannot be applied to the entire subduction zone, as it would require element sizes of a few kilometers over thousands of kilometers, leading to prohibitively high computation times. Instead, high-resolution local meshes designed for specific portions of the subduction zone could be developed. These refined models could provide detailed representations of slab structures and enhanced coastline descriptions, allowing, for example, a study of coastal long-term vertical motions (e.g. Jolivet et al., 2020; Oryan et al., 2024). Such high-resolution local models could be constrained at their boundaries by applying displacement fields derived from large-scale model simulations (Chiaruttini and Vattré, 2022) as presented in this paper.

## 7.4 Conclusion

We have developed a new finite-element mesh of the Chilean subduction zone, allowing the study of deformation associated with major earthquakes that have ruptured the Chilean subduction interface. This mesh accounts for regional variations in slab geometry (e.g., flat-slabs) and allows the study of their impact on postseismic deformation as well. It gives, in particular, a very valuable and non-biased framework to constrain the rheology of the asthenosphere during the seismic cycle. *Chile_Mesh_v1.0* is well suited for performing seismic inversions, postseismic deformation modeling, and full seismic cycle simulations. Specifically, coseismic slip inversions can be carried out using a consistent set of fault plane nodes along the entire Chilean subduction zone. This enables, for example, the direct reuse of slip distributions without degradation through interpolation onto a different fault grid when investigating the viscoelastic relaxation associated with a given earthquake. Additionally, *Chile_Mesh_v1.0* is compatible with different finite-element software used by the community (e.g. Zset/Zebulon, PyLith, Abaqus). Comparing coseismic and postseismic deformation generated by these different finite-element solvers can help assess their respective performances and improve consistency in subduction zone modeling. Finally, the adaptability of the mesh-building process described in this study is a significant advantage, enabling its use for generating finite-element meshes of other subduction zones with modifications, such as changes in the trench's path, slab topographic model, and vertical plane insertion.

## Acknowledgements

This work was supported by ANR grant number ANR-19-CE31-0003. Maps are made with Generic Mapping Tools GMT (Wessel et al., 2019) and its Python version (PyGMT). We thank the two reviewers and the Editor, Pr. Wenbin Xu, for their insightful and constructive comments.

## Data and code availability

To build this mesh, we used *Zset/Zebulon* (development version, http://zset-software.com/, last access: 15 June 2025) and *MMG* (version 5.7.3, https://www.mmgtools.org/, last access: 15 June 2025) software. *Zset/Zebulon* is a material and structure analysis software, co-developed by the 'Centre des Matériaux' Mines ParisTech (France), ONERA - the French Aerospace Lab (France), and distributed by Transvalor SA. The scripts and the mesh in various formats are available here: https://doi.org/10.5281/zenodo.16927951
(last access: 22 August 2025). All files needed to reproduce the reference case (Sect. 4) are provided in the git.

## Competing interests

The authors declare that they have no competing interests.

## References

Agata, R., Barbot, S., Fujita, K., Hyodo, M., linuma, T., Nakata, R., Ichimura, T., and Hori, T. Rapid mantle flow with power-law creep explains deformation after the 2011 Tohoku mega-quake. *Nature Communications*, 10(1):1-11, 2019. doi: 10.1038/s41467-019-08984-7.

Argus, D. F., Peltier, W. R., Blewitt, G., and Kreemer, C. The viscosity of the top third of the lower mantle estimated using GPS, GRACE, and relative sea level measurements of glacial isostatic adjustment. *Journal of Geophysical Research: Solid Earth*, 126 (5):e2020JB021537, 2021. doi: 10.1029/2020JB021537.

Backus, G. E. Converting Vector and Tensor Equations to Scalar Equations in Spherical Coordinates. *Geophysical Journal International*, 13(1-3):71-101, 07 1967. doi: 10.1111/j.1365-246X.1967.tb02147.x.

Barrientos, S. E. and Ward, S. N. The 1960 Chile earthquake: inversion for slip distribution from surface deformation. *Geophysical Journal International*, 103(3):589-598, 1990. doi: 10.1111/j.1365-246X.1990.tb05673.x.

Bedford, J., Moreno, M., Li, S., Oncken, O., Baez, J. C., Bevis, M., Heidbach, O., and Lange, D. Separating rapid relocking, afterslip, and viscoelastic relaxation: An application of the postseismic straightening method to the Maule 2010 cGPS. *Journal of Geophysical Research: Solid Earth*, 121(10):7618-7638, 2016. doi: 10.1002/2016JB013093.

Boulze, H., Fleitout, L., Klein, E., and Vigny, C. Post-seismic motion after 3 Chilean megathrust earthquakes: a clue for a linear asthenospheric viscosity. *Geophysical Journal International*, 231 (3):1471-1478, 07 2022. doi: 10.1093/gji/ggac255.

Broerse, T., Riva, R., Simons, W., Govers, R., and Vermeersen, B. Postseismic GRACE and GPS observations indicate a rheology contrast above and below the Sumatra slab. *Journal of Geophysical Research: Solid Earth*, 120(7):5343-5361, 2015. doi: 10.1002/2015JB011951.

Cathles, L. *Viscosity of the Earth's Mantle*. Princeton Legacy Library. Princeton University Press, 2015.

Celli, N., Lebedev, S., Schaeffer, A., Ravenna, M., and Gaina, C. The upper mantle beneath the South Atlantic Ocean, South America and Africa from waveform tomography with massive data sets. *Geophysical Journal International*, 221:178-204, 01 2020. doi: 10.1093/gji/ggz574.

Chiaruttini, V. and Vattré, A. Approche monolithique globale/locale à interface diffuse par intersection de maillage. In *CSMA 2022 15ème Colloque National en Calcul des Structures*, 2022. https://hal.science/hal-03687489v1/file/chiaruttini_vattré.pdf.

Cifuentes, I. L. The 1960 Chilean earthquakes. *Journal of Geophysical Research: Solid Earth*, 94(B1):665-680, 1989. doi: 10.1029/JB094iB1p00665.

Contreras-Reyes, E., Flueh, E. R., and Gretemeyer, I. Tectonic control on sediment accretion and subduction off south central Chile: Implications for coseismic rupture processes of the 1960 and 2010 megathrust earthquakes. *Tectonics*, 29(6), 2010. doi: 10.1029/2010TC002734.

Contreras-Reyes, E., Carvajal, M., and González, F. Offshore geometry of the South America subduction zone plate boundary. *Earth and Planetary Science Letters*, 651:119175, 2025. doi: 10.1016/j.epsl.2024.119175.

Dalton, C. A., Ekström, G., and Dziewoński, A. M. The global attenuation structure of the upper mantle. *Journal of Geophysical Research: Solid Earth*, 113(B9), 2008. doi: 10.1029/2007JB005429.

Deuss, A., Redfern, S. A. T., Chambers, K., and Woodhouse, J. H. The Nature of the 660-Kilometer Discontinuity in Earth's Mantle from Global Seismic Observations of <i>PP</i> Precursors. *Science*, 311(5758):198-201, 2006. doi: 10.1126/science.1120020.

Dumoulin, C., Doin, M.-P., and Fleitout, L. Heat transport in stagnant lid convection with temperature-and pressure-dependent Newtonian or non-Newtonian rheology. *Journal of Geophysical Research: Solid Earth*, 104(B6):12759-12777, 1999. doi: 10.1029/1999JB00110.

Freed, A. M. and Lin, J. Delayed triggering of the 1999 Hector Mine earthquake by viscoelastic stress transfer. *Nature*, 411(6834): 180-183, 2001. doi: 10.1038/35075548.

Freed, A. M., Bürgmann, R., Calais, E., Freymueller, J., and Hreinsdóttir, S. Implications of deformation following the 2002 Denali, Alaska, earthquake for postseismic relaxation processes and lithospheric rheology. *Journal of Geophysical Research: Solid Earth*, 111(B1), 2006. doi: 10.1029/2005JB003894.

Freed, A. M., Hashima, A., Becker, T. W., Okaya, D. A., Sato, H., and Hatanaka, Y. Resolving depth-dependent subduction zone viscosity and afterslip from postseismic displacements following the 2011 Tohoku-oki, Japan earthquake. *Earth and Planetary Science Letters*, 459:279-290, 2017. doi: 10.1016/j.epsl.2016.11.040.

Garaud, J.-D., Fleitout, L., and Cailletaud, G. Simulation parallelle de la relaxation post-sismique dans la région de Sumatra. *In Neuvième colloque en calcul des structures*, pages 585-590, 2009.

Glodny, J., Echtler, H., Figueroa, O., Franz, G., Gräfe, K., Kemnitz, H., Kramer, W., Krawczyk, C., Lohrmann, J., Lucassen, F., Melnick, D., Rosenau, M., and Seifert, W. *Long-Term Geological Evolution and Mass-Flow Balance of the South-Central Andes*, pages 401-428. Springer Berlin Heidelberg, Berlin, Heidelberg, 2006. doi: 10.1007/978-3-540-48684-8_19.

Hayes, G., Moore, G., Portner, D., Hearne, M., Flamme, H., Furtney, M., and Smoczyk, G. Slab2, a comprehensive subduction zone geometry model. *Science*, 362:eaat4723, 08 2018. doi: 10.1126/science.aat4723.

Hormazábal, J., Moreno, M., Ortega-Culaciati, F., Báez, J. C., Peña, C., Sippl, C., González-Vidal, D., Ruiz, J., Metzger, S., and Yoshioka, S. Fast relocking and afterslip-seismicity evolution following the 2015 Mw 8.3 Illapel earthquake in Chile. *Scientific Reports*, 13(1):19511, 2023. doi: 10.1038/s41598-023-45369-9.

Hu, Y. and Wang, K. Spherical-Earth finite element model of short-term postseismic deformation following the 2004 Sumatra earthquake. *Journal of Geophysical Research: Solid Earth*, 117(B5), 2012. doi: 10.1029/2012JB009153.

Hu, Y., Wang, K., He, J., Klotz, J., and Khazaradze, G. Three-dimensional viscoelastic finite element model for postseismic deformation of the great 1960 Chile earthquake. *Journal of Geophysical Research: Solid Earth*, 109(B12), 2004. doi: 10.1029/2004JB003163.

Hu, Y., Bürgmann, R., Freymueller, J. T., Banerjee, P., and Wang, K. Contributions of poroelastic rebound and a weak volcanic arc to the postseismic deformation of the 2011 Tohoku earthquake. *Earth, Planets and Space*, 66:1-10, 2014. doi: 10.1186/1880-5981-66-10.

Jolivet, R., Simons, M., Duputel, Z., Olive, J.-A., Bhat, H. S.,
and Bletery, Q. Interseismic Loading of Subduction Megathrust Drives Long-Term Uplift in Northern Chile. *Geophysical Research Letters*, 47(8):e2019GL085377, 2020. doi: 10.1029/2019GL085377.

Khazaradze, G., Wang, K., Klotz, J., Hu, Y., and He, J. Prolonged post-seismic deformation of the 1960 great Chile earthquake and implications for mantle rheology. *Geophysical Research Letters*, 29(22):7–1–7–4, 2002. doi: 10.1029/2002GL015986.

Klein, E., Fleitout, L., Vigny, C., and Garaud, J. Afterslip and viscoelastic relaxation model inferred from the large-scale post-seismic deformation following the 2010 Mw 8.8 Maule earthquake (Chile). *Geophysical Journal International*, 205(3): 1455–1472, 03 2016. doi: 10.1093/gji/ggw086.

Klein, E., Vigny, C., Fleitout, L., Grandin, R., Jolivet, R., Rivera, E., and Métois, M. A comprehensive analysis of the Illapel 2015 Mw8.3 earthquake from GPS and INSAR data. *Earth and Planetary Science Letters*, 469:123–134, July 2017. doi: 10.1016/j.epsl.2017.04.010.

Klein, E., Vigny, C., Nocquet, J.-M., and Boulze, H. A 20 year-long GNSS solution across South-America with focus in Chile. *Bulletin de la Société géologique de France, in press*, 2022. doi: 10.1051/bsgf/2022005.

Klotz, J., Khazaradze, G., Angermann, D., Reigber, C., Perdomo, R., and Cifuentes, O. Earthquake cycle dominates contemporary crustal deformation in Central and Southern Andes. *Earth and Planetary Science Letters*, 193(3-4):437–446, 2001. doi: 10.1016/S0012-821X(01)00532-5.

Li, S., Moreno, M., Bedford, J., Rosenau, M., and Oncken, O. Revisiting viscoelastic effects on interseismic deformation and locking degree: A case study of the Peru-North Chile subduction zone. *Journal of Geophysical Research: Solid Earth*, 120(6):4522–4538, 2015. doi: 10.1002/2015JB011903.

Li, S., Moreno, M., Bedford, J., Rosenau, M., Heidbach, O., Melnick, D., and Oncken, O. Postseismic uplift of the Andes following the 2010 Maule earthquake: Implications for mantle rheology. *Geophysical Research Letters*, 44(4):1768–1776, 2017. doi: 10.1002/2016GL071995.

Liu, T., Fu, G., She, Y., Meng, G., Zou, Z., Wu, W., Shestakov, N. V., Gerasimenko, M. D., Bykov, V. G., and Pupatenko, V. V. Postseismic deformation following the 2011 Mw9.0 Tohoku–Oki earthquake and its impact on Northeast Asia. *Geophysical Journal International*, 235(2):1479–1492, 08 2023. doi: 10.1093/gji/ggad314.

Lovery, B., Radiguet, M., Chlieh, M., Norabuena, E., Villegas-Lanza, J. C., Cresseaux, J., Ragon, T., Tsapong-Tsague, A., Tavera, H., and Socquet, A. Viscoelastic Relaxation Following the 2001 Mw 8.4 Arequipa Earthquake and Its Impact on the Interseismic Coupling of the South Peru Megathrust. *Geophysical Research Letters*, 52(12):e2024GL113879, 2025. doi: 10.1029/2024GL113879.

Maksymowicz, A., Ruiz, J., Vera, E., Contreras-Reyes, E., Ruiz, S., Arraigada, C., Bonvalot, S., and Bascuñan, S. Heterogeneous structure of the Northern Chile marine forearc and its implications for megathrust earthquakes. *Geophysical Journal International*, 215(2):1080–1097, 08 2018. doi: 10.1093/gji/ggy325.

Marsman, C. P., Vossepoel, F. C., D'Acquisto, M., van Dinther, Y., van de Wiel, L., and Govers, R. Unraveling Processes and Rheology of the Tohoku Earthquake Cycle Using Bayesian Inference. *Journal of Geophysical Research: Solid Earth*, 130(5): e2024JB029665, 2025. doi: 10.1029/2024JB029665.

Melnick, D., Moreno, M., Quinteros, J., Baez, J. C., Deng, Z., Li, S., and Oncken, O. The super-interseismic phase of the megathrust earthquake cycle in Chile. *Geophysical Research Letters*, 44(2): 784–791, 2017. doi: 10.1002/2016GL071845.

Melosh, H. J. and Raefsky, A. A simple and efficient method for introducing faults into finite element computations. *Bulletin of the Seismological Society of America*, 71(5):1391–1400, 10 1981. doi: 10.1785/BSSA0710051391.

Moreno, M., Melnick, D., Rosenau, M., Bolte, J., Klotz, J., Echtler, H., Baez, J., Bataille, K., Chen, J., Bevis, M., Hase, H., and Oncken, O. Heterogeneous plate locking in the South–Central Chile subduction zone: Building up the next great earthquake. *Earth and Planetary Science Letters*, 305(3):413–424, 2011. doi: 10.1016/j.epsl.2011.03.025.

Moreno, M., Melnick, D., Rosenau, M., Baez, J., Klotz, J., Oncken, O., Tassara, A., Chen, J., Bataille, K., Bevis, M., Socquet, A., Bolte, J., Vigny, C., Brooks, B., Ryder, I., Grund, V., Smalley, B., Carrizo, D., Bartsch, M., and Hase, H. Toward understanding tectonic control on the Mw 8.8 2010 Maule Chile earthquake. *Earth and Planetary Science Letters*, 321-322:152–165, 2012. doi: 10.1016/j.epsl.2012.01.006.

Nield, G. A., King, M. A., Steffen, R., and Blank, B. A global, spherical finite-element model for post-seismic deformation using *Abaqus*. *Geoscientific Model Development*, 15(6):2489–2503, 2022. doi: 10.5194/gmd-15-2489-2022.

Oryan, B., Olive, J.-A., Jolivet, R., Malatesta, L. C., Gailleton, B., and Bruhat, L. Megathrust locking encoded in subduction landscapes. *Science Advances*, 10(17):eadl4286, 2024. doi: 10.1126/sciadv.adl4286.

Peña, C., Heidbach, O., Moreno, M., Bedford, J., Ziegler, M., Tassara, A., and Oncken, O. Role of Lower Crust in the Postseismic Deformation of the 2010 Maule Earthquake: Insights from a Model with Power-Law Rheology. *Pure and Applied Geophysics*, 176(9):3913–3928, Sept. 2019. doi: 10.1007/s00024-018-02090-3.

Peña, C., Heidbach, O., Moreno, M., Bedford, J., Ziegler, M., Tassara, A., and Oncken, O. Impact of power-law rheology on the viscoelastic relaxation pattern and afterslip distribution following the 2010 Mw 8.8 Maule earthquake. *Earth and Planetary Science Letters*, 542:116292, 2020. doi: 10.1016/j.epsl.2020.116292.

Pollitz, F., Banerjee, P., Grijalva, K., Nagarajan, B., and Bürgmann, R. Effect of 3-D viscoelastic structure on post-seismic relaxation from the 2004 M= 9.2 Sumatra earthquake. *Geophysical Journal International*, 173(1):189–204, 04 2008. doi: 10.1111/j.1365-246X.2007.03666.x.

Pollitz, F. F. Gravitational viscoelastic postseismic relaxation on a layered spherical Earth. *Journal of Geophysical Research: Solid Earth*, 102(B8):17921–17941, 1997. doi: 10.1029/97JB01277.

Pollitz, F. F. Post-seismic relaxation theory on a laterally heterogeneous viscoelastic model. *Geophysical Journal International*, 155(1):57–78, 10 2003. doi: 10.1046/j.1365-246X.2003.01980.x.

Pollitz, F. F., Bürgmann, R., and Banerjee, P. Post-seismic relaxation following the great 2004 Sumatra-Andaman earthquake on a compressible self-gravitating Earth. *Geophysical Journal International*, 167(1):397–420, 10 2006. doi: 10.1111/j.1365-246X.2006.03018.x.

Resovsky, J., Trampert, J., and Van der Hilst, R. Error bars for the global seismic Q profile. *Earth and Planetary Science Letters*, 230 (3):413–423, 2005. doi: 10.1016/j.epsl.2004.12.008.

Ringwood, A. and Irifune, T. Nature of the 650–km seismic discontinuity: implications for mantle dynamics and differentiation. *Nature*, 331(6152):131–136, 1988. doi: 10.1038/331131a0.

Savage, J. C. A dislocation model of strain accumulation and release at a subduction zone. *Journal of Geophysical Research: Solid Earth*, 88(B6):4984–4996, 1983. doi: 10.1029/JB088iB06p04984.

Sodoudi, F., Yuan, X., Asch, G., and Kind, R. High-resolution im-
age of the geometry and thickness of the subducting Nazca
lithosphere beneath northern Chile. *Journal of Geophysical Re-
search: Solid Earth*, 116(B4), 2011. doi: 10.1029/2010JB007829.

Suito, H. Importance of rheological heterogeneity for interpreting viscoelastic relaxation caused by the 2011 Tohoku-Oki earthquake. *Earth, Planets and Space*, 69(1):1-12, 2017. doi: 10.1186/s40623-017-0611-9.

Suito, H. and Freymueller, J. T. A viscoelastic and afterslip postseis-
mic deformation model for the 1964 Alaska earthquake. *Journal*
*nal of Geophysical Research: Solid Earth*, 114(B11), 2009. doi:
*10.1029/2008JB005954*.

Sun, T., Wang, K., linuma, T., Hino, R., He, J., Fujimoto, H., Kido,
M., Osada, Y., Miura, S., Ohta, Y., et al. Prevalence of viscoelastic
relaxation after the 2011 Tohoku-oki earthquake. *Nature*, 514
(7520):84–87, 2014. doi: 10.1038/nature13778.

Sun, T., Wang, K., and He, J. Crustal Deformation Following Great Subduction Earthquakes Controlled by Earthquake Size and Mantle Rheology. *Journal of Geophysical Research: Solid Earth*, 123(6):5323–5345, 2018. doi: 10.1029/2017JB015242.

Trubienko, O., Fleitout, L., Garaud, J.-D., and Vigny, C. Interpretation of interseismic deformations and the seismic cycle associated with large subduction earthquakes. *Tectonophysics*, 589: 126-141, Mar. 2013. doi: 10.1016/j.tecto.2012.12.027.

Trubienko, O., Garaud, J.-D., and Fleitout, L. Models of postseismic deformation after megaearthquakes: the role of various rheological and geometrical parameters of the subduction zone. *Solid Earth Discussions*, 6:427-466, 2014. doi: 10.5194/sed-6-427-2014.

Vigny, C., Rudloff, A., Ruegg, J.-C., Madariaga, R., Campos, J., and Alvarez, M. Upper plate deformation measured by GPS in the Coquimbo Gap, Chile. *Physics of the Earth and Planetary Interiors*, 175(1):86-95, 2009. doi: 10.1016/j.pepi.2008.02.013.

von Huene, R., Corvalán, J., Flueh, E. R., Hinz, K., Korstgard, J., Ranero, C. R., and Weinrebe, W. Tectonic control of the subducting Juan Fernández Ridge on the Andean margin near Valparaiso, Chile. *Tectonics*, 16(3):474-488, 1997. doi: 10.1029/96TC03703.

Wang, K., Hu, Y., Bevis, M., Kendrick, E., Smalley Jr., R., Vargas, R. B., and Lauría, E. Crustal motion in the zone of the 1960 Chile earthquake: Detangling earthquake-cycle deformation and forearc-sliver translation. *Geochemistry, Geophysics, Geosystems*, 8(10), 2007. doi: 10.1029/2007GC001721.

Wang, K., Hu, Y., and He, J. Deformation cycles of subduction earthquakes in a viscoelastic Earth. *Nature*, 484(7394):327-332, Apr. 2012. doi: 10.1038/nature11032.

Wessel, P., Luis, J. F., Uieda, L., Scharroo, R., Wobbe, F., Smith, W. H. F., and Tian, D. The Generic Mapping Tools Version 6. *Geochemistry, Geophysics, Geosystems*, 20(11):5556–5564, 2019. doi: 10.1029/2019GC008515.

Wu, P. and Peltier, W. R. Viscous gravitational relaxation. *Geo-
physical Journal International*, 70(2):435–485, 08 1982. doi:
10.1111/j.1365-246X.1982.tb04976.x.

The article A 3D finite-element mesh for modeling large-
scale surface deformation induced by subduction megath-
rust earthquakes: Application to Chile © 2025 by Hugo
Boulze is licensed under CC BY 4.0.