Mechanism for Deep Crustal Seismicity: Insight From Modeling of Deformation Processes at the Main Ethiopian Rift

We combine numerical modeling of lithospheric extension with analysis of seismic moment release and earthquake b‐value in order to elucidate the mechanism for deep crustal seismicity and seismic swarms in the Main Ethiopian Rift (MER). We run 2‐D numerical simulations of lithospheric deformation calibrated by appropriate rheology and extensional history of the MER to simulate migration of deformation from mid‐Miocene border faults to ∼30 km wide zone of Pliocene to recent rift floor faults. While currently the highest strain rate is localized in a narrow zone within the rift axis, brittle strain has been accumulated in a wide region of the rift. The magnitude of deviatoric stress shows strong variation with depth. The uppermost crust deforms with maximum stress of 80 MPa, at 8–14 km depth stress sharply decreases to 10 MPa and then increases to a maximum of 160 MPa at ∼18 km depth. These two peaks at which the crust deforms with maximum stress of 80 MPa or above correspond to peaks in the seismic moment release. Correspondingly, the drop in stress at 8–14 km correlates to a low in seismic moment release. At this depth range, the crust is weaker and deformation is mainly accommodated in a ductile manner. We therefore see a good correlation between depths at which the crust is strong and elevated seismic deformation, while regions where the crust is weaker deform more aseismically. Overall, the bimodal depth distribution of seismic moment release is best explained by the rheology of the deforming crust.


Introduction
The depth distribution of seismicity in the East African Rift System (EARS) generally shows earthquakes only in the upper crust or alternatively a clear bimodal pattern with peaks in the upper crust and then in either the lower crust or upper mantle (e.g., Craig et al., 2011;Yang & Chen, 2010). While the upper crustal earthquakes are consistent with deformation of crust dominated by a quartz-rich composition typical of continents, a range of mechanisms are proposed to explain the deeper earthquakes. For example, Seno and Saito (1994) proposed that lower crustal earthquakes occur due to high pore fluid pressure from fluids migrating from the upper mantle. This mechanism is recently used to explain most of deep crustal seismicity in the Tanzanian and Kenyan rifts (e.g., Lindenfeld et al., 2012). Alternatively, in the Tanganyika and Main Ethiopian Rift (MER) (Figure 1), others propose that deep crustal seismicity is driven by brittle faults penetrating the entire crust (Lavayssiére et al., 2019;Lloyd et al., 2018). The tectonic faulting at these depths is enabled by the presence of strong crust, either because the entire lithosphere is anomalously thick and cold or because the lower crust has a mafic composition and is therefore anomalously strong. The bimodal distribution of earthquakes also suggests that the deviatoric stress in parts of the middle crust must be lower than that of the upper crust and of the deep crust/upper mantle (Yang & Chen, 2010).
In this manuscript we explore the deformation of lithosphere and crustal stress patterns to explain variations in the amount of seismicity in the crust. To this end, we use 2-D high-resolution viscoplastic numerical experiments. The model results are then compared with seismic moment release and earthquake b-values.

Tectonics of the Northern MER
Geodetic (Bendick et al., 2006;Bilham et al., 1999), structural (Kurz et al., 2007), and modeling (Corti, 2008) studies in MER indicate the migration of deformation from border faults to ∼60 km long, 20-30 km wide Quaternary to recent "magmatic segments" within the central rift floor (e.g., Ebinger & Casey, 2001). These magmatic segments are now the locus of active faulting, seismicity, and volcanism (Keir et al., 2006;Greenfield et al., 2019). The crust beneath the magmatic segments has higher than normal P wave velocities (V p ), interpreted as evidence for cooled mafic intrusions (e.g., Daly et al., 2008;Keranen et al., 2004;Mackenzie et al., 2005). The V p is particularly high at lower crustal depths suggesting that the lower  Hofstetter and Beyth (2003) and Keir et al. (2006). GPS velocity vector shows the motion of Somalia with respect to Africa (Nubia) (Iaffaldano et al., 2014, similar to Saria et al., 2014, for the rigid zones). The red triangles are active volcanoes from the Smithsonian Institution, Global Volcanism Program (https: //volcano.si.edu). The open black box shows the location of Figure 2a, and the open white box on the inset map shows the location of the main map. Elevation data are taken from GEBCO database (https://www.gebco.net).
crust is heavily intruded (Mackenzie et al.,2005). Thermal metamorphism due to this cooled intrusion makes the surrounding crust strong (Lavecchia et al., 2016;Muluneh et al., 2018), and such crustal strength is required to explain the observed melt chemistry in the MER (Armitage et al., 2018). In addition, the mafic rock type of the intrusions is stronger than the surrounding more felsic continental crust (Beutel et al., 2010).
Analysis of earthquake catalog puts constraints on the depth distribution of seismicity and seismogenic nature of the crust and mantle. The EAGLE catalog for 2001-2003 recorded ∼2,000 earthquakes ( Figure 2a) with magnitudes M L less than 4.0 (Keir et al., 2006). This catalog is the best available so far for the region in terms of accuracy of depth and location of earthquakes. The catalog has an error of ∼2 and ∼0.6 km in hypocenter and x and y directions, respectively. Seismicity has a maximum depth of ∼32 km, but most of the earthquakes nucleate at depth of <18 km (Figures 2b and 2c).
The crustal thickness models from the EAGLE active source experiment using P wave seismic velocity model show significant variations along and across the MER (e.g., Mackenzie et al., 2005;Maguire et al., 2006). The crust is ∼45 km thick beneath the northwestern plateau and 38 km thick beneath the eastern plateau. In the rift, the crust is > 35 km thick in the central MER and thins northward to ∼28 km in the northern MER. The crust consists of a number of distinct layers, most of which are traceable from the plateaus into the rift. The upper crust is ∼28 km thick beneath the northwestern and eastern plateaus. This layer thins to ∼18 km beneath the northern MER. The lower crust ranges in thickness from 14 km beneath the eastern plateau to ∼10 km beneath the northwestern plateau (e.g., Keranen et al., 2004). Beneath the northern MER, the lower crust has a thickness of <10 km (Keranen et al., 2009;Mackenzie et al., 2005). The base of the crust beneath the northwestern plateau and the rift is underlain by anomalously high velocity layer with a thickness of ∼10 km, which is interpreted as heavily intruded lower crust (e.g., Mackenzie et al., 2005).

Numerical Model
Previous 3-D numerical models of the MER successfully explain magmatic segmentations (Beutel et al., 2010), rift propagation and linkage (Brune et al., 2017;Corti et al., 2019), and kinematic consequences during oblique rifting (Duclaux et al., 2019). In the present study, our main interest is to combine numerical modeling of lithospheric extension with the depth distribution of seismicity.
We follow previous modeling approaches (Brune et al., 2017;Corti et al., 2019) and use rheological and thermal parameters to model the evolution of deformation in the northern MER. We also estimate the deviatoric stress available to drive extension at different crustal depths.

Governing Equations
We construct a 2-D box setup using thermomechanical finite element code ASPECT v2.0.0-pre (Advanced Solver for Problems in the Earth's ConvecTion, Bangerth et al., 2018;Glerum et al., 2018;Heister et al., 2017;Kronbichler et al., 2012;Rose et al., 2017) to model extension of the MER. Our model is based on previous ASPECT setups aimed at modeling continental rift dynamics (Corti et al., 2019;Glerum et al., 2020;Naliboff et al., 2020). We solve the incompressible flow equations for conservation of momentum (Equation 1), mass (Equation 2), and energy (Equation 3) assuming an infinite Prandtl number: (1) where η is viscosity, _ ε is strain rate tensor, u is velocity vector, P is pressure, ρ is density, ρ is adiabatic reference density, g is gravitational acceleration, κ is thermal diffusivity, ν h is artificial diffusivity, C p is specific heat capacity, and H is heat production. Density, ρ, is given as ρ 0 (1− α(T−T 0 )), where ρ 0 is the reference density, α is thermal diffusivity, T is temperature, and T 0 is the reference temperature.
For each compositional field c i , an additional advection equation (Equation 4) is introduced to Equations 1-3. As these equations contain no natural diffusion, artificial diffusivity ν h is introduced to stabilize advection (Kronbichler et al., 2012) 10.1029/2020GC008935 Geochemistry, Geophysics, Geosystems

Model Setup and Boundary Conditions
The model geometry comprises a domain of 500 × 160 km in x (cross-rift) and y (depth) directions, respectively ( Figure 3). The resolution of our model varies between 1 km in the region of interest (at depths shallower than 50 km and between 125 and 375 km in x direction) and 2 km outside this area above 130 km depth. The remaining model region has a resolution of 4 km.
Our modeling approach uses constitutive relationships for viscous and plastic rheology. Viscous flow follows a power law model for diffusion and dislocation creep (Karato & Wu, 1993) (Equation 5): where A is preexponent, n is the power law index, m is the grain size exponent, d is grain size, Q is the activation energy, _ ε ef f is the effective deviatoric strain rate, V is activation volume, and R is the gas constant. In case of diffusion creep, n = 1 and m > 0, while for dislocation creep n > 1 and m=0.
We simultaneously apply both the dislocation (η disl ) and diffusion (η diff ) creeps (van der Berg et al., 1993) for the viscous rheology (Equation 6) where η comp is composite of both viscous creep mechanisms.
Brittle/plastic rheology is implemented by rescaling the effective viscosity η pl ef f in such a way that the stress does not exceed the yield stress σ y (Equation 7) derived by the Drucker-Prager yield criterion (Equation 8).
where σ y is given by where P is the total pressure, ϕ is angle of internal friction, and C is the cohesion. When the viscous stress exceeds the plastic yield stress, the effective viscoplastic viscosity will be chosen as the plastic viscosity. Otherwise, the composite viscosity is used. The effective viscoplastic viscosity is fixed between lower and upper cutoff values of 10 19 and 10 24 Pa s, respectively.
The upper crust is modeled using a wet quartzite flow law (Rutter & Brodie, 2004) and the mafic lower crust and weak seed are represented by wet anorthite (Rybacki et al., 2006). Flow laws of dry and wet olivine (Hirth & Kohlstedt, 2003) represent lithospheric mantle and asthenospheric mantle, respectively. The undeformed crust has a thickness of 38,m, with 25 and 13 km thick upper and lower crust, respectively. The thickness of these layers is based on geophysical observations from the southeastern plateau (e.g., Keranen et al., 2009). Since the southeastern plateau has not been modified by magmatic underplating, unlike the northwestern margin, the crustal thickness represents the thickness prior to the opening of the MER. The initial lithosphere has a thickness of 120 km based on lithospheric thickness from the undeformed part of Africa (e.g., Fishwick, 2004). We set the horizontal component of velocity on the left and right boundaries to  (Keir et al., 2006) scaled to their magnitudes and colored with hypocentral depth. The polygons bound earthquakes that occur in the plateau and rift floor. Thin gray lines indicate faults in the rift (e.g., Boccaletti et al., 1998;Kazmin & Berhe, 1981). AA′ and BB′ profile lines show earthquake hypocenters depicted in the middle and bottom figures, respectively. The profiles project earthquakes within 20 km of the line in the northern and central MER. (b) seismicity continues to ∼32 km whereas in (c) seismicity is concentrated in the upper ∼12 km (i.e., by considering the hypocentral uncertainty of ±2 km). V ext /2 where V ext is the full opening rate. This material outflow is compensated by normal inflow V b through the model base. The tangential stress is 0 along the vertical and bottom boundaries, allowing for tangential motion. The top surface is a true free surface.
The surface temperature is kept constant at 0°C. The temperature at the bottom boundary is also fixed, at 1345°C. Lateral boundaries are thermally isolated. We prescribe an initial steady-state continental geotherm in the lithosphere, with an LAB temperature of 1300°and upper and lower crustal heat productions of 1 and 0.1 μW/m 3 , respectively. In the asthenosphere an initial adiabatic temperature profile is prescribed based on an adiabatic surface temperature of ∼1284°C. Boundary condition for composition is fixed to initial composition along the top and bottom boundary.
A Gaussian-shaped lithosphere-asthenosphere thermal and compositional perturbation with an amplitude and standard deviation of 5 and 10 km, respectively, helps localize deformation at the center of the model. We randomly distributed heterogeneities around the rift axis throughout the whole depth of the model representing the preweakened lithosphere (Dyksterhuis et al., 2007). The magnitudes of these random perturbations of initial plastic strain follow a broad Gaussian envelope of 100 km standard deviation and a strain amplitude of 0.4. Strain weakening is implemented as a linear decrease in friction angle from 30°to 9°for brittle strain between 0 and 1. For plastic strains larger than 1, it remains constant at 9°. The choice of friction angle is based on the observations of the strength of the crust and faults in the MER (Muluneh et al., 2018). Plastic strain, defined as the second invariant of the deviatoric strain rate times the time step, is tracked if the viscous stress exceeds the yield stress. For simplicity, cohesion is kept constant at a value of 20 MPa. All the rheological and thermal parameters as well as the ASPECT input file are given as supporting information.
Numerical modeling of lithospheric extension can be conducted either through stress or kinematic boundary conditions (Brune et al., 2016). Since the opening of the MER is well documented by a number of geodetic and plate kinematic observations showing a relatively constant extension rate, we use the kinematic boundary condition. We run the model using a full opening rate of 6 mm/yr (Iaffaldano et al., 2014), leading to 66 km opening during the last 11 Myr. In order to assess the effect of different velocity boundary conditions, we also run the model using a 4 mm/yr constant opening rate (DeMets & Merkouriev, 2016) corresponding to the long-term average during the last 16.5 Myr. Both models show similar deformation style at the end of model run. This implies that despite the difference in velocity boundary conditions, the deformation style is the same.

Model Robustness and Limitations
In order to validate the robustness of our results, we conducted additional models with viscosity cutoffs of 10 18 and 10 24 Pa s, which lead to almost indistinguishable results. We use time step of 5,000 years when applying the velocity of 6 mm/yr. Increasing or decreasing the time step by a factor of 2 does not change the results. Including linear cohesion softening did not significantly change the results.
Similar to previous modeling approaches, our experiment does not include the effect of a mantle plume, because the size of the studied plate boundary is in any case small compared to the extent of the African superplume (e.g., Nyblade et al., 2000). In addition, erosion, sedimentation, elasticity, and magmatic . Model geometry and boundary conditions. V ext is horizontal, full spreading rate of opening that drives extension and V b at the model base balances material outflow. A mechanically weak heterogeneity at the base of the lithosphere helps localize deformation at the center. We run the model for 11 Myr leading to 66 km extension, which is within the range of predicted extension for the northern MER. The dashed black box shows the subset (200 km × 60 km) of total domain of the model discussed in Figures 4 and 5. 10.1029/2020GC008935 Geochemistry, Geophysics, Geosystems underplating have not been implemented. Deformation in the MER is driven by constant, oblique kinematics since the onset of rifting at 11 Ma (e.g., Corti, 2008;DeMets & Merkouriev, 2016;Iaffaldano et al., 2014). We cannot investigate the role of obliquity (Brune, 2014;Brune et al., 2012;Duclaux et al., 2019) or the impact of along-strike mantle flow (Mondy et al., 2018) due to the 2-D nature of our experiment, but first-order deformation aspects are nevertheless expected to be very well represented.

Model Results
Our numerical experiment provides insight into deformation processes since the onset of rifting in MER. We discuss several outputs that are directly related to the seismogenic nature and strength of the crust.  (Ebinger & Casey, 2001). In general, our model evolves from asymmetric extension at the beginning of rifting to mostly symmetric rifting, generating a so-called asymmetric-symmetric pattern (Huismans & Beaumont, 2003). Figure 4a shows that the highest brittle strain rate (∼2 × 10 −14 s −1 ) is accommodated by narrow conjugate faults distributed within a 30 km wide zone in the upper crust. These narrow, high strain rate zones accommodate the ongoing extension in the rift. Toward the base of the upper crust, high strain rate is widely distributed without forming narrow shear zones.
Outside the rift zone, high strain rate is observed at the base of the upper crust enhanced by shear deformation due to the rheological contrast between the upper and lower crust. The thickness of the upper crust remains constant throughout model evolution, while the lower crust has thinned significantly within the width of the rift zone leading to broadly distributed strain. Figure 4b indicates that time-integrated strain concentrates on the major fault on the left side since model initiation. As extension proceeds, the high-angle fault rotates to very low dip angle in the lower crust where extension is broadly distributed. High-angle shear zone penetrates to the base of lower crust and accumulates strain at depth by forming ramp-flat-ramp structure (Figure 4b). In contrast, no significant strain accumulates on the right side of the rift until ∼6 Myr model time. Soon after 6 Myr model time, faults on the right side start accommodating brittle extension together with faults on the left side until 11 Myr model time. In the rift floor and right side of the model, deformation is mainly taken up by multiple, high-angle faults with relatively small topographic offset.

Viscosity
Viscosity plots show the low-viscosity asthenospheric material and the presence of competent layers (Figure 4c). Active shear zones in the upper crust and the whole lower crust are characterized by low effective viscosity (<10 21 Pa s). There is a slight deflection of the asthenospheric upwelling to the right guided by a high strain shear zone from the crust that forms a ramp-flat-ramp structure. The pattern of deflected asthenospheric upwelling is apparently similar to tomographic (Bastow et al., 2005) and magnetotelluric (e.g., Hübert et al., 2018) studies that show the offset of low-velocity anomalies away from the magmatic segments, although the clear pattern of the offset is difficult to explain due to the 2-D nature of our numerical experiment. Similar deflection of asthenospheric upwelling and rift axis jump is also observed in previous numerical modeling experiments (Huismans & Beaumont, 2003;Tetreault & Buiter, 2018).
Available estimates of the vertically averaged viscosity in the EARS range from 10 19.6 to 10 23 Pa s . Bendick et al. (2006) computed the viscoelastic relaxation effect following the 1993 dyking event in the northern MER. They found a best fit to observed displacement using a 15 km thick elastic/brittle crust over a viscous lower crust with viscosity of 1.125 × 10 18 Pa s. Our lower crustal viscosity estimate is higher 10.1029/2020GC008935 Geochemistry, Geophysics, Geosystems than that of Bendick et al. (2006). The difference could be due to the timescale of rapid dike injection and our modeling approach that considers the long-term deformation. High viscosity in the lowermost part of the lower crust hints at possible brittle deformation and hence might explain lower crustal seismicity in MER (Lloyd et al., 2018).

Brittle and Ductile Layers
The occurrence of brittle deformation is assessed at each time step based on the Drucker-Prager yield criterion. Figure 4d shows the presence of brittle layers in the upper and lower crust and upper mantle. At the end of model run (11 Myr), the brittle layer in the lower crust is consumed and only ∼10 km thick brittle upper crust and ∼10 km thick upper mantle deform in brittle manner. Figure 4d also shows the occurrence of a brittle layer below 600°C isotherm due to olivine dominated rheology of the upper mantle (McKenzie et al., 2005). Our numerical results at 11 Myr model time (Figure 4d) also show that the crust is seismogenic both in the uppermost crust and in the deeper crust. The lower part of the upper crust and lower crust deforms in ductile manner.

Deviatoric Stress
We explore the stress state arising from rift evolution and estimate the deviatoric stress that includes both the far-field and local components. Our model result offers a direct access to the stress tensor in 2-D. The second invariant of the deviatoric stress in 2-D is calculated using Equation 9: where τ xx , τ yy , and τ xy are the components of the deviatoric stress tensor.

10.1029/2020GC008935
Geochemistry, Geophysics, Geosystems section does not affect the discussion on the variation of stress with depth. The magnitude of stress increases from 20 to 80 MPa in the upper most crust. From 8 to 14 km depth, the stress decreases sharply to ∼10 MPa and then increases and reaches 160 MPa at ∼18 km before going back to 10 MPa at depth ≥30 km ( Figure 5).

Moment Magnitude and Earthquake b-Value
We used the earthquakes from the EAGLE catalog (Keir et al., 2006) to estimate the seismic moment release and b-values at different crustal depths. Seismic moment release (M o ) is the energy released by an individual earthquake and is a function of the earthquake's magnitude (Kanamori, 1983). The b-value describes the magnitude-frequency distribution of earthquakes, whereby a smaller value indicates a higher proportion of large earthquakes with respect to small earthquakes, and vice versa. Variations in the observed b-value have been attributed to changes in stress conditions and/or rock heterogeneity, and fluid diffusion (Marzocchi et al., 2020).
The calculation of seismic moment release (Figures 6a and 6b) indicates that the majority of seismic moment release occurs in the upper crust above 8 km. A small seismic moment is released at depths of 8-14 km in the rift (Figure 6b). We see a second peak in seismic moment release in the lower part of the upper crust at ∼16 km depth. Below this depth, seismic moment release is minimal. In order to obtain insight into the deformation style at different crustal depths and provide additional constrains to our modeling results, we estimate the b-value at depths ≤8, 8-14, and ≥14 km. These depth ranges are selected based on the depth distribution of seismicity, occurrence of low-magnitude seismic swarms (Keir et al., 2006), and sharp increase in shear strength of faults (Muluneh et al., 2018). We calculate the magnitude of completeness (M C ) using the maximum curvature method (Wiemer & Wyss, 2000) that takes the magnitude "bin" with the highest frequency and adds 0.2 magnitude units. We then use the maximum-likelihood calculation (Aki, 1965) and bootstrap analysis (Pickering et al., 1997) to calculate b-values and associated errors, respectively. The bootstrap analysis creates 10,000 data sets by randomly sampling the earthquake catalog and allowing for duplication (supporting information Figure S1), a b-value is calculated for each data set and we report the mean and standard deviation of these values.
At depths shallower than 8 km, the b-value is 0.75 ± 0.09 (466 earthquakes) and it increases to 0.91 ± 0.06 (1059 earthquakes) at 8-14 km depth. The b-value decreases to 0.82 ± 0.11 (396 earthquakes) below 14 km. The higher b-value at 8-14 km coincides with an increased number of earthquakes at ∼12 km (Keir et al., 2006), but due to the vast majority of earthquakes being relatively small in magnitude, there is a small seismic moment release.

Comparison and Discussion of Results
We compare our numerical model results with observed deformation styles in the northern MER. Then, the brittleness of layers and stress from model results are compared with depth distribution of seismicity, bvalue, and crustal strength. We assume that no seismicity occurs where the applied stress is less than the yield strength (Scholz, 1988) and that geologically reasonable viscoplastic deformation represents the earthquake cycle that happens on a much shorter time scale.
We note that the crustal thickness from our model output differs from observations in the northern MER (Figures 4 and 5). If we include the 5 km thick magmatic underplating (Maguire et al., 2006), the high deviatroic stress and brittle region will be within the crust, not the upper mantle; thus, our model does not necessarily predict upper-mantle earthquakes. The thickness of the lower crust in our model represents only the upper ∼5 km of the lower crustal thickness inferred from geophysical observations (e.g., Maguire et al., 2006). Since most of the seismicity concentrates in the upper crust, the difference does not affect the discussion on the depth distribution of seismicity and model results.

Deformation Style
Our 2-D model successfully describes temporal deformation pattern in MER: including migration of high strain rate from border faults to 20-30 km wide segments, broadly distributed ductile extension in the lower part of the upper crust and in the lower crust and high accumulation of brittle strain on the border fault. The migration of high strain rate from the border faults to rift axis faults and magnitude of strain rate are comparable to GPS (e.g., Kogan et al., 2012), structural (Agostini et al., 2011), and modeling (Corti, 2008) observations in the northern MER. Strain migration occurs without the input from magma (Figure 4a). Analog modeling experiment shows that migration of strain from the border faults to magmatic segments does not necessarily require the input from magma but can instead be caused by lithospheric thinning due to constant extensional kinematics (Corti, 2008).
The model results also indicate that brittle deformation occurs in two regions beneath the rift, in the uppermost crust and at the deeper crust. Deformation in the lowermost crust occurs in ∼70 km wide zone (Figure 4a). Although our model result shows broadly distributed deformation in the lowermost crust, similar to field geophysical observations (Keranen et al., 2009;Kogan et al., 2012), the width of the deforming zone differs quite significantly. For example, Keranen et al. (2009) suggest that deformation in the lower part of the upper crust and lower crust is distributed in ∼400 km wide zone. On the other hand, Kogan et al. (2012) argue that deformation is localized in ∼85 km wide zone with some deformation extending outside the structural MER. In order to explain the widespread deformation, Keranen et al. (2009) argue that the lower crust must be weak. This contradicts with more recent modeling study suggesting that lower crust in MER must be strong in order to explain melt chemistry (Armitage et al., 2018).

Deviatoric Stress and Seismic Moment Release
In the EARS the magnitude of deviatoric stress is in a range of 10-20 MPa (Mahatsente & Coblentz, 2015;Stamps et al., 2010), and hence, it appears to be too small to cause brittle failure on seismogenic faults (Craig et al., 2011;Scholz, 2002). Mahatsente and Coblentz (2015) argue that the ridge-push force from the oceanic part of Nubia (Africa)-Somalia plate is less than the integrated strength of the African plate and hence additional forces are required to deform the plate. The above studies report magnitude of vertically averaged deviatoric stress over the thickness of the lithosphere. Since both magnitude of deviatoric stress and strength vary with depth, detailed analysis on those parameters is crucial to understand the mechanism for the peak in midcrustal seismicity. Figure 6. The sum of seismic moment release in 1 km depth bins estimated using the relationship M 0 ¼ 10 3 2 Mw þ 9:1 Nm for plateau (a) and rift (b) shown by regions bounded by thick polygon. We make assumption that M L =M w . The red bars indicate the seismic moment released by low-magnitude earthquake swarms (Keir et al., 2006).
The second invariant of deviatoric stress in our models shows strong variation with depth. Figure 7 offers a clear representation of deviatoric stress and magnitude of seismic moment release with depth. The upper crust deforms with stress ranging from 20-80 MPa. This stress is higher than the shear strength of the upper crust, which is ∼20 MPa (Muluneh et al., 2018). A maximum of ∼1 × 10 15 Nm moment magnitude is released at 7 km depth. Between 8 and 14 km depths, the magnitude of stress decreases sharply to ∼10 MPa. The reduction in deviatoric stress is consistent with the observed occurrence of low-magnitude earthquake swarms with small seismic moment release and higher b-value ( Figure 7). Alternatively, we propose that earthquakes in the 8-14 km depth range, characterized by large numbers of small earthquakes, may be triggered by fluid release, potentially from cooling mafic intrusion (Keir et al., 2009). We observe a good Figure 7. Seismic moment release, variation of magnitude of deviatoric stress with depth at the center of the rift ( Figure  5) and earthquake b-value. Low stress region coincides with swarms of low-magnitude earthquake in the rift zone, small seismic moment release (red bars), and relatively higher b-value. Note that this layer has a shear strength ranging from ∼40 to ∼120 MPa (Muluneh et al., 2018, similar to global compilation of shear stresses from major boundary faults by Behr & Platt, 2014). The green lines indicate the depth ranges (≤8, 8-14, ≥14 km) at which b-values are estimated.

10.1029/2020GC008935
Geochemistry, Geophysics, Geosystems correlation between the peak in seismic moment release at ∼16 km depth and a peak to 160 MPa in the deviatoric stress in the lower part of the upper crust. Figure 8 summarizes the distribution of seismicity and faults in the MER. The combination of high accumulated brittle strain and stress together with high seismic activity points out that the Ankober border fault and surrounding regions are still active and accommodate the present-day opening of northern MER. Deformation at the rift center is accommodated by active fault zones.
Previous studies from other sectors of the EARS (e.g., Albaric et al., 2009) show that the depth distribution of seismicity can be fitted to different yield strength envelopes depending on tectonic settings. Several lines of evidences suggest that deformation style and seismicity in the MER vary along (e.g., Déprez et al., 2013;Muluneh et al., 2017) and across the rift (Keranen et al., 2009;Kogan et al., 2012). Future studies should address the role of contrasting rheologies and thermal properties between the plateau and the rift in controlling the depth distribution of seismicity (e.g., Chambers et al., 2019).

Conclusions
We present a detailed numerical modeling study of lithospheric extension and deviatoric stress state in the MER. Model results are compared with depth distribution of seismicity and seismic swarms, b-value, and seismic moment release in order to propose a mechanism for deep crustal earthquakes in the MER. We use a high-resolution, 2-D numerical experiment to model the evolution of deformation and stress using the most appropriate rheology for the lithosphere. Our model results successfully show the migration of deformation from border to rift center, similar to GPS and geophysical observations. Analysis of the deviatoric stress based on model results show that stress significantly varies with depth in the MER floor. The uppermost crust (i.e., ≤8 km) deforms with maximum stress of 80 MPa. The stress drops to 10 MPa at depth of 8-14 km and then increases to 160 MPa at ∼18 km. The peaks in stress correspond to peaks in seismic moment release in the MER crust. The low-stress depth range (i.e., 8-14 km) is characterized by ductile deformation, small seismic moment release, concentration of swarms of low-magnitude earthquakes, and higher b-value. We conclude that the bimodal depth distribution of seismic moment release in the MER is controlled by the rheology of the deforming crust.

Data Availability Statement
The EAGLE earthquake data used in this paper can be found in Keir et al. (2006). Figures were drafted using Generic Mapping Tools (Wessel & Smith, 1998). The ASPECT plugins used in the manuscript is available here (https://zenodo.org/record/3758145#.Xp28CFDRZTY; the DOI is 10.5281/zenodo.3758145.)