Earth and Planetary Science Letters Evidence for cross rift structural controls on deformation and seismicity at a continental rift caldera

In continental rifts structural heterogeneities, such as pre-existing faults and foliations, are thought to inﬂuence shallow crustal processes, particularly the formation of rift faults, magma reservoirs and surface volcanism. We focus on the Corbetti caldera, in the southern central Main Ethiopian Rift. We measure the surface deformation between 22nd June 2007 and 25th March 2009 using ALOS and ENVISAT SAR interferograms and observe a semi-circular pattern of deformation bounded by a sharp linear feature cross-cutting the caldera, coincident with the caldera long axis. The signal reverses in sign but is not seasonal: from June to December 2007 the region south of this structure moves upwards 3 cm relative to the north, while from December 2007 until November 2008 it subsides by 2 cm. Comparison of data taken from two different satellite look directions show that the displacement is primarily vertical. We discuss potential mechanisms and conclude that this deformation is associated with pressure changes within a shallow ( < 1 km) fault-bounded hydrothermal reservoir prior to the onset of a phase of caldera- wide uplift. Analysis of the distribution of post-caldera vents and cones inside the caldera shows their locations are statistically consistent with this fault structure, indicating that the fault has also controlled the migration of magma from a reservoir to the surface over tens of thousands of years. Spatial patterns of seismicity are consistent with a cross-rift structure that extents outside the caldera and to a depth of ∼ 30 km, and patterns of seismic anisotropy suggests stress partitioning occurs across the structure. We discuss the possible nature of this structure, and conclude that it is most likely associated with the Goba– Bonga lineament, which cross-cuts and pre-dates the current rift. Our observations show that pre-rift structures play an important role in magma transport and shallow hydrothermal processes, and therefore they should not be neglected when discussing these processes.


Introduction
The pathway taken by rising magma is influenced by local and regional stresses (e.g., Maccaferri et al., 2014) and lithological and/or rheological boundaries (e.g., Taisne and Tait, 2011). In old continental crust in particular, heterogeneities and structures, such as faults and lithological contrasts, are widespread, and they can strongly influence the location and geometry of magma reservoirs (e.g., Le Corvec et al., 2013). The roles of these competing factors have been demonstrated at many volcanoes. For example, Fig. 1. a) The Main Ethiopian Rift, with East Africa as inset. Red ellipses: rift calderas, scaled according to size and orientation of caldera rim fault after Wadge et al. (2016). Red triangles: non-caldera volcanoes. Black lines: intra-rift faults (Agostini et al., 2011). Spreading direction from Stamps et al. (2008). Yellow star: Addis Ababa b) Corbetti Caldera and surrounding region, showing the Wonji faults (red), Wendo Genet Scarp, Werensa Ridge and hypothesised Hawassa Caldera (purple line). The seismicity associated with the Wendo Genet Scarp is also shown (71 events with lateral uncertainty <10 km, recorded between January 2012 and January 2014) (Wilks, 2016). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) between this structure, the caldera long-axis, and the alignment of volcanic vents. Our observations are further supported by magnetotelluric interpretations of a large structure that cross-cuts the caldera and seismic data which identifies the structure to the east, where it extends down to ∼30 km through the crust. Seismic anisotropy measurements also indicate a concentration of stress along the structure cross-cutting the caldera, and show that stress perturbations are largest where the greatest surface deformation is observed. We discuss the potential sources which may have caused the deformation (magma, hydrothermal fluids, or meteoric water), and possible interpretations for this structure, including a pre-rift fault system, the edge of older solidified intrusions, or the rim of a Pleistocene caldera.

The East African Rift
The East African Rift system (EARS) is a ∼4,000 km long continental rift which defines the boundary between the Somalian and Nubian tectonic plates. Rifting occurs through a combination of magmatic and tectonic processes, and inherited structures and fabrics influence where and how this extension is accommodated (e.g., McConnell, 1972).
The Main Ethiopian Rift (MER) is the northernmost part of the EARS, and is an example of mature continental rifting (Chorow-icz, 2005). There are several Quaternary major silicic volcanic system in the MER, some of which have been observed deforming in recent decades: Corbetti, Bora, Haledebi and Aluto (Biggs et al., 2009;Hutchison et al., 2016). The most recent MER eruptions were at Tullu Moje (syn. Bora) in 1900, at Kone (syn. Gariboldi) in ∼1820, and in ∼1810 at Fantale (Wadge et al., 2016 and references therein). It has been suggested that the large silicic centres lie at the termini of magmatic segments (Ebinger and Casey, 2001;Keranen et al., 2004), where reduced extensional stresses facilitate longer residence times, favouring the development of silicic bodies through fractional crystallisation (Peccerillo, 2003;Hutchison, 2015). However, several centres lie along pre-rift faults, for example the elliptical calderas Kone, Gedemsa and Fentale (Acocella et al., 2002). An alternative mechanism for the formation of silicic centres is that reactivation of transtensional faults create regions of localised extension, focussing rising magma, and promoting magma reservoir formation (Acocella et al., 2002;Holohan et al., 2005). For example, near the Aluto caldera rhombic faulting of the border faults is associated with a pre-existing lithospheric weakness and sinistral oblique crustal shear (Boccaletti et al., 1998;Corti, 2009).

Corbetti caldera
The Corbetti caldera, in the southern MER, is the southernmost silicic centre along the MER (Fig. 1); further south the rift transitions to diffuse faulting and magmatism (Corti, 2009;Philippon et al., 2014). The caldera formed at 182 ± 18 ka, and is one of the largest in the EARS, measuring ∼10 by 15 km (Hutchison, 2015). The caldera scarp height is greatest to the west (∼200 m), and diminishes in height to the east, where it is unidentifiable. Corbetti is surrounded by agricultural land and is located within 15 km of two major population centres: Hawassa and Shashemene (∼400,000 people within 25 km 2 ).
There are two major centres of resurgent volcanism within the caldera: Urji (syn. Wendo Koshe) and Chabi, which have both erupted aphyric pantellerites. Urji, the western peak, has fall and flow pumice deposits, illustrative of high explosivity. Chabi, by contrast, is composed of obsidian flows, indicative of more effusive eruptions (Rapprich et al., 2016). Numerous thick fall deposits and obsidian flows indicate several eruptions have occurred at Corbetti since caldera formation. Tephra from the most recent Plinian eruption at Urji has been dated at 396 ± 38 BC (Rapprich et al., 2016), and at least four obsidian flows post-date this.
Magnetotelluric (MT) measurements and transient electromagnetic method soundings show a sharp resistivity gradient along the caldera long-axis (Gíslason et al., 2015); the north being more resistive than the south. This resistivity contrast is present in the upper 2 km, and >8 km below sea level. The shallow component of the feature is interpreted to reflect the vertical migration of hydrothermal fluids from depth, which then become entrained in the northwards local groundwater gradient. The deeper feature is attributed to a magma body (Gíslason et al., 2015).

Interferogram processing
InSAR (Interferometric Synthetic Aperture Radar) is a spacebased remote sensing technique, used to measure deformation of the Earth's surface (Massonnet and Feigl, 1998). InSAR uses the difference in the phase component of two radar images, acquired from approximately the same location but at different times, to produce an interferogram. We use SAR data from two satellites, ENVISAT (Image and Wide Swath modes) and ALOS, from between 2007 and 2009 to produce 75 interferograms (Table 1). ENVISAT Wide Swath interferograms were processed using the GAMMA software package (Werner et al., 2000). We used 19 scenes from ascending track 386 between October 2006 and August 2008. Interferogram selection was based on image pairs with perpendicular baselines less than 150 m and temporal baselines less than 200 days (Fig. A1). Interferograms with insufficient coherence and those unconnected to the network were then removed, leaving 20 interferograms (Table A1). We removed topographic phase contributions using the SRTM 30 m DEM (Farr and Kobrick, 2000), and filtered the interferograms using a Goldstein and Werner non-linear spectral filter (Goldstein and Werner, 1998), once with strength 0.8, and again with strength 0.6. Unwrapping was then done using the SNAPHU Minimum Cost Flow (MCF) algorithm with pixels with coherence less than 0.6 masked out (Chen and Zebker, 2001).
We processed data from ALOS, an L-band (23.6 cm wavelength) JAXA (Japanese Aerospace Exploration Agency) SAR satellite, using ISCE (InSAR Scientific Computing Environment) (Rosen et al., 2012). We used 5 acquisitions between June 2007 and December 2008 from ascending ScanSAR track 605 to produce 10 interferograms. Interferogram selection was made based on perpendicular baselines less than 500 m and temporal baselines less than 730 days (Fig. A1). Topographic phase contributions were removed using the SRTM 30 m DEM, and interferograms were resampled to 90 m. Each interferogram was then filtered twice (strength 0.4), before being unwrapped using the SNAPHU MCF algorithm, with unwrapping threshold of 0.1 (Chen and Zebker, 2001).
ENVISAT Image Mode data were processed using ISCE (Rosen et al., 2012). We used 10 Image Mode acquisitions from descending track 321 between October 2007 and March 2009 to produce 45 interferograms. Interferogram selection was based on image pairs with perpendicular baselines less than 800 m and temporal baselines less than 600 days (Fig. A1). Interferograms whose coherence does not extend across the caldera were excluded, to leave 28 (Table A1). Topographic phase contributions were removed using the SRTM 30 m DEM (Farr and Kobrick, 2000), after which pixels were multilooked to 120 m to increase coherence and reduce noise. We filter the interferograms using a Goldstein and Werner non-linear spectral filter with strength 0.8 (Goldstein and Werner, 1998). Interferograms were then unwrapped using the SNAPHU MCF algorithm with an unwrapping threshold of 0.1 (Chen and Zebker, 2001). Interferograms from all sensors were de-ramped where necessary.
We choose not to apply atmospheric corrections to our dataset as Corbetti has gentle, low topography (∼200 m), there are very few input data for weather models for East Africa, and the available atmospheric corrections are unable to match any turbulent signals in the data (Doin et al., 2009).

Surface deformation results
The ENVISAT IM data provides the clearest measure of the extent of the signal in the southern portion of Corbetti; it clearly shows a circular minor segment shape, with a sharp northern boundary ( Fig. 2g-l). The ENVISAT IM interferograms start in October 2007, and measure ∼1 cm of range increase between December 2007 and August 2008. The coherence was limited to the Chabi obsidian flows and an area west of the caldera that is not farmed.
The ALOS interferograms from June 2007 to December 2008 confirm the observation from the Envisat IM data, but span a longer time period, showing two distinct periods of deformation ( Fig. 2c-f). The first, between June 2007 to December 2007, is a period of range decrease of ∼2.5 cm in the southern portion of the caldera, with the same spatial pattern as in the ENVISAT WS data ( Fig. 2b-f). Interferograms post-December 2007 show range increase in the same region, totalling ∼1.5 cm by December 2008 ( Fig. 2e-l). The ALOS interferograms are the most coherent, illustrating the value of L-band InSAR for investigating ground deformation of arable land. However, ALOS interferograms with acquisitions on certain dates were affected by strong phase ramps, possibly due to orbital errors or, ionospheric or atmospheric delays, but de-ramping was able to account for this well on a local scale.
The ENVISAT WS interferograms represent the earliest measures of ground deformation in our analysis and show that prior to July 2007 there was no significant deformation at Corbetti (Fig. 2a). After July 2007, a north-south range change gradient of ∼2 cm over 2 km can be seen across the caldera, consistent with other observations (Fig. 2b). However, for the C-band WS data, coherence is limited to the Chabi obsidian flows in most interferograms.

Displacement direction
Using interferograms from ascending and descending tracks, which have different LOS vectors but measure the same signal, the vertical (U u ) and east-west (U e ) components of deformation can be estimated, but the sensitivity to deformation N-S is poor (e.g., Wright et al., 2004). We therefore formulate our equations making the assumption that this component of range change is zero. Mathematically, the range change observed by a satellite can be described by r =p.u, where p is the unit vector (p e , p n , p u ), pointing from the satellite to the ground in the local east, north and up directions, and u is the column vector of the components of displacement in the same reference frame, (u e , u n , u u ) T . This approach relies on the assumption that the ascending and descending images measure the same signal, which in the case of time-varying signals require them to have been acquired contemporaneously. To reflect this we selected interferograms with acquisitions close in time; the 26/12/2007-10/12/2008 descending ENVISAT and 23/12/2007-25/12/2008 ascending ALOS interferograms.
Decomposition of ascending and descending InSAR images into vertical and east-west components indicates the deformation in the south of the caldera between December 2007 and December 2008 is roughly vertical and ∼2 cm in magnitude (Fig. 3). The east-west component of this deformation shows some features of a similar magnitude, but their spatial extent is inconsistent with the signal seen elsewhere, and so are likely to be atmospheric artefacts. The deforming area in this ALOS interferogram includes a region of atmospheric delay, and contains some short wavelength features north of the caldera. Since the area that deforms in this period is the same as during June 2007 to December 2007, we assume that the uplift is also vertical.

Displacement time series
An individual interferogram records the deformation between two acquisitions, separated in time. To construct a time-series of deformation we used the SBAS (short baseline subset) method of Berardino et al. (2002). The displacement during incremental time steps is found from the vector of line-of-sight (LOS) displacements for each given pixel through a design matrix constructed to take into account the timespan of each interferogram. The cumulative displacement between any two dates can then be found by summing the relevant incremental displacements. This linear discrete inverse problem is under-constrained, and following Berardino et al. (2002) we use singular value decomposition, normalised using the L2 norm constraint. We applied a bootstrapping test and found our signal is not dependant on any particular satellite image, and therefore robust (e.g., Ebmeier et al., 2013).
Each interferogram provides relative phase changes, and for our time-series analysis we considered the relative displacement between a northern (fixed) and southern portion of the caldera, averaged over ∼125 km 2 and 38 km 2 respectively. The uncertainty in our time-series analysis is based on the mean variance in each dataset within 10 km 2 away from the caldera, where possible (∼0.79-0.84 cm). This is comparable to theoretical estimates of the variance of atmospheric noise over short length scales (∼10 km) from Emardson et al. (2003), which is 0.8 cm for any individual acquisition. Fig. 3 shows the result of time-series analysis from all three In-SAR datasets. Following our component analysis, we projected the LOS displacement into the vertical to better compare datasets with different LOS vectors. Individual WS interferograms indicate the absence of signal at Corbetti between October 2006 and July 2007 (Fig. 2a). The cumulative range change derived from the ALOS data The spatial extent of the deformation can be seen in Fig. 2b-l. Profiles though the region indicate that the northern edge of the signal is sharp; occurring over ∼2 km, and is co-incident with the caldera long axis (Fig. 2m-n). Fig. 2g-l shows that the signal is contained within the caldera to the west, but extends outside to the east. For some time-steps (e.g., ENVISAT IM displace- Fig. 3. a) Time-series of deformation from EMVISAT Image mode, ENVISAT Wide Swath mode and ALOS data. A map of Corbetti is included as an inset to show the regions between which this relative motion is calculated (southern relative to northern rectangle). b-c) ALOS and ENVISAT interferograms, between December 2007 and December 2008, used to determine the east-west (positive east) and vertical (positive up) components of deformation. d-e) east-west and vertical components of deformation. The dashed line shows the extent of the signal, identified in Fig. 2h. f) Profile along X-X through east-west and vertical deformation. Dashed line corresponds to location between triangle markers on X-X , indicating the deformation gradient. ment 27/08/2008-05/11/2008) the signal appears to extend into the northern portion of the caldera, but we attribute this to atmospheric artefacts.

Source modelling
To estimate the depth of a source that could cause a signal with such a sharp northern boundary, we forward model the deformation using a horizontal Okada rectangular dislocation model (Okada, 1985). The gradient of the deformation is a first-order indication of source depth, as shallow sources produce sharp edges, while deformation from deeper sources is smoothed by the elastic crust. We do not attempt to eliminate the possibility of an additional deeper source, but use the rectangular dislocation to produce a step-function at depth and hence estimate the maximum depth to the shallowest part of the deformation source. We plot profiles of deformation caused by a sill at depths of 5, 2, 1, 0.75, 0.5 and 0.25 km, normalised based on peak deformation. The modelled sill has a length (100 km), much greater than the kilometre length scale we are interested in, to ensure the shape of the deformation we observe is a result of a single sill edge only. We find that for peak-to-peak deformation to be contained within 2 km, a rectangular dislocation must be less than 0.75 km below the surface (Fig. 2o). For a 1 km depth, 90% of the deformation is contained within 2 km. A source with a tapered slip distribution would need to be shallower to produce the same deformation gradient. In this experiment we do not consider the interaction between the source and the fault, but we may expect source-fault interaction to amplify the deformation gradient, resulting in an underestimation of the source depth (e.g., Folch and Gottsmann, 2006).

Deformation mechanisms
Surface deformation occurs as a result of changes in volume or pressure in the subsurface. There are three possible sources that may have caused the deformation at Corbetti: magma, hydrothermal fluids, or meteoric water.
The sharpness of the northern boundary of the deformation implies that the source is shallow (∼1 km) (see Section 6) (Finnegan et al., 2008). The presence of magma at such shallow depths would result in surface manifestations, such as changes in fumarolic behaviour, possible phreatic or phreatomagmatic eruptions from interaction with meteoric water, changes in groundwater chemistry or felt seismicity, which have not been reported (e.g., Wicks, 2002;Jay et al., 2013).
Perturbations in pore fluid pressure associated with hydrothermal circulation can also result in surface deformation (e.g., Finnegan et al., 2008;Chaussard et al., 2014), thus it is important to consider coupled hydrothermal and magmatic systems when interpreting deformation at volcanoes with well developed hydrothermal systems. Pore fluid pressure changes can be driven by increased heat or fluid entering the system, or changes to fracture networks in response to subsurface stress changes (e.g., Bonafede, 1991;Rowland and Sibson, 2004). The response to these pressure changes can be either elastic, caused by seasonal rainfall variability, or inelastic, such as non-recoverable aquifer compaction (e.g., Lanari et al., 2004).
The MER and adjacent highlands have highly variable rainfall, which could cause seasonal deformation of shallow aquifers (Birhanu and Bendick, 2015). However, the ENVISAT Wide Swath data shows no deformation between October 2006 and March 2007, indicating that the signal is not part of ongoing seasonal variations, and so unlikely to be hydrological in nature as this time period covers the April-May rainy season.
We therefore propose that the deformation is related to perturbations in the hydrothermal system, caused by an increased flux of water or heat in response to deeper magmatic processes. An increased flux would result in a pressure and/or temperature increase in the hydrothermal reservoir, causing a volume increase and therefore uplift (Fig. 4a) (e.g., Miller et al., 2017). The overpressure then diffused away, possibly by the breaching of a barrier that previously confined the water, such as flow through newly formed cracks, causing subsidence (Fig. 4b) (Ali et al., 2015). The process has been observed at Campi Flegrei in the 1980s, and there are similarities between the temporal evolution of the deformation there and our observations at Corbetti (e.g., Bonafede, 1991;Battaglia et al., 2006). Surface deformation caused by coupled magmatic-hydrothermal systems has also been observed elsewhere in the East African Rift: at Aluto  The deformation is contained within the Corbetti caldera to the west and south, but extends outside the caldera to the east. At calderas which collapse via piecemeal or piston style mechanisms, sharp offsets occur around the caldera rim, and the caldera floor can become highly fractured (Lipman, 1997;Walter and Troll, 2001;Holohan et al., 2005). The variation in caldera scarp height at Corbetti suggests possible asymmetric collapse, and the continuity of the deformation outside the caldera could indicate an absence of a bounding caldera ring fault to the east. These observations are consistent with the collapse of Corbetti caldera via a trapdoor mechanism, where most of the collapse is in the west and the eastern rim acts as the hinge (Girard and Vries, 2005;Acocella, 2007). Alternatively, if slip occurred on the entire circumferential caldera fault during collapse, any structures that crossed this fault would be offset.

Seismicity
The seismicity in the Corbetti region has been studied using an array of broadband seismometers. These seismometers were operational for 2 yr (2012-2014) at 7 stations, with a maximum of 5 stations working at any given time (Fig. 1b), described in more detail in Wilks (2016). P-and S-wave first-breaks were manually picked, where coherent at three or more stations, and attributed weightings based on their quality. The software package NONLIN-LOC (Lomax et al., 2000) was used to locate earthquakes, using Pand S-wave arrival times and a one dimensional velocity model from Daly et al. (2008). Where possible, additional constraints on seismicity in the region came from stations deployed at the nearby Aluto volcano (Wilks et al., 2017).

Earthquake locations
Over the 2-yr deployment period 780 earthquakes were located, 224 of which were within 15 km of Corbetti caldera centre. In contrast, the Aluto volcano, which lies roughly 75 km north of Corbetti, experienced over 2000 similar sized earthquakes in the same time period (Wilks et al., 2017). At Corbetti, most of these events were located between Urji and Chabi down to a depth of 9 km. These earthquakes are associated with volcanic deformation that occurs after our InSAR observations, and do not give any further information on the subsurface structure.
However, during the period January 2012 to January 2014, 71 earthquakes were recorded associated with the ∼E-W trending Wendo Genet scarp and Werensa Ridge, ∼650 and ∼350 m high respectively and located ∼10 km to the east of Corbetti (Fig. 1).
The presence of slickensides on these features suggests there has been strike-slip motion on these faults, although no offsets have been reported (Mohr, 1968;Korme et al., 2004;Rapprich, 2013). Left-lateral strike-slip displacement here is consistent with models of MER kinematics, derived from structural data, focal mechanisms and GPS velocities (e.g., Muluneh et al., 2014). Most of the seismicity occurred between 7 and 15 km, with the shallow subsurface comparatively aseismic (9 events <7 km).
Between 15-20 km no events occur, but 20-32 km, there are 25 earthquakes (with maximum lateral uncertainty <10 km). These deepest events depict a linear structure that dips towards the southwest and are unlikely to be associated with transient volcanic processes at Corbetti given their distance from the caldera. The depth and distribution of the seismicity, on a linear plane down to 30 km, indicates that this structure extends down throughout the crust. The occurrence of seismicity here, along-strike of the deformation we identify within the caldera, suggests that the structure which cross-cuts Corbetti continues outside the caldera, and cuts across the border faults.

Seismic anisotropy
We use shear-wave splitting analysis to evaluate seismic anisotropy in the upper crust in the Corbetti region. The vertical alignment of sub-seismic length-scale cracks and fractures in the crust leads to variations in seismic velocities with direction and polarisation. The propagation of two independent shear waves with orthogonal polarisations is perhaps the most unambiguous indicator of seismic anisotropy. The alignment of fracture reveals the orientation and the anisotropy of the stress field (e.g., Verdon et al., 2008). Fractures will align in the direction of maximum horizontal compressive stress, which is revealed by the polarisation of the fast shear-wave (φ). The delay time (δt) between the fast and slow shear-waves is proportional to the fracture density or difference in maximum and minimum horizontal stress. It is also sensitive to the compliance of the fractures, which is related to properties such as permeability and fluid content. For the purposes of this work, we neglect the influence of any intrinsic horizontal crystal alignment in the rocks of the shallow crust or any fine scale horizontal layering, as vertically propagating shear waves will be less sensitive to such anisotropy (see Verdon et al., 2009, for more discussion of this).
Shear-wave splitting analysis is performed on ∼1200 sourcereceiver paths to the Corbetti seismic stations. For details of the methodology, the reader is referred to Wuestefeld et al. (2010). 28 measurements produce acceptable splitting results. We neglect measurements with errors >20 • in φ and 0.02 s in δt, and any source-receiver paths inclined >4 • . The range of δt is up to 0.31 s, which corresponds to a shear wave anisotropy of up to 9.2%. Shear wave splitting is accrued along the waves' travel path, so in subsequent figures we plot it at the event-station midpoint. Fig. 5a shows polar histograms of orientations of the fast shear waves at each station, as well as the overall trend of φ. It is apparent that most measurements show a fast shear-wave polarisation that is parallel to the Werensa ridge, and is consistent with the orientations expected for the principle stresses on an east-west striking strike-slip fault. The orientation of the fast shear-wave polarisation is not the same as current plate motion or extension direction implied by the Wonji Fault belt. However, there is a secondary cluster of rift parallel orientations, which occur in paths outside the caldera to CO3E and CO7E, the stations which also lie furthest from the centre. The magnitude of the splitting is highest in the southern half of the caldera. These observations suggest that the most intense fracturing is in the deforming region, as might be expected, and that the regional stress field is most strongly perturbed within the caldera. Outside the edifice, the regional stress field appears to dominate.
Observations of splitting from deep events (down to ∼35 km) below the Wendo Genet scarp and Werensa ridge further support the hypothesis that the cross-rift structure plays a role in modifying the stress field. Fig. 5b shows the dichotomy in φ between splitting in ray-paths travelling north of this zone to station C03E, and those within in it (to C02E). Splitting is cross-rift-parallel for the southern paths and rift-parallel for the northern ones, which also show more anisotropy, though both are relatively weak compared to paths within the caldera.

Caldera geometry and locations of post-caldera volcanism
We hypothesise the subsurface structure that limits the extent of the deformation also influenced the magma plumbing system at Corbetti, specifically the formation of the pre-caldera magma reservoir and the location of post-caldera volcanism. We test this hypothesis by considering the geometry of the surface features (e.g., Acocella et al., 2002;Le Corvec et al., 2013). The shape of a caldera at the surface is thought to reflect the shape of the precaldera magma reservoir and we therefore test whether the caldera rim is elliptical in shape, and whether the long axis of this ellipse is consistent with the deformation boundary. Pre-existing structures can influence the geometry of magma reservoirs, causing them to be elongate in the orientation of the structure (Holohan et al., 2005;Robertson et al., 2016). Secondly, since faults can act as a pathway for magma migration, we considered whether there is a preferential alignment of post-caldera vents (Mazzarini et al., 2013;Muirhead et al., 2015).

Caldera rim geometry
We digitised 80 points that describe the location of the exposed caldera rim, identified using high resolution optical imagery (Google Earth, 2016) and published maps, using QGIS (Rapprich, 2013;Gíslason et al., 2015). We then inverted for the long and  Wilks et al. (2017). Polar histograms (dark blue) are shown at each station which recorded at least three observations. The inset (light blue) histogram shows fast orientations for the whole data set. Blue, black, red and purple lines are the same as Fig. 6, which show the orientation perpendicular to border faults, of the caldera long axis, plate motion and deformation boundary, respectively. Circles show earthquake locations giving shear wave splitting observations, with colour indicating the strength of shear wave anisotropy. Lines connect events (black dots) with stations (white triangles). Note that anisotropy is largest for paths closest to the centre of the caldera, and that the dominant trend is sub-parallel to the trend of the Wendo Genet scarp and Werensa ridge. b) Variation of shear wave splitting with path from events near the Wendo Genet scarp and Werensa ridge. Coloured bars show the orientation of the fast shear wave measured at stations C02E (cyan) and C03E (blue), respectively, with the length of the bar proportional to the shear wave anisotropy, with a minimum bar length of 2% anisotropy. Coloured circles also show anisotropy as per the scale shown (note the scale is different between subplots). Circles and bars are plotted at the event-station midpoint. Black dots show earthquake locations and lines connect events and stations. Note that paths which spend longer to the north of the scarp show a rift-parallel φ trend, whilst those to the south have φ closer to the scarp strike. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) short axis lengths, caldera orientation, and centre point from these points using the method of Szpak et al. (2014). This method uses an approximate maximum likelihood approach which combines the accuracy of orthogonal methods and the speed of algebraic methods to find the solution to the equation of a conic that is non-degenerate.
The method seeks to minimise the Sampson distance, an algebraic approximation of the orthogonal distance between points and a candidate ellipse (Szpak et al., 2014). The use of the algebraic Sampson distance allows the mathematical equation of the conic to be expressed in terms of geometric parameters of an ellipse (orientation, length of long and short axes and centre point location), and give a quantitative measure of the uncertainties in the form of a covariance matrix. This method makes the assumption that the noise associated with the location of each point on the caldera is independent and Gaussian.
The ellipse that best describes the surface expression of the Corbetti caldera rim is centred at 38.381 • E 7.192 • N, and has a 13.8 ± 0.4 km long axis orientated 097 ± 3 • , and a 11.3 ± 0.2 km short axis (Fig. 6). The ratio of long axis to short axis length de- Fig. 6. Geophysical and structural data indicating there is a pre-existing fault trending ∼097 • through Corbetti. a) The black line shows the caldera long axis, black cross the caldera centre and black line the exposed caldera rim. Red circles identify vents or craters decametre in scale or larger, as identified by satellite optical imagery. Black triangles demarcate the apexes of Urji (west) and Chabi (east). The purple dashed line shows the location of the pre-existing fault as observed using InSAR (c). b) Histogram shows the mean vent-structure distances for 10,000 simulations of 16 vent locations inside the caldera for a E-W (red), N0-S (blue), NW-SE (green) and NE-SW (yellow). Dashed vertical lines indicate the mean vent-structure distance for each orientation of the corresponding colour. c) InSAR data from 23/12/2007-25/12/2008 showing the deforming region. d-e) Magnetotelluric data after Gíslason et al. (2012). Inset are the orientations of the caldera long axis, InSAR gradient, current plate motion direction (Stamps et al., 2008), and the perpendicular to the border faults. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) fines the caldera ellipticity, which at Corbetti is 0.82 ± 0.03: we therefore consider the caldera to be elliptical. The caldera long axis is co-incident with the boundary of the deformation region (Fig. 6). Since the shape of the caldera rim is taken to reflect the geometry of the pre-collapse magma storage region, it implies that the structure influenced magma migration prior to the caldera collapse at 182 ± 18 ka.

Post-caldera volcanism
Crustal structures can also act as pathways for migrating magma, influencing the location of small vents (e.g., Korme et al., 1997Korme et al., , 2004Le Corvec et al., 2013). To test whether the structure influenced magma migration since the caldera collapse at Corbetti, we quantify the locations of post-caldera volcanism in relation to the caldera geometry. We digitised the location of 16 post-caldera vents greater than 10 m in diameter, identified using high resolution optical imagery (Google Earth, 2016) and published maps (Rapprich, 2013;Gíslason et al., 2015;Rapprich et al., 2016), using QGIS (Fig. 6).
Our hypothesis is that post-caldera volcanism is influenced by a subsurface structure, taken to be coincident with the caldera long axis. This predicts that vent locations will be closer to the caldera long axis than a random distribution of vents within the caldera and can be statistically tested by comparing the mean distance between the mapped vents and the caldera long axis, and the same measurement for a synthetic, random dataset of vent locations. We only used vents within the caldera, which we defined as the area inside the caldera rim, where exposed, and within the best fitting ellipse that describes the caldera where there is no clear rim. To find the probability that the vent locations inside Corbetti are closer to the caldera long axis than if generated at random, we simulated 10,000 other 16 vent locations and found the mean distance for each simulation. The proportion of simulated mean distances less than the observed mean distance gives the probability that randomly formed vents would be located closer to the caldera long axis than the observed vents. Fig. 6 shows the distribution of random mean vent-long-axis distances for 10,000 simulations. The proportion of simulations with mean distances less than the measured mean distance, 1390 m, defines the probability that the actual distribution of vents are located at random. For a structure aligned with the caldera long-axis, this proportion is 3.0%, which demonstrates statistical significance. Furthermore the major centres of resurgent volcanism, Urji and Chabi, both lie along the ellipse long axis.
However, clustering of vents in the centre of the caldera would also produce a small mean vent-structure distance, but the vents will lie equally close to a line of any orientation through the caldera. As such, we find the probability that the vents are distributed closer to a 'structure' orientated N-S, NE-SW and NW-SE than the observations. These probabilities are 0.46, 0.24 and 0.25 respectively, much higher than for the hypothesised E-W structure (0.03). However, even in a homogeneous medium vent locations are unlikely to be random because established magma pathways are often reused and surface topography exerts stresses that may influence magma migration (e.g., Pinel and Jaupart, 2003;Roman and Jaupart, 2014;Xu and Jónsson, 2014).
Nonetheless, we conclude the post-caldera volcanism is located closer to a subsurface planar structure than randomly distributed vents would be. This suggests that such a structure influenced magma migration over the timescales of post-caldera vent formation, which is on the order of tens of thousands of years (Rapprich, 2013;Hutchison, 2015).

Nature of subsurface structure
In this section we discuss the candidates for the subsurface structure, of which several are plausible in rift settings: for instance the stock of an earlier volcano, the rim of a preceding caldera or a pre-caldera fault. The candidate must be able to explain the orientation and ellipticity of the caldera, the distribution of post-caldera vents, and the horizontal and vertical extents inferred from seismic and InSAR observations.
Beneath the caldera complexes in the MER solidified magmatic intrusions have been identified with gravity and seismic surveys (Cornwell et al., 2006;Maguire et al., 2006). An intrusion beneath the northern half of the Corbetti caldera, with its edge along the caldera long axis, would enhance deformation in the southern portion of the caldera relative to the north, as observed. This is because crystallised silicic material has different material properties (e.g., rigidity and permeability) to partially molten intrusions (e.g. Hickey et al., 2013). Material properties that are spatially variable in this way would also explain the shear wave splitting measurements that suggest differences in fracture density between the north and the south. While we cannot discount this explanation, it is unable to explain the caldera orientation or the location of post-caldera vents, which occur in both the northern and southern parts of the caldera.
An alternative explanation is that the cross-cutting structure is related to the ring fault of the Hawassa Pleistocene collapse caldera (Woldegabriel et al., 1990;Rapprich, 2013). There is a change in strike from east-west within the caldera, to north-west-south-east at the Werensa ridge and Wendo Genet scarp, which would be consistent with a caldera ring fault. However, caldera faults typically extend to depths of less than 10 km (see Saunders, 2001;Cole et al., 2005;Holohan et al., 2005), and magnetotelluric data within the caldera shows the structure extends to at least 10 km (Gíslason et al., 2015) (Fig. 6), while outside the caldera, seismicity on the Wendo Genet scarp and Werensa Ridge indicate the structure extends to 30 km depth (Wilks, 2016) (Fig. 1).
Our preferred candidate for the structure that cross-cuts Corbetti is a pre-existing fault. Pre-rift faults cross-cut the Precambrian basement throughout Ethiopia, and many have been reactivated by the current phase of rifting (Korme et al., 2004;Chorowicz, 2005;Corti, 2009). For example, the pre-Jurassic Ambo Fault Zone (AFZ) and the Yerer-Gugu lineament controlled the Pleistocene development of the transfer zone between the northern and central MER around 8 • N (Bonini et al., 2005) (Fig. 1). The AFZ can be identified perpendicular to the rift as an off-rift low velocity structure (Bastow et al., 2005). Another pre-rift structure, the Goba-Bonga lineament, is thought to have impeded the northwards propagation of the Kenyan Rift at ∼11 Ma and stalled MER rifting during the Miocene (Bonini et al., 2005). Oligio-miocene (20-30 Ma) volcanism is aligned along the AFZ and the Goba-Bonga lineament on the western rift flank (Korme et al., 2004;Corti, 2009), suggesting such structures remain important magma pathways for tens of millions of years.
Corbetti caldera lies on the cross-rift projection of the Goba-Bonga lineament. The Werensa ridge and Wendo Genet scarp are thought to be the surface expression of the Goba-Bonga lineament, where it obliquely intersects the rift border fault (Žáček et al., 2014). The structures do have a different strike to the structure identified within the Corbetti caldera, but a curvature of the shallow resistive anomaly (Fig. 6) can be seen connecting the two: from east-west inside the caldera to north-west-south-east to the east. Furthermore, it is not uncommon for fault systems to change strike along their lengths (e.g., Sengör et al., 2005). At the Werensa ridge and Wendo Genet scarp there is seismicity (71 earthquakes, magnitudes between 0.65 and 4.10, January 2012-January 2014) down to ∼30 km ( Fig. 1) (Wilks, 2016). The Moho here is ∼38 km deep (Keranen and Klemperer, 2008) and so the depth extent of the seismicity is evidence that this fault cross-cuts almost the entire crust and has the potential to influence magma storage and transportation on a crustal scale. The occurrence of earthquakes at lower crustal depths is not atypical in the MER, specifically seismicity has been observed ∼35 km deep at similar crustal-scale pre-rift structures such as the Yerer-Tullu Wellel (Fig. 1) (Keir et al., 2009).
This interpretation of the structure provides an explanation for all of the observations, including the pre-caldera elliptical magma body, post-caldera magmatism, deformation, resistivity and seismicity. This major crustal structure acts to guide the vertical transport of hydrothermal fluids into the shallow subsurface, and as a barrier to horizontal fluid flow, defining the lateral extent of the hydrothermal reservoir. The presence and location of the structure is also consistent with the hypothesis that some MER silicic centres formed where transtensional faults cross-cut the rift, creating regions of localised extension or weaker material that promotes magma reservoir formation (Acocella et al., 2002).

Conclusions and comparisons
We show that pre-rift structures influence magmatic and hydrothermal processes over a range of timescales at the Corbetti caldera, Ethiopia. A rift-oblique structure influenced (1) surface deformation associated with a fault bounded reservoir, which shows the influence of faults on hydrothermal circulation over annual timescales; (2) the location of post-caldera volcanism, which highlights the influence on shallow magma transportation pathways over tens of thousands of years; and (3) the caldera geometry, indicative of a control on crustal magma storage over hundreds of thousands of years.
This work raises questions regarding the influence of pre-rift oblique structures in continental rifting, both in terms of their influence on magmatic processes but also their role in strain accommodation. A comparison can be drawn to the Taupo Rift System, New Zealand, where pre-existing oblique structures align with volcanic domes (e.g., the Tarawera Dome Complex, Cole et al., 2010) and control hydrothermal circulation (Rowland and Sibson, 2004). In contrast, in Icelandic rift zones, there is no basement continental crust and to date there appears to be little evidence of oblique structures influencing volcanism, suggesting that a primary control on magma pathways depends on crustal structure.
This study demonstrates the importance of combining different techniques and datasets that give observations of multiple processes over multiple timescales to understand magmatichydrothermal systems. Our observations have implications in understanding the relative importance of heterogeneities that affect the development of magmatic systems, especially in active tectonic regimes such as continental rifts. It provides insight into how these same processes influence hydrothermal reservoir formation, and also how they respond to external influences. LICS and a NERC New Investigators grant (NEI00816/1). MW was funded by an EPSRC studentship. AN was funded by a Leverhulme Early Career Fellowship (ECF-2014-496) and Bristol University Microseismicity Projects. We thank the Aluto Research and Geophysical Observations (ARGOS) project for seismic data, and ESA for ENVISAT data. ALOS data were provided through ESA third party mission. Seismic data (2012-2014) is available from the Incorporated Research Institutions for Seismology (http :/ /iris .edu) under the temporary network code XM. This work is a contribution to the Natural Environment Research Council (NERC) funded RiftVolc project (NE/L013932/1, Rift volcanism: past, present, and future).