Climatic control on the location of the Southern Andes volcanic arc

Volcanic arcs at convergent plate margins are primary surface expressions of plate tectonics. Although climate affects many of the manifestations of plate tectonics via erosion, the upwelling of magmas and location of volcanic arcs are considered insensitive to climate. In the Southern Andes, subduction of the Nazca oceanic plate below the South American continent generates the Southern Andes Volcanic zone. Orographic interactions with Pacific westerlies lead to high precipitation and erosion on the western slopes of the belt between 42-46°S. At these latitudes, the topographic water divide and the volcanic arc are respectively farther and closer to the subduction trench than at lower latitudes, despite a constant subduction dip angle along strike. Here, we use thermomechanical numerical modeling to investigate how magma upwelling is affected by topographic changes due to orography. We show that a leeward topographic shift may entail a windward asymmetry of crustal structures accommodating the magma upwelling, consistent with the observed trench-ward migration of the Southern Andes Volcanic Zone. A climatic control on the location of volcanic arcs via orography and erosion is thus revealed.


Southern Andes volcanic arc and climate
Dehydration of a subducting oceanic plate commonly occurs at ~100 km depth and leads to partial melting of mantle rocks 1,2 . Upwelling buoyant magmas then feed the volcanic arc on the upper continental plate [1][2][3][4][5] . The slab dip angle, thus, exerts a primary control on the distance between the volcanic arc and the subduction trench 2,6 , while crustal and lithospheric structures affect the magma ascent toward the surface [3][4][5] . The topography generated by compressional strain and magma upwelling during ocean-continent subduction is subject to climate-driven erosion [7][8][9][10][11] . The redistribution of the surface masses by erosion, in turn, affects the lithospheric deformation, partial rock melting, and magma transfer, thereby influencing the structural and magmatic evolution of a compressional orogen [7][8][9][10][11][12][13] . If the topography acts as an orographic barrier, in particular, enhanced precipitation and erosion on the upwind side of the orogen force a leeward migration of the topographic range 9 . While the role of orography and associated erosion gradients on the structure and topography of an orogen was addressed before 7-10 , its forcing on the location of a volcanic arc has been overlooked.
In the Southern Andes between ~34-46°S, subduction of the Nazca oceanic plate beneath South America bounded to the south by the northward migrating Chile Triple Junction (CTJ) generates the Southern Andes Volcanic Zone (SVZ) and an orogenic wedge with mean elevation between 1-3 km above sea level [14][15][16][17][18][19][20] (Fig. 1a). Throughout the Plio-Quaternary, steady oblique subduction along strike occurs at a rate of ~7 cm/yr (Refs. [21][22][23] with slab dip angle of ~25° (Refs. [24][25][26]. Three latitudinal zones account for statistically significant changes in the relative distances between the subduction trench, the topographic water divide, and Quaternary volcanoes ( Fig. 1, Supplementary Table 1, Methods). The average distance between each volcano and the subduction trench decreases from ~300 km north of 42°S (zones 1 and 2), where volcanoes are located above the 100 and 120 km slab isodepths, to ~260 km south of ~42°S (zone 3), where volcanoes are located above the 80 and 100 km slab isodepths (Fig. 1a,b). This ~40 km trenchward migration of the volcanic arc between zones 1-2 and zone 3 cannot be explained by a southward increase in the slab dip angle [24][25][26] . Instead, south of ~40°S the amount of precipitation increases on the western slopes due to orographic interactions with westerly winds 27 . The increase in erosion rates on the upwind side of the range 28,29 leads to a ~70 km eastward migration of the topographic water divide (Fig. 1c) and contributes to a ~1000 m decrease in the mean altitude 15,18,20 of the orogenic wedge from zone 1 to 3. Volcanoes are found west and east of the topographic water divide north of ~40°S, but they are systematically west of the water divide south of ~40°S (Fig. 1a,d). Quaternary erosion rates 29 . Convergence plates velocities 21 . Depth of the top of the Nazca slab 25,26 . Late Miocene arc front adapted from 19 . b-d, Scatter and box plots of the relative distances between the volcanoes, the subduction trench and the topographic water divide (Supplementary Table 1). Negative values indicate volcanoes to the west of the water divide. Whiskers (w) extend to the most extreme data points (without considering outliers, plotted with '+' (Supplementary Table 1), and correspond to ~2.7σ and 99.3% coverage of the data. Outliers, defined as values greater than 3 + w( 3 − 1 ) or smaller than 1 − w( 3 − 1 ) , where q1 and q3 are the first and third quartiles respectively, correspond to within-plate volcanoes 19 located more than 400 km east of the subduction trench.

Numerical analysis
We use a lithosphere-asthenosphere thermomechanical visco-elasto-plastic numerical geodynamic model (see Methods) to test the hypothesis of an upwind (i.e., trench-ward) migration of the SVZ forced by orography. A mantle-melting zone (MMZ) with 20 km diameter is imposed in the center of the model domain. The upward transfer and emplacement of magma into the crust occurs through a 3 km wide magmatic channel, which allows for a simplified representation of magmatic percolation by hydrofractures [3][4][5] , diffusion 32 , porous flow 33 , and reactive flow 34 through the rheologically stronger mantle lithosphere 35 (Supplementary Table 2). The parametric study (Supplementary Table 3 Table 3.

Tectonics, orography and volcanic arcs
The topography of an orogen commonly adjusts to lateral erosional gradients over timescales between 0.01-0.1 Ma (Ref. 36), which are significantly shorter than the Pliocene to present interval, during which orography and the onset of glaciation conditioned the orogen growth 28,29,37 in absence of major tectonic changes 22,23,38,39 . We can thus interpret our results in the light of the multi-Ma timescale orogenic and volcanic arc evolution.
The break-up of the Farallon Plate into the Nazca and Cocos oceanic plates 22 The climatic control on magmatism, thus, occurs not only through a forcing from ice building-melting 12,13 , erosion 11,44 and/or sea level changes 45 on the magma production, but also through an orographic forcing on the magma upwelling path, affecting the location of volcanic arcs. We expect this recognition to be broadly applicable to other volcanic arcs worldwide and particularly relevant for volcanic arc ore deposits, the exploitation of which is primarily conditioned by the regional topography itself. Because volcanic arcs provide a substantial contribution to the evolution of climate across timescales, this recognition provides an additional evidence of the tight coupling between climate and plate tectonics 46,47 .

Online Content
Methods, Additional References, and Supplementary Tables and Figures are available for this paper.

Authors Contributions
VAPM contributed to the development of the initial scientific question, performed the study and wrote the manuscript. PS conceived the initial scientific question, helped with the data interpretation and numerical investigation, and contributed the manuscript writing. CS contributed to the development of the initial scientific questions, helped with the interpretation of the geological data and numerical models, and contributed to manuscript writing.

Competing interest declaration
The authors declare no competing interests.

Corresponding Author line
Correspondence and request for materials should be addressed to VAPM.

Measurements of volcano-trench-water divide distances and zonal analysis
The relative distances between the volcanoes, the subduction trench, and the topographic water divide were measured based on volcanoes coordinates provided by Ref. 31. The location of the subduction trench and the topographic water divide was traced using Google Earth at a scale ~1:4,000,000. Measurements (Supplementary Table 1) account for the closer point along the trench and the closer point along the water divide to each volcano. The subdivision of the SVZ in zones 1, 2 and 3 is based on the petrochemical and structural arrangement of the arc [16][17][18][19] and help us to highlight meaningful correlations between measured distances.

Thermomechanical model
The numerical model is based on conservative finite-differences with marker-in-cell technique and incorporates temperature-dependent rheologies for solid and partially molten rocks 35,48 . The model accounts for a (1) mechanical, (2) visco-elasto-plastic rheological, (3) thermal, and (4) partial rock melting component. A short description of the models' components is provided hereafter. Additional details and full description of the numerical techniques are provided elsewhere 11,35,49 .

Mechanical Component
The continuity equation (1) (1) where is the local density, is time, ⃗ is the velocity vector, and ∇ is the divergence operator, allows for the conservation of mass during the displacement of a geological continuum.
The momentum equation (2) (2) where is the stress tensor, and are spatial coordinates, and is the i-th component of the gravity vector, describes the changes in velocity of an object in the gravity field due to internal and external forces. 2) are the shear viscosity for diffusion and dislocation creep, respectively, 0 is the material static viscosity, is the diffusion-dislocation transition critical stress, n is the stress exponent, is the activation energy, is the activation volume, P is pressure, R is the gas constant, T is the temperature, and ̇ is the second invariant of the strain rate tensor. Then, the viscous deviatoric strain rate tensor, ′ ( ) (equation 4), is computed as: where, ′ is the deviatoric stress tensor, is the Kronecker delta, ̇ is the volumetric strain rate (e.g., related to phase transformations), is the bulk viscosity.
Elastic deformation is reversible and assumes proportionality of stress and strain. The elastic deviatoric strain rate tensor, ′ ( ) (equation 5), is computed as: where, is the shear modulus, and ̆′ is the objective co-rotational time derivative of the deviatoric stress tensor.
Plastic (or brittle) localised deformation occurs at low temperature in the upper part of the lithosphere, after reaching the absolute shear stress limit, (equation 6), is defined as: where, is cohesion, is the effective internal friction angle. The plastic strain rate tensor, , is then computed as: where is the plastic multiplier which satisfies the plastic yielding condition = .
At the lithospheric scale, all deformation mechanisms occur jointly and the overall viscoelasto-plastic rock strain rate tensor, ′ ( ) (equation 8), is defined as:

Thermal Component
Heat conservation during advective and conductive heat transfer in the continuum is computed by the energy equation (9): where is specific heat capacity at a constant P, is the thermal conductivity, + + + are the volumetric heat productions by radiogenic, shear, adiabatic and latent heat, respectively. ∝ , = ′′ ( ) , and and are defined in Supplementary Table 2.

Partial Melting Component
Partial melting occurs between the wet solidus, , and dry liquidus, , of the considered lithologies 50-52 (Supplementary Table 2 where is compressibility, is thermal expansion and , 0 , and 0 are density, pressure and temperature of rocks at surface conditions.

Reference Model Setup
The initial setup measures 200 × 120 km in the x and y dimensions (Fig. 2