Joint flexural-density modeling of the Taltal, Copiapó, and Iquique hotspot ridges and the surrounding oceanic plate, offshore Chile

Oceanic hotspot tracks are thought to be formed due to the interaction of a moving overlying oceanic lithosphere and a mostly stationary mantle plume (Morgan, 1971; O’Connor and Duncan, 1990). Depending on the thermal and elastic properties of the oceanic lithosphere, the resultant magmatism can be either (1) extrusive, forming seamounts in the so-called shield stage (usually associated with rough adjacent seafloor), or (2) intrusive (typically with smooth and broad seafloor; Richards et al., 2013; Orellana-Rovirosa and Richards, 2017). Large seamounts represent major loads on top of the lithosphere, which consequently deflects to compensate the effect. Similarly, intrusive magmatic bodies represent subsurface loads and mass additions that modify the structure of the lithosphere as well (Watts, 2001; Watts et al., 2010; Contreras-Reyes et al., 2019, 2021a). Thus, the study of the physical properties of the oceanic lithosphere around and beneath hotspot provinces provides important clues for better understanding not only the present-day but also the long-term evolving thermal and mechanical properties of the oceanic lithosphere.

Geodynamic models and geophysical constraints show that hotspot tracks characterized by broad regions of shallow and smooth seafloor (swell topography) preferentially occur in young oceanic plates at the time of hotspot-lithosphere interaction (Orellana-Rovirosa and Richards, 2017). Examples include the Marquesas Islands (McNutt and Bonneville, 2000) and the Nazca and Iquique Ridges (IR; Contreras-Reyes et al., 2021a, 2021b). In contrast, rough topography accompanied by relatively high and isolated extrusive hotspot volcanoes occurs in old and more rigid oceanic lithosphere (Orellana-Rovirosa and Richards, 2017), such as the Juan Fernández (Kopp et al., 2004; Manríquez et al., 2014) and Hawaii (Watts et al., 2021) hotspot tracks. In addition, Watts et al. (2010) pointed out that seamounts formed near a spreading center on weak and young oceanic plate are locally more compensated and more buoyant, being characterized by an overthickened crust. Examples are the Nazca Ridge (Hampel et al., 2004; Contreras-Reyes et al., 2021a) and Cocos Ridge (Walther, 2003; Sallarès et al., 2003). In contrast, seamounts formed on old oceanic plate are more regionally compensated and present little or no anomalously thick crust (Watts et al., 2010), such as in the case of the Juan Fernández Ridge (Kopp et al., 2004; Manríquez et al., 2014) and Canary Islands (Watts et al., 1997).

Geometrical characteristics of the oceanic plates have been modeled by flexural lithospheric deflection in response to variations in elastic thickness T_e (a proxy for the strength of the lithosphere) and vertical loading (Watts et al., 1980; Tassara et al., 2007; Kalnins and Watts, 2009; Manríquez et al., 2014; Guan et al., 2019; Contreras-Reyes et al., 2021a, 2021b). Variations in T_e can be theoretically related to variations of the lithospheric age, composition, temperature, strain rate, and strength (Burov and Diament, 1995). However, diverse studies have suggested that the variations in T_e are primarily associated with fracturing and hydration degree and less dependent on slab age (Contreras-Reyes and Osses, 2010; Craig and Copley, 2014; Manríquez et al., 2014; Hunter and Watts, 2016; Zhang et al., 2018; Garcia et al., 2019, and references therein).

In general, the physical properties of the oceanic plates are modeled with a relatively homogenous density of the crust and upper mantle and/or with constant T_e values at the regional or local scale. However, independent geophysical studies show important variations in P-wave velocity (V_p) and densities within the oceanic crust and upper oceanic mantle around the trench–outer rise region, seamounts, and oceanic fracture zones (Kopp et al., 2004; Contreras-Reyes et al., 2008, 2010; Wei et al., 2017; Cai et al., 2018; Shillington et al., 2015; Shillington, 2018; Obana et al., 2019; Watts et al., 2021; Xu et al., 2022; MacGregor et al., 2023; Dunn et al., 2024), which should impact the vertical stress balance of the flexure modeling. To tackle this issue, we propose a two-dimensional (2-D) forward modeling schema to calculate jointly the flexure of the oceanic lithosphere in response to bathymetric load and its internal density structure.

As a case study, the portion of the oceanic Nazca plate westward from the central Chilean margin (~24°S to 30°S) has been selected. This area (Fig. 1) is characterized by the presence of the Taltal and Copiapó hotspot ridges (TR and CR), which currently collide roughly orthogonally to the margin and subduct beneath the South American plate. On the other hand, the IR subducts obliquely to the northern Chilean Trench (Bello-González et al., 2018; Contreras-Reyes et al., 2021b). Thus, this geodynamic setting provides the opportunity to study the interaction between active hotspots, aseismic ridges of different size and geometry, deformation and fracturing of the crust, and the physical structure of the oceanic crust and mantle. In this work, we study the physical structure of the oceanic Nazca plate by means of gravity models with a novel 2-D coupled flexural-density modeling schema, using satellite gravity and bathymetry data complemented with local new multibeam bathymetry data. Particularly, we aim to gain insights into how magmatic material is transferred from a deep mantle plume to the overlying lithosphere and how the latter deforms beneath the volcanic edifices and/or in response to magmatic material pounded above or below the oceanic Moho.

The Nazca plate host several hotspot tracks formed by different hotspot mantle plumes (Fig. 1), including: (1) the Nazca Ridge formed at the Easter–Salas y Gómez hotpot mantle plume system (Ray et al., 2012), (2) the Juan Fernández Ridge formed at the Juan Fernández hotspot (Yáñez et al., 2001; Manríquez et al., 2014; Lara et al., 2018), and (3) the IR formed at the Foundation hotspot (Bello-González et al., 2018; Contreras-Reyes et al., 2021b). Paleo-kinematic studies suggest that the TR was formed >1500 km eastward from the East Pacific Rise on an off-ridge (hotspot mantle plume located far away from the spreading center) setting by the interaction of the San Félix hotspot and the relatively old oceanic Nazca plate (35–40 Ma; Fig. 1; Bello-González et al., 2018). A similar genesis was proposed for the CR located ~300 km southward (Fig. 1). Likely, the CR was formed by the interaction of the Caldera hotspot and the 35–40 Ma oceanic Nazca plate on an off-ridge setting (Fig. 1). Bello-González et al. (2018) already noted the parallelism between the tracks of the TR, CR, and Juan Fernández Ridge where their respective hotspots are located 900–1000 km westward of the Chilean trench and ~1700 km eastward of the East Pacific Rise (Fig. 1).

Plate kinematic models predicted that the collision point of the TR and CR with the Chilean Trench remained quasi-stationary from 40 to 30–25 Ma (Bello-González et al., 2018). A southward migration of the collision points occurred from 23 to 15 Ma (reaching a maximum velocity of 300 km/m.y. along the trench). From 15 Ma to present, the collision point has remained quasi-stationary again (Bello-González et al., 2018). The seismic structure beneath the TR and CR has not yet been studied; however, because this hotspot track was formed in an off-ridge setting, one might expect the absence of an overthickened crust. Gravimetric data can be used to study the structure of the TR and CR and test this hypothesis.

As was mentioned above, the study zone is crossed by the broad swell of the IR (Figs. 1 and 2). This long-wavelength ridge includes numerous seamounts and presents a continuous trace northward from ~25°S while to the south seeming to be bifurcated in two branches. The northern branch (IRnb in Fig. 2A) is extended to about ~30°S and the southern branch (IRsb in Fig. 2A) continues further south. Seismic and gravity models (to the north of ~25°S) show that the IR is associated with a wide crustal root that reaches ~18 km depth, which is equivalent to ~5–8 km of thickening in comparison to the surrounding crust (Contreras-Reyes et al., 2021b; Myers et al., 2022). The three hotspot ridges of the study zone (IR, TR, and CR) are developed in a region of the Nazca plate that presents crustal ages from 47 to 40 Ma (Müller et al., 2008; see Fig. 2B) with two discontinuities, which are the Mejillones fracture zone (northward from 26°S; MFZ in Fig. 2B) and the Challenger fracture zone (southward from 30°S; ChFZ in Fig. 2B).

To generate regional density and flexural models in the study zone, we use the satellite free-air gravity grid published by Sandwell and Smith (2009) and Sandwell et al. (2013, 2014), which includes the global topography grid of Smith and Sandwell (1997). Both grids have pixels of 1 min each and are available at https://topex.ucsd.edu/cgi-bin/get_data.cgi.

At lower scale, most of the high-resolution multibeam bathymetry surveys are concentrated in the trench region, north and south of the study area, and do not cover the bathymetric features associated with the TR and CR. In this work, we present primary multibeam bathymetry above two seamounts of the CR (Fig. 2). These data were acquired in February 2018 by the Chilean project AUB 170003–CONICYT/ANID using the Kongsberg EM-122 multibeam of the R/V Cabo de Hornos. The velocity of the vessel was kept between 9 and 13 km/h during the acquisition. Data processing was performed using MB-System seafloor mapping software (https://www.mbari.org/technology/mb-system/), and the final bathymetric grid was generated using Generic Mapping Tools (GMT) software (https://www.generic-mapping-tools.org/) with pixel size of 100 × 100 m.

We developed a 2-D forward algorithm to model jointly the flexure of the oceanic lithosphere in response to bathymetric load and its internal density structure. This technique allows study of the role of the crustal and mantle density structure and lateral variations of elastic thickness T_e(x) in the lithospheric flexure (i.e., as a function of the horizontal distance x along the profile). The general 2-D equation for flexural deflection w(x), under a thin-layer approximation (Turcotte and Schubert, 1982), is given by:

where the bending moment (M_b) can be expressed as:

which gives the following differential equation:

where D(x) is the flexural rigidity, which depends on T_e, Young’s modulus E, and Poisson’s ratio ν via:

F₀ is the horizontal stress, and G corresponds to a gravitational unbalanced stress (vertical load), which can be expressed as a function of w(x) and density contrast Δρ(x,z) with respect to a reference density model. For instance, in the simplified case of a one-layer crustal model with constant density ρ_crust overlying a mantle with constant density ρ_mantle and topography Topo(x), G is given by the widely used expression:

where g is the acceleration due to Earth’s gravity.

It is important to note that T_e(x) can be interpreted as an “effective” elastic thickness (McNutt et al., 1988; Burov and Diament, 1995). Considering the case of a set of k decoupled layers with different elastic thickness T_e(x)_j, the effective elastic thickness of the total set of layers is:

To include 2-D variations of the density structure, our approach considers a stack of F horizontal layers with arbitrary variable thicknesses [h(x)_{n = 1…F}] and variable density at the top and bottom of each layer [ρ_top(x)_{n = 1…F0}, ρ_bottom(x)_{n = 1…F0}], as is schematized in Figure 3. Our algorithm solves Equation 3 considering G(x) as gravitational unbalanced stress (or uncompensated load) generated by all layers at the base of the model (flat bottom of layer F at Z_max; see Fig. 3).

At any specific horizontal coordinate x, the internal density of the layer n is defined as a linear function of depth (z), i.e.: ,

then, by simple integration in z [from Z(x)_n to from Z(x)_n+1], the vertical load [L(x)_n] due to the weight of layer n is:

Finally, the gravitational unbalanced stress generated by all layers at the base of the model is:

where L_ref is the vertical load at the base of the model due to the reference density model.

The first deflected interface (which could be the Moho) is defined at the base of layer M as:

where ZMR is the non-deformed (flat) reference level for this interface (see Fig. 3). Then, according to Equations 7 and 8, the internal density gradient and the vertical load of the layer M is also a function of w(x). All layers below layer M (from M + 1 to F – 1; Fig. 3) have a deflected geometry but preserve their thicknesses, i.e.:

Finally, the base of the deepest layer is defined as a flat surface at the depth Z_max. These definitions of depths as functions of w(x) are included in Equation 9 to determine the gravitational unbalanced stress at the base of the model G(x), which in turn gives the deflection w(x) after the numerical solution of Equation 3. The algorithm solves the flexure equation, Equation 3, with finite differences in a similar form as proposed by Contreras-Reyes and Osses (2010).

The variable elastic thickness [T_e(x)] is represented by a continuous and derivable function generated by a sum of hyperbolic tangents (see the Supplemental Material¹), which allows an exact evaluation of D(x) and its derivatives along the modeled profile. Numerical tests show that this analytical representation of T_e(x) improves the stability of the finite differences solution, keeping the flexibility to model strong changes of T_e(x).

At each stage of the forward modeling procedure, the gravity response of the modeled density structure [which includes the solution of w(x)] is evaluated by the algorithm GGrad (Maksymowicz et al., 2015) and compared with real gravity data. GGrad can be directly used to calculate the gravity response of a complex 2-D layered structure. This forward model strategy allows a simultaneous modeling of the gravity signal (free-air or complete Bouguer anomalies) and crustal geometry, with independent constraints for Moho depth and bathymetry and/or elevation data. The modeling algorithm is implemented in FORTRAN90, using OCTAVE (https://octave.org/) and GMT routines to build graphic and interactive interfaces.

Theoretically, this 2-D coupled flexural-density algorithm could be used to model the flexural response of 2-D loads (i.e., bathymetric loads and density anomalies infinitely extended in the direction perpendicular to the profile). Accordingly, in the case of bathymetric elevations with three-dimensional (3-D) geometries, as is the case of individual seamounts, the flexural effect of the bathymetric load is overestimated. However, the 2-D modeling can precisely represent the structure of the regional bathymetric and gravity anomalies elongated in the direction perpendicular to the profiles. As will be discussed in the following section, the 2-D modeling algorithm was applied to bathymetry and gravity profiles generated by the average of bathymetry and gravity data in a thick band around the profile tracks.

The parameter D in Equation 3 encompasses all elastic properties of the media into the flexural analysis and formally depends on T_e, E, and ν (Equation 4). However, in order to provide a more intuitive parameter for interpretation, we model the flexure varying the elastic parameter T_e and keeping constant E and ν at the values of 7 × 10¹⁰ Pa and 0.25, respectively, which has been previously used for flexural model along the Chilean margin (Contreras-Reyes and Osses, 2010; Contreras-Reyes et al., 2021b).

The reference (undeformed) models of the oceanic crust for flexural modeling were based on available seismic tomography models along active wide-angle refraction seismic profiles located in the northern Chilean margin (at 22°S and 23.5°S; Contreras-Reyes et al., 2012; Sallarès and Ranero, 2005; see locations in Fig. 2A). These V_p models were converted to densities through the empirical Nafe-Drake transformation curve (Brocher, 2005). Table 1 summarizes the one-dimensional density structure considered to generate the reference models for flexural modeling.

Considering the general seafloor elevation of the profiles (far from seamount features), the Moho discontinuity (bottom of layer 3, Table 1) has a reference depth of ~11.5 km (equivalent to ZMR value in Equations 10 and 11; see also Fig. 2). On the other hand, the base of the model (flat bottom of layer 4, Table 1; Z_max in Fig. 2) is defined in all profiles at 100 km depth with a density of 3.33 g/cm³; no density variation is modeled below 100 km depth. This value (3.33 g/cm³) is commonly assumed as a reference average of the upper mantle density (e.g., Bott, 1991; Hunter and Watts, 2016). According to this definition, we start with a density model with mantle densities of 3.24 g/cm³ just below the Moho (based on the P-wave velocities of seismic topographies at 22°S and 23.5°S) increasing linearly with depth to 3.33 g/cm³ at 100 km. This smooth reference mantle density gradient in depth (Table 1) and the modeled anomalies (introduced in section 4.3) determine small density variations below 30–35 km depth.

The processed swath bathymetry track acquired by the Chilean project AUB 170003 aboard the Chilean vessel R/V Cabo de Hornos shows the first-ever detailed view of two seamounts along the CR (Figs. 2A and 4). These seamounts, named here as Copayapu I and Copayapu II, are immediately seaward of the Chilean trench. As is expected, the horizontal width of bathymetry track is progressively thinner in the zone of the highest elevations, corresponding in this case to the summit of Copayapu II seamount. According to processed data, Copayapu II has a maximum elevation of ~−1015 m and presents a stellate shape with prominent radial flank rift zones, some of them extending ~15 km. Its upper section is coronated by a partially empty crater and summit scarps (Fig. 4). At least two parasitic flank cones are observed on its western flank. To the east, Copayapu I seamount (~2000 m height) has a lower basal diameter (~20 km) and exhibits a more regular profile with less-developed flank rift zones. Its rounded summit reaches a lower elevation compared to Copayapu II, although considering the bending of the Nazca plate in the outer-rise zone, both seamounts would have similar elevation above seafloor. In fact, analyzing a bathymetric profile through the seamounts (Fig. 4A) and using the seafloor around these features to generate an approximate geometry of the bent oceanic plate by linear interpolation (dotted red line at the upper panel in Fig. 4C), we observed summits >2000 m above the seafloor in both seamounts (2606 m in Copayapu II and 2145 m in Copayapu I; see lower panel in Fig. 4C).

It is worth noting that near the eastern limit of the profile, a bathymetric high of >500 m above the surrounding seafloor is observed. This feature corresponds to a horst structure associated with the outer-rise deformation in the trench zone, which is representative of the deformation style and roughness of the Nazca plate observed near the trench axis along the study zone (Fig. 2A). These data (complemented by lower-resolution bathymetry; Figs. 1 and 2) show that both seamounts have a general conical symmetrical shape, with eroded flanks characterized by radial incisions and a rough surface radial-chutes-and-ridges morphology, which is expected for relatively old seamounts formed far from weak oceanic structures such as fracture zones (Chaytor et al., 2007). This morphology is also observed in some seamounts along the Juan Fernández Ridge to the south (e.g., Lara et al., 2018) although without a clear linear relation between erosive patterns and age.

To highlight the short-wavelength morphology of the seamounts and the structures of the oceanic plate, Figure 4B shows the slope gradient map of the observed bathymetry (formally | (elevation)|). The seafloor shows a clear deformation pattern where the outer-rise flexure is associated with structures with near north-south direction. These outer-rise structures reach larger vertical offsets at the horst-and-graben sequence observed near the trench, which is also observed as a pervasive deformation on the Copayapu I seamount. A second set of large structures is observed with northwest direction, probably associated with the original fabric of the oceanic plate at this latitude as is evidenced by iso-contours of Nazca plate age (Fig. 2B) and magnetic anomalies (see the Supplemental Material [^{footnote 1}]). Similar features have been observed by Ranero et al. (2005) and Poli et al., (2017) northward and southward of the study region. One of these northwest lineaments seems to limit the northeastern flank of the Copayapu I seamount, suggesting that this original oceanic fabric could have been reactivated during the outer-rise flexure process, which is also suggested by the lineament (or scarp) observed in the southwestern flank of Copayapu II seamount (Fig. 4B). The deformation pattern observed with these new bathymetric data confirms activity of the outer-rise bending faults and suggests a possible reactivation of the primary seafloor structures along this portion of the margin, as previously observed to the north and south (Ranero et al., 2005; Poli et al., 2017; Geersen et al., 2018a, 2018b). These reactivated structures could favor the hydration of the oceanic crust and upper mantle, as was proposed by several authors for this and other convergent margins worldwide (Ranero et al., 2005; Contreras-Reyes et al., 2008; Contreras-Reyes and Osses, 2010; Shillington et al., 2015, 2018; Kita and Ferrand, 2018; Myers et al., 2022).

In order to study the first-order relation between oceanic elastic thicknesses, bathymetry, and gravity signal, we generated a simplified set of 2-D joint flexural-density models along nine gravity profiles (Fig. 2). These simple models were obtained using constant T_e values ranging from 1 to 20 km (specifically 1, 5, 10, 15, and 20 km) and keeping constant the reference density structure below the bathymetry (Table 1), except for the two lower layers which depend on the Moho depths (i.e., these layers depend on the calculated flexure; see section 3.2). Figure 5 shows two examples of this exercise for profiles L1a and L5, which are representative for the areas near the trench and far from the subduction zone, respectively (the other profiles can be seen in the Supplemental Material).

As was mentioned above (section 3.2), the 2-D joint flexural-density algorithm was applied to bathymetry and gravity profiles generated by the averaging (perpendicular to the profile) of bathymetry and gravity data in a band around the profile tracks. In this case, we use a 60-km-thick band centered on the profile tracks shown in Figure 2. This procedure aims to avoid the overestimation of the flexural effect of the loading associated with 3-D bathymetric features (such as individual seamounts), highlighting the regional bathymetric and gravity trends elongated in the direction perpendicular to the profiles. Thereby, in all 2-D models presented in Figures 5–8 (and in the Supplemental Material), the bathymetry and gravity curves correspond to these nine averaged profiles. With this procedure, the obtained models across the axis of the ridges should be considered as an approximate representation of the density and T_e structure.

Observing the simple models obtained for constant T_e values, the general shape and amplitude of the gravity anomalies cannot be explained exclusively by a pattern of T_e values. In some cases, the gravity signal is better explained by high T_e values (>5 km), as in the northern portion of profile L1a (Fig. 5). However, in other regions, the low T_e values seem to generate more realistic signals (as in profile L5, Fig 5). In general, near the TR and CR, the shape of the gravity signal can be approximated by these models, but not the amplitude of the gravity signal. On the other hand, far from the ridge traces, several gravity anomalies are completely misfitted (see, e.g., the middle and southern portion profile L1a). In a global analysis, the best root-mean-square (RMS) values are obtained in general for high T_e values (usually ≥5 km) in the eastern profiles (L1 to L3a) and for lower T_e values (usually ≤5 km) in the western profiles (L4 to L5). The latter suggests a pattern of westward decrease of T_e, but the associated RMS values are high (between ~5 to ~12 mGal), which reflect a large misfit of the main anomalies, even considering the highest nominal data error reported for the gravity database (±3 mGal; Sandwell et al., 2013).

Considering these first results, we generate a second set of models with variable T_e curves along the profiles to study the effect of T_e variation on the gravity anomalies. After testing numerous T_e curves, we reach models that reduce the original RMS values to between 4% and 45% (16% on average, from RMS values ~7.6 mGal to values ~6.4 mGal). Figure 6 shows four examples of the obtained models (along profiles L1a, L3, L4, and L5; see other profiles in the Supplemental Material). This better fit is obtained in general with reductions of T_e values around the axes of the TR and CR; however, broad sections of the profiles remain poorly represented by the density models. This is clear in profile L1a (Fig. 6) and in the other eastern profiles (L1, L2, and L2a; see the Supplemental Material), where even under extreme variations in T_e, the gravity signal is not well represented, suggesting large density variations along the profiles. Westward, the gravity signals of the profiles L3 to L4a (Fig. 6; Supplemental Material) are better explained by models with variable T_e (in comparison to models of lines L1 to L2a), but several local gravity anomalies differ in amplitude and shape with the modeled gravity. In particular, the model of profile L5 (Fig. 6) shows that a significant decrease in T_e around the active San Félix–San Ambrosio hotspot leads to a huge crustal root below this feature (~10 km thicker than the surrounding crust), but the shape of this large gravity anomaly is not well represented. The latter implies a large (and probably unrealistic) amount of extruded material in this active hotspot to explain a large deflection of the crust, which in turn suggests that a variable density structure is needed to explain the gravity signal associated with the San Félix–San Ambrosio islands.

Starting with the models presented in the previous subsection, we construct a set of 2-D coupled flexural-density models that includes density variations to fit, as much as possible, the observed average gravity data along the nine profiles. Figure 7 shows examples of these density models and the corresponding gravity fits along profiles L1a, L3, L4, and L5 (see the other profiles in the Supplemental Material). We observe that all these models explain the data with RMS values <1 mGal (most of them ~0.4 mGal) and the shape and amplitude of every gravity anomaly is now properly adjusted.

Due to the scarcity of independent constraints for the Moho geometry and crustal structure in the area, we do not define the densities of the specific geological or lithological units, but the obtained density variations should be indicative of the regional density variations and/or thickness changes of the geological or lithological structure. Then, the resulting models can be seen as a family of possible solutions, which should be updated according to future geophysical studies. However, some general characteristics emerged naturally during the forward modeling process. The short-wavelength gravity anomalies are explained mainly by density variations in the upper portion of the crust, while the long-wavelength anomalies are modeled by variations in the lower crust and upper mantle density structure. On the other hand, there is an obvious trade-off between crustal (or shallow) and upper mantle density (or deep density structure), which implies that for modeling a given gravity anomaly it is possible to include decrease of the crustal density with an increment of the mantle density or vice versa. In addition, the resulting Moho depth (amplitude of the deflection) responds to the general density of the crust and upper mantle, where an increase of the crustal and/or upper mantle density increases the deflection of the crust, and vice versa.

These trade-offs imply that the fitting of a low gravity anomaly (for instance) is possible with a decrease of the crustal density, the upper mantle density, or both. However, this decrease in the density tends to decrease the deflection (shallowing the Moho discontinuity), which in turn tends to increase the modeled gravity; i.e., the reduction of the density could have a contrary effect in the gravity model than initially expected. Thus, a low gravity anomaly could be obtained balancing the decrease of crustal densities with an increase in the upper mantle density and/or adding a relatively high value of T_e to minimize the deflection of the model.

An example of this complex trade-off is observed at profile L1a (Fig. 6), where to account for the large reduction of the regional gravity to the south of the profile, it is necessary to include a decrease in the upper crustal densities and a low-density layer below the Moho discontinuity. To reduce the upward flexure generated by the general reduction of the densities in this profile, a high T_e value must be employed. Otherwise, the reduction in density leads to an unrealistically thin crust (<3–2 km), which in turn increases the modeled gravity, making it impossible to fit the observed data. This is a clear example of a type of global physical constraint that is obtained by the joint flexural-density modeling approach.

Figure 8 summarizes the density and T_e distribution obtained along the nine profiles by the 2-D coupled flexural-density modeling with density variation. In broad segments of the eastern profiles (L1, L1a, L2, and L2a), the models show a layer with low density below the Moho discontinuity. This layer (UL in Fig. 8) is necessary to fit the amplitude and shape of the regional gravity anomalies around the TR and CR. We assume a constant density of 3.1 g/cm³ for this layer (lower than the 3.24 g/cm³ of the reference model; see Table 1), which is an arbitrary value but consistent with the low V_p values observed below the Juan Fernández Ridge (Kopp et al., 2004) and with those of the underplated material modeled to explain the swell seafloor bathymetry associated with the Nazca Ridge (Contreras-Reyes et al., 2019, 2021a).

In general, the upper crust in these profiles (L1 to L2a) shows lower densities around TR and CR (and to the north and south of these ridges) in comparison to the segments between TR and CR. Westward, profiles L3, L3a, L4, and L4a systematically show a thicker crust between the TR and CR and to the south but lower densities in the upper crust mainly around the TR and CR and to the north and south of these ridges. These thick crustal zones are correlated with the presence of the IR, which intersects the TR and CR at different latitudes (see Fig. 2). In contrast to the models for the eastern profiles, these models do not include the UL unit, and the deflection due to topography load is enough to explain the main features of the gravity anomaly. This is also observed in the models with a simple density structure (Fig. 6; Supplemental Material) where the models for profiles L3 to L4a have lower RMS values than the other models.

Finally, the model along profile L5 (Figs. 7 and 8) shows that the inclusion of a UL unit below the San Félix–San Ambrosio islands explains the shape of the gravity anomaly of this active hotspot. Because of the presence of the upper mantle unit and the large decrease in T_e below the TR axis, the Moho discontinuity is slightly upward deflected, showing a small crustal root below the San Félix–San Ambrosio islands. In contrast, the CR only shows a small crustal root in response to the local bathymetric load associated with this ridge.

The results of 2-D joint flexural-density models for a simple density structure (Fig. 6; Supplemental Material [see ^{footnote 1}]) are similar to those expected for 2-D flexural models with constant crustal and mantle densities (e.g., Bodine et al., 1981; Wessel, 1993; Contreras-Reyes et al., 2021b). In these models with a simple density structure, the vertical load generated by bathymetric highs is compensated by crustal roots with different wavelengths depending on the considered T_e. Worldwide, the general shape of the gravity signal associated with oceanic ridges is roughly explained by this kind of simple model, but clearly, there are important features of the gravity anomalies that cannot be modeled with this simple approximation (Hwang and Kim, 2016; Watts et al., 2021). In our results, most of the small gravity features and some regional trends are not correctly explained by the simple density models, even considering a variable-T_e model along the profiles. Examples of this are the smooth gravity high of >20 mGal, observed between the TR and CR in profile L1a (between horizontal distances of ~250 km and ~350 km), the regional gravity decrease observed southward the CR position in the same profile, and the wide signal associated with the TR at profile L5 (Fig. 6). This observation highlights the relevance of density inhomogeneities of the oceanic crust and upper mantle, which is somehow expected, considering the variability of the seismic velocity structure observed around seamounts, fractures zones, spreading centers, and other oceanic features in diverse regions (e.g., Caress et al., 1995; Tsuji et al., 2007; Bai et al., 2019). On the other hand, the density inhomogeneities control changes in the vertical load on the lithosphere, which must be considered to model its variable elastic properties by the application of general nonlinear flexure (Equation 3).

During the forward modeling process, we note that the gravity models have moderate sensitivity to the amplitude of T_e variations. Hence, we can include increases and decreases in T_e along the profiles to improve the data fit, but several T_e curves lead to a similar RMS value. In those cases, we select T_e curves favoring models with low RMS values and smooth lateral variations of the Moho depth and crustal and/or mantle density and discard models with strong thickening and thinning in zones with little change in bathymetric loads. It is important to remark that different T_e curves in the same profile generate models with similar density and geometric structures (see examples in Fig. S6.1 in the Supplemental Material). This behavior indicates that the inclusion of a variable density structure along the profiles is a key strategy to explain the gravity signal, at least in the study zone. Observing the simple models (Fig. 5; Supplemental Material), we note that T_e variations determine important changes in the Moho geometry and gravity response, but several regional gravity features remain unexplained without a more complex density structure. In fact, the large improvement in data fit from constant-T_e models to final models is mostly due to the inclusion of density variation (Figs. 5–7; Supplemental Material). This phenomenon could be characteristic of the study zone and not characteristic of zones where large seamounts and others bathymetric features control the gravity response at the regional scale. The uncertainty of the models leaves the final T_e models to be considered as representative of possible solutions that could be modified after the availability of new independent geophysical constraints.

To complement the discussion, Figure 9 summarizes the results obtained in the nine profiles (Fig. 8) by the interpolation of the Moho depth (Fig. 9A), average crustal density (from the seafloor to the Moho; Fig. 9B), base of the modeled UL unit (Fig. 9C), and the average of the upper mantle density from the Moho to 30 km depth (Fig. 9D). The regions with high T_e (>15 km) and low T_e (<5 km) are highlighted with red and blue dotted lines, respectively (according to T_e models along the profiles).

In a regional view (Fig. 9), most of the high-T_e regions are located to the east, which was also suggested by our constant-T_e models (Fig. 5). This regional pattern can be associated with a trenchward age increase of the Nazca plate (see Fig. 2B), although at local scale the models show a more complex behavior characterized by a T_e reduction beneath the TR and CR. Reduction in T_e above seamounts has been related to lithospheric reheating in hotspots and/or by mechanical weakening due to fracturing and loading processes (Watts and Ribe, 1984; Wessel, 1993; Lee et al., 2009). While thermal reheating (rejuvenation of the lithosphere) is a possible explanation for the T_e decrease above the San Félix–San Ambrosio active hotspot (profile L5), local T_e decreases in other sectors of the TR and CR could be associated with pervasive fracturing and lithospheric deformation near the ridge axes (probably increased by tectonic reactivation in the outer-rise zone; Fig. 4).

As is observed in Figure 9A, the wide IR is associated with zones of deep Moho discontinuity (and high crustal thickness), while the TR and CR are more related to relatively small and local crustal roots (see also Fig. 8). This variability could be associated with the formation of these oceanic ridges at different hotspots. The IR is formed at the Foundation hotspot (Bello-González et al., 2018), located in a thin, hot, and weak lithosphere near a spreading center (East Pacific Rise). In contrast, the TR and CR are formed at hotspots located at relatively old, cold (and rigid) zones of the Nazca plate (Figs. 1 and 2), where the mantle plume (and the associated extruded material) should weaken the overlying thick and old oceanic lithosphere. These models for the TR and CR are consistent with previous observations at the Juan Fernández Ridge, where active seismic studies have shown an almost undeflected Moho discontinuity below the O’Higgins Seamount Group (Kopp et al., 2004). These similarities in the Moho geometry suggests a similar seamount formation process at the San Félix–San Ambrosio, Caldera, and Juan Fernández mantle plumes.

Analyzing the density structure, we observed a complex distribution of the average crustal density with zones of high densities mainly located to the south of ~25°S (Fig. 9B). This change can be interpreted as a latitudinal segmentation of the oceanic crust delineated by the presence of the TR, and may be related to the age transition associated with the Mejillones fracture zone. As is also observed in Figure 8, the TR and CR tracks are characterized (in general) by local decreases in crustal density, which is mainly due to the thickening of the zone with densities <2.8 g/cm³ in the upper crust. These local decreases in average crustal density may be related to the presence of extruded material and/or the increase of fracturing due to the deformation and flexure of the ridges. Similar crustal structure is observed in several aseismic ridges worldwide (e.g., Caress et al., 1995; Kopp et al., 2004; Nishizawa et al., 2009; Stratford et al., 2015).

As was described above (section 4.3), in broad segments of the eastern profiles L1, L1a, L2, and L2a, the presence of a low-density upper mantle unit (UL) is necessary to explain the gravity anomalies. In Figure 9C, we show the base of the UL unit (which is equal to Moho depth in zones without UL). Compared to the Moho depth (Fig. 9A), UL shows a local thickening below the TR and CR traces and to the north and the south of these features (see also Fig. 8). Farther to west, profile L5 shows a local thickening of UL below the San Félix–San Ambrosio islands. Due to the inherent trade-off between density and size of the anomalies, the thickness of the modeled UL is controlled by its density (3.1 g/cm³, based on the V_p values observed by Kopp et al. [2004] below the Juan Fernández Ridge), but its general shape (location of high- and low-thickness zones) remains unaltered under variation of its density.

The UL layer with the low V_p values observed by Kopp et al. (2004) below the O’Higgins Seamount Group (in the Juan Fernández Ridge) could be interpreted as a hydrated zone of the upper mantle. These authors related the mantle hydration process with the fracturing of the crust reactivated in the context of the outer-rise zone, which seems to be also characteristic of the area, at least for the CR (Fig. 4). In comparison with our study zone, the outer-rise fracturing seems to be more pervasive (and with a local northeastern strike) around the Juan Fernández Ridge (Poli et al., 2017), but the broad extension of the UL to the south of the CR in profiles L1 and L1a and to the north of the TR along profiles L2 and L2a (Fig. 8) supports an association with hydration processes linked to outer-rise tectonics far from the ridge axis. A similar characteristic was observed by Kopp et al. (2004) beneath the Juan Fernández Ridge; however, in all our eastern profiles, UL is thicker below the ridge axis, which also could suggest the presence of underplated material below the seamounts and/or a deeper hydration process favored by the load and flexure around the seamounts. Beyond the presence of possible underplated and/or hydrated material, the average density of the crustal mantle (Fig. 9D) shows a smooth eastward decrease, which can be interpreted as a general increase of fracturing and mantle hydration related to the outer-rise tectonics. On the other hand, the thick UL unit observed below the San Félix–San Ambrosio islands is consistent with a process of thermal reheating in the active hotspot, as also is suggested by the low T_e values modeled in the area.

The gravity anomalies associated with the IR are well explained by crustal roots in absence of a UL unit. Therefore, hydrated upper mantle and underplated material are not explicitly modeled beneath this ridge. The obtained Moho deepening below the IR is consistent with the previous flexural and gravity models of the IR presented by Contreras-Reyes et al. (2021b) and with the unique wide-angle seismic profile that crosses the IR (Myers et al., 2022) to the north of the study zone. Particularly, Myers et al. (2022) showed that the upper mantle below the IR does not present a reduction in V_p, which is consistent with our model beneath the IR. In relation to the Moho geometry, we observed some variability in the interaction between the large IR crustal roots and the TR and CR ridges (Figs. 8 and 9). While a thick crustal root is observed in the intersection between the IR and TR at profile L3, the IR root sems to be unperturbed by the TR along profile L3a. Similarly, the CR does not show an evident root in the zone of the southern branch of the IR (profile L4a). On the other hand, profile L4 shows crustal roots below the TR and CR immediately to the north and south of the IR. This variability suggests that the thick crust associated with the IR (and its southern branches) does not favor later extrusion of volcanic material in the hotspot nor subsequent local flexure.

One of the most salient features observed in this study is a clear variability of the oceanic plate structure along the ridges and around them. In the western portion of the study zone (profile L5), the San Félix–San Ambrosio hotspot (TR hotspot) seems to be much more active than the Caldera hotspot (CR hotspot), forming a larger bathymetric high, probably associated with a significant thermal reheating and underplating process. In other places, the crustal roots of the CR are larger, having the same size as observed in the TR (e.g., profile L4), and in the eastern profiles (profiles L1 to L2a) the presence of hydrated upper mantle and/or underplated material seems to be associated with the CR and/or TR with variable thicknesses. While fracturing, fault hydration and fault reactivation can be associated with changes in tectonic setting far from the hotspots (e.g., in the context of the outer-rise zone), the magnitude of underplated units and extruded material should reflect mostly the activity of the hotspot. Then, the variable magnitude of UL unit and Moho deflection observed in the models can be interpreted as evidence of the variation in time of the hotspot activity, favoring different rates of underplating and extrusion. This genetic variability should generate significant changes in the subsequent processes of crustal fracturing and mantle hydration, determining a mechanical segmentation of the oceanic plate along the subduction margin. Accordingly, the results show the importance of the acquisition of more geophysical and geological data along the ridges to reconstruct the temporal evolution of the activity of the San Félix–San Ambrosio and Caldera hotspots.

In the study area, the Nazca–South America subduction zone shows segmentation in the large-earthquake sequence, in interplate coupling, and in the tectonics and morphology of the continental wedge, which has been related to the subduction of the bathymetric highs of the ridges and the frictional properties in the long term (geological time) and short term (Contreras-Reyes and Carrizo, 2011; Maksymowicz, 2015; Molina et al., 2021; Cubas et al., 2022). Numerous authors have shown that the hydration of the interplate boundary and upper mantle is a key factor to explain the effective basal friction coefficient, coupling, interpolate and intraplate seismicity, and slow-slip events (Dahlen, 1984; Chesley et al., 2021; Wang and Lin, 2022, and references therein). Then, the fracturing type, plate hydration, and underplating variability observed seaward from the trench should be expected for the subducted portion of the oceanic plate, controlling (at least in part) the seismotectonic inhomogeneities of the subduction zone. Figure 9 shows clusters of seismicity (U.S. Geological Survey, 2023) at the latitude of the TR and CR subduction but also high seismicity to south of ~27°S. The same features are observed in the detailed seismological study of González-Vidal et al. (2023), who also showed a decrease in the current locking degree around and to the south of the CR. This pattern could be related to the pervasive fracturing observed in the CR seaward from the trench (Fig. 4) and to the broad hydrated and/or underplated unit (UL) associated with the CR along profiles L1 and L1a.

In 1922, the southern portion of the margin (~26.5°S to 30°S) ruptured during a large (~M_s 8.5–8.8) megathrust earthquake. Nonetheless, historical earthquakes of only moderate magnitude (~M_s 8) ruptured to the north and south of ~26°S and between the latitudes of the TR and the CR trench collision zones (Ruiz and Madariaga, 2018; Kanamori et al., 2019). In addition, a slow-slip event has been detected below the onshore forearc (~27°S–27.5°S; Klein et al., 2018), which has been related to a low-V_p, low- V_s (S-wave velocity), and relatively low V_p/V_s anomaly interpreted as a subducted bathymetric feature of the CR (Pastén-Araya et al., 2022). All these seismological features show a mechanical inhomogeneity of the interplate boundary, which might be a result of the interaction between the inhomogeneous oceanic plate (observed in this study seaward from the trench) and the temperature and fluid conditions of the upper plate in the subduction zone. In this sense, our results show that the oceanic plate has significant physical heterogeneity in density, fracturing, deformation, and underplating and/or mantle hydration, which are not exclusively confined to the aseismic ridge. This variability should be considered in explaining the seismotectonic properties of subduction margins.

It is worth noting that the poor data fit associated with simple models, with constant and variable T_e (Figs. 5 and 6; Supplemental Material), shows that the inhomogeneities of the density structure are necessary to explain the gravity signal in detail. For instance, bathymetry and flexure of a homogenous plate can explain the rough trends and shape of the gravity signal but are not enough to account for numerous gravity anomalies that can be associated with, for example, variations in fracturing or mantle hydration. In our view, this opens a new research line in flexure modeling to explore the oceanic plate structure. According to this, our final models are possible solutions for the gravity signal of the study zone, but several aspects of the models could be modified after the availability of new independent constraint of the Nazca plate structure. Trade-offs are present in the model, mainly related to the geometry of the Moho discontinuity and the deep extension of the UL unit. However, the scarce seismic studies (in relation with the margin length) and the smooth north-south variation of the Nazca plate around the study area (Kopp et al., 2004; Sallarès and Ranero, 2005; Müller et al., 2008; Contreras-Reyes et al., 2012, 2014) indicate that the Moho discontinuity should not present significant variations outside of the ridge axis (i.e., in the “normal” oceanic plate). This confirms the presence of the modeled low-density upper mantle units, because otherwise, unrealistic and rapid Moho deepening (even far from high bathymetric loads) are needed to explain the regional decreases in gravity.

In a global view, our 2-D joint flexural-density approach shows that the effect of density variations on lithospheric flexure cannot be avoided in a detailed modeling of plate structure. Also, regional and local variation in T_e is expected, and therefore, its derivative should be properly included in the modeling process. In this sense, our result motivates a more extensive application of this methodology to extend and refine the knowledge of the oceanic lithosphere in different geodynamical contexts. Future methodology advances along this line should be directed to include the outer-rise flexure effect of trench bending (similarly to Contreras-Reyes and Osses, 2010; Manríquez et al., 2014; Contreras-Reyes et al., 2021a, 2021b) in a schema of 2-D joint flexural-density modeling and the development of a 3-D version of the algorithm to study the 3-D structure of the oceanic plate at the scale of individual seamounts.

Our analysis reveals a complex physical structure of the Nazca plate in the study zone, related to the diverse genetic process of the Taltal, Copiapó, and Iquique hotspot ridges. The Iquique Ridge, formed in young oceanic lithosphere, is characterized by a wide and deep crustal root, while the Taltal and Copiapó Ridges, both formed in old oceanic lithosphere, show relatively small and local crustal roots generated during the weakening process (lithospheric reheating) in the hotspot mantle plumes.

The obtained 2-D joint flexural-density models show evidence of this weakening process and its temporal and spatial variability. Currently, the San Félix–San Ambrosio hotspot is characterized by high activity reflected in the shallowing of the modeled Moho discontinuity accompanied by a thick underplated unit below. On the other hand, the zone near the Caldera hotspot shows a small perturbation of the oceanic lithosphere, interpreted as evidence of current low hotspot activity. Trenchward, the Taltal and Copiapó Ridges (as well as the Iquique Ridge) present high variability in crustal root size, probably related to temporal variability in the activity of the hotspots.

The results show variability of crustal density at different scales. High average crustal densities are observed to the south of ~25°S, and at small scale, reductions in crustal densities are observed around the ridge tracks. These density anomalies may be related to variations in the oceanic plate fabric; age; and fracturing, flexure, and tectonic activity along aseismic ridges and the outer-rise zone. At the eastern limit of the Copiapó Ridge, near the trench, the high-resolution swath bathymetry around Copayapu I and Copayapu II seamounts confirms the activity of the trench outer-rise bending faults and suggests a possible reactivation of the oceanic fabric structures along this portion of the margin. This pervasive fracturing is consistent with the presence of broad low-density zones below the Moho, interpreted as hydrated zones of the upper mantle, which also could be reflecting the variability of underplated material below the axis of the ridges.

The observed physical inhomogeneities of the oceanic plate Moho depth, T_e, underplating, fracturing, and mantle hydration determine a spatial segmentation that should be considered to explain the mechanical behavior of the subducted oceanic plate, including the interplate seismotectonic features of this and other subduction margins.

This work shows that the modeling of the flexural effect of density inhomogeneities is necessary to explain the gravity signal of the oceanic plate. The coupled modeling of complex density structure and variable elastic thickness and its derivatives should be properly considered to obtain a more detailed view of oceanic plate structure. Our results motivate future methodological developments in 2-D and 3-D joint flexural-density models and a more complete coverage of active seismic and swath bathymetry in oceanic plates to study more deeply their mechanical inhomogeneities.

This work was funded by the Agencia Nacional de Investigación y Desarrollo (ANID) under the program Fondo Nacional de Desarrollo Científico y Tecnológico (FONDECYT), grant 1210101. Maksymowicz and Lara acknowledge the support of Millennium Institute on Volcanic Risk Research–Ckelar Volcanoes ANID/MILENIO grant ICN2021_038. Discussion undertaken in this contribution is also tributary to FONDECYT grant 1211792. Our work benefited greatly from thoughtful and probing reviews from Science Editor Andrea Hampel, from reviewer Jörg Ebbing, and from an anonymous reviewer.