Transient Deformation and Stress Patterns Induced by the 2010 Maule Earthquake in the Illapel Segment

Evaluating the transfer of stresses from megathrust earthquakes to adjacent segments is fundamental to assess seismic hazard. Here, we use a 3D forward model as well as GPS and seismic data to investigate the transient deformation and Coulomb Failure Stresses (CFS) changes induced by the 2010 Maule earthquake in its northern segment, where the Mw 8.3 Illapel earthquake occurred in 2015. The 3D model incorporates the coseismically instantaneous, elastic response, and time-dependent afterslip and viscoelastic relaxation processes in the postseismic period. We particularly examine the impact of linear and power-law rheology on the resulting postseismic deformation and CFS changes that may have triggered the Illapel earthquake. At the Illapel hypocenter, our model results in CFS changes of ∼0.06 bar due to the coseismic and postseismic deformation, where the coseismic deformation accounts for ∼85% of the total CFS changes. This is below the assumed triggering threshold of 0.1 bar and, compared to the annual loading rate of the plate interface, represents a clock advance of approximately only 2 months. However, we find that sixteen events with Mw ≥ 5 in the southern region occurred in regions of CFS changes > 0.1 bar, indicating a potential triggering by the Maule event. Interestingly, while the power-law rheology model increases the positive coseismic CFS changes, the linear rheology reduces them. This is due to the opposite polarity of the postseismic displacements resulting from the rheology model choice. The power-law rheology model generates surface displacements that fit better to the GPS-observed landward displacement pattern.


INTRODUCTION
Megathrust earthquakes induce local and large-scale deformations (e.g., Hu et al., 2004;Vigny et al., 2011;Moreno et al., 2012), which may last several decades due to viscoelastic relaxation (Khazaradze et al., 2002;Hu et al., 2004;Wang et al., 2012). The crust continues to move after the coseismic instantaneous response, exhibiting different deformation patterns that are observed by GPS measurements (Wang et al., 2012;Tomita et al., 2015;Klein et al., 2016). Inland GPS data show that the crust moves in the sense of the main shock in front of the rupture area, which in most subduction zones represents trenchward movement. At neighboring segments of the rupture area, however, the GPS horizontal displacements exhibit landward motion (Heki and Mitsui, 2013;Tomita et al., 2015;Klein et al., 2016;Loveless, 2017;Melnick et al., 2017), as envisioned by Anderson (1975). While the trenchward motion patterns are mostly attributed to afterslip on the fault interface and viscous relaxation in the lower crust and upper mantle (e.g., Hu et al., 2004;Hergert and Heidbach, 2006;Wang et al., 2012;Peña et al., 2019Peña et al., , 2020, the driving mechanisms controlling landward patterns following megathrust earthquakes are currently a matter of controversial debate. For instance,  and Tomita et al. (2015) showed that landward acceleration of the oceanic plate following the 2011 Tohoku-oki earthquake in Japan may be attributed to upper mantle viscoelastic relaxation. Alternatively, other authors have proposed that the slab pull balance of forces after large earthquakes is the driving mechanism of the GPSobserved postseismic landward patterns (Heki and Mitsui, 2013;Yuzariyadi and Heki, 2020).
The deformation patterns induced by the coseismic and postseismic deformation are of vital importance, as they can potentially modify the state of stress at adjacent segments, which may be mature segments critically stressed, and ultimately triggering large-magnitude earthquakes. Examples of earthquake triggering because of coseismic stress release and its subsequent postseismic deformation can be found in the Sumatran-Andaman region, Indonesia (McCloskey et al., 2005;Pollitz et al., 2006;Hughes et al., 2010;Qiu and Chan, 2019) and earthquakes in strike-slip fault systems such as the 1999 Hector Mine earthquake in Southern California, United States (Freed and Lin, 2001) and the 2016 Visso earthquake in Italy (Tung and Masterlark, 2018).
On the 27th of February 2010, the Maule earthquake occurred in Central Chile with a magnitude of Mw 8.8, breaking a 500 km long plate boundary segment along strike (Figure 1 and Vigny et al., 2011;Moreno et al., 2012). About 5 years later, a segment approximately 200 km north of the Maule rupture area ruptured with the Mw 8.3 Illapel earthquake that occurred on September 16th, 2015 (Figure 1). This earthquake broke a presumably mature segment, where the last event occurred in 1943 with a magnitude of Mw 7.9 (Beck et al., 1998). In the Illapel segment (enclosed area in Figure 2), an enhancement of landward displacement was observed following the 2010 Maule event (Figure 2 and Klein et al., 2016;Melnick et al., 2017;Ruiz et al., 2016). Using GPS data, Ruiz et al. (2016) revealed an increase of the annual interseismic velocity in this segment of approximately 15% between the 2010 event and the end of 2014, but they did not measure the locking degree directly or the stress changes due to this increase. Melnick et al. (2017) found an increase of 20% in the apparent locking degree by using GPS observations in the period within the 2010 Maule and the 2015 Illapel events. In the study of Melnick et al. (2017), however, the effect of exponential transient components, e.g., viscous relaxation, were not included, as they only focused on the change in the linear trend component of the observations between the 2010 and 2015 events to calculate the locking changes. By contrast, Klein et al. (2016) showed that a possible mechanism to explain such a pattern may be viscoelastic relaxation in the asthenosphere, as well as in a deep, low-viscosity channel along the fault interface, using a linear Burgers rheology model. Despite these advances in understanding the kinematics of these patterns, it remains unclear how the 3D power-law rheology affects the landward postseismic displacements and how stresses were built-up on the fault interface in the Illapel segment due to the coseismic and postseismic deformation of the 2010 Maule event.
Here, we use the 3D geomechanical-numerical model of Peña et al. (2020) and cumulative GPS observations between the 2010 Maule and 2015 Illapel earthquakes to examine the transient deformation and stress changes in the Illapel segment induced by the coseismic and postseismic deformation of the 2010 Maule earthquake. We first investigate the differences in the resulting surface postseismic displacements due to the choice of model rheology and how this affects the calculation of the Coulomb Failure Stress (CFS) changes. Finally, we evaluate if the CFS changes induced by the co-and postseismic effects of the 2010 Maule earthquake could have triggered the 2015 Illapel event.

MODEL SETUP AND GPS OBSERVATIONS
We use the 3D geomechanical-numerical model of Peña et al. (2020) to investigate both the deformation and stress transient patterns induced in the Illapel segment due to the 2010 Maule event. The model described in Peña et al. (2020) considers the slab and Moho geometries obtained from Hayes et al. (2012) and Tassara et al. (2006), respectively. The model extends 4000 km in the WE direction, 2000 km in the NS direction, and 400 km in the vertical direction (see Peña et al. (2020) for further details). The postseismic deformation following the 2010 Maule event is simulated by combining afterslip and viscous relaxation jointly through a forward modeling approach. To do so, we employ the best-fit simulation of Peña et al. (2020), which considers afterslip on the fault interface, and power-law rheology with dislocation creep processes in the crust and upper mantle (see elastic and creep parameters in Supplementary Tables 1, 2). In this simulation, most of the afterslip occurs at depths <60 km (Figure 2A), while viscous relaxation occurs in the continental lower crust and upper mantle beneath the volcanic arc (see Figure 8A in Peña et al., 2020). As discussed by Peña et al. (2020), this model was favored not only because of the best fit to cumulative GPS observations and GPS time-series but also because of the correlation between the modeled afterslip and the aftershocks. For comparison, we also consider a simulation with linear rheology, which consists of a whole elastic crust and linear viscosity in the upper mantle. The afterslip distribution from the linear rheology case is displayed in Figure 2B.
In this study, we compute the CFS changes on the megathrust fault interface in the Illapel segment due to the coseismic and postseismic deformation of the 2010 main shock. We calculate the CFS changes employing the linear relationship between shear and normal stresses given by, where τ is the change in shear stresses, µ the effective coefficient of friction and σ the change in normal stresses (e.g., King et al., 1994;Stein, 1999). Following the approach of Moreno et al. (2018), we also consider a heterogeneous effective coefficient of friction on the fault interface. The study of Moreno et al. (2018) showed that this approach could explain the occurrence of moderate-size earthquakes in southern Chile better. Moreno et al. (2018) estimated µ = 0.04 for regions where the locking degree is ≥ 0.8 and µ = 0.01 for regions with a locking degree < 0.8. We employ a similar approach, testing for three different locking degree threshold values (0.8, 0.5, and 0.2), where µ is set to 0.04 and 0.01 above and below the locking threshold, respectively. These low values of µ are supported by a recent study by Dielforder et al. (2020) and by the inversion of stress orientations of moment tensors from subduction earthquakes as well. This suggests a frictional zonation with a higher effective coefficient of friction in highly locked domains as compared to creeping ones (Hardebeck and Loveless, 2018). We tested end-member µ values as well (Supplementary Figure 1), but they do not result in any major CFS changes. This agrees with Pollitz et al. (2006), who found similar stress changes by considering a range , which represents 5.5 years. Therefore, the resulting modeled surface displacements displayed in the Illapel segments account for this period, although the differences considering 6 years are negligible. Thus, green (trenchward) and yellow arrows (landward) represent the horizontal cumulative 5.5-years surface displacements from the GPS data used in this study, while the red arrows the ones from the model in the same period.
of µ between 0 and 0.8 to investigate the stress changes in the Sumatra-Andaman region. The locking coefficient is taken from Tilmann et al. (2016) who use an inversion of the interseismic GPS observations using an inversion method based on FEMgenerated Green's functions. Their implemented fault geometry and elastic parameters follows the approach of Moreno et al. (2012), which are also consistent with those used in the 3D model of Peña et al. (2020). The resulting numerical problem is resolved with the commercial finite element code Abaqus, version 6.14 and the CFS changes are computed on the fault interface using GeoStress Add-on (Stromeyer et al., 2020) for Tecplot 360 EX, version 2019 R1.
To investigate the differences in deformation of the Illapel segment due to the choices of model rheology, we use the cumulative postseismic GPS displacements from Melnick et al. (2017) between the 2010 Maule and 2015 Illapel earthquakes. In this data set, the effect of aftershocks and seasonal components have been removed (further information can be found in supporting information in Melnick et al., 2017). Following the approach from previous studies (e.g., Klein et al., 2016;Weiss et al., 2019;Peña et al., 2020), we have also removed the interseismic trend measured prior to the 2010 Maule earthquake to thus account only for the postseismic effects following the 2010 Maule event (Figure 2).

RESULTS
The surface displacements are presented in Figure 2 and the coseismic and postseismic CFS changes are presented in Figures 3, 4. Finally, we compare the cumulative CFS changes   Figure 2 shows the cumulative GPS displacements between the 2010 Maule and 2015 Illapel earthquakes versus the modeled ones. In the Maule rupture zone, the predicted horizontal displacements from the modeled afterslip and viscoelastic relaxation do not exhibit large differences between the linear or power-law model rheology. However, this is not the case for the patterns in the Illapel segment. There, the observations display landward displacement patterns with a strong north component and a maximum cumulative displacement of 3.9 cm after ∼5.5 years. By comparing the GPS data with model predictions, we find that the model with power-law rheology can better explain those observed postseismic landward patterns than the linear rheology case. The power-law rheology simulation results in landward motion with a mostly eastwards azimuth and with a similar average magnitude compared to the observations (dashed-line rectangle in Figure 2A). This simulation also shows a maximum displacement at approximately 30.5 • South, decreasing to the north and to the south from this latitude. This is also in agreement with the GPS data. Conversely, the linear rheology case results in trenchward motion, with a dominant northwest direction and slightly larger magnitude than the observed horizontal displacement (Figure 2B). Given that the afterslip distributions from both rheology model cases show similar magnitude and patterns in the closest region to the Illapel segment (Figure 2), the main deviations in their resulting horizontal displacement patterns are due to differences in location of viscous relaxation from the rheology model choice.

Surface Displacements Patterns
In the power-law rheology case, most of the viscous relaxation occurs in the continental lower crust, beneath the volcanic arc. Conversely, in the linear rheology case in the continental mantle wedge, beneath the forearc, as well as in the oceanic mantle (see Figure 8 in Peña et al., 2020).

CFS Changes
We also estimate the CFS changes in the Illapel segment to evaluate to what extent the Maule earthquake may have triggered the seismicity in the Illapel segment and subsequently the 2015 Illapel main shock (Figures 3, 4). The CFS changes are calculated on the megathrust fault interface for regions shallower than 60 km depth (enclosed pink area in Figure 1). Additionally, since we are interested in evaluating the CFS changes on the fault interface, we use the focal mechanism parameter of the 2015 Illapel earthquake from the USGS National Earthquake Information Center (NEIC), which reported a strike of 353 • , a dip of 19 • , and a rake of 83 • . We also use the same receiver fault parameters for inferring the CFS changes at the hypocenters of the seismicity between the 2010 Maule and 2015 Illapel events, since most of the events share the same focal mechanisms as the Illapel main shock, especially for events > M w 5.7 (Supplementary Figure 2). We also investigated the potential role of variations of the receiver fault parameters by exploring plausible variations of the focal mechanisms of the seismic events, but we found only minor effects (Supplementary Figures 3-7). Thus, we first present the impact due to the choice of model rheology for the calculation of the CFS changes at the hypocenter of the 2015 Illapel main shock. Figure 4 shows the CFS changes from the preferred model alongside the seismicity recorded between the Maule and Illapel events. This seismicity catalog is also acquired from the NEIC. Figure 3 shows the resulting CFS changes induced by coseismic and postseismic deformation of the Maule earthquake in the Illapel segment. As described in the model setup section, we consider three cases with friction zonation on the fault interface by assigning a higher effective coefficient of friction (µ = 0.04) to regions where the locking degree is ≥0.8, 0.5, and 0.2 than its surroundings (µ = 0.01). These we hereafter label as scenario I, II, and III, respectively. Although the choice of the friction zonation has an impact on the resulting CFS changes (Figure 3), it has a negligible effect on the resulting displacement patterns. The largest CFS changes are found in regions closest to the 2010 Maule rupture region, with a maximum magnitude close to 1 bar (0.1 MPa) that decreases northwards to values smaller than 0.001 bar. For the postseismic period, the model with power-law rheology results in positive CFS changes, while the one with linear rheology leads to negative values. We find that the coseismic stress transfer dominates the CFS changes in all three scenarios. It accounts for approximately 85% of the total CFS changes at the hypocenter of Illapel main shock (Figures 3A,D,G). At the Illapel main shock hypocenter, the coseismic CFS changes are 0.052, 0.051, and 0.006 bar for scenarios I, II, and III, respectively (Figures 3A,D,G). For the postseismic period using the powerlaw rheology model, the postseismic CFS changes are 0.008, 0.008, and 0.001 bar for scenarios I, II, and III, respectively (Figures 3B,E,H). On the contrary, negative CFS changes are found from the linear rheology model during the postseismic period. Here, we find CFS changes of −0.004, −0.011, and −0.003 bar for scenarios I, II, and III, respectively (Figures 3C,F,I).
Since the model with power-law rheology can better explain the landward GPS displacements, we continue our analysis using the CFS results from this model and compare to seismic events with M w ≥ 5.0 between the 2010 Maule and Illapel events. The seismicity is mostly located in the surroundings of highly locked regions (≥0.8) as shown in Figure 4A. These events also fall well inside regions of higher CFS changes compared to scenarios II and III (Figures 4A-C). In particular, sixteen events in the southern part of the segment, as shown in region R1 (black dashed rectangle in Figure 4A), occurred in regions of CFS changes larger than 0.1 bar (Figures 4A,D). On the other hand, in the same region R1, we find five and four events for scenarios II and III, respectively, with CFS changes ≥ 0.1 bar (Figures 4B,C,E,F). Additionally, in all scenarios, the CFS changes at the hypocenter of the largest event with magnitude Mw 6.7 (April 17, 2012) is higher than 0.1 bar (Figures 4D,E,F). Supplementary Figure 1 show that soon after the 2010 Maule earthquake, two events with magnitude Mw 5.2 occurred 30 min (E2, ∼31.7 • S) and 6 days after the Maule event (E1, ∼29.9 • S), while a seismicity cluster C1 (∼32.6 • S) composed of 6 events with magnitudes Mw 5.0-5.3 occurred 9 days after the Maule main shock. Our results indicate that the coseismically positive CFS changes calculated at E1, E2, and C1 prevail in scenario I with 0.003 bar, 0.06 bar, and 0.13 bar, respectively; however, the CFS for E1 in scenario II is slightly larger than in scenario I (with 0.004 bar).

DISCUSSION
Landward deformation patterns following megathrust earthquakes have been observed by modern GPS data (e.g., Heki and Mitsui, 2013;Melnick et al., 2017), but their controlling mechanism is not fully understood. Our results show that the choice of rheology has a large effect on the azimuth of deformation at segments adjacent to the mainshock. Here, an upper mantle with homogenous and linear viscoelastic material properties cannot predict the observed landward pattern. In contrast, a model considering second-order, power-law rheology in the crust and upper mantle can better explain this pattern (Figure 2A). An alternative approach to explain this pattern is the one of Klein et al. (2016), who considered an upper mantle with Burgers rheology, as well as a deep subduction channel with low viscosities (∼10 17 Pa s). On the other hand, Yuzariyadi and Heki (2020) showed that a generic 3D model with linear Maxwell rheology in the upper mantle exhibited trenchward patterns. Because of this discrepancy with the GPS data, Yuzariyadi and Heki (2020) hypothesized that deep slab acceleration may also play a key role, in agreement with the study of Heki and Mitsui (2013). Nevertheless, our findings indicate that the surface displacement patterns from the 3D forward model of Yuzariyadi and Heki (2020) are produced by the simplicity of the rheological structure of their model, which considers a linear viscoelastic mantle, whose viscosity is homogeneously distributed in the whole upper mantle, and a fully elastic crust.
Even though the fit of the postseismic displacement with the forward model that uses non-linear rheology is much better in comparison with linear viscoelastic implementation, it cannot perfectly predict the north component of the GPS signal. However, this has hardly any impact on the estimates of the CFS changes as the stress traction vector is resolved in the rake vector direction of 83 • from the Illapel earthquake, which means that there is marginally any contribution from the north component of the displacement to the stress changes.
The CFS changes produced by the coseismic and postseismic deformation of the Maule event are relatively small. Nevertheless, these small values are sufficient to trigger seismicity, as statically shown in other subduction and strike-system settings with a threshold magnitude as small as 0.1 bar (e.g., Stein, 1999, and references therein;Ziv and Rubin, 2000). Moreover, laboratoryderived, rate-and-state friction laws predict that small CFS changes can promote the increase of seismicity (e.g., Stein, 1999, and references therein). Notably, in the southern region R1 (Figure 4A), five events Mw ≥ 5.5 occurred in regions with a change of CFS > 0.1 bar ( Figure 4D). Our findings reveal that at the hypocenter of the two larger events (Mw 6.4 and 6.7) the CFS increased by ∼0.63 and 0.14 bar, respectively. Furthermore, our results indicate a potential triggering of the events in region R1 due to the joint contribution of coseismic and postseismic deformation ( Figure 4A).
Yet, no strong increase in seismicity rate changes has been reported after the 2010 Maule event in the whole Illapel segment (Supplementary Figure 7 and Heidarzadeh et al., 2016;Poli et al., 2017;Ruiz et al., 2016), but so far only events ≥ Mw 4.0 have been considered. Huang and Meng (2018) used data provided by the Centro Sismológico Nacional de Chile (CSN) to extract repeating earthquakes (Mw ≥ 2.5), employing a matched-filterenhanced method. They considered ∼4 years before the Illapel earthquake and found an increase of seismic activity in the Illapel segment 100 days before the Illapel main shock but it remains unclear how seismicity evolved in the first ∼2 years following the 2010 Maule event, a period in which the postseismic deformation associated to the Maule event prevails (e.g., Bedford et al., 2016;Klein et al., 2016;Peña et al., 2020). As shown by Toda et al. (2011) for the case of the enhanced seismicity rates following the 2011 Tohoku-Oki earthquake in Japan, an enhanced seismic catalog considering M ≤ 4.0 events may provide key information to ultimately evaluate the triggering of the events in the southern region under the CFS criteria.
On the other hand, at the hypocenter of the Illapel main shock and events further north, the calculated CFS changes (<0.06 bar) are below the assumed threshold CFS of 0.1 bar. Furthermore, if we use the shear stress drop of 30-60 bar of the Illapel earthquake estimated by Yin et al. (2016) with a 3D kinematic model and assuming that this shear stress has been accumulated since the previous Mw 7.9 event on the same segment in 1944 (Beck et al., 1998), the annual shear stress loading rate would be between 0.42-0.84 bar/yr. Changes in the normal stress due to overburden can be assumed to be very small. Consequently, the CFS changes of ∼0.06 bar that is transferred co-and postseismically to the Illapel hypocenter represent only 7-14% of the annual loading rate. Thus, given that the CFS changes are below the threshold of 0.1 bar and the small clock-advance of approximately a few months only, we have two options how to interpret the results. Either there is no triggering and earthquake occurrence is a Poisson process in subduction zones, or another mechanism is acting that is not represented in the model used for the stress transfer calculations.
Additional processes transferring CFS more efficiently may be related to dynamic triggering due to the passage of seismic waves (e.g., Gomberg et al., 2001;Freed, 2005;van der Elst and Brodsky, 2010), or changes in pore pressure due to changes in the permeability in the fault interface (e.g., Brodsky et al., 2003;Elkhoury et al., 2006;Manga et al., 2012). The first one may change the state of the stresses of the medium due to the passage of the seismic waves, both at near and remote distances, but is a very short-term mechanism (e.g., Freed, 2005;van der Elst and Brodsky, 2010). Given the small coseismically CFS values inferred at the hypocenter of very short-term events E1 (0.003 bar) and E2 (0.07 bar, Figure 4A), 9 days and 30 min after the Maule event, respectively (Supplementary Figure 1), dynamic wave processes are likely candidates for the occurrence of these events. On the other hand, pore-fluid processes last longer and produce surface deformation patterns confined close to the rupture zone, i.e., mostly offshore at subduction zones, rather than at further distance (e.g., Masterlark, 2003;Hu et al., 2014;McCormack et al., 2020). Particularly, pore pressure changes would significantly affect the effective normal stresses (e.g., Brodsky et al., 2003;Elkhoury et al., 2006;Hughes et al., 2010;Manga et al., 2012), and would also explain the delay in response due to the slow diffusion of the fluids. In addition, its resulting surface displacement field is generally smaller than that from viscous relaxation, although the CFS changes are larger than those from viscous relaxation only (Masterlark et al., 2001;Hughes et al., 2010). Therefore, our CFS calculations in the Illapel segment during the postseismic period may be underestimated by not considering pore-pressure processes. Its integration into the 3D model would require further constraints from geodetic observations, potentially from InSAR data using its high spatial resolution.

CONCLUSION
We examine the co-and postseismic deformation and stress change pattern on the Illapel segment induced by the 2010 Maule earthquake using a 3D geomechanical-numerical model with linear and power-law rheology. Our results show that the postseismic horizontal deformation patterns in the Illapel segment are particularly sensitive to the choice of the rheology model. Here, the power-law rheology model can better explain both the landward directions and magnitude of the cumulative GPS displacement vectors. The opposite surface displacement patterns from the choice of rheology also have a major impact on the resulting postseismic CFS changes, but compared to the co-seismic stress transfer, only 15% of the CFS changes are transferred postseismically to the hypocenter location of the 2015 Illapel earthquake. From magnitudes of CFS changes it seems quite unlikely that the Illapel earthquake was triggered by the Maule event -such a small CFS change of <0.06 bar at the hypocenter converts to a clock-advance of approximately 2 months. This finding would not change with a better fit of the modeled north component of the GPS observation. Hence, either no triggering occurred or a local secondary process such a pore-pressure changes, which is not considered here, potentially contributed to a remote triggering mechanism.
We also compared the cumulative CFS changes to the seismicity that occurred between the Maule and Illapel events, using the model with the power-law rheology. Here, particularly in the southern region of the Illapel segment, sixteen seismic events occurred in regions of positive CFS changes above the triggering threshold of 0.1 bar, indicating a potential triggering because of the Maule event. Therefore, our study emphasizes the importance of using power-law rheology models to also investigate surface landward patterns and stress changes induced by megathrust earthquakes in neighboring segments. Future investigation should especially focus on assessing the transfer of stresses from the 2014 Iquique in its southern segment (Loa segment) and from the 2015 Illapel earthquake in its northern segment (Atacama segment), gaps remaining with potential magnitudes exceeding M 8.0.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: Melnick D. et al. (2017). "The superinterseismic phase of the megathrust earthquake cycle in Chile." Geophysical Research Letters 44 (2)

ACKNOWLEDGMENTS
CP would like to thank the Initiative and Networking Fund of the Helmholtz Association through the project "Advanced Earth System Modelling Capacity (ESM)." The authors would also thank Jonathan Bedford for his helpful comments on the manuscript. The authors are grateful to the editor, QQ, and the three reviewers, QZ, C-HC, and LG, for their constructive comments. All data used are properly cited in the reference list, figures, and tables.