Abstract

Volcanic activity continued to occur along the length of the Baja California Peninsula (northwestern Mexico) even after the cessation of subduction during the middle Miocene. This volcanism occurred mainly in monogenetic volcanic fields, erupting lavas with a wide variety of compositions, including: adakites, niobium-enriched basalts, high-niobium basalts, and high-magnesian andesites. The chemical compositions of these magmas suggest an origin in partially melted basaltic oceanic crust that was subsequently subducted below the peninsula.

Several attempts have been made to explain the origin and compositional diversity of post-subduction volcanism in Baja California. Many of these attempts rely on the hypothesis that the magmas were formed through adiabatic decompression of upwelling asthenosphere in direct response to the formation of a window or tear in the subducted slab. This process, however, cannot offer a satisfactory explanation for all existing observations, particularly the lithospheric structure, of Baja California.

Here, we present a physical model that sheds light onto the origin of the post-subduction volcanism in Baja California. The model calls upon viscous dissipation as the causative agent of volcanism. Our starting conjecture is that shearing along a low-viscosity channel confined between the stalled oceanic slab and continental crust of Baja California peninsula caused partial melting at moderate depths following cessation of subduction. Our modeling results show that it is indeed possible for rocks to reach their solidus temperature by means of this mechanism. Numerical results indicate that shear heating could lead to a temperature increase of close to 200 °C at a depth of 30 km, sufficient to produce more than 30% melt by volume.

INTRODUCTION

Volcanic activity has occurred along the Baja California Peninsula (northwestern Mexico) since Oligocene time. The volume and composition of volcanic activity, however, have varied over time, in conjunction with changes in the tectonic setting of the peninsula. From the late Oligocene to early Miocene, a continuous belt of calc-alkalic continental arc volcanism occurred along what is now the eastern coast of peninsular Baja California (Hausback, 1984; Sawlan and Smith, 1984). Then, during middle Miocene time, the western margin of the North America plate shifted from a convergent to a transform plate margin (Mammerickx and Klitgord, 1982). The transform fault system developed between the Pacific and North America plates in the vicinity of the old trench forming the San Benito–Tosco–Abreojos fault system (TAFS) (Spencer and Normark, 1979). At ca. 6 Ma, a new plate boundary was created along the axis of the Gulf of California (Stock and Hodges, 1989; Umhoefer et al., 1994). Geological and geophysical data suggest that since then, Baja California has behaved as a microplate, with motion distinct from those of the Pacific and North America plates (Spencer and Normark, 1979; Legg et al., 1991; Dixon et al., 2000; Fletcher and Munguía, 2000; Michaud, 2004).

Since the middle Miocene, and despite the cessation of subduction, volcanism has continued up to Quaternary times and is manifested in several monogenetic volcanic fields located west of the extinct volcanic arc (Fig. 1). With the exception of the Santa Clara volcanic field, which is located near the Pacific shore of the peninsula, the loci of these volcanic fields are aligned along a narrow belt near the center of the peninsula. Moreover, post-subduction lavas display significant compositional variability. These magmas include adakites, niobium-enriched basalts (NEBs), high-niobium basalts (HNBs), and high-magnesian andesites that are locally termed bajaites (Aguillon-Robles et al., 2001; Benoit et al., 2002; Calmus et al., 2003; Calmus et al., 2011; Castillo, 2008, 2009; Negrete-Aranda and Cañon-Tapia, 2008; Pallares et al., 2008). All of these rocks fall within the family of adakitic rocks, which have been proposed to range from pristine slab melts, to adakite-peridotite hybrid melt, to melt derived from peridotite metasomatized by slab melt (Castillo, 2012). What is important, however, is that these magmas exhibit both mantle and slab-derived compositions.

It has been conjectured that right after the cessation of subduction during the middle Miocene, the stalled slab of the Farallon plate, which fragmented into the Magdalena and Guadalupe microplates, tore or broke off (Dickinson, 1997; Benoit et al., 2002; Calmus et al., 2003; Ferrari, 2004; Pallares et al., 2008; Calmus et al., 2011). The asthenospheric mantle ascended through the tension gash in the stalled slab (Fig. 2), decompressed, and partially melted the edges of the subducted slab to produce post-subduction volcanism. These ideas, however, are being re-examined in light of recent seismic observations; these observations show that several of the post-subduction volcanic fields sit atop a positive shear velocity anomaly (Fig. 3). This anomaly has been interpreted as a slab remnant still attached to the oceanic plate (Persaud et al., 2007; Obrebski and Castro, 2008; Zhang et al., 2007, 2009; Wang et al., 2009; Long, 2010).

Here, we present a simple physical model that explores a possible genetic relationship between the recently imaged slab and post-subduction volcanism in south-central Baja California. In its most fundamental form, the model consists of a low-viscosity channel confined between the stalled slab and continental crust of the peninsula that heats up by viscous dissipation. We use this model to investigate if the differential coupling between the slab and the overriding Baja California Peninsula can generate enough heat to produce partial melting of the slab remnant and of the lower continental crust of Baja California.

The model is sensitive to initial and boundary conditions, but for reasonable viscosity values, rates of plate motions, and thickness of the low-viscosity channel, we will show that it produces temperatures in excess of the solidus of wet basalt and hydrated peridotite, the possible sources of magma in the study area. Additionally, from the model, we are able to infer the chemical composition of the magmas in the simulations and compare them with observations of the temporal and spatial distribution of post-subduction volcanism.

CRUST AND MANTLE STRUCTURE OF CENTRAL BAJA CALIFORNIA

Recent seismic studies of the shallow mantle underneath the Baja California Peninsula (Wang et al., 2009, 2013; Long, 2010; Zhang et al., 2007, 2009; Obrebski and Castro, 2008; Persaud et al., 2007) cause us to question the formation mechanics of shallow slab gaps and/or slab windows (Dickinson, 1997; Benoit et al., 2002; Calmus et al., 2003, 2011; Ferrari, 2004; Pallares et al., 2008). For example, seismic tomography in the Gulf of California clearly shows that there is no tensional gash beneath south-central Baja California (Wang et al., 2009, 2013; Zhang et al., 2007, 2009). In fact, Rayleigh-wave seismic tomography studies have revealed a high-velocity anomaly in the shallow mantle underneath the continental crust of south-central Baja California. This anomaly appears to be a fossil slab still attached to the Guadalupe and Magdalena microplates (Fig. 3A; Wang et al., 2009, 2013; Zhang et al., 2007, 2009).

Even though the Rayleigh-wave tomography offers only a first-order picture of the lithospheric structure of Baja California, the remnants of the slab imaged by Wang et al. (2009) and Zhang et al. (2007, 2009) have been resolved at much higher resolution by Obrebski and Castro (2008) and Persaud et al. (2007). Further information about the anisotropic structure of the lithosphere can be obtained through receiver functions (RFs) analysis, and by forward modeling the anisotropic features observed on radial and transverse ground motion components (Levin and Park, 1998). Possible evidence of anisotropy in receiver functions obtained at several stations of the NARS-Baja Seismic Network were reported earlier by Persaud et al. (2007), and later modeled by Obresbki and Castro (2008). The former authors presented possible forward models for the anisotropy of the shallow lithosphere beneath three stations located in the different tectonic provinces that compose the northern-central Gulf of California region. The models show reasonable consistency with recent tectonic processes that occurred in this area. For instance, modeling of station NE75 of the NARS-Baja array in the central peninsula (Fig. 3B) suggests the presence of either deformed oceanic crust or metamorphosed mantle wedge. The metamorphosed layer is situated ∼5 km beneath the continental crust, which has a local thickness of 30 km. The lowermost 4 km of the continental crust beneath NE75 also seems to have registered the effects of shear associated with the subduction of the Farallon slab (Obrebski and Castro, 2008) (Fig 3B).

Observations by Persaud et al. (2007) reveal a similar structure for central Baja California. Their analysis of station NE75 is consistent with that of Obrebski and Castro (2008) and points to the existence of a slab beneath this station at a depth of ∼30 km (Fig. 3B). Thus, it is clear that different authors performing independent studies and using different methodologies produced similar models for station NE75. It is also noteworthy that the presence of a slab beneath south-central Baja California is in good agreement with mechanical models developed by Bohannon and Parsons (1995) and Plattner et al. (2010). Those authors proposed that, because of the buoyancy of the Farallon slab, most of the fragments of that slab may have stalled beneath the North America borderland after the cessation of subduction. Furthermore, a resistivity model of the continental crust, developed by Romo-Jones (2002) for an area near parallel 28° N, revealed a flat, conductive anomaly across the peninsula at a depth of ∼30 km. Romo-Jones interpreted the conductive body as a contact between Baja Peninsula crust and the stalled slab of the Guadalupe microplate.

NUMERICAL MODEL OF POST-SUBDUCTION VOLCANISM

We here develop a simple thermomechanical model of the lithosphere of Baja California. The model is motivated by observations in xenoliths recovered from lavas of the San Quintín volcanic field in northern Baja California (not shown on map in Fig. 1). Those xenoliths sampled an active shear zone in the mantle and lower crust of the peninsula (Cabanes and Mercier, 1988; Palasse et al., 2012). The shear zone sampled is at a depth of <35 km (1 GPa) and has a temperature between 800 and 950 °C. The model is also motivated by recent GPS data, which show that the TAFS still accommodates ∼1 cm/yr of slip with respect to the Pacific plate (Fig. 1; Dixon et al., 2000; Plattner et al., 2007, 2009, 2010). Estimates of frictional heating along strike-slip faults show that they can produce enough heat to increase the temperature of adjacent rocks by more than 200 °C (Leloup et al., 1999). Numerical experimentation shows that in the extreme situation of a slip rate of 10 cm/yr and a stiff lithosphere, shear heating at the strike-slip shear zone could raise the temperature by as much as 590 °C at the Moho and by as much as 475 °C at a depth 20 km (Leloup et al., 1999).

Our model is framed for the period starting when subduction had stopped and ending when the transcurrent regime had settled in. Our simulation starts at 14.5 Ma, during the middle Miocene. It was during the period between 14.5 and 12 Ma that the TAFS became a transform boundary (Spencer and Normark, 1979; Michaud et al., 2006). With the model, we investigate whether the differential motion of the peninsula causes partial melting in the mantle, continental crust, and the stalled slab imaged by the tomography.

The model is specific for the south-central part of the Baja California Peninsula (Fig. 1). This region is where seismic imaging shows that the fossil slab is still attached to the oceanic plate (Persaud et al., 2007; Obrebski and Castro, 2008; Zhang et al., 2009; Wang et al., 2009; Long, 2010; Figs. 3 and 4). It is also in this region where observations constrain, within a few kilometers vertically, the structure of the continental crust and upper mantle (e.g., Persaud et al., 2007; Obrebski and Castro, 2008; Fig. 3B).

Figure 4 presents the conceptual tectonic model for Baja California on which we base our simple physical construct. Following the ideas of Nicholson et al. (1994) and Plattner et al. (2010) and the observations of Persaud et al. (2007) and Obrebski and Castro (2008), the model has an oceanic slab that tectonically underlies the base of the crust for a distance of ∼200 km from the trench. The slab is attached to the Pacific plate (Wang et al., 2013) and is being dragged below the North American crust parallel to the TAFS. Also, as suggested by Obrieski and Castro (2008), the model incorporates a mantle sliver a few kilometers thick between the oceanic slab and the overlying Baja California continental crust. We assigned to this layer a material behavior similar to that of the asthenosphere; we extrapolated estimates of the rheological behavior of the shallow mantle derived from xenoliths in the San Quintín volcanic field further to the north (Palasse et al., 2012).

On the problem of whether it is dynamically feasible for a slab to be dragged beneath the overriding plate, Pikser et al. (2012) concluded that along-strike translation of a remnant slab is a viable process. Results of three-dimensional finite element models demonstrate that the viscosity contrast between slab and surrounding mantle is greater than or equal to the requirements needed for the edgewise translation of the slab.

Our physical construct now incorporates the general features of the complex structure presented in Figure 4. Similar to other models, however, it is highly idealized. To begin with, it is two-dimensional (Fig. 5). The modeled cross-section is oriented parallel to the axis of the Baja California Peninsula and passes through the Santa Rosalía, San Ignacio, and La Purísima volcanic fields (Figs. 1 and 5). This orientation coincides with the Pacific plate flow lines; this, in turn, allows us to use a plane strain formulation in our model. Additionally, the modeled region contains three horizontal layers of variable mechanical and thermal properties, representing the continental crust, mantle wedge, and stalled slab.

The equations governing the behavior of the layers are in an Eulerian frame where material does not physically move, but the internal flow velocity and temperature fields are computed based on boundary conditions. Such a frame is often adopted in fluid mechanics and geodynamical problems (Gerya, 2010). The mechanical and thermal behavior of the layers is discussed next.

Given the long-lived activity of the post-subduction volcanic fields (∼10 m.y.), we chose a simple linear Newtonian rheology for the layers. This means that at time scales larger than a few million years, the lithosphere behaves as an incompressible viscous fluid (Ranalli, 1995). The equations describing the flow of viscous substances are the Navier-Stokes equation and the continuity equation (Turcotte and Schubert, 2002): 
graphic
 
graphic
where ρ is the density of the rocks, v is the velocity of the flow field, forumla is the symmetric part of the strain-rate tensor, p is the internal pressure, and η is the viscosity. The continental crust and slab are stiff layers whereas the mantle sliver forms a low-viscosity channel.
The change in temperature in the layers is modeled by means of the convection-diffusion equation: 
graphic
where T is the temperature field, Cp is the heat capacity of the rocks, q is the conductive heat flow, Hr is the rate at which heat is released by radioactive decay, and Hv is rate of heat produced by viscous dissipation. Radioactive sources are confined to the upper continental crust and decay exponentially with depth (Turcotte and Schubert, 2002): 
graphic
In Equation 4, θ is the radioactive heat generation of the rocks and hr is the characteristic thickness within the continental crust over which the radioactive materials are concentrated (Table 1).
Now, viscous dissipation is the heat produced during irreversible non-elastic deformation and is given by (Turcotte and Schubert, 2002): 
graphic
The model also deals with the production of partial melt. As a first-order approximation of this complex phenomenon, it is assumed that the degree of melting increases linearly with temperature according to the following relations (Gerya, 2010): 
graphic
where m is the volumetric fraction of melt, Ts is the solidus temperature and Tl is the liquidus temperature of magma. Observe that these two last critical temperatures strongly depend on pressure and chemical composition of rocks. Here we use the empirical relationship presented in Gerya (2010) and references therein for dry granite for the continental crust (see Table 1). This assumption is justified by a long line of granitic intrusions along Baja California that may have left anhydrous residues at depth. For the mantle wedge and upper basaltic section of the slab, we use the parameterization by Katz et al. (2003); that parameterization models melt formation from hydrated peridotite and basalts (see also Hess, 1989; Hirschmann, 2000; Johannes, 1985; Poli and Schmidt, 2002). More than 100 m.y. of subduction history in the area supports this last assumption.
It is noteworthy that when rocks undergo melting, heat is absorbed by the liquid phase in the form of latent heat. This, in turn, is accompanied by an increase of the heat capacity of the liquid phase. Because our model is continuous, we address this effect by using an effective heat capacity, Cpeff, where the model is undergoing partial melting (Burg and Gerya, 2005): 
graphic
In Equation 7, Hl is the latent heat of melting. The model also incorporates empirical laws for thermal conductivity and thermal diffusivity, the physical properties governing conduction of heat. Experimental data indicate that these properties strongly decrease with increasing temperature (Whittington et al., 2009). Thus, the middle to lower crust acts as a thermal insulator. This is an important effect, because it provides a positive feedback between heating, thermal insulation, and partial melting (Whittington et al., 2009). Here, we use the following equations for the temperature dependence of thermal diffusivity (κ) in the continental crust (Whittington et al., 2009): 
graphic
For the mantle wedge and slab, we use the following relations for the thermal conductivity (λ) (Gerya, 2010): 
graphic
where λ0 = 1.18 W m–1 K–1 and a = 474 W m–1. λ and κ are related to one another by: 
graphic
It should be noted, however, that κ and λ are not completely proportional, owing to an increase in Cp with increasing temperature. For a typical crust composition, the following equations describe how Cp changes with temperature (Whittington et al., 2009): 
graphic
In this last expression, Cp is in joules per mole per Kelvin.

Table 1 summarizes the values of the material parameters used in the model and other empirical relations that govern their behavior. Figures 5 and 6 illustrate the material structure of our model. The thickness of the first two layers is reasonably well constrained by seismic observations (Persaud et al., 2007) and the thickness of the slab is estimated by means of the classical solution of the cooling of a half-space using a slab age based on seafloor isochrons along the extinct subduction zone (Nicholson et al., 1994). The viscosity of the mantle sliver is 5 × 1019 Pa-s (Palasse et al., 2012). Notice that due to limitations in our modeling code, low-viscosity blocks (η = 1017 Pa-s) with very high thermal conductivity (λ = 500 W/mK) flank our jelly-sandwich model (i.e., a low-viscosity channel confined between the stalled slab and continental crust of the peninsula). These side blocks are useful modeling artifacts and do not represent geological entities. We introduced the side blocks to minimize boundary effects by corner flow. As we will show, the low-viscosity zones do not affect the results in the center of the model, which is the area of interest.

Further, note in Figure 6 the abrupt changes in mechanical and thermodynamical properties among the layers. We will see that these provide the conditions necessary for slab melting (producing adakite) and production of slab melt–related magmas (e.g., bajaite and NEB). Now, in our model, the bulk of melt production is confined to the low-viscosity channel and the first 10 km of the simulated slab, which corresponds to the oceanic crust. For that reason, we decided not to include a sub-slab mantle layer in our modeled region (lowermost 10 km of slab). That is, the slab has a uniform chemical composition of oceanic crust.

The following boundary conditions are imposed. Both the top and bottom boundaries are isothermal surfaces. The former has a temperature of 300 K (∼25 °C) whereas the latter has a temperature of 950 K (∼680 °C). The temperature at the bottom is constrained by the solidus temperature of the modeled slab (Fig. 6C). These temperatures correspond to a thermal gradient of 12 K/km, which is compatible with the low heat flow observed in the forearc of subduction zones (∼30 mW/m2). Along the lateral boundaries, we impose insulating boundary conditions (i.e., the conductive heat flow, q, is zero) and reflective boundary conditions for the flow field. As an initial condition, we assume a linear distribution of temperature given by a temperature gradient of 12 K/km.

The top boundary of the model is stress free, meaning that topography can develop on this boundary. By contrast, the slab is set in forced motion along its bottom boundary (Fig. 5): 
graphic
where vs(t) is the velocity of the slab. For middle and late Miocene times, the synthetic slab slides under the continental crust at a rate of 40 mm/yr. This rate of movement is based on tectonic reconstructions, which show that at that time the TAFS accommodated most of the plate boundary transform motion (Spencer and Normark, 1979; Atwater and Stock, 1998; Michaud et al., 2006). For Pliocene to present times, we use a rate of 7 mm/yr. In the model, the shift occurs instantaneously at ca. 6.5 Ma. At 6.5 Ma the Gulf of California became the new plate boundary, leaving only a small fraction of the deformation to be accommodated by the TAFS. Recent GPS measurements of residual velocity at the TAFS are in good agreement with the parameters used in the model for recent time simulations (Plattner et al., 2010). More importantly, geodynamic models call for the slab to remain attached to the oceanic crust (Spencer and Normark, 1979; Nicholson et al., 1994; Michaud et al., 2006).

RESULTS OF NUMERICAL EXPERIMENTATION

Model Equations 1–12 were solved numerically by means of the finite difference method. We used the code developed by Gerya (2010); that code was extensively adapted for the boundary conditions and material structure of the problem addressed here. The node spacing used in the finite difference grid is 1500 m in the x direction and 500 m in the y direction. Our simulation starts at 14.5 Ma, during the middle Miocene, and uses discrete time steps of 500 k.y. Figure 7 presents the thermal evolution for the Baja California lithosphere as predicted by Equation 3, and Figure 8 illustrates the production of melt as per Equation 6.

Figure 7 shows that a thermal anomaly starts to grow in the modeled mantle wedge within a few million years of initiation of the differential motion between the simulated slab and continental crust. By 9 Ma the thermal anomaly is large enough to reach the solidus temperature and partial melting starts (Fig. 7). The higher thermal conductivity of the crust (Fig. 6B) allows for a rapid cooling of the generated heat, whereas the slab’s lower thermal conductivity dissipates the generated heat slowly, allowing for a heat buildup. This effect is enhanced by the presence of melt, which further depresses the thermal diffusivity (see Equation 7). Observe that melt is being generated within the oceanic slab, near the interface between the slab and mantle wedge (Figs. 8 and 9A). Further, note that our model lacks a mechanism for magma transport (e.g., Aharonov et al., 1997; Vasilyev et al., 1998); therefore, the melt stays in situ. By 6 Ma, the time at which strike-slip motion shifted to the axis of the gulf, the temperature in the shear zone reaches 1000 K; this temperature elevation is sufficient to cause 40% of partial melting of the underlying slab (Fig. 9A). Also notice that the thermal anomaly crosses the wet peridotite solidus line, producing partial melting of the mantle sliver. After this time, the thermal anomaly loses expression; note that most of the heat subsequently generated in the mantle wedge is slowly dissipated by the oceanic slab, while the crust rapidly cools. Toward the end of the simulation, the thermal anomaly barely crosses the solidus line of wet basalt/gabbro. Consequently, partial melting in the slab is only a few percent (Fig. 9B).

Our model successfully explains the occurrence of bajaites and the association of adakite and NEB in post-subduction volcanism. As discussed below, lavas of such chemical composition are systematically linked to melting of the mantle wedge that had been previously metasomatized by slab melt through the formation of amphibole (Rogers et al., 1985; Saunders et al., 1987; Rogers and Saunders, 1989; Calmus et al., 2003). The numerical experiment produces melting in the required time frame: the late Miocene. Perhaps more importantly, the modeled conditions are capable of producing partial melting even after the velocity of the slab is reduced by a factor of four at 6 Ma, that is, to the slip rate currently observed along the TAFS (Fig. 9B).

DISCUSSION

From studies by Rogers and Saunders (1989) and Rogers et al. (1995) to recently published reports, the slab-melting model has been considered as the most viable mechanism to explain the atypical post-subduction volcanism in Baja California (Bourgois and Michaud, 2002; Benoit et al., 2002; Calmus et al., 2003, 2011; Ferrari, 2004: Pallares et al., 2008; Castillo, 2008). However, the initial, basic premise of the slab-melting model has undergone several modifications over the past decade. For instance, Benoit et al. (2002) suggested that the subduction of an active ridge would have opened a slab window through which the fertile sub-oceanic mantle would have upwelled, partially melting the slab. The altered basaltic crust from the edges of the tear would have also melted, generating slab melts that either reached the surface and formed high-silica adakites or metasomatized the mantle and generated high-field-strength-element- (HFSE) enriched amphiboles (Defant et al., 1992; Benoit et al., 2002; Castillo, 2008). Melting of the amphibole-metasomatized lithospheric mantle would have produced HNBs and NEBs. At the same time, a relatively higher degree of melting of the metasomatized mantle would have generated tholeiitic magmas. Finally, melting of the amphibole-metasomatized mantle, where >5% garnet was stable, produced bajaites.

From the previous arguments, it can be seen that the slab-melting model offers a reasonable explanation for the peculiar post-subduction volcanism of Baja California. The model, in fact, is supported by comparable rock associations elsewhere. However, as discussed by Negrete-Aranda and Cañón-Tapia (2008) and Castillo (2008), the essential hypothesis of the slab-melting model, that is, the subduction of an active spreading center and the subsequent formation of a slab window, does not relate well with observed evidence. Further, the model is at odds with the more recent observations of the seismic structure and thermal state of the peninsula, and with numerical models (Wang et al., 2009; Zhang et al., 2009; Burkett-Andrews and Billen, 2009a, 2009b). A further problem is that recent studies have shown that adakitic rocks can be produced through a number of different processes. These processes include melting of the mafic lower crust, slab anatexis and high-pressure crystal fractionation of basaltic magma, and low-pressure crystal fractionation of basaltic magma plus magma mixing processes (Castillo, 2006, 2008, 2009; Thorkelson and Breitsprecher, 2005; Petrone and Ferrari, 2008).

To achieve an adakitic composition, amphibole is required as a reactant and garnet is required as a residual phase (Drummond et al., 1996). According to Thorkelson and Breitsprecher (2005), adakitic signatures can be explained in terms of pressure-temperature paths, along which slab anatexis is possible at depths of ∼5–150 km. Within this broad range of slab melting, there is a smaller field, within which both garnet and hornblende are stable in the slab and adakitic melts are predicted in their pressure-temperature and phase stability diagrams. The region of possible adakite formation lies in a field bordered by the wet basalt solidus, the garnet-in line, and the amphibole-out line. Pressure-temperature paths, which pass through the field, indicate that adakite melt generation is constrained to ∼25–90 km. These general observations are in good agreement with our model predictions because results presented in this paper locate the melting of the slab within this depth range.

The simulations presented here recreate pressure-temperature conditions required by the atypical post-subduction volcanism in Baja California. They do so without the upwelling of asthenospheric material through a slab window or other complex geodynamic processes. More importantly, our model conforms to the lithospheric structure of Baja California, as imaged by seismological methods, and the relative motions of plates.

We now contrast our model results with petrogenetic models of post-subduction volcanism. We favor no particular petrogenetic model in our analysis. However, we find that our model predicts the necessary thermal conditions needed to explain the observed geochemical compositions. We will focus on the bajaite and the adakite-NEB association, which are by far the most voluminous volcanic products scattered throughout the peninsula.

Bajaite is of particular interest, because it is a type of high-magnesian andesite, which is compositionally akin to bulk continental crust (Kelemen, 1995; Rudnick and Gao, 2003). In the case of Baja California, numerous authors have considered bajaite as derived from the melting of mantle peridotites previously metasomatized by slab-derived fluids in a high-thermal regime (Rogers et al., 1985; Saunders et al., 1987; Rogers and Saunders, 1989; Calmus et al., 2003, 2011; Pallares et al., 2008; Castillo, 2006, 2008).

In this regard, our model is consistent with the latter ideas, as the thermal anomaly predicted by the model causes partial melting of both the simulated mantle sliver and the underlying slab (Figs. 8 and 9). However, it is noteworthy that there are still a few inconsistencies in this model. The mantle wedge in our model is only 5 km thick (Figs. 4 and 5); thus, it is uncertain whether there were enough pressure and temperature differences in such a thin mantle to produce the chemically and isotopically different bajaites, HNBs/NEBs, and tholeiites. Field relationships indicate that adakites, bajaites, and HNBs/NEBs were erupted almost simultaneously; this implies that slab-melt metasomatism and partial melting of the mantle wedge occurred almost at the same time.

Figure 10 shows that a slight compositional change in the upper crust leads to numerical solutions consistent with some of the alternate models for the production of the previously mentioned adakites. For example, let us assume that the lower crust is composed of basic residues left by the granitic intrusions found at shallower crustal levels. As discussed below, a heterogeneous layered continental crust with a wet basaltic bottom section has lower solidus temperatures than that of granitic rocks (Fig. 10). Consequently, the lower crust would generate partial melt under shearing at its interface with the mantle. This may explain the generation of adakitic rocks according to the mechanism proposed by Castillo (2008, 2009).

The results presented in this paper clearly show that the rather complex slab-window tectonic settings formerly invoked are not necessary to create the thermal conditions needed to produce the particular geochemically diverse post-subduction volcanism in Baja California. The thermal anomaly in the numerical model, however, is sensitive to the prescribed width and depth of the shear zone, and strain rate. For example, if the width of the channelized shear flow zone were thicker than 10 km, the shear zone in the model would produce a faint anomaly, and no melt would be produced at all. The reason for this result is that the rate of heat generation in the low-viscosity channel scales inversely with the square of the thickness of the shear zone (Turcotte and Schubert, 2002).

The model is also highly sensitive to the depth of the shear zone. For example, detachments located at depths shallower than 24 km do not result in a temperature anomaly. This result occurs because the thermal conductivity of the continental crust changes roughly inversely with depth (Equation 8 and Fig. 6B). A thin continental crust is more conductive than a thick one, provided that the thermal gradient is the same; as a result, the formation of a thermal anomaly is inhibited. Similarly, simulations in which the shear zones are located at depths in excess of 28 km result in modest thermal anomalies. The model, therefore, is sensitive to the thickness of the modeled slab; this effect is likely caused by the rather strict boundary condition, T = 950 °C, imposed at its bottom.

The thermal history of our model is also sensitive to the initial velocity. In our derivation, we assumed 40 mm/yr as the velocity of the slab during late Miocene time and 7 mm/yr for Quaternary time. Numerical experimentation shows that the slower the initial velocity, (1) the weaker the thermal anomaly becomes, (2) the longer it takes to develop, and (3) the sooner melt production shuts down. Basically, we were unable to fit the volcanic history of the peninsula for vs(t = 0) ≤ 20 mm/yr.

As discussed in the preceding paragraphs, only within a narrow parameter space does the model develop growing thermal anomalies leading to partial melting of the oceanic slab. The limited viable parameter space of the model suggests a simple explanation as to why volcanic fields are aligned into such a narrow belt (Fig. 1). To see why this is the case, consider a scenario in which the volcanoes are situated above the optimal melt production window, the slab gently dips toward the east, and both the mantle sliver and continental crust are wedge shaped (see east-west cross section in Fig. 4). The model tells us that as the continental crust thins (by even a few kilometers) west of the melt production window, the lower crust loses its thermal insulating properties. This loss hinders the generation of melt. As the low-viscosity channel widens to the east, it becomes too thick to generate any substantial heat; this loss of heat generation causes melt production to shut down.

The main simplifying assumption in the model is that rocks have constant viscosity. In this regard, numerous studies have shown that the long-term viscosity of rocks strongly depends on temperature following an Arrhenius-like viscosity law (Ranalli, 1995). In spite of this limitation, geodetic and seismic results show that the viscosity structure used here bears resemblance to the viscous layering of western North America. For example, modeling of geodetic data in the Salton Trough (southern California and northern Baja California), an area of high heat flow being deformed at high strain rate, indicates the continental crust is stiff (Pollitz et al., 2001). Moreover, viscoelastic modeling done by Fay and Humphreys (2005) shows that the southern California crust has a long-term mean viscosity of 1023 Pa-s under near-surface conditions, and that this value decreases gradually with depth to a value of 1021 Pa-s at the mantle-crust interface. Further support comes from the long history of magmatic intrusions in Baja California. These pervasive granitic batholiths may have left an anhydrous lower continental crust of increased strength (e.g., Maggi et al., 2000).

On the other hand, studies on the upper-mantle viscosity under western North America indicate it has a near-solidus temperature and a viscosity in the order of 1019 Pa-s (Freed and Bürgmann, 2004; Dixon et al., 2004). Further evidence that the mantle wedge beneath western North America may have been weak comes from seismic observations and numerical simulations in low-angle subduction zones (e.g., Manea and Gurnis, 2007; Pérez-Campos et al., 2008; Li et al., 2012). These studies show that low-angle subduction involves a rheologically weak zone between the subducting and upper plates. Consequently, the low-viscosity mantle wedge suggested by our model could have been inherited from the previous history of low-angle subduction along the western margin of Baja California (Bohannon and Parsons, 1995).

With respect to the rheological behavior of the oceanic slab, there are two important effects that were left out in our model. The first one is a strong negative feedback between production of magma and shear heating. As Mg-silicate fluids are many orders of magnitude less viscous than mantle rocks over the same pressure regime (e.g., Karki and Stixrude, 2010), it is expected that the low-viscosity channel would expand into the basaltic slab. (Note that the melt was initially produced at the interface between slab and mantle wedge; Fig. 8.) Because the rate of heating strongly depends on the thickness of the shear zone, its expansion will rapidly weaken the slab and self-limit melt production. The second effect is a positive feedback between extraction of melt and shear heating. That is, as magma is being produced, it percolates upwards, leaving behind a thinner residual slab that is prone to further melting.

The preceding discussion invites the obvious question: Which of these two processes governs the long-term rheological response of the slab? In trying to answer this question, we need to look into the overall permeability of the slab, mantle, and crust; the density contrast between melt and rock; pore volume compressibility; and viscous compaction length (Aharonov et al., 1997; Vasilyev et al., 1998). Nevertheless, as a first-order observation, it can be established that for magmas to preserve their geochemical signatures (e.g., NEBs, HNBs, and bajaites), they must rise from depths of at least 30 km to the surface without extensive re-equilibration with surrounding mantle and crust. These lines of evidence indicate that melt flow in Baja California was focused through high-permeability structures. This, in turn, suggests that once the adequate conditions for magma transport develop, the positive feedback between extraction of melt and shear heating should dominate over the self-limiting process of production of melt.

Finally, we would like to point out that the model requires the slab to be dragged a distance of ∼350 km to the north of the main study area. Unfortunately, existing tomographic models show considerable uncertainty in the north of the peninsula with overall average to low velocities in this region. Indeed, areas of low velocities in the north seem to be consistent with proposed slab-window theories (Zhang et al., 2007). In any case, the slab is unambiguously rendered by tomography up to 29° N. This explains ∼230 km of the required ∼350 km of lateral translation. The remaining ∼30% discrepancy can be explained by a combination of deformation, erosion, and thermal equilibration along the lateral edge of the slab, which could have been caused by compression heating or viscous dissipation. A more detailed plate reconstruction history, high-resolution seismic tomography, and more robust geodynamic models could probably resolve this issue.

SUMMARY

No consensus yet exists regarding the specific source regions and processes involved in producing magmas parental to the anomalous post-subduction volcanism in Baja California. However, taken at face value, our model supports a common origin or generation process for adakites, NEBs, and bajaites from the subducted slab and metasomatized mantle due to the high thermal regime. This work reveals that shear heating might have created such thermal anomaly and the conditions that triggered the production of post-subduction volcanism. Model results are consistent with first-order observations in the sense that it predicts the formation of bajaites, adakites, and HNBs/NEBs in a narrow melt-production window. Moreover, when the model is modified to incorporate a more realistic crustal composition, we are able to accommodate alternate models of the production of post-subduction lavas as well.

According to our model, after ∼14 Ma, the stalled slab under the Baja California crust became heated by a growing thermal anomaly located at a shear zone between them. Syn-subduction magmatism was extinguished, and eventually was replaced by magmas of slab composition. This response was driven largely by physical, thermal, and chemical changes to both the underlying mantle and the subducted slab. During active subduction times, the Farallon-Pacific subducted slab beneath Baja California was young and therefore buoyant, and the mantle wedge that formed the shear zone was hydrated and of low viscosity and was produced as a result of the long history of subduction in the area. The thermal structure modeled in this paper supports the general idea that post-subduction volcanism of Baja California resulted from melting of young oceanic crust. The model presented here differs from other models, however, by invoking and quantifying significant amounts of shear heating as the cause of slab melting. Our analysis underlines the importance of considering the effects of shearing in young oceanic crust during the transition from a convergent to a strike-slip boundary.

We are grateful to the Exploration Office of PEMEX-PEP for the access to borehole temperatures and other confidential material, which helped in the development of ideas presented in this paper. The authors thank Pat Castillo, Arturo Gomez-Tuena, David Hilton, Stephen Smith, and an anonymous reviewer for their very constructive comments and suggestions. We would also like to acknowledge the technical expertise of Victor Manuel Frias-Camacho, Sergio Arregui-Ojeda, Jessica Salas, and Martin Francisco Pacheco-Romero. This research was funded by the Mexican Council of Science (CONACYT) grant no. 60647. R.N. is also grateful to CONACYT for providing her with post-doctoral fellowship 165662 and the National Science Foundation grant no. GRD 8806. Additional funds were provided by Programa para el mejoramiento del profesorado de la Secretaría de Educación Publica (SEP-PROMEP) grant no. /103.5/12/335 to R.M.S.