Forearc density structure of the overriding plate in the northern area of the giant 1960 Valdivia earthquake

The objective of this work is to analyse the density structure of the continental forearc in the northern segment of the 1960 M_w 9.6 Valdivia earthquake. Regional 2D and local 3D density models have been obtained from available gravity data in the area, complemented by new gravimetric stations. Models are constrained by independent geophysical and geological information and new TEM and MT soundings. The results show a segmentation of the continental wedge along and perpendicular to the margin, highlighting a high-density anomaly, below the onshore forearc basin, that limits the late Paleozoic–early Mesozoic metamorphic basement in the region where Chaitenia terrane has been proposed. A progressive landward shift of this anomaly correlates with the high slip patch of the giant 1960 M_w 9.6 Valdivia earthquake. Based on these results, we propose that the horizontal extension of the less rigid basement units conforming the marine wedge and Coastal Cordillera domain could modify the process of stress loading during the interseismic periods, and also that changes in position and extension of the late Paleozoic–early Mesozoic accretionary complex could be linked with the frictional properties of the interplate boundary. This analysis provides new evidence of the role of the overriding plate structure in the seismotectonic process in subduction zones.

1 Introduction

The physical structure of the oceanic and continental plates have had an important role in the long- and short-term deformation process of the subduction margins. On the other hand, the tectonic activity has modified the internal structure and geometry of the tectonic plates (i.e. Bilek et al., 2003; Hackney et al., 2006; Hicks et al., 2014; Contreras-Reyes and Carrizo, 2011; Bassett and Watts, 2015; Poli et al., 2017). This geodynamical feedback is evinced by spatial correlations between the physical segmentation of the continental wedge, and ruptures of large megathrust earthquakes (i.e. Contreras-Reyes et al., 2010; Li and Liu, 2017; Martínez-Loriente et al., 2019; Molina et al., 2021). Examples of this are the spatial correlation between gravity (density) anomalies in the continental wedge and the location of high slip patches in large earthquakes (> 7.5–8 M_w, Song and Simons, 2003; Wells et al., 2003; Álvarez et al., 2014; Bassett and Watts, 2015; Bassett et al., 2016; Schurr et al., 2020), which suggests that changes in normal stresses on the seismogenic zone have a role on the seismic rate and slip propagation during large earthquakes (Tassara, 2010; Maksymowicz et al., 2015, 2018; Molina et al., 2021). On the other hand, changes of the continental wedge geometry have been associated with variations of the interplate boundary friction at the maximum slip patches of the large 2011 M_w 9.0 Tohoku-Oki, 2010 M_w 8.8 Maule and 1960 M_w 9.6 Valdivia earthquakes (Cubas et al., 2013a, b; Maksymowicz, 2015; Contreras-Reyes et al., 2017; Molina et al., 2021).

Diverse works have highlighted the importance of the transition between accretionary prisms (or highly fractured frontal units) and the more rigid rocks of the continental basement as a tectonic limit, controlling, at least partially, the upward propagation of coseismic slip, foreshocks and aftershocks during large megathrust earthquakes (Scholz, 1998; Contreras-Reyes et al., 2010; Moscoso et al., 2011; Kodaira et al., 2012; León-Ríos et al., 2016; Maksymowicz et al., 2017, 2018; Tsuji et al., 2017). At the same time, the downdip limit of the megathrust earthquakes has been related (among other factors) to physical properties of the mantle wedge and deep interplate boundary (Peacock and Hyndman, 1999; Seno, 2005; Wang et al., 2020), which are modified by fluid subduction, slab dehydration and the presence of basal accretionary complexes (Moreno et al., 2018; Menant et al., 2019). Nevertheless, less attention has been paid to the internal physical structure (and lithology) of the continental crust above the downward limit of the megathrust, even considering that all forearc units above the fragile–ductile limit should work as a part of the same mechanically coupled system (van Dinther et al., 2012; Comte et al., 2019).

In this context, we have explored the continental forearc density structure of the Nazca–South America subduction zone in a segment where the high slip patch of the giant 1960 M_w 9.6 Valdivia earthquake ruptured (Fig. 1). As mentioned before, this slip patch correlates not only with a low gravity anomaly above the marine forearc (Wells et al., 2003) and low continental slope angles (Maksymowicz, 2015), but also, with a landward extension of Paleozoic metamorphic outcrops on the shore (Fig. 1a). Furthermore, ages and petrological data of continental basement rocks (metamorphic and plutonic rocks) suggest a complex ancient history of accreted terranes (Ramos et al., 1986; Rapalini, 2005) that constitutes the present continental crust in the area (Fig. 1b). Particularly, the recent proposal of an oceanic terrane accreted against the Gondwana margin during Devonian times (Chaitenia, Hervé et al., 2016, 2018, Ct in Fig. 1b) could determine changes in the internal structure of the continental forearc, southward of ∼ 40^∘ S. However, the exact limits of this basement configuration remain poorly constrained. In order to reveal the crustal structure of this active portion of the Chilean margin, this work presents the results and interpretation obtained from regional 2D density models, extended from Nazca plate to the Andes Cordillera (Fig. 3), and a local 3D density inversion of the continental forearc (red rectangle in Figs. 2 and 3). The models include magnetotellurics (MT) and transient electromagnetic (TEM) measurements, as well as available independent geophysical and geological data to constrain forward modelling and 3D inversion.

Figure 1 Geotectonic settings of the studied zone. (a) Tectonic and morphostructural features above the coloured bathymetric/topographic elevation grid. Metamorphic and plutonic/intrusives outcrops are shown. Oceanic fracture zones (Mocha, Valdivia and Chiloé) are indicated in black. The deformation front (continental wedge toe) at the Nazca–South America trench is indicated by a bold black line. The blue lines correspond to continental structures identified at the surface (SERNAGEOMIN, 2003; Melnick and Echtler, 2006), including the Liquiñe–Ofqui fault system (LOFS), Mocha–Villarrica fault zone (MVFZ) and Lanalhue fault zone (LFZ). The red triangles illustrate active volcanoes. The red arrow indicates the direction of Nazca–South America convergence, and the green lines represent the iso-slip contours of the Valdivia earthquake according to Moreno et al. (2009). (b) Schematic map of basement units after Hervé et al. (2018) and other elements as in (a).

2 Geotectonic settings

The study zone, located between 38.5 and 42.5^∘ S (Fig. 1), is part of the south-central Chilean margin, where the oceanic Nazca plate subducts beneath the continental South American plate. The current rapid convergence rate (∼ 6.6 cm yr⁻¹, Kendrick et al., 2003; Vigny et al., 2009) determines high seismotectonic activity, including the occurrence of mega-earthquakes such as the giant 1960 M_w 9.6 Valdivia earthquake (the largest instrumentally registered worldwide). In the long term, this subduction process has been continuously active since Jurassic times (Charrier et al., 2007), being superimposed to ancient tectonic processes of Gondwana structuration (Ramos et al., 1986), and generating the current configuration of the continental South American plate western border. Marine seismic studies, to the north and south of the study area, indicate that the structure of the continental wedge shows physical and tectonic segmentation from the trench to the coast, characterized by active accretionary prisms along the lower slope regions, compressional geometries, and the development of confined slope basins inside the middle and upper slope, while the shelf region exhibits forearc basins with a complex deformation style structured by normal and inverted faults (Bangs and Cande, 1997; Geersen et al., 2011; Becerra et al., 2013; Bangs et al., 2020). Consistently, V_p models derived from wide-angle seismic refraction, at 38^∘ S (Contreras-Reyes et al., 2008) and south of 43^∘ S (Contreras-Reyes et al., 2010), present changes in the deep structure of the continental wedge that can be interpreted as transitions between the accretionary prism, paleo-accretionary rocks and continental basement. Regarding the geometry of the marine forearc, the continental wedge shows a narrow continental slope (defined between the deformation front and shelf break; see Fig. 1) to the south of ∼ 41^∘ S. This morphological change corresponds to a decrease in the slope angle at the northern region of the Valdivia earthquake rupture (38.5–41^∘ S), which in turn can be interpreted as a decrease in the effective friction coefficient (${\mathit{\mu }}_{\text{b}}^{\ast }$) at the interplate boundary (Dahlen, 1984; Cubas et al., 2013b; Maksymowicz, 2015).

Onshore, three major trench-parallel morphostructural units from west to east can be observed (Fig. 1a): (1) the Coastal Cordillera (CC), where old rocks of a paired metamorphic belt are exposed (Hervé, 1988); (2) the Central Depression (CD) characterized by the presence of unconsolidated Quaternary sediments overlaying Cenozoic deposits (Jordan et al., 2001); and (3) the Principal Cordillera (PC), where the active volcanic arc is currently located. In a close spatial relation with the volcanic arc, the prominent Liquiñe–Ofqui Fault System (LOFS, Fig. 1) stretches along more than 1000 km between 37 and 46^∘ S (Cembrano et al., 1996). This continental structure has been interpreted as right-lateral strike-slip system that currently concentrates most of the crustal intraplate seismic activity in response to oblique Nazca–South America convergence (Lange et al., 2008; Orts et al., 2012) and exhumation at these latitudes (Adriasola et al., 2005; Glodny et al., 2008). Moreover, numerous tectonic lineaments and fault zones have been described (SERNAGEOMIN, 2003; Melnick and Echtler, 2006), generally showing north-west and north-east orientations. According to Melnick et al. (2009), the kinematics of LOFS generates intense deformation in its northern limit, explaining the deformation associated with large north-west strike continental faults (as LFZ) and the eastward bending of the CC.

Accretionary metamorphic complexes, associated with late Paleozoic–early Mesozoic subduction, are exhumed along the study zone (Hervé, 1988; Duhart et al., 2001; Willner et al., 2004; Hervé et al., 2013). These units correspond to a paired metamorphic belt, which includes the Western and Eastern Series (WS/ES) formed under high $P/T$ and low $P/T$ conditions, respectively. WS has been interpreted as a basal accretionary complex while ES is interpreted as a frontal accretionary prism and/or as the shallow sedimentary units deformed by the basal underplating of WS units (Glodny et al., 2005; Willner et al., 2005). This paired metamorphic belt is observed continuously at the CC, but the width of their outcrops varies along the margin (see Fig. 1a). Between ∼ 38 and 40^∘ S, and southward of ∼ 41.5^∘ S, outcrops of WS are observed eastward, near the western limit of PC. Thus, between ∼ 40 and ∼ 41.5^∘ S, the eastern limit of these units is not defined due to the presence of the CD deposits and could form most of the forearc basement, or it could be confined near the coast. Westward of accretionary metamorphic complexes and north of 38^∘ S, the Coastal Batholith (late Paleozoic intrusive rocks) is observed along CC, but southward (in the study zone) the outcrops of this ancient volcanic arc bends to the southeast and becomes part of the PC. Younger Plutonic and intrusive rocks, related to magmatic arcs from Mesozoic to Cenozoic times (Andean tectonic cycle), are observed along the PC near the position of the active volcanic arc and the LOFS, forming the North Patagonian Batholith (Charrier et al., 2007; Hervé et al., 2018; SERNAGEOMIN, 2003; SEGMAR, 1997; see Fig. 1a).

The continental crust of the western border of South America was configured, during Paleozoic times, by collisions of allochthonous terranes against Gondwana (Rapalini, 2005). To the north of the study zone, Chilenia terrane (Ch in Fig. 1b) collided during Devonian times (Ramos et al., 1986; Hyppolito et al., 2014, and references therein), but its southern extension is roughly defined and could be present in the northern area of the study zone. Southward, the geodynamic evolution of the margin during Devonian to Triassic times has been explained with a double subduction system (Hervé et al., 2016). These authors proposed the development of an island arc (named as Chaitenia, Ct in Fig. 1b) parallel to the margin colliding with Gondwana during Carboniferous times (Hervé et al., 2016, 2018; Rapela et al., 2021, Ct in Fig. 1b). If this hypothesis is correct, the continental crust of the current forearc corresponds to Chaitenia, south of ∼ 40^∘ S. However, it is important to point out that the limits between all these terranes are poorly constrained in the study zone owing to the scarcity of basement outcrops.

3 Data and methods

3.1 Gravimetric database and processing

We compile a gravimetric database (see Fig. 2), including public databases and new measurements in the studied area. The resulting merged database includes (1) onshore gravimetric data acquired by Chilean and European institutions in the Central Andes from 1982 to 2006, originally compiled by Schmidt and Götze (2006); (2) 167 new gravimetric stations acquired by our group in 2019; (3) marine gravimetric profiles available in the GEODAS database data (NOAA); and (4) a satellite gravimetric grid from Sandwell and Smith (https://topex.ucsd.edu/cgi-bin/get_data.cgi, last access: 20 December 2021, Sandwell and Smith, 2009; Sandwell et al., 2014) to cover marine gaps and regions to the south of 42^∘ S. Bathymetric and topographic databases merge an onshore elevation grid (SRTM elevation grid, Jarvis et al., 2008) and swath bathymetry data of the studied zone (Flueh and Grevemeyer, 2005), complemented by Global Topography V18.1 (Smith and Sandwell, 1997).

The new gravimetric data were distributed to fill in some observed gaps in onshore studies, and to complement and validate gravity and topographic information from old stations. The gravity acquisition was made using a Lacoste & Romberg G-411 gravimeter with a digital upgrade (http://www.gravitymeter-repair.com, last access: 20 December 2021) funded by ANID-FONDECYT project no. 11170047. Elevation was obtained by differential GPS using Topcon HiperV instruments of the University of Chile (DGF). GPS data were processed with the permanent GPS bases of the Chilean national seismological network (Centro Sismólogico Nacional, http://www.csn.uchile.cl/red-sismologica-nacional/red-gps/, last access: 20 December 2021), and the new gravity measurements were tied to the absolute gravity stations available in the study zone (International Gravimetric Bureau (BGI), https://bgi.obs-mip.fr/, last access: 20 December 2021). Estimated precision of new gravity measurements is under ±0.01 mGal (10⁻⁵ m s⁻²) and obtained elevation errors of differential GPS data are under ±0.5 m. The data were corrected to obtain the complete Bouguer anomaly (CBA) using standard correction processes (Blakely, 1995; Lowrie, 2007): tide correction, instrumental drift correction using daily repetitions at base stations, normal gravity correction, free-air, Bouguer, and terrain corrections. These processes were conducted considering a 2.67 g cm⁻³ reduction density. Earth tide correction was removed from the new data according to the Longman (1959) algorithm. Normal gravity correction of new data considered the subtraction of the theoretical gravity of the WGS-84 ellipsoid. Free-air correction of all onshore data was calculated as 0.3086h (mGal), where h is ellipsoidal height in metres (Lowrie, 2007). Due to the inhomogeneities in the elevation measurement techniques used in old onshore data acquisition (registered between 1982 and 2002, Schmidt and Götze, 2006), we prefer to use SRTM elevation data to perform the free-air and Bouguer corrections of these old gravity data. The terrain correction of all data was calculated following a combination of the algorithms proposed by Kane (1962) and Nagy (1966) and with a high-resolution SRTM elevation grid. The terrain correction includes topographic data located up to ∼ 300 km around each station. GPS data processing, gravity data processing and all figures presented here consider geographic coordinates in the datum WGS84 and WGS84-18S for UTM coordinates.

The spatial coverages of different gravity databases (satellite, marine and onshore) present areas of interception (Fig. 2) where they can be compared to determine the average gravity differences (constant average shifts). These shifts were used to generate a merged database levelling all data to the values observed in the new acquired data. The free-air values of the onshore stations were used to move the free-air anomaly of Sandwell and Smith (satellite data) to the same level (adding a shift of −17.3 mGal to Sandwell and Smith data), and finally GEODAS free-air data (marine lines) were levelled and merged with the other data (adding a shift of −24.78 mGal to GEODAS data) to calculate the CBA (Fig. 3).

3.2 Density modelling

3.2.1 2D regional forward gravity models

In order to study the regional structure of the continental wedge and subduction zone, we modelled five profiles (P1_Toltén, P2_Unión, P3_Osorno, P4_Llanquihue and P5_Chepu; see Fig.2), which run perpendicular to the trench at latitudes of 39.25, 40.2, 40.5, 41 and 42^∘ S, respectively. These profiles were extracted from the regional complete Bouguer anomaly grid obtained from the merged gravity database (Fig. 3). 2D forward modelling was performed by using the GravGrad modelling scheme (Maksymowicz et al., 2015), allowing the calculation of the gravimetric response of a stack of layers with arbitrary shape. The densities inside each layer can be varied along the vertical and horizontal directions. As the gravity anomalies are not exclusively dependent on the density structure below each data (i.e. it should be modelled considering masses around the profile), the elevation (bathymetry and topography) in 2D modelling is an averaged elevation profile which includes data inside a ∼ 40 km wide band around the profile (i.e. averaging the elevation to 20 km on each side of the gravity profile). A 40 km wide band is a reasonable assumption considering that the wavelengths of the CBA gravity anomalies along the profiles are mostly larger than ∼ 40–50 km. On the other hand, this parameter is not critical for the obtained 2D model. In fact, a completely different value (e.g. 10 km wide) can be considered with minor modification in the resulting density model (see more details in the Supplement).

As a forward modelling procedure, GravGrad allows the user to iteratively modify the geometry and densities of all layers in the model to fit the observed gravity anomalies (CBA in this case). Section 3.3 and Fig. 2 describe the independent information used to constrain the slab geometry, continental Moho depths and sedimentary thickness at the CD basin. The general density structures of the slab, continental plate and mantle were based on seismic and seismological V_p models available at the zone, converted to density by the empirical Nafe–Drake transformation curve (Brocher, 2005).

Figure 2 Geophysical database in the study zone and the location of 2D and 3D models. The blue lines indicate the location of five 2D regional density model, and the red rectangle is the zone in which local 3D inversion was obtained. The black dots designate gravity stations compiled by Schmidt and Götze (2006), onshore, and in GEODAS database (NOAA) offshore. The green dots illustrate the complementary gravity stations acquired by our group under the ANID-FONDECYT project no. 11170047. The blue and red dots correspond to TEM soundings and MT stations acquired under the aforementioned FONDECYT project. Magenta dots indicate the MT stations presented by Segovia et al. (2021) and cyan dots correspond to TEM soundings published by DGA (2012). Receiver function profiles obtained by Dzierma et al. (2012a) are shown with cyan lines. The location of seismic V_p–depth models (Contreras-Reyes et al., 2008, 2010; Bangs et al., 2020; Maksymowicz et al., 2021) and the seismic reflection lines presented by Jordan et al. (2001) and González et al. (1989) are indicated with magenta lines. The yellow dots correspond to the location of boreholes (McDonough et al., 1997).

3.2.2 3D gravity inversion

Regarding a more detailed analysis of the continental density structure onshore, a 3D inversion was performed in a central patch of the study zone, where a large gravity maximum is observed parallel to the Coastal Range (red polygon in Fig. 3). This onshore 3D density model was obtained using the UBC-GIF GRAV3D v3.0 software (Li and Oldenburg, 1998). The algorithm inverts the residual Bouguer anomaly (RBA) to derive a 3D density anomaly model of the crust. The residual Bouguer anomaly was generated by subtracting a first-order polynomial trend from complete Bouguer anomaly data (see Supplement). The 3D mesh has $\mathrm{67}\times \mathrm{80}\times \mathrm{102}$ blocks (in X, Y, Z directions, respectively). The horizontal mesh size is 3000 m × 3000 m. Due to the progressive sensitivity decrease in the gravity inversion to sources at depth, UBC-GIF developers recommend using smaller cells near the surface and increase the cell thickness with the model depth (https://www.eoas.ubc.ca/ubcgif/iag/index.htm, last access: 20 December 2021). Accordingly, the cell size gradually grows from 100 to 1500 m in our model, reaching 70 km in depth.

A total of 3514 onshore gravity data points were used for the 3D inversion, generating a regular spaced grid (with a grid size of 3000 m × 3000 m). After numerous preliminary experiments, we set the length scale parameters of the UBC-GIF algorithm (Li and Oldenburg, 1998) as 6000, 6000 and 3000 m in X, Y and Z directions, respectively. These length scale parameters define the horizontal and vertical smoothness of the solution, and preferred values are the double of the horizontal and vertical cell size used to discretize the media. This criterion is one of those recommended in the software manual (https://www.eoas.ubc.ca/ubcgif/iag/index.htm, last access: 20 December 2021), but also it is important to highlight that sensitivity tests show that under a strong variation of these parameters the obtained solution of density structure shows similar features (see details in the Supplement). The 3D inversion was constrained by information from the geological map of 1:1 000 000 scale (SERNAGEOMIN, 2003), MT stations, TEM stations and onshore seismic lines (Fig. 2). Accordingly, we performed the inversion with the following definitions for surface sediments and the basement: a minimum homogeneous thickness of 500 m was assigned to the entire area of the model where Quaternary sediments are exposed. These cells can take densities between 1.9 and 2.1 g cm⁻³. The next 500 m (five cells) in depth correspond to a transition zone, where the blocks could be sediments or rock and can vary between 1.9 and 2.7 g cm⁻³. The next 300 m (three cells) in depth corresponds to a second transitory zone, where the blocks could be fractured rock or consolidated rock and can vary between 2.4 and 3 g cm⁻³. Below, the blocks corresponding to the basement can take values between 2.5 and 3 g cm⁻³. Finally, below 7500 m depth we constrained the model to have greater densities than the background (2.67 g cm⁻³). Then, those deep cells can take values between 2.67 and 3 g cm⁻³ in order to ensure more realistic vertical gradients in the lower constrained deep portion of the model.

To include the presence of the main lakes in the zone, the model is forced to be water in the blocks that correspond to lake, assigning them a density of 1 g cm⁻³. The bathymetry of the two first lakes was obtained from Chilean National Oceanographic Service (http://www.shoa.cl/php/inicio, last access: 20 December 2021), and in the case of Llanqihue Lake, a mean of 200 m of depth was considered. Similarly, to consider the gravimetric effect of the sea, the model was forced to be water in the blocks above bathymetry.

3.3 Geophysical constraints

3.3.1 Available geophysical information

Independent geophysical data were used to constrain the 2D density models and onshore 3D inversion (Fig. 2). This information includes (1) the available 2D velocity–depth models at different latitudes (Contreras-Reyes et al., 2008, 2010; Bangs et al., 2020; Maksymowicz et al., 2021), used as a reference for the general structure of the oceanic plate and marine continental wedge after V_p-to-density conversion according to the empirical Nafe–Drake transformation curve (Brocher, 2005); (2) interpretation of reflection seismic profile (in depth) at ∼ 42^∘ S (González et al., 1989); (3) the Quaternary sedimentary thickness and the top of the Paleozoic basement observed in the onshore ENAP seismic lines Z5B-010A and ZDO-001 (McDonough et al., 1997; Jordan et al., 2001) and ENAP boreholes (McDonough et al., 1997; Honores et al., 2015); (4) the SLAB 2.0 model (Hayes et al., 2018) to constrain the deep slab geometry; (5) Moho depth along the profiles presented by Dzierma et al. (2012a); and (6) 1D electrical resistivity models using magnetotelluric measurements obtained by Segovia et al. (2021), 1D electrical resistivity models from TEM measurements presented by DGA (2012), and 1D electrical resistivity models from new MT and TEM measurements. These resistivity models were used to constrain the thickness of young sedimentary fill at CD. At the MT/TEM stations where 1D resistivity models do not reach the base of the young sedimentary fill (by cultural electromagnetic noise or limited penetration in thick sedimentary fill areas), we define values of minimum sedimentary thickness, aiming to decrease uncertainties in the density modelling.

3.3.2 Electromagnetic methods to constrain gravity measurements

The new magnetotelluric data (red dots in Fig. 2) were collected using Metronix ADU-08 data loggers and MFS-07 induction coil magnetometers along with Pb-PbCl electrodes. Time series data were recorded between 12 and 24 h. All sites were processed using the robust method based on Egbert and Booker (1989). The 1D resistivity models of new and previously measured data (from Segovia et al., 2021) were obtained using occam (Constable et al., 1987 and Bostick, 1977) algorithms implemented in WinGLink (Schlumberger, version 2.21). See data and models in the Supplement.

The transient electromagnetics measurements (blue dots in Fig. 2) were carried out utilizing the ABEM WalkTEM (ABEM, 2016). In general, a central loop set-up was used with a transmitter loop size of 100 m × 100 m or 40 m × 40 m. TEM stations were modelled by using Interpex-IX1D TEM software, generating 1D resistivity depth models using the ridge regression algorithm (see data and models in the Supplement).

4 Results

4.1 Complete Bouguer gravity anomaly (CBA)

It is necessary to describe the main characteristics observed in the CBA at the study zone (Fig. 3) before analysing density models. The general aspect of the CBA is a sequence of bands with high and low gravity, roughly parallel to the margin. Offshore, we observe the low CBA associated with the deep trench seafloor and its sedimentary fill. It is important to notice that this low CBA extends several kilometres landward from the deformation front (toe of the continental wedge), which implies the presence of low-density units at the lower slope of the continental wedge. The main feature observed in the slope and shelf area is the low CBA zone (L1 in Fig. 3a) extended from ∼ 38.5 to ∼ 41^∘ S, correlating with a decrease in general slope angles at the same latitudes (Fig. 1). This morphologic and gravimetric anomaly is also correlated with the maximum slip patch of the giant Valdivia earthquake, as highlighted by several authors (Wells et al., 2003; Maksymowicz, 2015; Contreras-Reyes et al., 2017; see Fig. 3b).

Onshore, the regional aspect of CBA is an eastward trend of gravity decreases from the coast to PC, mostly related to the presence of continental roots below the Andes (Tassara et al., 2006; Tašárová, 2007). Therefore, high CBA anomalies are observed along the coast (Fig. 3a), but their amplitude decreases between ∼ 38.7 and ∼ 40^∘ S, where a relatively low CBA anomaly is observed (L2 in Fig. 3a). It is necessary to consider that L2 is spatially correlated with a zone of landward extension of the CC and metamorphic complex outcrops (WS/ES).

To the east, a sequence of gravity lows, with sparse gravity maximums, is correlated with the eastern part of the CD basin, the current volcanic arc and LOFS, suggesting a complex density structure at the PC zone. Between ∼ 40 and ∼ 41.5^∘ S (and probably southward) a prominent positive anomaly can be seen above the western portion of the CD basin (H1 in Fig. 3a), indicating the presence of a high-density body elongated to the north-east and covered by sedimentary fill of CD. This interesting forearc gravity maximum was observed by Hackney et al. (2006) based on the same onshore data, and has been confirmed by our new complementary stations.

Figure 3 Complete Bouguer anomaly (CBA) at the study zone. (a) The main gravimetric forearc features observed in the area (see main text for details). The blue lines illustrate the tracks of 2D forward density models, and the red rectangle indicates the zone where 3D density inversion was carried out. Dotted black ellipses show the approximated extension of anomalies L1, L2 and H1 along with other elements seen in Fig. 1a. (b) Crustal tectonic structures and CBA. The figure includes the contours of the metamorphic Western and Eastern Series and late Paleozoic batholith outcrops. Dotted black lines are interpreted gravimetric lineaments, e.g. Valdivia–Futrono lineament (VFL), Bahía Mansa–Choshuenco lineament (BMCHL) and Osorno lineament (OL). The other elements are as in Fig. 1a.

By a visual inspection of the CBA map, we interpret a set of gravity lineaments (Fig. 3b) with north-west and north-east directions. Derivative filters (directional derivatives, slope gradient and analytical signal) applied to CBA helped to identify these regional trends (see Supplement). This qualitative interpretation confirms the location of fault zones previously identified at the surface (SERNAGEOMIN, 2003; Melnick and Echtler, 2006), suggesting their continuity through the forearc and, in some cases, their seaward extension (e.g. Valdivia–Futrono lineament, VFL in Fig. 3b). Additionally, new gravimetric lineaments are identified in CBA, suggesting the presence of large structures affecting the basement units (e.g. Bahía Mansa–Choshuenco lineament, BMCHL in Fig. 3b). The H1 anomaly is limited to the north-west by the Osorno lineament (OL in Fig. 3b), which presents continuity with an identified west-dipping reverse fault in the south-west (SERNAGEOMIN, 2003; Melnick and Echtler, 2006; Hackney et al., 2006; Encinas et al., 2021), indicating that the geometry of H1 has a structural and tectonic control.

4.2 2D density profiles

Figure 4 shows the results of the 2D forward gravity models obtained through the five studied profiles. As observed in Fig. 4a, the modelling process allowed a good fit to be attained between observed and calculated gravity data, associated with low RMS values (≤4.0 mGal) in comparison with the total amplitude of the gravity anomaly. According to these results, the marine structure of the overriding plate can be described as a general landward increase in density between the deformation front (DF at the trench) and the coastal area, where it is possible to define at least two internal units (Fig. 4b to f). The first one corresponds to a frontal low-density unit ($\mathit{\rho }<\sim \mathrm{2.5}$ g cm⁻³) of about 25–35 km width, with a rapid landward horizontal density gradient. This frontal unit is roughly correlated with the lower slope of the continental wedge. The second unit is characterized by a lower landward horizontal density gradient and shows densities between ∼ 2.5 and ∼ 2.8–2.9 g cm⁻³. This middle wedge unit is extended from the lower slope to the coast by ∼ 70 km at profiles P2 to P5 and is slightly wider (∼ 90 km) at the northernmost profile P1_Toltén. Immediately below the seafloor, all profiles present marine forearc basins with variable thicknesses (< 5 km) and densities lower than ∼ 2.3 g cm⁻³. A few kilometres westward from the coast, a continental wedge shows a transition to higher densities landward (higher than ∼ 2.9–3.0 g cm⁻³ in the deep portion of the crust). This transition can be described as a landward limit of the middle wedge unit and seems to have a west-dipping geometry.

Onshore, the upper portion of the continental forearc (the upper ∼ 10–15 km) displays a sequence of low- and high-density zones. Below CC (and metamorphic complex outcrops), the shallow densities are generally higher than 2.5 g cm⁻³ and downward. However, this region is not particularly dense. In fact, below the sedimentary fill of the CD basin we observe a high-density maximum in the five 2D profiles (H1 in Fig. 4b to f). Then, the results confirm the presence of a high-density zone associated with the high CBA anomaly already described (D1 in Fig. 3a), suggesting its prolongation to the north-east and south-west. Comparing the 2D profiles, we notice that D1 is progressively closer to the coast, southward from profile P2_Unión (∼ 40^∘ S); i.e. D1 presents a north-east trend, as suggested before in the CBA description. It is important to note that the presence of D1 is clear in all profiles except for P2_Unión, in which this density anomaly is slightly raised from a more homogenous model of the upper continental crust.

To the east of D1, all profiles show another high-density zone (D2 in Fig. 4b to f). It is important to consider that the large LOFS approximately correlates with the western limit of D2 in profile P1_Toltén, with the eastern limit of D1 in profile P2_Unión and with the eastern border of D2 at profiles P4_LLanquihue and P5_Chepu, which suggests a structural relation between the deep geometry of the high-density anomalies (D1 and D2) and LOFS.

Figure 4 2D regional forward models (a) complete Bouguer anomaly (CBA) curves along five modelled profiles. Coloured curves are the modelled gravity signals, and the corresponding observed data are presented in grey. (b) Density–depth model along the profile P1_Toltén, morphological limits of the continental wedge (as lower slope, LS) are indicated. Dotted black lines illustrate the approximate limits of frontal low-density unit and the middle wedge/shelf unit of the marine continental wedge. The thin vertical line indicates the limit between CC and CD. Red triangles correspond to the active volcanoes located near the profile. Pink dots correspond to slab geometry according to the SLAB2.0 model (Hayes et al., 2018). Green dots depict the continental Moho depths obtained by receiver function analysis (Dzierma et al., 2012a). Black dots indicate the base of the poorly compacted shallow sedimentary layer according to MT and TEM soundings, and cyan dots correspond to the minimum thickness of this sedimentary layer according to MT and TEM soundings. White lines limit a zone of high V_p and low $V_{\text{p}}/V_{\text{s}}$ obtained by Dzierma et al. (2012b). Thin dotted yellow and cyan lines limit electrically conductive and resistive zones (C and R) according to Kapinos et al. (2016). (c) Density–depth model along the profile P2_Unión. (d) Density–depth model along the profile P3_Osorno. Grey dots show the base of the poorly compacted shallow sedimentary layer according to onshore seismic profiles and ENAP boreholes (McDonough et al., 1997; Jordan et al., 2001; Honores et al., 2015). Other elements as in (a). (e) Density–depth model along the profile P4_Llanquigue. Thin dotted yellow and cyan lines limit electrically conductive and resistive zones (C and R) according to Segovia et al. (2021). Other elements as in (d). (f) Density–depth model along the profile P5_Chepu. Grey dots represent the base of shallow sedimentary layers according to seismic reflection data (González, 1989) and other elements as in (a). Individual figures of each profile are presented in the Supplement.

4.3 3D local density model

As explained before, a large latitudinal change of onshore forearc continental structure is observed in the central profiles (P2_Unión, P3_Osorno and P4_Lannquihue), where D1 seems to have a north-east trend and where D2 is observed near the LOFS and the arc. This motivates the development of a local 3D density inversion in the forearc area, to derive the detailed structure of upper continental crust with an independent model strategy.

The 3D inversion modelled the input residual Bouguer anomaly (Fig. 5a) with high precision, as is observed in Fig. 5b, where differences between modelled and observed data are in general lower than ±1 mGal. The results show density contrast anomalies to about 20 km depth (Figs. 6, 7 and in the Supplement), which means that deeper anomalies are mostly contributing to a regional linear trend of the CBA at the scale of 3D local inversion.

Figure 5 (a) Residual Bouguer anomaly (RBA) used as input for 3D inversion. (b) Difference between final modelled and input RBA data.

Four constant depth slices through the final 3D model are shown in Fig. 6. In the slice at 49 m b.s.l. (Fig. 6a) it is possible to observe the geometry of the CD basin as yellow zones that correspond to density contrasts between −0.77 and −0.57 g cm⁻³. Also, it is possible to observe the areas that correspond to seawater as blue zones with density contrasts of −1.67 g cm⁻³ and sedimentary areas below the lakes with density contrasts of about −0.73 g cm⁻³. Similar density structures are observed in Fig. 6b (at 1049 m b.s.l.).

Figure 6 Density slices from the 3D density model at different depths. (a) Slice at 49 m b.s.l. Blue lines correspond to continental structures identified at the surface (SERNAGEOMIN, 2003; Melnick and Echtler, 2006); VFL highlights the Valdivia–Futrono lineament, and OL corresponds to Osorno lineament (see Fig. 3b). Segmented green line illustrates the border of CD, and the external limits of the lakes are highlighted with a cyan line. Active volcanoes are marked by black triangles. Segmented black lines indicate the approximate borders of H1 and H2 density anomalies (see main text for details) defined at the deeper slice shown in (d). Partially hidden zones are outside of onshore data considered for 3D inversion. (b) Slice at 1049 m b.s.l. Note the change of colour scale in relation to (a). Other elements as in (a). (c) Slice at 5149 m b.s.l. Other elements as in (b). (d) Slice at 7149 m b.s.l. and other elements as in (b) .

Figure 6c and d show slices at 5149 and 7149 m b.s.l., respectively. In this figure it is feasible to observe the prominent high-density zone D1 under the CD, being consistent with the observed one in CBA and 2D models (Figs. 3 and 4). In the 3D model, D1 covers ∼ 230 km along the strike and ∼ 80 km in the horizontal direction, being oriented ∼ 25^∘ E in the same direction as the western border of CC. The density contrast of this structure is higher than 0 g cm⁻³ in most of the areas, reaching 0.3 g cm⁻³ at denser zones. South of ∼ 40.75^∘ S, D1 is limited with WS outcrops to the west, while to the north-west it is bordered by the OL lineament (see Fig. 3b). The northern limit of D1 correlates with the presence of an VFL lineament. To the east, a low-density lineament in the same D1 direction, with density contrasts from −0.17 to 0 g cm⁻³, can be seen. This low-density band is about 10 to 15 km wide and is limited to the east with the high-density zone D2. The 3D density model shows variations inside D1 and D2, which are formed by high-density zones (density contrast ≥ 0.3 g cm⁻³) merged with lower-density areas.

To analyse the density variation with depth, Fig. 7 presents six W–E vertical cross sections of the 3D model at different latitudes (c1 to c6 in Fig. 6a). All cross sections show that D1 is below the CD and correlated to the eastern limit of the CC and WS outcrops, confirming the results obtained by 2D regional modelling. The top of the D1 anomaly is obtained around 5 km b.s.l., displaying a trend of deepening to the north (also suggested by 2D modelling). At a lower scale, the geometry of D1 is characterized by two lobes, and its western and eastern borders seem to be tilted to the east and west, respectively (segmented black lines in Fig. 7). Eastward, D2 is modelled in the southern region at profiles c4, c5 and c6 (Fig. 7d to f), and it is also characterized by two lobes. The western limit of this anomaly suggests an inclination to the east.

Figure 7 Density–depth profiles extracted from the 3D inversion model (see location in Fig. 6a). (a) Profile c1 at UTM-North = 5 560 500 m (WGS84-18S). Thin vertical lines indicate the limits of CC, CD. Black dots mark the base of shallow sedimentary unit, according to TEM/MT measurements, and cyan dots correspond to the minimum sedimentary fill according to the TEM/MT soundings that do not reach the basement. Grey and blue dots show respectively the base of Quaternary sediment and the top of the Paleozoic basement according to onshore seismic profiles and ENAP boreholes (McDonough et al., 1997; Jordan et al., 2001; Honores et al., 2015). The interpretation of the approximate borders of the D1 anomaly is highlighted with dotted black lines. The partially hidden zones are outside of onshore data considered for 3D inversion. (b) Profile c2 at UTM-North = 5 530 500 m (WGS84-18S). (c) Profile c3 at UTM-North = 5 500 500 m (WGS84-18S). Grey and blue lines show the shallow sedimentary unit and the top of the metamorphic basement, according to seismic reflection profile ZDO-001 (see main text for details) and other elements as in (b). (d) Profile c4 at UTM-North = 5 467 500 m (WGS84-18S). The approximate borders of the high-density anomalies D1 and D2 are highlighted with dotted black lines. (e) Profile c5 at UTM-North = 5 440 500 m (WGS84-18S). Grey and blues lines show the shallow sedimentary unit and top of the metamorphic basement, according to seismic reflection profile Z5B-010A (see main text for details) and other elements as in (d). (e) Profile c5 at UTM-North = 5 410 500 m (WGS84-18S).

5 Interpretations and discussions

These obtained results exhibit a landward segmentation of the continental wedge density structure observed from trench to arc (see an interpretative schema at Fig. 8a). Offshore, the frontal portion of the continental wedge (to ∼ 25–35 km landward from the deformation front) presents low densities with a rapid horizontal increment of densities, interpreted as a compaction process in the active accretionary prism along the south-central Chilean margin (Maksymowicz et al., 2015). It is also evidenced by seismic studies in the region (Moscoso et al., 2011; Tréhu et al., 2019; Contreras-Reyes et al., 2008, 2010; Bangs and Cande, 1997; Bangs et al., 2020). To the east, below the sedimentary fill of slope and shelf basins, the continental wedge is characterized by a second unit of higher density and lower horizontal density gradient (middle wedge unit, MWU). This unit can be associated with fractured basement rocks and/or more compacted units of a paleo-accretionary prism. In this sense, Contreras-Reyes et al. (2008) at ∼ 38^∘ S and Contreras-Reyes et al. (2010) at ∼ 43^∘ S interpret this unit (in V_p-depth profiles) as a paleo-accretionary prisms of an undetermined age between Mesozoic and Tertiary. On the other hand, at ∼ 39 and ∼ 40.5^∘ S, Bangs et al. (2020) suggest that the Paleozoic–early Mesozoic accretionary complex (WS/ES) can extend further seaward to the eastern limit of the active accretionary prism (seismic backstop), in accordance with the interpretation of marine seismic data (and boreholes) of González (1989) at ∼ 42^∘ S. However, the exploration boreholes presented by González (1989) were drilled in the shelf basin area and therefore do not provide direct information about the age of the continental basement in the western portion of MWU.

Landward from MWU, the next segment correlates on the surface with the morphostructural domain of CC and shows a density increase with respect to the marine wedge, but lower densities compared to continental crust below the CD and PC. Therefore, this CC domain is clearly related to the Paleozoic–early Mesozoic accretionary complexes (WS/ES) and their continuity to depth. Gravity modelling techniques do not define a downward limit of WS/ES (without independent deep constraints). Nevertheless, interpretations of seismic reflection data at ∼ 38.25^∘ S (Krawczyk et al., 2006; Ramos et al., 2018) showed the downward prolongation of WS/ES reaching deep levels near the continental Moho interface (∼ 30 km depth). As previously mentioned, the seaward limit of WS/ES is not defined by direct lithological observations; their presence beneath the shelf basin is confirmed by exploration boreholes (González, 1989). Thereafter, the relative rapid change in velocity associated with the transition between MWU and CC domain (dotted grey line in Fig. 8) is interpreted as a structural limit (rather than a lithological change of the basement). This structural limit is probably associated with the development of the shelf basin and a general seaward increase in fracturing within the continental wedge. This structural interpretation seems to be confirmed by Contreras-Reyes et al. (2008) at ∼ 38.25^∘ S, where continental intraplate seismicity (located by Haberland et al., 2006) is aligned with this limit, as well as the intraplate seismicity located at 39.5^∘ S by Dzierma et al. (2012c).

Figure 8 Schematic interpretation of continental forearc structure at regional scale. (a) Interpreted profile based on the 2D density model along P3_Osorno line. The main geological/physical units are hatched with different colours: light green for active accretionary prims, grey for Paleozoic–early Mesozoic metamorphic accretionary complex (WS/ES) and light red for high-density continental crust interpreted as evidence of Chaitenia terrane in the forearc region. Dark red areas correspond to anomalies with high-density contrast (≥ 0.1 g cm⁻³) according to the 3D model along the P3_Osorno line. Light blue indicates Mesozoic to Cenozoic deposits of slope and shelf basins and Cenozoic deposits in the CD basin. Segmented green, grey and red lines indicate the landward limits of the active accretionary prism, middle wedge unit and CC domain (see main text for details). (b) Interpretation map of continental forearc structure. Light red areas highlight the high-density zones (H1 and H2) defined at ∼ 7 km depth in the 3D density model (Fig. 6), and cyan areas correspond to high V_p–low $V_{\text{p}}/V_{\text{s}}$ zones identified at 10 km depth by Dzierma et al. (2012b). As in (a) segmented green, grey and red lines indicate the landward limits of the active accretionary prism, middle wedge unit and CC domain and other elements as in Figs. 1a and 3b.

The CC domain extends landward to the contact with the D1 anomaly (dotted red line in Fig. 8). The eastern border of CC range at the surface correlates almost exactly with the western border of D1 in the 3D model (Figs. 6 and 7), which is also observed in 2D regional profiles (Fig. 4). Accordingly, we understand that the continental crust of the CC domain is deformed against a denser (and probably more rigid) block of the continental crust observed here as the D1 anomaly. As mentioned before, the lineament OL (Figs. 1b, 6 and 8) is continued to the south by a west-dipping reverse fault that limits CC and CD (SERNAGEOMIN, 2003; Melnick and Echtler, 2006; Hackney et al., 2006; Encinas et al., 2021). This is an example of the contractional deformation styles that could be generated in the eastern border of CC by the depth contact between the CC domain and D1. It is interesting to note that the onshore refraction seismic profile ZDO-001, located to the west of D1 (in the CC along the profile P3_Osorno; see location in Fig. 2), shows the inversion of an Oligo-Miocene normal fault, while the seismic profile Z5B-010A (located to the south of P4_Llanquihue profile) presents a minor contractional deformation in the CD sequences, above the D1 (Jordan et al., 2001). Deep below the CC, Maksymowicz et al. (2021) shows seismic reflectors at deep crustal levels with east-dipping angles, which is consistent with the geometry of the western border of D1 (Fig. 7) and supports a structural relation between the metamorphic complexes (WS/ES) and D1. However, the resolution of the density model is not enough to calculate the precise inclination of the western border of D1, and this structural relation is only suggested by the results.

In the northern profile P1_Toltén, D1 correlates with the high V_p and low $V_{\text{p}}/V_{\text{s}}$ anomalies obtained by Dzierma et al. (2012b). Inside the region highlighted by white contours in Fig. 4b, these authors show $V_{\text{p}}/V_{\text{s}}$ values lower than 1.74, contrasting with values higher than 1.78, eastward and westward. In same region V_p and V_s models reach values ∼ 4 % and ∼ 8 % higher than surrounding regions, respectively. Then, at least at the profile P1_Toltén, the correlation between D1 anomaly and the change in elastic properties is clearly observed. Considering this V_s velocity anomaly and an increase in density of about 0.05 g cm⁻³ (associated with the D1 anomaly, Figs. 4a and 7), we estimated an increase in shear modulus of the order of 20 % in comparison to the surrounding regions (at the same depth). To the south, this seismic velocity anomaly shows a clear continuity with the D1 geometry noticed in the 3D density model (Fig. 8b). This continuity indicates that D1 is a primary characteristic of the continental crust, southward from ∼ 39^∘ S, and this supports the interpretation of D1 as a dense–rigid zone. The latitudinal analysis of these independent geophysical models establishes that those basement units associated with D1 are progressively shifted to the east (and taken away from the trench), northward from 41.5^∘ S (Fig. 8b). In other words, the portion of the continental wedge formed by the MWU and CC domain is ∼ 50 km wider at 39.5^∘ S compared to that observed at 42^∘ S.

Outcrops of the late Paleozoic batholith near ∼ 40^∘ S are observed in the western border of WS/ES and D1 anomaly, possible implying that D1 is the southward continuation of the late Paleozoic batholith. However, outcrops of this Batholith are described at ∼ 40.3^∘ S (Deckart et al., 2014) and ∼ 42.5^∘ S (SERNAGEOMIN, 2003) near to a volcanic arc, indicating a possible association of this unit with the D2 anomaly. In this case, south of 39^∘ S, D1 should be a high-density basement unit located westward from the late Paleozoic (Pennsylvanian) batholith. An interesting candidate to fit these conditions is Chaitenia terrane (Ct in Fig. 1b), which is described as an island arc, accreted to the Gondwana margin during late Devonian times (Hervé et al., 2016, 2018). The northward limit of D1 (high V_p and Low $V_{\text{p}}/V_{\text{s}}$ anomaly, Dzierma et al., 2012b) can be roughly defined by the MVFZ. This structure could be interpreted as a limit between Chaitenia and Chilenia terrane to the north (Fig. 8b). This interpretation raises an interesting question about the role of the Chaitenia–Chilenia limit in the observed westward shift of the late Paleozoic batholith southward of 38^∘ S and its relationship with the continental deformation generated by the kinematics of LOFS (Cembrano et al., 1996; Melnick et al., 2009; Geersen et al., 2011). On the other hand, Plissard et al. (2019) observed that outcrops of mafic and ultramafic (serpentinites) rocks associated with WS (south of 39^∘ S) show P–T patterns (and structural characteristics) that allow these units to be interpreted as rock located below an incipient back-arc basin during Devonian times (380–370 Ma), which were incorporated into the subduction channel, reaching depths of about 60 km downward in the interpolate boundary (during Carboniferous time), to be finally exhumed in the eastern border of an accretionary wedge during Permian times. Under this interpretation the eastern border of the Devonian island arc (associated with Chaitenia terrane by Hervé et al., 2016, 2018) corresponds to an incipient back-arc, rather than a subduction zone, but the process finally ends in the accretion of the Devonian island arc to the Gondwana margin. Again, these accreted units could be related to the D1 anomaly in the region.

Beyond the lithology and age of D1 and D2 anomalies it is necessary to highlight the spatial association between the active volcanism and the main lineaments of the LOFS (Lara and Folguera, 2006; Sánchez et al., 2013; Díaz et al., 2020). Figure 6 shows that most of the quaternary volcanoes are located above the local regions of relatively low-density contrast (in general < 0.0 g cm⁻³ below 5 km depth) inside D2. This local 3D pattern of density anomalies is not easy to interpret in 2D regional models, because they are averaging the density structure around the profiles, and shows that the 3D local inversion is a relevant methodology (complementary to 2D regional analysis) to observe medium-depth and shallow-density structures in the upper crust. These local regions of relatively low-density contrast could respond to more fractured regions of the upper crust as a response of deep structures associated with branches of LOFS and other continental structures presents below the CC and CD. In fact, relatively low-density zones may be related to active volcanic processes observed along this fault system in the Araucanía, Los Ríos and Los Lagos districts. As shown by Díaz et al. (2020), relatively low electric resistivity is found at depths between 7 and 15 km below the local trace of the LOFS, east of Osorno volcano, associated in this case with a zone of partial melt related to a deeper ascent of basaltic magmas enhanced by the LOFS, and therefore a lower density compared to its surroundings.

The upward migration of magmas should generate local weakening zones in the overriding plate, and consequently, the continental crust in the active volcanic zone should present pervasive fracturing, fluid migration and lower density. Hence, the basement extended to the east of D1 could correspond to a similar lithology but be affected by the pervasive fracturing and fluid migration processes associated with an active volcanic arc and LOFS. This interpretation is supported by the increase in $V_{\text{p}}/V_{\text{s}}$ values to the east of D1 (Dzierma et al., 2012b), at least at the profile P1_Toltén (Fig. 4a). On the other hand, Kapinos et al. (2016) and Segovia el al. (2021) describe electrically conductive anomalies eastward from LOFS and high resistivity values beneath the CD basin (Fig. 4b and e). To the west of D1, these authors also establish conductive anomalies associated with the CC domain, reinforcing the interpretation of a transition from a highly deformed and fractured basement (related to deep units of WS/ES) to a denser/rigid basement below CD.

It is already known that the rupture propagation during large earthquakes, the interseismic deformation (including aftershocks and foreshocks), and the interplate locking are complex processes that depend primary on the frictional properties at the interplate boundary (subduction channel) and the stress field evolution (Scholz, 1998; Perfettini and Avouac, 2004; Tassara, 2010; Moreno et al., 2018; Im et al., 2020). As mentioned before, the segment of the continental wedge that includes MWU and CC domains, i.e. fractured and/or metamorphic basement units, is progressively wider to the north of 42^∘ S. This structural change correlates with the patch of high coseismic slip of the 1960 M_w 9.6 Valdivia earthquake (Fig. 8b), which added to the correlation with gravity anomaly L1 (Fig. 3a) and with changes in slope morphology, suggesting a link between the megathrust seismotectonics and physical properties of the overriding plate. In this regard, we propose that the MWU and CC domains correspond to a portion of the continental plate displaying a higher elastic and permanent deformation compared to the rigid basement landward (Chaitenia/Chilenia). Consequently, the change in the horizontal extension of this unit should modify the process of stress loading during the interseismic periods.

Due to the relatively scarce seismological data in the area and the long recurrence time for large events, it is difficult to conceptualize the complete seismotectonic story of the study zone. Nevertheless, some observations seem to support our hypothesis. Firstly, the rupture zone and aftershocks (including continental intraplate events) of the M_w 7.6 earthquake occurred in 2016 at 43.5^∘ S. This event was the largest since 1960 in the rupture area of the Valdivia earthquake (Moreno et al., 2018; Lange et al., 2018) was located at the base and within the CC domain, in the western border of a high V_p-low $V_{\text{p}}/V_{\text{s}}$ anomaly (Lange, 2008). This velocity anomaly is a clear continuation of D1 to the south of the studied area. On the other hand, historical (not instrumentally recorded) megathrusts events activated this segment of the margin in 1737 (M_w∼7.5) and 1837 (M_w∼8). They have been associated with ruptures extended to the south of ∼ 39^∘ S (Kelleher, 1974; Lomnitz, 2004), indicating that the northern portion of 1960 Valdivia earthquake could have different mechanical properties.

Lithology and internal deformation style inside and at the base of MWU and CC domains can play a different but complementary physical role in the seismotectonic segmentation of the margin. The high slip patch of Valdivia earthquake also correlates with the segment where the geometry of the marine continental wedge (seaward from shelf break) is consistent with a decrease in the effective friction coefficient (${\mathit{\mu }}_{\text{b}}^{\ast }$) at the interplate boundary (Maksymowicz et al., 2015). This suggests oversaturated fluid conditions in the subduction channel, at least in the western portion MWU at the study zone. At the same time, according to Menant et al. (2019), the deformation style of basal accretionary complexes (typically an antiformal stack of duplexes) favoured upward fluid fluxes from the interplate boundary, generating dewatering and the increase in ${\mathit{\mu }}_{\text{b}}^{\ast }$ in some adjacent regions of the subduction channel (mainly downward from basal accretionary complex). Several authors have suggested the presence of this deformation style in the deep zone of WS units (Krawczyk et al., 2006; Ramos et al., 2018; Moreno et al., 2018; Maksymowicz et al., 2021). Under this interpretation, the widening of MWU and CC domains to the north of ∼ 42^∘ S could favour high friction in the deep region of the interplate boundary (below CC domain) and a relatively low friction in the seaward portion of MWU. Therefore, the position and horizontal extension of the WS could be linked to changes in the frictional properties along the megathrust. However, more studies should be done to explore the seaward limit of WS/ES and the internal structure of MWU and CC domain.

6 Conclusions

2D and 3D density models of the forearc show a landward and latitudinal segmentation of the continental wedge in the studied zone. Offshore, the active accretionary prism limit, with a more competent basement below the middle wedge and shelf, exhibits a landward increase in density, probably associated with a progressive decrease in fracturing. To the east, the Coastal Cordillera domain presents an increase in the upper crust densities but reaches lower values than those observed in the high-density anomaly below the Central Depression. Northward from ∼ 42^∘ S, this high-density anomaly is seen progressively further from the trench, determining a northward widening of the middle wedge and Coastal Cordillera. This feature correlates with the high slip patch of the giant 1960 M_w 9.6 Valdivia earthquake.

Based on geological information, we associate the middle wedge unit (at least its eastern portion) and Coastal Cordillera domain with the late Paleozoic–early Mesozoic accretionary complex, and the high-density anomaly below the Central Depression as geophysical evidence of Chaitenia terrane. The deformation style at the eastern border of the Coastal Cordillera and seismological studies support the hypothesis of a more rigid behaviour of the continental crust below the Central Depression. Accordingly, we propose that changes in the horizontal extension of the middle wedge unit and Coastal Cordillera domain should have modified the process of stress loading during the interseismic periods, and that changes in position and extension of the late Paleozoic–early Mesozoic accretionary complex could be linked to the frictional properties of the interplate boundary.

Our results highlight the role of the overriding plate structure in the seismotectonics process in subduction zones, but more studies are necessary to understand the changes in physical properties (elasticity, temperature, among others) associated with the geological story of the margin. This work motivates similar analysis of the continental basement in other subduction margins, as in the 2010 M_w 8.8 Maule earthquake and M_w 9.0 Tohoku-Oki rupture zones.