The spatial clustering of fracture networks and vents in basaltic volcanic fields has been analyzed in three sectors of the East African Rift System, the classical example of an active continental rift. Fracture trace maps and monogenetic basaltic vents have been thus collected in the Afar Depression, in the Main Ethiopian Rift, and in the Virunga Belt (Western Rift). The mapped vents are generally younger than 2 Ma, and most are of Holocene age.

All the analyzed fracture networks have self-similar clustering with fractal exponents (Df) varying in the 1.54−1.85 range. Also, vents show a self-similar clustering with fractal exponents (Dv) in the 1.17−1.50 range. For all the studied sectors, the relationship Df > Dv has been observed. The fractal exponents for vents (Dv) of power-law distributions are computed in a range of lengths with a lower and an upper cutoff. The upper cutoff (Uco) for the fractal clustering of vents in the studied sectors of the East African Rift System are compared with the respective crust thickness derived by independent geophysical data. The computed Ucos for the studied sectors well match the crust thickness in the volcanic fields. A preliminary conceptual model to explain the relationships between the upper cutoffs of the fractal distribution of vents and the thickness of the crust in the volcanic fields is thus proposed in the light of the percolation theory.

Fluids commonly migrate through the Earth's crust in hydro-fractures such as mineral veins or dykes, stopping at various depths in the crust (e.g., Watanabe et al., 1999; Dahm, 2000a, 2000b; Gudmundsson and Brenner, 2004). It is widely accepted that the transport of magma from lower crustal-upper mantle source/reservoir regions to its final level of emplacement mainly occurs through dykes and sills (Lister and Kerr, 1991; Petford et al., 1993; Rubin, 1995; Petford et al., 2000). These fluid-driven fractures or hydro-fractures are generally opening-mode fractures (e.g., Gudmundsson, 2002; Gudmundsson and Brenner, 2004). The magma ascent rate depends on magma properties (including viscosity, density, temperature, and heat content), country rock properties (including temperature, density, thermal conductivity, and permeability), and the stress field. Several lines of evidence, e.g., high-density xenolith settling in basalt magmas (Basu, 1977; Spera, 1980; Petford et al., 2000), numerical analysis (Dahm, 2000a, 2000b), and magma cooling rates (Maaloe, 1973), indicate velocities of magma ascent of 10−2 ms−1 up to 1 ms−1, which imply high bulk permeability of the crust. Rock-fracturing processes enhance the bulk permeability of the crust and allow the ascent of magma at rates that are akin to the time-scale characterizing magmatic activity (Rubin, 1993; Petford et al., 1993; Petford et al., 2000; Canon-Tapia and Walker, 2004). Fluid-driven fracturing (fractures filled by magma, i.e., dyke formation) is thus the viable mechanism for emplacing magma within the crust (e.g., Turcotte, 1982; Hutton, 1996; Petford et al., 2000). Depending on magma buoyancy relative magnitude, crust's fracture toughness, crustal mechanical discontinuities, and magma availability, dykes may stop at some levels in the crust or, eventually, they may also construct magmatic chambers or reach directly the surface generating eruptions (e.g., Ida, 1999; Dahm, 2000a; Gudmundsson, 2002; Taisne and Tait, 2009).

In particular, eruptions of basaltic magmas imply the transfer of magmas from deep reservoirs up to intermediate magma chambers in middle upper crust or directly to the Earth's surface. Polygenetic volcanoes repeatedly erupt from the same general vent (summit or crater) and are formed if basaltic magmas stop at intermediate to shallow crustal magma chambers, whereas monogenetic volcanoes erupt only once (e.g., MacDonald, 1972) and are constructed when magma directly erupts from feeders. According to Connor and Conway (2000) volcanic fields are dominantly basaltic in composition and are formed by monogenetic vents each produced by a single episode of volcanic activity and associated to feeder dykes. The occurrence and spatial distribution of monogenetic eruptive structures within volcanic areas are linked to fracture systems and associated stress fields (Takada, 1994a). Moreover, morphometric parameters of monogenetic cones, such as cone elongation, breaching direction, and cone alignment, indicate the direction of fractures acting as magma feeders (Tibaldi, 1995).

It has been proposed that fractures filled by magma (i.e., dikes) tend to coalesce during their ascent to the surface, thereby controlling the final level of magma emplacement. The actual distribution of volcanic vents at the surface, i.e., the formation of monogenetic and/or polygenetic volcanoes, is mainly controlled by the magma input rate and the crustal strain rate (e.g., Fedotov, 1981; Takada, 1994a, 1994b). High-strain rate or small magma input rate promote the formation of monogenetic volcanoes whereas low-strain rate and high magma input rate mainly generate polygenetic volcanoes.

It is important to emphasize that basaltic monogenetic vents testify to the presence of deep crustal or subcrustal magma reservoirs directly connected via fracture network to the surface, involving a hydraulic connection through the whole crust or a large portion of it between source and surface. Moreover, the correlation between vent distribution and fracture network properties is such that the spatial distribution of vents may be studied in terms of self-similar (fractal) clustering (Mazzarini and D'Orazio, 2003; Mazzarini, 2004, 2007; Mazzarini et al., 2008), as in the case of fracture networks (Bonnet et al., 2001). Findings based on this approach suggest that, for basaltic volcanic fields in a deformed continental crust, the distribution of monogenetic vents is linked to the mechanical layering of the crust. Vents tend to cluster according to a power-law distribution defined over a range of lengths bounded between a lower limit (Lower cutoff, Lco) and an upper limit (Upper cutoff, Uco); the upper cutoff approximates the thickness of the fractured medium (crust). This correlation has been studied in volcanic fields within extensional continental settings in back-arcs, such as in southernmost Patagonia (Mazzarini and D'Orazio, 2003) and in continental rifts such as the Main Ethiopian Rift (Mazzarini, 2004) and Afar (Mazzarini, 2007), as well as in contractional continental settings as the Andean foreland (Mazzarini et al., 2008).

In this work we present a qualitative model to explain the observed relationships between the upper limit of the fractal distributions of basaltic monogenetic vents and the crust's thickness on the base of the percolation theory (Stauffer and Aharony, 1992). In this model the monogenetic vents represent the intersection of the connected part of the actual fracture network, which effectively ties the source of magma to the surface (i.e., the backbone in the percolation theory, see below), with the surface and their self-similar clustering linked to the fracture self-similar clustering.

Successively, the model will be tested in three different sectors of the East African Rift System (EARS) in a context of continental extension and new data on vent and fracture distributions will be discussed in the light of the proposed model and critically compared with available independent geophysical data.

The spatial distribution of volcanism in terms of volcano spacing and its relationship with fracture patterns and lithospheric structure has been studied since the early 1970s in oceanic and continental settings (e.g., Vogt, 1974; ten Brink, 1991). In particular, the spacing of central volcanoes within continental rift settings has been linked to the elastic thickness of the lithosphere (Mohr and Wood, 1976). Several authors have investigated the spatial and temporal distribution of volcanoes and monogenetic cones, focusing on the close relationship between tectonics and volcanism, for example in the Michoacan-Guanajuato Volcanic Field (Hasenaka and Carmichael, 1985; Connor, 1987, 1990), in the Camargo Volcanic Field in Mexico (Aranda-Gomez et al., 2003), and in the Springerville Volcanic Field in Arizona (Connor et al., 1992; Condit and Connor, 1996). Vent alignment has often been used to infer the direction of the minimum horizontal principal stress (Nakamura, 1977; Lutz, 1986; Wadge and Cross, 1988), and vent distribution has been used as evidence for structural control on vent location (Connor, 1990; Connor et al., 1992) and to outline the importance of strain rate in the style of volcanism (Takada, 1994a; Alaniz-Alvarez et al., 1998; Mazzarini et al., 2004).

Volcanoes tend to form clusters, showing statistical self-similarity in both space and time as documented in hot spots on oceanic crust (Shaw and Chouet, 1991) and in continental magmatic arcs (Pelletier, 1999). Vent clustering in basaltic volcanic fields has been also described in the Newer Volcanic Province, southeast Australia (Lesti et al., 2008) and in Afar, east Africa (Mazzarini, 2007).

As previously discussed, the occurrence of basaltic monogenetic volcanism implies the existence of a portion of the fracture network allowing a direct connection between the magma reservoir at depth and the surface (vents). Clearly this feature points to the geometric and hydraulic characteristics of the actual fracture network. On this basis, Mazzarini (2004) developed a simple model for visualizing the relationship between vents and the geometric properties of fractures. This model assumes that the aperture of a fracture is greatest at its barycenter and that volcanic vents erupt at the point of maximum fracture aperture; the resulting vent distribution is thus closely linked to the hydraulic properties of the actual fracture network.

The main geometric features of a fracture network (fracture attitude, aperture, spacing, intersections, length, and density) are generally measured and mapped at different scales of observation, showing scale invariance spanning several orders of magnitude (e.g., Marrett et al., 1999; Bonnet et al., 2001; Bour et al., 2002). The way in which fractures fill space (i.e., the spatial distribution of fractures) strictly depends on the spacing of fractures, that in turn is correlated with the thickness of the fractured medium calculated on the basis of the stress saturation model (Wu and Pollard, 1995; Gross et al., 1995; Ackermann and Schlische, 1997).

Hydraulic features of fractures such as fracture connectivity and aperture are scale invariant (Bonnet et al., 2001). In particular, the spatial clustering of a fracture network, represented as the fracture barycenter, has been directly linked to the hydraulic properties of the fracture network (Renshaw, 1999; Bour and Davy, 1999; Darcel et al., 2003).

The connectivity of fractures defines the portion of the existing fracture network that hydraulically connects the system boundaries, allowing fluids to flow (Margolin et al., 1998; Darcel et al., 2003). Connectivity mainly depends on fracture size (length), density, orientation, and on the spatial correlation among fractures (e.g., Renshaw, 1999; Berkowitz et al., 2000; Darcel et al., 2003).

Fracture lengths in nature often display a power-law distribution in the form
graphic
where N(l) is the number of fractures whose length is longer than l and a is the fractal exponent (e.g., Bonnet et al., 2001, and references therein).
A robust way to define how fractures fill space is to analyze their self-similar clustering (Bonnet et al., 2001). This is accomplished by computation of the correlation sum that accounts for all the points at a distance of less than a given length l (Bonnet et al., 2001, and references therein) following the equation
graphic
where C(l) is the correlation sum and D is the fractal exponent.

These two parameters strongly control the connectivity of the fracture network, high values of the a exponent in Equation (1) imply high frequency of short fractures, whereas low values of the a exponent mean that the system is mainly controlled by long fractures (Fig. 1A). Moreover, the higher the values of the D exponent in Equation (2), the more the fractures are homogenously distributed, on the other hand the lower the D exponent, the higher the fracture clustering (Fig. 1A).

Generally fracture networks exhibit fractal length distributions with exponent a < 3 (Renshaw, 1999) or a < 4 (Bonnet et al., 2001). The case where exponent a < 1 is very infrequent implying that system connectivity is controlled by very long fractures spanning through the whole system. For a > 3 the network connectivity is ruled by faults smaller than the system size and classical laws for percolation theory apply (Bour and Davy, 1997).

Fracture density in terms of fracture length distribution and fracture spatial distribution mainly control the overall permeability of the system. A critical fracture density for which the network is connected allowing a direct connection between the network's boundaries must exist. The percolation theory quantifies the existence of the critical fracture density for which the network is connected and defines some important parameters that describe the behavior of a network as it becomes a percolating network (Orbach, 1986; Stauffer and Aharony, 1992; Song et al., 2005). If the density of fractures in a network (ρ) is greater or equal to a critical density (ρc), the network becomes connected and a cluster of connected fractures spanning the whole sample will form (the infinite cluster; Stauffer and Aharony, 1992); ρc is defined as the percolation threshold, and for a two-dimensional (2D) system it is 0.59 (Stauffer and Aharony, 1992). The connected portion of the network is a subset of the existing fracture network (e.g., Roberts et al., 1998, 1999).

A simple example of a percolating network is a grid; the link between nodes of the grid is defined as a bond. The grid on the left side of Figure 1B has a bond density of ~44% and there is no connection between the lower and the upper side of the square. The grid on the right side of Figure 1B has a bond density of ~60% and it is at the percolation threshold, in this case a cluster of connected fractures (the infinite cluster) connects the lower and the upper side of the square (blue portion of the grid).

During the transfer of magma within the crust, as simplified in Figure 1B, the pressure gradient is essentially vertical and the connected network spans the whole path between source and the final location of the magma. No assumption is done about the isotropy of the medium. Anyway, in layered sequences lateral connectivity of fractures is a function of the fracture spacing that in turn is a function of the layer thickness (e.g., Wu and Pollard, 1995).

The portion of the infinite clusters that actually allows the connection, for example, the portion of the fracture network in which fluid actually flows, is named backbone. Above the percolation threshold (i.e., ρ ≥ ρc) the network and the backbone are statistically self-similar (Orbach, 1986; Stauffer and Aharony, 1992). A percolating network is self-similar in a well-defined range of lengths comprised between a lower and an upper cutoff, that is the percolating network is fractal for length scales r between α ≤r ≤ ξ; the lower cutoff α could be seen as the elementary grid, that is there are no bonds connecting pairs of nodes shorter than the elementary grid. The upper cutoff is often referred to as the percolation correlation length (ξ) (see Orbach, 1986). The percolation correlation length (ξ) is the pair connectedness length for percolation or, in other words, the maximum dimension of voids in the infinite clusters, in the case of r ≥ ξ the percolating network appears homogeneous (Orbach, 1986; Stauffer and Aharony, 1992). The definition of self-similar clustering for the analyzed spatial correlation of fractures (i.e., computation of the fractal exponent, see Bonnet et al., 2001, and references therein) is thus performed for a range of lengths (the size range) between lower and upper cutoffs. According to the definitions described above for a percolating fracture network, the backbone in the infinite clusters is the actual fluid pathway thus ξ represents the distance between the system's boundaries (i.e., the source and the surface for vents).

The upper cutoff is here considered to be directly linked to the mechanical layering of the medium. Mandelbrot (1982) suggested that there are upper and lower cutoffs for the scale-invariant characteristics of fractures (e.g., spacing, length, and density) and that these are a function of mechanical layers and rock properties. Experimental studies on normal fault populations suggest the presence of upper and lower cutoffs in the power law describing the distribution of the geometric properties of fractures (Ackermann et al., 2001). Moreover, the thickness of both sedimentary beds and the crust controls the scaling law of fractures and earthquakes (Pacheco et al., 1992; Davy, 1993; Ouillon et al., 1996).

The dependence of the fracture network spatial distribution on the rheologic layering of the medium (i.e., the crust) can thus be inferred from the backbone of the network (i.e., the vents). The connected fracture network allows basaltic magma to rise to the surface from deep crustal or subcrustal reservoirs, passing through most or the whole of the crust. Analysis of the fractal character of the spatial distribution of vents can thus reveal the mechanical layering of the crust. Assuming a direct genetic and spatial link between fracture and vent (e.g., Connor and Conway, 2000; Mazzarini, 2004), scale invariance in vent distribution reflects the fractal properties of the backbone. The backbone of the percolating network should be more clustered than the network (Stauffer and Aharony, 1992), i.e., the fractal D values estimated through Equation (2) for the backbone (Db) should be lower than the value for the network (Dn).

The proposed correlation between the distributions of vent and fractures and the crustal thickness is here explored in the East African Rift System (EARS). The EARS is a classical seismically and volcanically active continental rift extending several thousands of kilometers in a N-S direction accommodating the extension between the Nubian (Africa) and Somalian plates (Fig. 2; e.g., Rosendahl, 1987; Braile et al., 1995; Morley et al., 1999; Chorowicz, 2005). The boundary between the Nubia and Somalia plates is complex and the position of the pole of rotation is still a matter of debate (e.g., Horner-Johnson et al., 2007, and references therein). Anyway, a strong consistency exists between geologic and geodetic data, suggesting that the ESE-WNW Nubia–Somalia relative motion may have remained steady over the past 3 Myr (Corti, 2009, and references therein).

Bulging of the lithosphere, crustal extension, and consequent widespread volcanism have been ascribed to the impinging of one or two plumes on the base of the East African lithosphere in late Eocene–early Oligocene (e.g., Ebinger and Sleep, 1998; Rogers et al., 2000) or, more recently, to multiple plume branches rising from a deep-seated mantle upwelling (the African superplume, e.g., Furman et al., 2006).

Crustal stretching (final width/initial width; β), crustal thickness, and intensity of volcanism greatly vary along the EARS moving from south, western branch Virunga area, to the north, Afar Depression (Ebinger, 1989a, 1989b; Braile et al., 1995; Ebinger and Furman, 2002; Beyene and Abdelsalam, 2005; Dugda et al., 2005, 2007; Corti, 2009).

Three sectors of the EARS are here investigated in terms of their fracture and vent spatial clustering: the Afar Depression (AD), the Main Ethiopian Rift (MER), and the Virunga Belt (VB; in the EARS's western branch). These sectors are characterized by a different degree of crust stretching, crust thickness, and intensity of volcanism (Table 1 and Fig. 2).

Afar Depression

The Afar Depression is located at the confluence of the Main Ethiopian Rift, the western Gulf of Aden, and the southern Red Sea (Fig. 2). The Afar Depression is mainly floored by Pliocene and younger silicic volcanic rocks and basalts and by Quaternary sediments (Beyene and Abdelsalam, 2005; Mazzarini, 2007, and references therein). Interaction between the southern Red Sea and Aden oceanic ridges and the Afar stretched continental crust led to the formation of rifts and associated volcanism (Manighetti et al., 1998; Lahitte et al., 2003a, 2003b, and references therein). Basaltic volcanism is Pleistocene to Holocene and has produced several scoria cones and eruptive fissures (Lahitte et al., 2003a, 2003b, and references therein). The AD is characterized by strong crustal attenuation: Bouguer gravity data indicate crustal thinning (Makris and Ginzburg, 1987; Woldetinsae and Gotze, 2005) with an average thickness of ~25 km. Inverse modeling of gravity data shows a crustal thickness of 23–24 km in the AD (Tiberi et al., 2005). Seismic refraction data imaged a crustal thicknessof 23–25 km in the southern Afar Depression, and only 15 km in the northern part of the Afar Depression (Berckhemer et al., 1975; Prodehl and Mechie, 1991; Prodehl et al., 1997). Analysis of receiver functions from broadband seismic data (Dugda et al., 2005, 2007) reveals that the crust is 25 km thick in the Afar Depression. Both gravity and seismic data reveal a crustal thickness of 35–40 km in the shoulders of the rifts. Crustal stretching in AD varies from β = 1.7 in its southern part to β > 2 (Berckhemer et al., 1975; Eagles et al., 2002; Tiberi et al., 2005; Beyene and Abdelsalam, 2005).

Main Ethiopian Rift

The Main Ethiopian Rift (MER) connects the Afar depression, at the Red Sea–Gulf of Aden junction, with the Turkana depression and Kenya Rift to the south (Fig. 2). It is commonly divided into three sectors: north, central, and south (Hayward and Ebinger, 1996; Corti, 2009). The MER contains a large quantity of volcanic products characterized by a bimodal distribution of basic and acidic magma types (e.g., Trua et al., 1999). The MER is a late Miocene NE-SW– (in the north) to N-S– (in the south) trending fault bounded basin, filled by late Miocene to Holocene volcanic rocks and continental sedimentary deposits (e.g., Corti, 2009, and references therein). The rift valley is intersected by narrow NNE–trending structures (the Wonji Fault Belt; Mohr, 1987) consisting of fault bounded lava units aged less than 0.5–0.3 Ma (Morton et al., 1979; Chernet et al., 1998; Abebe et al., 2005). Exposed volcanic products consist of basalts, rhyolites, ignimbrites, and pyroclastic deposits. Monogenetic activity, consisting of spatter cones, scoria cones, maars, and lava domes (e.g., Mazzarini et al., 1999; Abebe et al., 2005; Rooney et al., 2007; Corti, 2009), is widespread on the rift floor as well as along the rift margins (e.g., Rooney et al., 2007). These vents are formed by evolved lavas (mainly rhyolites) in the case of domes, and by basalts in the case of cones. Vent ages range from ca. 5 Ma (mainly domes) to Holocene, with prevalent basaltic activity since late Pliocene–early Pleistocene (Morton et al., 1979; Chernet et al., 1998; Mazzarini et al., 1999; Abebe et al., 2005).

The crustal structure of MER has been imaged during gravity and seismic surveys. The crustal thickness estimated by analysis of seismic data and by gravity data inversion is 28 ± 2 km (Makris and Ginzburg, 1987). Gravity data inversion imaged a crustal section in the rift valley composed by 3–5 km of sedimentary infill and by upper crust ~20 km thick (Mahatsente et al., 1999). Seismic refraction data imaged a crustal thickness of 28–33 km in the northern Main Ethiopian Rift floor and ~40 km at the rift margins (Berckhemer et al., 1975; Prodehl and Mechie, 1991; Prodehl et al., 1997). Furthermore, analysis of receiver functions and joint inversion of Rayleigh wave group velocities and receiver functions from broadband seismic data reveal that the crust is 30–35 km thick in the Main Ethiopian Rift (Dugda et al., 2005, 2007).

A new seismic experiment has been conducted in MER in 2003 (EAGLE project; Maguire et al., 2003, 2006). EAGLE data show crust thickness along MER varying from 30 to 35 km in the north segment of the rift to ~35–40 km in the central segment (Keranen et al., 2004; Maguire et al., 2006; Keranen et al., 2009; Corti, 2009). At the junction of the north and central segments of MER (at the latitude of Nazret; Figs. 2 and 3D), the new seismic data highlights 30–35 km of crust thickness (Keranen et al., 2009). New seismic data by the EAGLE project (Keranen et al., 2004) also support that basaltic volcanic fields in the rift floor (e.g., Debre Zeyt; Mazzarini et al., 1999) are localized on right-stepping en echelon N-S–trending zones (e.g., Ebinger and Casey, 2001), where uppermost brittle crust is only 10–11 km thick lying above a lower crust that has been strongly modified by magmatic intrusions. Crustal stretching along MER varies from β = 1.1 in its south and central segments to β = 1.7 in its northern segment (Ebinger et al., 1993; Ebinger and Furman, 2002; Wolfenden et al., 2004; Tiberi et al., 2005).

Virunga Belt

The Virunga Belt (VB) is a volcanic province that lies near the northern end of the western branch of the East African Rift System: the Western rift (Fig. 2). In contrast to the widespread volcanism observed in the Eastern Rift, MER, and Afar, volcanism in the Western Rift is spatially and volumetrically small, leaving the vast majority of the Western Rift devoid of magmatism (Braile et al., 1995; Ebinger and Furman, 2002). The volcanic activity started since 11 Ma has produced basaltic lava flows ~500 m thick flooring the rift, successive activity focused in near E-W corridors, with a period of tholeiitic and alkalic volcanism beginning at 2 Ma and continuing to historic time. Erupted magmas are characterized by silica-undersaturated mafic products and the local development of evolved products from central volcanoes (Ebinger 1989a, 1989b; Ebinger and Furman, 2002). Analysis of broadband seismic data reveals that the crust is 30–35 km thick in the Westrn rift and the Virunga belt area (Dugda et al., 2005). Locally, crust thickness is supposed to be less than 30 km in basins located in transfer zones as the VB area (see Ebinger 1989a, 1989b; Corti et al., 2002). Analysis of earthquake catalog and focal depth distribution in the Virunga area (Albaric et al., 2009) shows that the number of crustal earthquakes is maximum at 10–15 km of depth and then decreases down to a depth of ~32 km. Recently, teleseismic P-wave receiver function analysis (Tuluka, 2009) imaged a crust mantle transition at a depth varying from 30 to 42 km, and a low velocity zone at depth varying from 18 to 30 km in the area of the large active Nyamuragira and Nyiragongo volcanoes. In the Virunga area we thus assume a crust thickness in the range 20–30 km (Table 1). Crustal stretching is low, β < 1.2 (Ebinger, 1989a, 1989b).

Data Collection

The analyzed data sets consisted of the locations of mongenetic vents and the fracture trace maps produced for the Afar Depression, Main Etiophian Rift, and Virunga Belt. In these areas, exposed volcanic products consist of lava flows, pyroclastic deposits, and ignimbrites associated with central volcanoes, spatter and scoria cones, maars, and lava domes. The mapped vents correspond to different types of monogentic volcanoes such as cinder or lava cone, domes, and maar. The ages of mapped vents range from ca. 3 Ma to Holocene (Table 1). More evolved lavas (mainly rhyolites) are generally associated with domes, tuff rings, and central volcanoes, whereas cones and maars are generally basaltic (Mazzarini et al., 1999; Corti et al., 2003; Lahitte et al., 2003a, 2003b; Mazzarini, 2004; Mazzarini et al., 2004; Abebe et al., 2005; Mazzarini, 2007).

Vent locations have been acquired by careful inspection of Landsat ETM+ images (e.g.,Goward et al., 2001; http://landsathandbook.gsfc.nasa.gov/handbook.html). Mosaics of ETM images were thus used to map a large portion of EARS (courtesy of the Maryland University Global Land Cover Facility). The images are geo-referenced to the UTM projection (zone 37 N –WGS84 for AD and MER data sets, and zone 35 S –WGS84 for VB data set) and displayed as RGB false color composites (with ETM+ band 7 in the red channel, ETM+ band 4 in the green channel, and ETM+ band 2 in the blue channel). The original spatial resolution of Landsat ETM+ images (pixel size) is 30 m. The 15 m spatial resolution of the Landsat ETM+ mosaic was obtained through a color transform using the 15 m geometric resolution of the Landsat ETM+ panchromatic band (Janza et al., 1975; Vrabel, 1996). Vent locations for the Afar Depression have been acquired by Mazzarini (2007). The vent locations for the Main Ethiopian Rift and in the Virunga Belt have been newly acquired through methods described in Mazzarini (2007) using Landsat ETM+ images in order to obtain accuracy in vent locations comparable to that for the Afar data set. The new data set for vents in the Main Ethiopian Rift contains basaltic vents that are located inside the rift valley in order to smooth out any influence of crustal structure beneath the rift margins, thus differing from the data set used by Mazzarini (2004).

Detection of vents by satellite image analysis provides thus a sample of the true vent population in the volcanic field. This sample clearly does not contain vents with a diameter smaller than a few pixels (i.e., less than 60 m), and vents that have been covered by younger volcanic products and continental deposits. To obtain a meaningful sample of the true vent population a few hundreds of vents should be detected (see below).

In order to check how good the sampling was, a random sample of 20% of the vents has been removed from the MER data set. The fractal dimension of this new data set remains practically constant (less than 0.01% of variation) and the error on the estimation of the upper and lower bounds of the distribution is 1%–2%. As discussed later, the sampling of hundreds of vents for the self-similar clustering gives a robust estimation of fractal exponent and cutoffs of the vent distribution.

Structural analysis of DEM and satellite images is based on the detection of rectilinear features (lineaments) of regional to continental scale. Lineaments consist of sharp tonal differences and alignments of morphological features (e.g., volcanic cones, triangular facets, portions of streams, aligned ridges, and crests). Within tectonically active areas, lineaments usually match the fracture network and the main fracture patterns (e.g., Wise et al., 1985; Abebe et al., 1998; Zakir et al., 1999). We assume that lineaments are the intersection of the fracture network with the Earth's surface; we thus consider lineaments as traces of fractures and will hereafter refer to lineaments as fractures. According to Pollard and Aydin (1988), we define as fracture any brittle rupture of rocks defined by two nearly planar parallel surfaces that meet at the fracture front and with the relative displacement small compared to the fracture length. In this scheme, joint and fault differs for the different type of the relative displacement of the two blocks (Mode I, II, or III).

Fracture trace maps of the studied sectors of the EARS have been derived by the analysis of shaded relief images derived from the digital elevation model and satellite images. For each studied area, images with different sun azimuth are derived from the SRTM digital elevation model with a cell size of 90 m (Farr and Kobrick, 2000; http://www2.jpl.nasa.gov/srtm/); in order to emphasize features a vertical exaggeration (3×) was used. The images have been rotated at random angles (in Microstation Bentley software environment), making fractures identification independent from image orientation and a 1:250,000 scale has been utilized.

Satellite images are analyzed using Google Earth software (http://earth.google.com/) where panchromatic SPOT images (http://www.spotimage.com) with pixel resolution ranging from 2.5 to 20 m are available. Surfing in Google Earth allows a complete rotation of images and the continuous variation of points of view thus providing a very useful tool for lineament mapping.

Geomorphologic features such as aligned ridges and valleys, straight drainage channel segments, and linear scarp faces have been mapped, identification of fracture traces in images is also based on color and tonal differences between different surface units (Figs. 3A–3C).

Clearly, in active volcanic fields magmatic products tend to cover previous volcanic constructs and topography (i.e., fault escarpment) and competition between magmatic and tectonic activity control the overall morphology of rift zones. The relationship between dikes, topography, and deformation in rifts zones has been a matter of debate. Dike intrusion focused subsidence and faulting (Rubin and Pollard, 1988), accommodated crust stretching, and prevented slip on faults (Parsons and Thompson, 1991), and the blade-like dike focuses deformation at the rift axis, leading to the formation of an axial rift valley (Behn et al., 2006). Examples of faulting and vents in Afar and MER are reported in Figure 3; it is apparent that cones as well as fault traces are well expressed.

We compared the azimuth distribution of detected lineaments with field data collected in the MER (Fig. 3D). Azimuth distribution of field data (joints and faults) has a main peak at N-S, N10°E and a secondary peak at N45°E-N50°E (Fig. 3E). Lineaments show a main peak at N20°E-N30°E and a secondary peak at N45°E-N50°E (Fig. 3F). There are ~10° between the main peaks of field and lineament data sets, while secondary peaks are similar. This difference could be easily explained tacking into account the different scale of observation, in fact: (i) the average length of detected lineaments is ~6 km (Table 2); and (ii) the most recent deformation (N-S fractures) is the most expressed at the outcrop scale.

Satellite and DEM derived images clearly show volcanic features as well as fractures and faults (Fig. 3), and vent positions were identified in a GIS environment (ArcView 3.3). The location of vents was identified on the satellite images with an accuracy of one pixel (i.e., 15 m).

Methods

The detected vents were analyzed in terms of their spacing and self-similar clustering. Vent spacing (or separation) is analyzed by computing the average minimum distance between vents. The coefficient of variation (CV) (Gillespie et al., 1999, and references therein) for the distribution of vent spacing describes the degree of vent clustering. CV > 1 indicates clustering of vents, CV = 1 indicates a random or Poisson distribution of vents, and CV < 1 indicates anticlustering (a regular distribution) of vents. CV is defined as
graphic
where s is the standard deviation and m is the mean.
The spatial distribution (self-similar clustering) of fractures and vents was analyzed by calculating the correlation exponent D applying the two-point correlation function method (Bonnet et al., 2001; Mazzarini 2004, 2007). For a population of N points (fracture barycenters or vent centers), the correlation integral C(l) is defined as the correlation sum that accounts for all the points at a distance of less than a given length l (Bonnet et al., 2001, and references therein). In this approach, the term C(l) is computed as
graphic
where N(l) is the number of pairs of points whose distance is less than l. If scaling holds, Equation (2) is valid, and the slope of the curve in a log(C(l)) versus log(l) diagram yields the D value.

Following Equation (2), the computed D value is valid for a defined range of distances (l). The distance interval over which Equation (2) is valid is defined by the size range. For each analysis, the size range of samples is in turn defined by a plateau in Δlog(C(l))/Δlog(l) (i.e., the local slope) versus log(l) diagram: the wider the range the better the computation of the power-law distribution (Walsh and Watterson, 1993). The derivation of the cutoffs is a crucial point and is generally not trivial, especially when the local slope does not show a regular and wide plateau. The choice of the zones where the plateau is well defined and the determination of the lower and upper cutoffs (Lco and Uco, respectively) are done according to Mazzarini (2004) by selecting the wider length range for which the correlation between log(l) and local slope is greatest.

A substantial database is considered to exceed the threshold of sampling lines (or vents) that is required in order to obtain stable statistical results; in fact truncation and censoring affect the computation of the fractal distribution (e.g., Bonnet et al., 2001). Indeed, this threshold has varied considerably among different authors (see André-Mayer and Sausse, 2007). About 200 values were recommended by Priest and Hudson (1976) and Bonnet et al. (2001). In this study, the used data sets ranging from ~350 to more than 5500 values are considered to be robust enough for statistical analyses.

The spacing of analyzed fracture networks and vents is listed in Table 2. In all the studied sites the average length and spacing of fractures are in the ~4 to ~6 km and ~1 to ~3 km ranges, respectively. Vent spacing for the analyzed data sets is very similar (~1 km) and it varies in the ~0.9 to ~1.3 km range; all the analyzed data sets show values of the coefficient of variation >1 (Table 2).

The spacing of fractures in AD varies in the <0.1 to ~28 km range with an average of 2.3 km; the fracture lengths have an average of 4.4 km in the 0.1 to ~37 km range (Table 2). Vents are clustered (CV = 1.1, Figure 4, Table 2) and show spacing values varying between 0.1 and ~18 km with an average of 1.2 km (see also Mazzarini, 2007). Both fractures and vents show self-similar clustering defined over an order of magnitude (Table 3). Fractures have D = 1.54, Lco = 2.5 km, and Uco = 23.6 km; vents have D = 1.42, Lco = 1.2 km, and Uco = 23.4 km. The upper cutoffs of distributions of both fractures and vents are very similar and the fractal exponent for the fractures is higher than that for vents.

In the MER site, spacing of fractures ranges between <0.1 and ~38 km with an average of 2.9 km; fracture length is in the ~0.2 to ~42 km range with an average of 5.7 km (Table 2). Vents are clustered (CV = 1.5, Figure 5, Table 2) and show spacing values varying between ~0.2 and ~21 km with an average of 1.3 km. Self-similar clustering of fractures is defined between Lco = 2.1 km and Uco = 41.7 km with a fractal exponent D = 1.61. Self-similar clustering of vents has exponent D = 1.17, Lco = 2.8 km, and Uco = 10.1 km; also in this case the fractal exponent of the fracture distribution is higher than that of the vent (Table 3). In this case the upper cutoff of the vent distribution is lower than the upper cutoff of the fracture distribution.

In the VB fracture spacing is in the ~0.2 to ~29 km range with an average of 1.5 km; fracture length ranges between ~0.4 and ~23 km and averages 4.3 km (Table 2). Vents are clustered (CV = 1.2, Figure 6, Table 2) and show spacing values varying between ~0.1 and ~10 km with an average of 0.9 km. Also for VB, self-similar clustering of fractures and vents is well defined; fracture has D = 1.85, Lco = 1.4 km, and Uco = 22.1 km. Vents in VB have distribution with fractal exponent D = 1.50, Lco = 1.0 km, and Uco = 21.6 km. In VB the upper cutoffs of fracture and vent distributions are similar and the fractal exponent of the fracture distribution is higher than that of the vent distribution (Table 3).

Summarizing, in the analyzed sectors of the East African Rift System both fracture networks and vents show well-defined self-similar clustering with the fractal exponent of fracture distribution higher than that of vent distribution (Table 3). The higher clustering of vents than the clustering of fractures is well defined considering some parameters, namely: the ratio (Sr) of spacing of fractures (Sf) to the spacing of vents (Sv), the ratio (Dr) of fractal exponent (Df) of fracture distribution to the fractal exponent (Dv) of vent distribution and the value of the fractal exponent (Dv) of the vent distribution.

The values of the Dr ratio are >1 for all the analyzed data sets (Table 4), consistently with the higher degree of vents’ clustering than fractures’ clustering. The Sr ratio shows values between 2 and 1.7 (Table 4) pointing to a higher closeness for vents than for fractures. Finally, we observe that, for the analyzed data sets, the higher the Sr ratio the lower the Dv value (Table 4).

Before we go further into the results, a few points have to be discussed. (1) Is a bi-dimensional fractal analysis valid to characterize a network that actually operates in three dimensions? (2) Is the isotropy of the crust a prerequisite for the model to work? (3) How does evolution of volcanic field account for crust thickness changes?

Most of the observations of fracture distribution (length, density, self-similar clustering) derive from the analysis of fracture trace maps, that is from analysis of 2D data. If size distribution of the fractures follows a power law, trace lengths in an intersecting plane are also power law with an exponent, a2D, equal to a3D 1 (e.g., Marrett, 1996; Piggott, 1997; Berkowitz and Adler, 1998; Bonnet et al., 2001). Similarly, the intersection of a 3D fractal by a plane results in a fractal with D2D equal to D3D − 1, according to fractal theory (Mandelbrot, 1982). Moreover, it has been observed that joints in thinly layered sedimentary rock typically span the mechanical thickness of the layer and the general linear relationship between spacing S and layer thickness H has been observed (e.g., Gross, 1993; Narr and Suppe, 1991; Renshaw, 1997). The general relationship S = f(H) results thus in fracture systems that can be idealized as two-dimensional with the fracture length defined along the layer (i.e., bedding plane, mechanical layer).

Self-similar clustering of fractures is performed by applying Equation (4) to points selected along the fracture traces, namely the fracture barycenters. As stated in Bour et al. (2002) the result of this analysis is valid whatever the point used to define fracture location (barycenter, fracture tips, or any point taken at random in the fracture). This is because the derivation of two-point correlation dimension is statistically dominated by the numerous small fractures, for which the error in the determination of the precise spatial location is relatively insignificant. For the same reason the censoring effects (fractures that intersect the system boundaries; Pickering et al., 1995) did not affect the derivation of the fractal dimension. Therefore vents could be assumed as points along the fracture traces and they are samples of the intersection of the backbone with the surface. Providing a large number of points (>200, see Bonnet et al., 2001), the sampling of the backbone is thus statistically robust.

Self-similar clustering does not require isotropy of the crust, magma exploits already formed fractures and creates by itself new ones to reach the surface. In the Lizard Ophiolite Complex, Cornwall (Jolly and Sanderson, 1997) and in the Mull swarm in Scotland (Jolly and Sanderson, 1995) dikes exploited favorably oriented preexisting fractures and created new ones. The actual fracture network will thus consist of all the structures that allow a hydraulic connection to operate in order to transfer magma from depth to surface (an overall vertical movement). Therefore, as a preliminary consideration, the anisotropies due to crust mechanical layering and occurrence of weakness planes (inherited structures) favor the extrapolation of 2D analysis to a 3D system. This is markedly different than the model of Takada (1994a, 1994b), where intruding magma creates the fractures (hydrofracturing) in an isotropic crust.

Volcanic fields evolve in time in terms of vent distribution, density, and age as, for example, in the Springerville (Arizona; Condit and Connor, 1996) and Michoacan-Guanajuato (Mexico; Hasenaka and Carmichael, 1985; Connor, 1987; Connor, 1990) volcanic fields.

When is a fully infinite cluster formed in a fracture network? Infinite cluster forms as the fracture network density is greater than the percolation threshold (see above). At the early stage there are only a few vents that clearly are connected to the magma reservoir at depth. System boundaries could be connected by an infinite cluster or by a unique bond that spans all of the system. The few vents formed in Afar during the volcano-tectonic crisis in September 2005 (Wright et al., 2006) are an example of a dike ~10 km in height that spanned the whole crust between the magma reservoir and the surface. In the proposed model the fully self-similarity and the strict relationships between the Uco and the distance between the source and the surface are achieved after some time, the time the network necessitates to reach the threshold density (ρc). We tentatively assume that the density at the percolation threshold is achieved after the crust in the volcanic field reached its fracture saturation. Indeed, after a certain amount of strain in layered sequences the ratio between fracture spacing and layer thickness becomes constant and does not change as strain increases (e.g., Wu and Pollard, 1995). In order to model works properly, a common depth for the magma reservoir is assumed for the analyzed vents. Clearly, the style of volcanism may change through time, generating more evolved magma and fissural magmatism and determining the formation of central volcanoes (e.g., Lahitte et al., 2003b; Mazzarini et al., 2004). In this case the distribution of volcanic vents may be controlled by strain rate and brittle/ductile layering of the crust.

The analyzed vents along the EARS show ages varying between 0 and 2 Ma (see above and Table 1), whereas available geophysical data image the current crustal structures. A question now rises, how does the vents self-similar clustering account for the continuous crustal stretching in a deforming crustal setting such as an active continental rift?

In order to answer this question as a first-order approximation, a pure shear deformation for the crustal thinning along the EARS is here assumed and an average Poisson's ratio ν~0.3 for the crust in east Africa (Zandt and Ammon, 1995) is used; the extensional rate times the duration of volcanism in the studied sectors of the EARS (i.e., 2 Myr) is considered the horizontal strain (εx). During the last 2 Myr extensional rate is ~20 mm/yr for the Afar Depression (AD; Eagles et al., 2002), and ~6 mm/yr and 2 mm/yr for the Main Ethiopian Rift (MER) and the Virunga Belt (VB), respectively (Stamps et al., 2008). The horizontal strain (εx) is related to the vertical strain (εy) by the simple relationship εx =−νεy, where ν is the crustal Poisson's ratio. Accordingly, since late Pliocene–early Pleistocene (ca. 2 Ma ago) crust stretching reduced the crust thickness by ~10 km in AD, ~4 km in MER, and ~1 km in VB. With the exception of Afar Depression, these values are of the same order of the errors in the computation of the Ucos (Table 3). An error of 15% for both upper cutoff and crust's thickness is thus assumed, accounting for the measurements’ precision (Tables 1 and 3) and of the role of the crustal stretching that acted in the EARS since 2 Myr.

As expected for a percolating cluster (i.e., the backbone) the clustering of the vents is higher than the clustering of the fracture network (Df > Dv; Tables 3 and 4). In the Afar Depression both fractures and vents have Uco values that fit in very well with the crust thickness (Table 3) as also reported for the vent distribution in the northern portion of the depression (n-AD; see Mazzarini, 2007).

In the Virunga Belt area (Table 3) the Uco value is less than 30 km as suggested by geophysical and geologic observations (Ebinger, 1989a, 1989b; Corti et al., 2002; Albaric et al., 2009). Anyway, also for the Virunga Belt area both fractures and vents have very similar Uco values (Table 3).

The analyzed vents in the Main Ethiopian Rift have a self-similar clustering bounded by aUco value (~10 km, Table 3) that is lower than the upper cutoff (~27 km) derived in the same area by Mazzarini (2004). This difference is due to the difference between the two analyzed data sets. Vents in Mazzarini (2004) included also cones in the western border area outside the rift valley and with ages older than 2 Ma (Mazzarini et al., 1999), whereas the data set used in this analysis is composed only by very young vents in the rift valley.

The EAGLE data show that the crust in MER is in the range 30–40 km (e.g., Keranen et al., 2004; Maguire et al., 2006; Keranen et al., 2009; Corti, 2009), well above the 10–11 km derived by the self-similar clustering of vents collected in the rift-floor. Moreover, EAGLE data also show that crust beneath magmatic segments in the MER (see Ebinger and Casey 2001) is strongly modified by huge magma injection generating strong rheological layering (Corti, 2009 and references therein). We infer that when the crust is strongly modified by localized magmatic processes, as in the case of northern MER, vent self-similar clustering record local crust layering induced by magma emplacement within the crust (i.e., ~10 km is the thickness of the upper competent layer).

Geophysical and geochemical data (Keir et al., 2006; 2009; Rooney et al., 2007; Bastow et al., 2005) suggest that MER represents the transition from a broad mechanical crustal extension (faulting and stretching) toward narrow zones of magma intrusion prior to the onset of sea-floor spreading (e.g., Ebinger and Casey, 2001; Buck, 2006). On the other hand, velocity structure of the uppermost mantle beneath MER shows no correlation with crustal magmatic segmentation (Bastow et al., 2005; 2008), suggesting that the different distribution of crustal deformation and magmatism are not related to mantle features. Corti (2008), on the base of analogue modeling, concluded that rift evolution and segmentation are controlled by rift obliquity independently by magmatic processes. Moroever, Agostini et al. (2009) suggested that oblique reactivation of pre-existing crustal weakness plays a major role in the rift evolution and architecture.

Volcanism clearly clusters along magmatic segments in the rift reflecting the localized strong modification of the crustal structure by magma intrusion. To define whether magmatic segmentation is controlled by magmatism or rift obliquity and crustal inherited structures is beyond the aim of this contribution. Anyway, in the Main Ethiopian Rift, since Pleistocene, deformation and volcanism focused in the rift floor (see Ebinger and Casey 2001; Corti, 2009 and references therein) and vent self-similar clustering appears to record such as tectono-magmatic setting (Table 3).

In the Main Ethiopian Rift the fracture network shows Uco (41.7 km) well matching the derived regional crust thickness of ~40 km (e.g., Keranen et al., 2004; Maguire et al., 2006), while Uco for vents is 10.1 km. The apparent discrepancy between fractures and vents Ucos in MER could be explained by differences in the extent of the two data sets: (1) fractures span the whole rift architecture from the border to the floor and encompasses a long history in the rift evolution, that is from ca. 10 Ma to the present (Corti, 2009, and references therein) thus accounting for an initial stage of broad mechanical crustal extension; and (2) vents are localized in the rift floor and account for the Pleistocene to Holocene evolution of the rift with focused magma intrusion in lower crust. The occurrence of weak and partially molten lower crust beneath the MER (Keranen et al., 2009) could thus result from recent rift evolution and is well recorded by volcanism distribution.

The classical segmentation of the MER in northern, central, and southern segments, based on different geophysical, geological, and deformation characters (Bonini et al., 2005; Keranen and Kemplerer, 2008; Corti, 2009), is not here considered being the vents mainly located in the central and northern MER segments. Future analysis will focus on the characterization of the MER segmentation.

In the Afar Depression although the average crust thickness is ~23 km, as calculated from vent self-similar clustering, local variations exist (see Mazzarini, 2007). Differently, in the MER vent self-similar clustering records very localized crust structure, and only vents outside the rift floor record the average crust thickness (e.g., Mazzarini, 2004). In fact, the rifting process in Afar is more evolved than MER and the magmatic modification of the lithosphere is widespread whereas magmatic modification of the crust is only localized in narrow zones in the MER (e.g., Ebinger and Casey, 2001).

The crust thickness derived from geophysical data is plot against the Uco values (Fig. 7) defined for the self-similar clustering of the vents in the sectors of the EARS together with data from the Pali Aike volcanic field in southernmost Patagonia (Mazzarini and D'Orazio, 2003), the Payen volcanic field in the Andean foreland (Mazzarini et al., 2008), and in the northern sector of the Afar Depression (Mazzarini, 2007). Despite the differences in the geodynamic settings of the analyzed volcanic fields, a linear correlation (Uco = 1.0481CT − 1.8367, R2 = 0.9856 where CT is the crust thickness in kilometers) is clearly defined, indicating a strong correlation between the upper cutoff of the fractal distribution of vents and the crust thickness in the volcanic fields.

Fracture networks and monogenetic vents have been analyzed in volcanic fields in three sectors along the East African Rift System: the Afar Depression, the Main Ethiopian Rift, and the Virunga Belt, respectively (Fig. 2). Despite differences between the three investigated sectors in crust's stretching, crust's thickness, and in the intensity of the volcanism, both fractures and vents show a clear self-similar clustering. In particular, vents in all three analyzed sectors in the EASR as well as in other volcanic fields in Afar and southern South America (Mazzarini and D'Orazio, 2003; Mazzarini, 2007; Mazzarini et al., 2008) are more clustered (higher self-similar clustering) than the fracture network, that is Df > Dv (Tables 3 and 5). This is expected in the Percolation Theory for a percolating cluster (i.e., the backbone) in respect to the network (Orbach, 1986; Stauffer and Aharony, 1992). Accordingly, for a percolating fracture network, the backbone in the infinite clusters is the actual fluid pathway, thus the percolation correlation length (ξ) represents the distance between the system's boundaries (i.e., the source and the surface for vents). As a consequence, in percolating networks, the percolation correlation length (ξ) represents the upper cutoff (see Orbach, 1986).

Notably, there is a clear correspondence between crustal thickness in the sites of analyzed volcanic fields (derived by independent geophysical data sets) and the upper limit of the computed size ranges (upper cutoffs) defining the fractal distributions of the vents (Table 3 and Fig. 7). These data suggest that the value of the upper cutoff derived by the analysis of the vents self-similar clustering defines the maximum distance in the backbone covered by fluids when rising to the surface from a reservoir at depth.

We thank G. Corti and an anonymous reviewer for valuable comments and suggestions that improved the manuscript. I.I. benefited from the project FIRB “Sviluppo di nuove tecnologie per la protezione e difesa del territorio dai rischi naturali (FUMO) - UR14 (Scientific Responsible Valensise G.)” funded by the Italian Ministero dell'Istruzione, dell'Università e della Ricerca.