A consensus is emerging from studies of continental rifts and rifted margins worldwide that significant extension can be accommodated by magma intrusion prior to the development of a new ocean basin. However, the influence of loading from magma intrusion, lava extrusion, and sedimentation on plate flexure and resultant subsidence of the basin is not well understood. We address this issue by using three-dimensional flexural models constrained by geological and geophysical data from the Main Ethiopian Rift and the Afar Depression in East Africa. Model results show that axial mafic intrusions in the crust are able to cause significant downward flexure of the opening rift and that the amount of subsidence increases with decreasing plate strength accompanying progressive plate thinning and heating during continental breakup. This process contributes to the tilting of basaltic flows toward the magma injection axis, forming the typical wedge-shaped seaward-dipping reflector sequences on either side of the eventual rupture site as the new ocean basin forms.

During continental rifting the extending continental lithosphere is progressively modified by magmatic processes and mechanical deformation until it breaks and gives rise to new oceanic lithosphere bordered by young rifted continental margins. Although Earth materials are weak in extension, cratonic lithosphere retains significant strength during rifting, as evidenced in rift flank uplift and flexural isostasy studies of active rifts (e.g., Weissel and Karner, 1989; Ebinger et al., 1989). Progressive extension and plate thinning cause decompression melting of the asthenosphere at progressively shallower depths, leading to an increase in magma production through time (e.g., White and McKenzie, 1989; Buck, 2006), primarily focused in axial volcanic segments that mark the new plate boundary (e.g., Barberi and Varet, 1977; Ebinger and Casey, 2001; Wright et al., 2006). The sustained emplacement of the new magma into a narrow zone of the lithosphere and eruption at the surface modify the plate density profile, causing significant internal and surface loading. These loads promote plate flexure along the long and narrow intrusion zones, and change the state of stress within the plate (Beutel et al., 2010). Despite the importance of magmatism during late-stage rifting and continental breakup, the influence of loading from magma intrusion, lava extrusion, and sedimentation on plate flexure and resultant subsidence of the basin is not well understood. In this paper we describe flexural models that quantify the effect of axial intrusion and surface volcanism on basin subsidence within volcanic rifts, and demonstrate the increasing role of magmatic loading with progressive plate thinning and weakening. Extensional stresses induced by the plate bending beneath the volcanic load may facilitate the rise of magma from the mantle. These models are parameterized from seismic imaging, gravity surveys, and studies of plate strength (elastic thickness of the lithosphere) of the Main Ethiopian Rift (MER) and the Afar Depression in East Africa (Fig. 1), which represent unique modern-day analogs for processes occurring during continental breakup. The East African Rift zone is one of the few locations on Earth where the rift inception to rupture progression can be analyzed (e.g., Hayward and Ebinger, 1996; Bastow and Keir, 2011).

The MER and Afar Depression formed as a result of the divergence of the Nubian, Somalian, and Arabian plates (Fig. 1) above an anomalously hot (e.g., Rooney et al., 2012; Ferguson et al., 2013), slow-seismic-wavespeed mantle (Montelli et al., 2004; Benoit et al., 2006; Bastow et al., 2005, 2008). Three plates interact within the Afar Depression, which comprises the subaerial Red Sea and Aden rifts and the northern termination of the MER. Whereas the MER opening is constrained to ∼6 mm/yr, ∼N100°E-directed motion between Nubia and Somalia (e.g., Bilham et al., 1999), both the Red Sea and Gulf of Aden systems show time-averaged extension rates of ∼15–20 mm/yr related to the ∼N35°E-oriented Africa-Arabia rifting (Fig. 1; e.g., Vigny et al., 2006; ArRajehi et al., 2010). These rates are augmented by discrete, large-volume dike intrusion events that accommodate decades to centuries of plate opening (e.g., Wright et al., 2012; Grandin et al., 2011). Geological, geodetic, and seismological measurements of active strain and magma intrusion across the Afro-Arabian rift zone show that south of ∼16°N the rift bifurcates into two branches (the main Red Sea and the subaerial Red Sea rift in Afar), with an along-strike partitioning of extension. North of ∼16°N, extension is within the main Red Sea Rift, whereas moving south of 16°N, extension is accommodated progressively on land within the Afar Depression (McClusky et al., 2010).

Geological and geophysical data show that the MER-Afar system opened by faulting and magmatism, with along-rift variations in style of extension interpreted as the expression of different stages in an evolutionary rift sequence (e.g., Hayward and Ebinger, 1996; Keir et al., 2013). The crustal thickness decreases from >35 km beneath the southern-central MER to <15 km in northern Afar (e.g., Makris and Ginzburg, 1987; Dugda et al., 2005; Stuart et al., 2006; Hammond et al., 2011), concomitant with a decrease in elastic thickness of the lithosphere from >15 km in the southern-central MER to ∼5 km in northern Afar (Fig. 1; Hayward and Ebinger, 1996; Pérez-Gussinyé et al., 2009). The seismogenic layer thicknesses show a similar decrease from ∼15 km in the MER (Keir et al., 2006) to ∼5–7 km in northern Afar (Ayele et al., 2007; Nobile et al., 2012). Coincident with crustal thinning, geophysical and geological data indicate a northward increase in the volume of mafic magmatic intrusions and extrusive volcanism, which is also associated with a reduction in the elevation of the rift floor from ∼1700 m above sea level in the southern-central MER to below sea level in the Danakil Depression of northern Afar (e.g., Bastow and Keir, 2011; Keir et al., 2013). Whereas in the southern-central MER extension is mostly accommodated through tectonic faulting at rift margins, in the northern MER and Afar strain migrated in the past 3–7 m.y. from Oligocene and Miocene border faults to ∼15-km-wide, 60-km-long axial zones of localized faulting and magmatism, where magma intrusion plays a larger role than faulting in strain accommodation (e.g., Ebinger and Casey, 2001; Wolfenden et al., 2005; Ebinger et al., 2013). Geophysical data suggest that as much as 50% of the crust beneath the axial zones of localized faulting and magma intrusion is new igneous material. Crustal tomography reveals the presence of anomalously fast, elongate bodies in the middle to lower crust extending along the rift axis as shallow as ∼10 km (Fig. 2; Keranen et al., 2004; Daly et al., 2008). These 20-km-wide and 50-km-long bodies are interpreted to represent cooled mafic intrusions (Keranen et al., 2004; Daly et al., 2008) containing at least 40% gabbro (Cornwell et al., 2006). The lateral variations in crustal properties and the presence of magmatic underplating make it difficult to unambiguously image the Moho across the rift. In this study we assume that the entire width of the rift extended by ∼50%, and that the narrow magmatic segments approach 100% extension, primarily through the intrusion of new igneous crust through a process analogous to slow-spreading mid-ocean ridges (e.g., Iceland). The vertical and horizontal dimensions and density contrast of the large gabbroic bodies intruded into the crust during late-stage rifting have been interpreted to influence the vertical motions, state of stress, and architecture of the continental lithosphere (Keir et al., 2006; Beutel et al., 2010).

The data have been interpreted to indicate an along-axis variation in rift evolution, from embryonic rifting in the southern MER to incipient oceanic spreading in northern Afar (e.g., Hayward and Ebinger, 1996; Pagli et al., 2012), where Pleistocene–Holocene narrow axial basalt ranges (e.g., Erta Ale) are interpreted to be subaerial equivalents of oceanic spreading centers (e.g., Barberi and Varet, 1977). In the Danakil Depression of northern Afar, basalt flows from the axial ranges typically flow away into the lower lying, evaporite-rich basin (e.g., Pagli et al., 2012; Field et al., 2012), creating a thick (up to ∼5 km) basin stratigraphy of thinly interbedded basalts and evaporites (Talbot, 2008), a geology typical of the seaward-dipping reflector sequences observed at volcanic rifted margins (e.g., Planke and Eldholm, 1994; White et al., 2008). Tilting of basaltic flows toward the rift axis has been observed in other areas of Afar (e.g., Tendaho Graben and Manda Hararo Rift, Acocella et al., 2008; southeastern Afar margin in Djibouti, Le Gall et al., 2011), with structural features displaying striking similarities with those of coastal flexures creating seaward-dipping reflector sequences at volcanic rifted margins (Le Gall et al., 2011). Our study therefore provides important constraints on the initiation of seaward-dipping reflector sequences at magmatic margins.

Model Setup

For the analysis of the vertical deflection due to a surface or subsurface load on an elastic plate overlying an inviscid substratum, we have used the numerical model by Li et al. (2004). The model is based on a high-order discretization with multigrid technique to solve the differential equation governing plate deflection, which is given by:
where D is the flexural rigidity of an elastic plate, w(x,y) is the vertical deflection of the plate, ρm is the density of mantle, ρfill is the density of the material occupying the space created by the flexural subsidence, and q(x,y) is the net force per unit area exerted by the sum of applied loads: q(x,y) = ρloadghload(x,y), where ρload and hload are the density and thickness of the applied loads, respectively, and g is gravitational acceleration. The flexural rigidity D defines the resistance to bending of a continuous elastic plate:
where E is Young’s modulus, ν is Poisson’s ratio, and Te is the effective elastic thickness of the plate; Te represents the cumulative strength of a multilayered lithosphere to loading over the 104–105 yr time scales of isostasy. The model simulates reversible and instantaneous (purely elastic) deformation once the load is applied.

Although both the Afar Depression and the northern MER are highly faulted, the spatial dimensions of the magma intrusion zones and the geological characteristics of rifting justify the use of a thin, continuous plate model. Current deformation is largely accommodated at the rift axis, where the thinned lithosphere is strongly modified by the extensive magma intrusion and extension occurs through a combination of magmatic intrusion and faulting. The zones of active deformation show small offset faults that may have formed above dike intrusions (Rowland et al., 2007), and there is no structural or stratigraphic information for large offset normal faults within the narrow magmatic segments (e.g., Wolfenden et al., 2005; Rowland et al., 2007). The time and length scales of active and time-averaged deformation indicate that localized dike intrusion is more important in accommodating crustal strain over the past ∼1 m.y. (e.g., Ebinger and Casey, 2001; Wolfenden et al., 2005; Nobile et al., 2012; Ebinger et al., 2013). The border faults, both in Afar and the (northern) MER, are relatively inactive on the basis of seismic and field evidence, so that the faults are effectively locked at rift margins (e.g., Keir et al., 2006).

The three-dimensional (3D) model represents a plate with 250 × 250 cells, with a mesh size of 5 × 5 km2; the load is applied on top of the continuous elastic plate (see Fig. 3). No initial topography of the bending plate (e.g., differences in elevation between the rift floor and the surrounding plateaus) was included in the models; this topography would simply be superimposed on the predicted plate bending. The following modeling results thus show the deformation as a result of the igneous load. A similar approach was used by Beutel et al. (2010) who studied deflection of two en echelon, rift segment–scale magma bodies (see discussion following).

The models are constrained by geophysical data of the tectonically active rift valleys of the MER zone. The values of Te were varied between 4 km and 15 km, according to the northward reduction in Te constrained by the gravity and topography coherence studies of Ebinger and Hayward (1996) (Fig. 1). Values of Young’s modulus (70 GPa) and Poisson’s ratio (0.30) were taken from Beutel et al. (2010). For the typical dimensions and densities of the axial intrusion we used constraints from controlled and passive source seismic tomography, models of magnetotelluric data, and 2D and 3D models of gravity data (Keranen et al., 2004; Tiberi et al., 2005; Cornwell et al., 2006; Whaler and Hautot, 2006; Daly et al., 2008) across axial magmatic segments in Ethiopia (Fig. 2). We used rectangular intrusions with 20 km width, 20–40 km length, and 10–16 km thickness, and a density contrast between the crust and the mafic intrusion of 210 kg/m3 (Fig. 2C; Cornwell et al., 2006). The influence of the additional load imposed by synrift sedimentation was also considered: we used a density of 2380 kg/m3 for the sedimentary rocks (Cornwell et al., 2006), which were considered to completely fill the depression resulting from the intrusion-related flexure. In addition, in some models we considered the external load imposed by the presence of volcanic edifices, whose dimensions were based on volcanoes from the northern MER (such as Boseti: basal diameter 20 km, height 750 m, and volcanic rock density 2700–2900 kg/m3).


We explored the deflection of the lithosphere induced by internal and external loading with varying plate thickness and size of magma intrusion both with and without synrift sedimentation (Figs. 4 and 5). The calculations show that the elastic plate bends in an ∼100–150-km-wide region around the mafic intrusion zone with the amplitude of plate bending at a maximum above the center of the intrusions (Figs. 4A and 5). At the margins of the subsiding region, minor uplift (a few meters maximum) is observed, giving rise to a small forebulge as shown in Figure 5. The magnitude of the bending is primarily dependent on Te (i.e., flexural thickness), the presence of synrift sediments, and the dimensions (i.e., volume) of the intrusions (Figs. 4B–4D).

The main control is exerted by variations in Te, the decrease of which induces a prominent, nonlinear increase in plate bending: from ∼200 m at Te of 15 km, to ∼400 m at Te of 10 km, and to ∼1000 m at Te of 5 km (Fig. 4B). At all values of Te, additional loading from synrift sediments approximately doubles the amount of subsidence (Fig. 4B). Increasing the length (Fig. 4C) or the thickness (Fig. 4D) of the intrusion also increases the subsidence. Our calculations demonstrate that the load of a volcanic edifice increases the amplitude of the plate bending, but its contribution beyond that created by the internal load is limited (typically <10%–20%), with the effect more pronounced at lower Te and when the size of the magma intrusion is smaller (i.e., decreasing the internal load; Figs. 4B–4D). In addition, the presence of a volcanic edifice has a more limited influence on plate bending than that of the sedimentary infill, despite volcanic rocks having a higher density (2700 kg/m3) than sediments (2380 kg/m3), a result related to larger volume of the subsiding region filled by sedimentary rocks than that of the volcanic edifice. Comparison of plate bending in the case illustrated here of a single mafic intrusion with the case of two en echelon bodies (Beutel et al., 2010) shows no significant variation in subsidence, with differences never exceeding ∼50 m (Fig. 6). Additional tests made with variations in density and height of the surface volcano suggest no significant influence of these parameters on plate bending.

The relation between the maximum subsidence w0 and Te is highly nonlinear and can be fit by a power function in the form w0 = A TeB, where A and B are parameters varying in the range 2.8–8.7 × 103 and 1.1–1.5, respectively, as a function of the variable load (presence of volcano, sediment infill, and length and thickness of the intrusion). This relation is of the same form as the dependence in the case of a line load (w0 = q0α3/8D; i.e., w0A TeB, where q0, α, and D are load, flexural parameter, and flexural rigidity, respectively, and B = 3/4; e.g., Turcotte and Schubert, 2002).

Data from continental rifts and rifted margins worldwide indicate that significant extension can be accommodated by magma intrusion prior to the development of a new ocean basin (e.g., Ebinger and Casey, 2001; Thybo and Nielsen, 2009; Bastow and Keir, 2011). The emplacement of large volumes of magma into the lithosphere (and eruption at the surface) modifies its density profile and transfers heat to the plate, affecting the patterns of long-term subsidence predicted by theoretical models of stretching and thinning of the continental lithosphere (e.g., Buck, 2004; Daniels et al., 2014). In particular, intrusion of large magma volumes in the lithosphere, as observed in the MER, may cause significant internal and surface loading resulting in plate flexure and rift subsidence (Beutel et al., 2010).

Despite some inherent limitations (e.g., constant flexural rigidity and Te throughout the model domain, simplified shape of the intrusions), our simple flexural models have allowed quantification of this process, significantly extending the previous state of stress studies by Beutel et al. (2010). Although not included in our continuous plate models of intruded loads and surface (volcano) loads, analytical and numerical models provide insights into the role of in-plane tensional stresses on the shapes of the plate deflection beneath the volcanic loads (e.g., ten Brink, 1991; Ebinger and Hayward, 1996). Assuming a tensional stress of 2 × 1012 Nm–1 (e.g., Kusznir et al., 1991) applied to an 800-m-high line load on a plate with Te of 10 km would reduce the maximum subsidence by ∼12% and broaden the depression by ∼25%. Our results, therefore, are probably overestimates of the maximum subsidence.

The results show a nonlinear increase in amount of subsidence when Te decreases, with the implication that as the plate weakens during progressive thinning and heating during rifting, loading-induced subsidence will become an increasingly important mechanism driving rift evolution. In the MER, where Te is typically >8–10 km, the effect of the large mafic intrusions imaged by geophysical data is relatively small (Beutel et al., 2010). However, thickening of the youngest sediments above the axial intrusions is suggested by the geometry of the low-velocity upper crustal layers imaged in seismic tomography by Keranen et al. (2004; see Fig. 2B): the uppermost units thicken dramatically directly above the intrusions, in contrast to deeper units that have a more uniform thickness across the rift. In addition, bending of the rift floor is visible in topographic profiles across magmatic segments, such as along a profile from the Asela-Langano boundary fault system into the Aluto-Gedemsa magmatic segment (Pizzi et al., 2006). Although the overall riftward bending may result from initial (Miocene) phases of basin downwarping (e.g., Zanettin and Justin-Visentin, 1975; Kazmin et al., 1980), the architecture of recent (Pliocene–Pleistocene) rift-related sediments (e.g., riftward tilting of Pleistocene–Holocene alluvial fans) at the eastern (Asela-Langano) margin (see Agostini et al., 2011) points toward a recent tilting likely induced by axial magma intrusion (see also Pizzi et al., 2006; Fig. 7A).

An overall bending of the rift toward axial zones of intensive dike intrusion and magmatic construction has been observed in the southern Red Sea rift (Fig. 7B; Wolfenden et al., 2005), as well as in other regions of central-southern Afar (e.g., Tendaho Graben and Manda Hararo Rift, Acocella et al., 2008; southeastern Afar margin in Djibouti, Le Gall et al., 2011), creating an axis-dipping wedge of primarily basaltic flows analogous to seaward-dipping reflector sequences imaged on volcanic margins worldwide (e.g., Eldholm et al., 1989; Planke and Eldholm, 1994). Intrusion-induced subsidence is expected to increase in northern Afar, where Te decreases significantly to ∼5 km in the Danakil Depression (Ebinger and Hayward, 1996). In this area, the basin floor reaches below sea level and is dominated by young evaporite-rich basins bisected by the basalt-rich Pleistocene–Holocene Erta Ale volcanic range. Here, an active axial intrusion shallower than ∼10 km depth is modeled by InSAR (interferometric synthetic aperture radar) data (e.g., Amelung et al., 2000; Pagli et al., 2012; Nobile et al., 2012). Limited borehole data suggest that the basin is filled by 3 km or more of interbedded Pliocene–Holocene evaporites and basalt flows (Hutchinson and Engels, 1972), suggesting that rapid subsidence of the region toward and below sea level accelerated during the past 5 m.y. Our models suggest that as much as 1 km of this subsidence (i.e., ∼33%) is likely to be caused by the intrusion-related loading and plate flexure, which also accounts for the thickness of the basin fill. The remainder is caused by the response of the lithosphere to faulting and lithospheric thinning, as well as localized subsidence above dikes as an elastic response directly above the mode I opening at depth.

Loading of the lithosphere has been shown to create flexural stresses capable of controlling processes such as magma ascent and emplacement during extension and the spacing of volcanoes in rift settings; the load imposed by volcanic edifices has been considered the major parameter contributing to these flexural stresses (ten Brink, 1991). Our results, however, suggest that the contribution to subsidence by the load of a volcanic edifice is relatively small compared to the intrusion-induced subsidence. Our analysis indicates that volumetrically large crustal intrusions, rather than volcanoes, may have a strong control on rift morphology over time. The base of the bent plate will be in extension, facilitating the rise of magma, and serving to maintain the along-axis magma intrusion zones.

Figure 8 summarizes the general consequences of flexure caused by internal loads. After the initial phases of tectonic rifting, which may create asymmetric basins with depocenters in correspondence to the main boundary faults, strain migrates to axial zones of localized faulting and magma intrusion, where significant subsidence is initiated when large volumes of magma intrude the lithosphere and accommodate extension. Intrusion-related subsidence becomes an increasingly important mechanism driving increasing rift subsidence as the plate weakens during progressive rifting leading to the formation of new columns of hot weak lithosphere at nascent mid-ocean ridges. In the final stages of breakup, when thinning of the heavily intruded, weakened plate induces a pulse of decompression melting and a significant increase in extrusive magmatism (see also Nielsen and Hopper, 2004; Bastow and Keir, 2011), the axial intrusion-induced subsidence serves to progressively reorient the basalt sequences and interbedded sediments to dip toward the future breakup boundary. Consequently, subsidence caused by high-density crust at zones of intrusion is a factor in the reorientation of originally subhorizontal flows to form the seaward-dipping reflector sequences at the newly formed magmatic rifted margins as the rift breaks the continent and the rifted margins continue to subside below sea level while the heat transferred to the plate during rifting gradually dissipates (e.g., Wolfenden et al., 2005; Beutel et al., 2010; Bastow and Keir, 2011; Keir et al., 2013). Thus, this process may represent a contributing factor to the development of the seaward-dipping reflector sequences, which are normally interpreted as either resulting from progressive oceanward flexuring of the lavas due to differential loading by the newly erupted flows (isostatic models; e.g., Palmason, 1980; Mutter et al., 1982; Mutter, 1985; Eldholm, 1991) or tectonic features associated with fault development (e.g., fault-related flexures controlled by major continentward-dipping normal faults; e.g., Nielsen, 1975; Tard et al., 1991; Barton and White, 1997; Geoffroy et al., 1998; Geoffroy, 2005; Quirk et al., 2014).

We thank Joel Ruch, an anonymous reviewer, and Associate Editor Eleonora Rivalta for the detailed comments that helped to improve the manuscript. We also thank Cindy Ebinger for discussions and comments on earlier versions of this work. Agostini thanks David Coblentz for his support during his stay at the Los Alamos National Laboratory. Agostini and Corti thank Giordano Montegrossi for discussions. Bastow acknowledges support from the Leverhulme Trust. Corti and Keir were supported by CNR (Consiglio Nazionale delle Ricerche) Short Term Mobility grants 2007 and 2013. Keir also acknowledges support from Natural Environment Research Council grant NE/L013932/1. Ranalli’s participation in this project was partly supported by Carleton University.