The distribution of volcanism in the Beta‐Atla‐Themis region of Venus: Its relationship to rifting and implications for global tectonic regimes

A new analysis of the spatial relationships between volcanic features and rifts on Venus provides new constraints on models of planetary evolution. We developed a new database of volcanic features for the Beta‐Atla‐Themis (BAT) region and used nearest neighbor measurements to determine relationships between different types of volcanic features and the rifts. Nearest neighbor analysis shows that all the dome‐type and corona‐type subpopulations tend to cluster. Rift associations were inferred from the deviation of a feature's population distribution (as a function of distance from rift) from that of a random population. Dome‐type features in general have no discernible relationship with rifts. Most corona‐type features have a strong association with rifts, with intermediate and large volcanoes also tending to occur close to or on rifts. Shield fields, on the other hand, tend to occur away from rifts. Our new evidence supports classifications of rifts on Venus into different types, possibly by age, with a shift from globally dispersed (more uniform) volcanism toward the more rift‐focused distribution, which suggests a shift in tectonic regime. Our observations are consistent with recent models proposing the evolution of Venus from a stagnant lid regime to a subcrustal spreading regime. We also present evidence for a failed rift on Venus and note that this process may be analogous, albeit on a larger scale, to a proposed model for the evolution of the East African rift system.


Introduction
The nature of Venus' long-term tectonic evolution is an enduring enigma, for which competing hypotheses have been proposed and extended [e.g., Basilevsky and Head, 1998;Guest and Stofan, 1999;Ivanov and Head, 2015]. The evolutionary models in the following summary have different implications for the likely distribution of geological features at the surface. Here we use spatial analysis to explore the relationships between volcanic and tectonic features and to assess the validity of these models. We discuss how the distribution of volcanic features on the surface of Venus informs us of the nature of internal driving volcano-tectonic processes, the nature of planetary resurfacing processes, and how these affect the long-term tectonic evolution of Venus.
On Earth, the lithosphere is continuously recycled via the process of plate tectonics. Lithospheric plates are created at mid-ocean ridges and are then transported to subduction zones where they are destroyed and recycled. Continental crust is less dense than oceanic crust, so continental materials commonly escape subduction. On Venus the relatively dry, more viscous mantle inhibits Earth-like plate tectonics McKenzie, 1996, 1998]; indeed, none of the characteristic features associated with plate tectonics have been observed on Venus, suggesting that these processes may not occur at all. Experimental work does indicate, however, that localized plume-induced subduction may potentially occur on Venus [Davaille et al., 2017], suggesting an ongoing resurfacing mechanism. Figure 1. Idealized sketch summaries of (a) terrestrial plate tectonics and (b) a stagnant lid regime on Venus [e.g., Solomatov and Moresi, 1996] as discussed in the text. The dimensions and layer thicknesses (for which Venus estimates vary widely) are not to scale and are intended to compare the overall structure rather than absolute quantities. elastic, "rigid" mantle lithosphere passes with increasing depth to high-viscosity convecting mantle, transitioning when a rheological boundary is crossed [Solomatov and Moresi, 1996]. Ghail [2015] recently proposed a subcrustal spreading model in which a petrological transition at the crust-mantle boundary within the "stagnant lid" lithosphere may allow these two lithospheric components to become separated at a detachment horizon. Regions of upwelling and downwelling mantle may then drive a form of plate tectonics beneath the crust, rejuvinating the lithosphere. On Earth, spreading at rifts facilitates decompression melting and migration of magma to shallow crustal depths resulting in volcanic activity at the surface coincident with, or very close to, the rift axis [e.g., McKenzie and Bickle, 1988;Langmuir et al., 1992;Perfit and Davidson, 2000]. Under a stagnant lid scenario one would not necessarily expect an Earth-like distribution of rift-related volcanic features. Venus requires further study in order to fully understand potential interactions between the underlying tectonics and volcanism at the surface.
The observation that the global distribution of impact craters is highly uniform [Phillips et al., 1992;Schaber et al., 1992;Strom et al., 1994], and by implication the surface is everywhere of a similar age, led to the view that this stagnant lid scenario is punctuated by periods of global resurfacing [Parmentier and Hess, 1992]. These are described in terms of global lithospheric recycling events as proposed in  and further developed in Turcotte [1995Turcotte [ , 1996, Turcotte et al. [1999], and Romeo and Turcotte [2008]. The crater distribution has also been found to be statistically consistent with equilibrium resurfacing, whereby resurfacing is ongoing and localized [Phillips et al., 1992;Bjonnes et al., 2012].
The "directional" model of Venus evolution developed by Basilevsky and Head [1998], and subsequently extended Head, 2000a, 2002;Ivanov and Head, 2011, proposes a series of global epochs characterized by distinct geological regimes, with each period representing a globally synchronous stratigraphic marker characterized by the dominant volcanic/tectonic activity. This directional model does, however, face a number of challenges. Hansen [2000] raised doubts about the methodology for determining a globally uniform stratigraphy. Further, Guest and Stofan [1999] and later studies [Addington, 2001; cite numerous examples evident in the Magellan data that suggest that the directional model markers are not necessarily globally synchronous and may have been formed over longer periods, and more dispersed through time in a "nondirectional" model. In the directional model, the relative ages of extensional structural features (rifts) are defined as either "young" or "old." These definitions are detailed in Krassilnikov et al. [2012] and are based on the assumption that the regional plains unit (defined in Basilevsky and Head [2000a]) is a globally synchronous marker. Rifts defined as "old" predate the regional plains and must demonstrate evidence of embayment by the plains material. Similarly, the rifts defined as "young" postdate the regional plains and must display evidence of having been formed after the plains emplacement, for example, cross cutting the unit. Full definitions of rift classification, and the occurrence and types of structures identified as rifts can be found in Krassilnikov et al. [2012]. Even assuming the plains unit as a global marker, inherent uncertainty in rift ages is introduced as the actual duration of plains emplacement is unknown. For example, a rift may cut early plains material leading to it being defined as "young". Later in its emplacement, the plains material might later embay another rift leading to it being defined as "old" although it may in fact be younger if it formed later in the plains emplacement episode. In addition, use of the regional plains marker as a global relative age boundary is open to debate. Without this assumption the age definitions are limited to geographically localized relationships. The terms "young" and "old" are therefore best treated as both relative and geographically localized terms, with potential global relevance.
While we use the terms "young" and "old" rifts in the rest of the paper, this is to test their relevance when describing the distribution of volcanic features on Venus and is not predicated on acceptance of a global directional model of Venus evolution. For example, in contrast to directional models, the purely nondirectional models (equilibrium resurfacing and stagnant lid) predict that there should be no systematic differences between "young" and "old" rifts. The subcrustal lid model predicts that some rifts (Parga Chasma and Hecate Chasma, Figures 2b and 2d) overlie regions of faster lid extension and have more voluminous magmatism over a greater distance from the rift axis than rifts overlying slower lid extension rates (Devana Chasma and Dali-Diana Chasmata). Inspection of Figures 2b and 2d may contain evidence in support of this as the former (predicted to be fast) rifts are characterized by broad regions of predominantly "young" rifting, in contrast to the latter (predicted to be slow) rifts, which appear to be relatively narrow "young" rifts (Devana Chasma) or contain higher proportion of "older" rift segments (Dali-Diana). These models can be tested by the spatial analyses presented in this paper.
To understand the link between surface rifting and volcanism on Venus, the nature of volcanism must first be characterized. The SAR (synthetic aperture radar) data returned from Venus by the NASA Magellan mission in the 1990s heralded the beginning of our detailed understanding of volcanic features on Venus [e.g., Guest et al., 1992;Head et al., 1992;McKenzie et al., 1992;Pavri et al., 1992;Stofan et al., 1992;Gregg and Greeley, 1993]. These observations built on the legacy of Earth-based observations as well as the Venera and Pioneer programs. Most notable among these studies in terms of characterizing the classification and global cataloguing of features was Head et al. [1992]. This work classified volcanic edifices ranging from a few kilometers in diameter, commonly occurring in clusters known as shield fields, to immense volcanoes hundreds of kilometers in diameter. The eruptive nature of volcanism on Venus is thought to be dominantly effusive due to the observed surface geomorphology , high surface pressure [Taylor, 2010], and largely mafic composition observed by the Venera and Vega landers [Kargel et al., 1993]. There have, however, been recent studies suggesting that explosive activity may be possible under certain circumstances based on evidence from theoretical modeling [Airey et al., 2015] and by scrutinizing the radar properties of volcanic deposits [Ghail and Wilson, 2013], although this is likely to be rare by comparison.
Many of the intermediately sized volcanoes show a characteristic wide, flat-topped morphology and were classified as steep-sided domes. Some of these display mass-wasting features on their flanks leading to their subclassification as fluted domes. Isolated calderas not associated with an obvious volcanic edifice were also identified and catalogued, as well as another subclassification of intermediately sized volcanoes known as anemones, which derive their name from the radiating petal-like radar bright lava flows. In addition to volcanoes, Head et al. [1992] also catalogue an apparently unique family of volcano-tectonic features, first described by Barsukov et al. [1984], known as coronae, novae, and arachnoids displaying, to a greater or lesser extent, concentric fractures (coronae), radial fractures (novae), or a combination of the two (arachnoids); these features are collectively termed "coronoids" henceforth in this study.
The wide array of shield and cone volcanoes on Venus occurs in a variety of environments, providing clues to their method of formation. The largest volcanoes may achieve their large size by virtue of the apparent lack of plate motion concentrating plume activity beneath a particular region of the crust for an extended AIREY ET AL.
VOLCANISM AND GLOBAL TECTONICS ON VENUS   [Stofan et al., 2001a]. They frequently coincide with regional topographic rises and rift junctions [Stofan et al., 1995]. The large range of volcanic edifice sizes suggests a correspondingly large range in the scale of mantle features impinging on the lithosphere, with the broad geographic distribution interpreted as evidence for mantle plumes more numerous than on Earth [Stofan and Smrekar, 2005]. The nature of the smallest volcanoes, and their tendency to occur in discrete fields, implies that they may be formed by distinct individual small-scale eruptions occurring over a magma source [Ivanov and Head, 2004]. Basilevsky and Head [1998] describe them as relatively old features according to their geological history model, predominantly predating the regional plains stratigraphic marker described in Ivanov and Head [2004].
Several theories regarding the formation mechanisms for the steep-sided and fluted domes have been proposed. Pavri et al. [1992] suggested two methods of formation, extrusion of an evolved, high-viscosity magma, or of a volatile enhanced basaltic foam. Later work developed a model describing a basaltic lava, with a cooling crust that is continually fracturing and annealing while spreading laterally into the distinctive pancake shape [Stofan et al., 2000]. Most recently, detailed modeling of dome emplacement in "constant volume" and "time-variable volume" regimes has been carried out [Quick et al., 2016], greatly improving our understanding of these features. This work suggests a composition likely to be comparable with terrestrial basaltic andesite, rather than rhyolite, and emplacement times of 1.6 to 16 years. Previous qualitative studies note that these features commonly occur in clusters [Pavri et al., 1992;Ivanov and Head, 1999].
The coronoid family of volcano-tectonic features has attracted a wealth of research into their formation mechanisms and distribution. It has long been suggested that these features occur due to the upwelling of magma diapirs or plumes Squyres et al., 1992;Stofan et al., 1992;Smrekar and Stofan, 1997]. It has been proposed that these may form a genetic continuum, with radial fracturing characteristic of novae associated with mantle upwellings and fracture annuli characteristic of coronae associated with subsidence . Stofan et al. [2001b] break down the corona population into two subdivisions: type 1 coronae with >50% fracture annuli and the much more poorly defined type 2 coronae with <50% fracture annuli, the former occurring predominantly on rifts and fracture belts and the latter occurring predominantly on the plains. Models of coronae formation include the delamination and deformation of the lithosphere via magma upwelling [Smrekar and Stofan, 1997], the deformational response of the lithosphere when loaded via magmatic intrusion [Dombard et al., 2007], upwelling-induced crustal convection [Gerya, 2014], and ring-like dripping via Rayleigh-Taylor instabilities of a dense layer at the base of the lithosphere at plume margins [Piskorz et al., 2014]. Arachnoids are thought to be formed by uplift and relaxation by magmatic diapirs formed over mantle plumes [Aittola and Kostama, 2000;Krassilnikov, 2002].
The global distribution of coronae on Venus is nonrandom Squyres et al., 1993] with high concentrations noted in the Beta-Atla-Themis (BAT) region (Figures 2b and 2d) [Squyres et al., 1993], and at midlatitudes . Martin et al. [2007], however, found their distribution to be indistinguishable from random within 1500 km of Parga Chasma. An apparent genetic link between coronae and rift systems is well documented Baer et al., 1994;Smrekar and Stofan, 1997;Martin et al., 2007;Piskorz et al., 2014] but has yet to be adequately explained. Arachnoids on the other hand have previously been noted to tend to cluster on the plains, in contrast to novae, which tend to follow sparse chains [Aittola and Kostama, 2000]. Analysis of a subset of ∼100 coronae carried out by Krassilnikov et al. [2012] showed that 97% of the sample predated the regional plains stratigraphic marker, but almost half displayed post regional plain emplacement volcanism, suggesting that coronae may be long-lived features. Assuming some directionality in Venus evolution, they also proposed type 2 coronae to be generally stratigraphically older than type 1 coronae. This late-stage volcanism is identified as a component of the corona formation model described by Gerya [2014].

The Beta-Atla-Themis (BAT) Region of Venus
This work focuses on the region of most intense rifting, the BAT region, which is a large, roughly triangular, area following three large rift zones connecting Beta, Atla, and Themis Regiones (Figure 2b). The emphasis here is on localized spatial relationships between key distinct feature types, with a view to helping define the broader tectonic context when considered together. The area is well mapped and the volcanic features have been catalogued Crumpler and Aubele, 2000;Stofan et al., 2001b] and includes a region of proposed currently active volcanism [Shalygin et al., 2015]. The defined area contains all feature types under investigation and, in most cases, in sufficient numbers with which to derive significant numerical results. Although the The BAT rifts are zones characterized by extension, as evidenced by densely packed linear arrangements of extensional features such as graben and normal faults, with individual features often attaining tens of kilometers in width and hundreds of kilometers in length Head, 2011, 2015]. These observations, along with the fact that they occupy topographic highs [Ivanov and Head, 2011], suggest that they formed in association with mantle upwelling beneath the lithosphere, which in turn drove volcanism and rifting at the surface [Ivanov and Head, 2015].

Data Sources and Mapping
The base geospatial environment onto which the volcanic features and rift zones are mapped is a composite of the Cycle 1 SAR radar maps retrieved by the NASA Magellan mission. The original spatial resolution of this product was 110 m along-track resolution and 101 to 250 m cross-track ground resolution (dependent on latitude). This study uses the standard full resolution data product, which was resampled from the original to 75 m gridded resolution as part of the postmission data processing (described in Saunders et al. [1992]). The study area is defined as 40 ∘ N to 50 ∘ S and 90 ∘ E to 310 ∘ E, to capture the full BAT region. ArcGIS 10 was used to process the mapping and for analysis.
The coordinates from existing volcano databases were used to map the volcanic features prior to the geospatial analysis. The volcano-tectonic feature catalogue of Crumpler and Aubele [2000], including all the features mapped in Head et al. [1992], was reduced down to the areal extent of the study region and added to the base map. Subsequently, the more comprehensive database of coronae outlined in Stofan et al. [2001b] was added, which resulted in numerous conflicting classifications. Resolving this followed the convention of allowing the latter definitions to supersede the former while retaining the former's definitions when they were absent from the latter database. In addition, other observed features absent from both databases were mapped on in order to expand the hybrid database. These additional features include previously unidentified, or omitted, structures of all types shown in Figure 2d with the exception of fluted domes and calderas.
Shield volcano size classifications are defined based on their diameter as large (≥100 km), intermediate (≥20 km to <100 km), and small (<20 km). Small volcanoes are so numerous that they are not tabulated in the database unless they occur within shield fields, which are included as single entries. As the original database is only accurate to the nearest 0.5 ∘ , each point was subsequently centered on the feature in question. The rift sections were mapped onto the base map using the classification in Krassilnikov et al. [2012] incorporating "old" (predating regional plains unit) and "young" (postdating regional plains unit) according to previously defined stratigraphic relationships [Basilevsky and Head, 2000a] and stored as polyline vector data sets. The rifts were mapped in this manner in order to represent rifts in the region in a way that could be used to test the concepts inherent to the directional model.

Analytical Techniques
Spatial analysis of the mapped volcanic features allows quantification of geospatial relationships to answer the following questions: (1) how do certain volcanic features occur relative to others of their own kind?, (2) how do they occur relative to features of other kinds?, and (3) how do they occur relative to rift axes?
The nearest neighbor index (NNI) was calculated for each type of volcanic feature. This quantity determines the ratio of the observed mean distance between points (D O ) to the expected mean distance between points in a random distribution (D E ). An NNI of 1 indicates a random distribution, while NNI>1 or <1 are indicative of, respectively, either more uniform (or dispersed) or more clustered distributions [Mitchell, 2005]. Previous applications of this method to analyze the clustering of geospatially referenced data include the distribution of craters on Venus [Hauck et al., 1998] and the analysis of earthquakes in the Red Sea [Al-Ahmadi et al., 2014]. Equations (1) to (3) describe how ArcGIS 10 calculates the NNI.
where d i is the distance between point i and its nearest neighbor, n is the number of points, and A is the total area of the study region. The feature distribution data used in this study is projected onto a geographic coordinate system for Venus provided in ArcGIS 10; therefore, the distances quoted are for that over the surface of a sphere as opposed to a straight (chordal) line through the surface of the planet. As a point of reference for each NNI calculated, 50 random distributions of points for each sample size, over the same area, were generated using the ArcGIS built in algorithm, along with their own calculated NNI values. The mean random NNI for each sample size was calculated for comparison with the real populations with the corresponding sample size. In addition to calculating the NNI for each feature type, the distance to nearest neighbor value for each point in both the observed and random data sets was also stored so that the population distribution histograms could also be produced. Of the 50 simulated distributions, the one with the NNI closest to NNI = 1 (true random) was selected for this purpose.
To determine the features' relationships with the rift axes, the distance of each feature point on the map from the nearest rift (of any type) was measured in order to identify any preferred distance from, or affinity for occurring near, a rift of any age. These distances were measured once more with just the distance to the "young" rifts considered in order to provide comparative data sets that could be analyzed separately based on relative rift formation age. This allows us to consider whether our data support a more directional or nondirectional model for Venus' evolution, which will become an important part of the discussion. Population distribution histograms were also produced for these data as well as the corresponding results recorded using the previously selected random data sets for comparison.
In order to display on the map where any clustering may be concentrated, feature kernel density maps were produced. This was achieved using the number of similar features within a given search radius (10 3 km) from each point, giving a unique value to each cell (gridded to 5 km) as a measure of feature density. When displayed using a color ramp, the density magnitude per cell (in features per km 2 ) for each feature type was then projected onto the base map in order to produce the feature density layer. The chosen search radius and shading intensity were selected in order to enable qualitative assessment of feature distribution and grouping tendencies between populations. They provide sufficient resolution to highlight regions where even relatively small clustered or linear groupings, and in particular where more closely distributed features, occur when considered on a scale comparable to the global rift distribution.

Results
The NNI values for the volcanic feature data sets described in section 2.2 are shown as the red crosses in Figure 3, sorted into their respective family groups as described in the key. Note that they fall predominantly in the lower half of the plot. With the exception of Group III, all the observed data sets fall, to a greater or lesser degree, comfortably outside the range of random distributions and within the "clustered" field, statistically significant at p = 0.05. The blue dots represent the NNI values from the 50 simulated random distributions; their associated mean values are shown as the black circles. All the mean simulated random NNI values fall between 1 and 1.24, indicating a tendency for truly random distributions with small sample sizes to shift toward a uniform distribution. This is due to edge effects resulting from the fact that the nearest neighbors of features close to the study region's boundary may be outside the domain within which it is calculated. For smaller sample sizes, there will be a commensurately larger effect of this as the likelihood of the nearest neighbor being outside the domain will increase with D E (the expected mean distance between points in a random distribution). The overall effect of this is for the NNI of a feature whose actual nearest neighbor is outside the domain to be larger than it should be. As the true populations are affected in exactly the same manner as the random populations of corresponding sample size against which they are compared, relative comparative analysis remains valid whereas absolute values must be considered with this in mind. Table 1 shows the sample sizes and NNI for all the feature classes.
The extent to which random populations with smaller values of n fall in the "dispersed" field can be shown to be dependent on the value of n as shown in Figure 4. The tendency for the volcanic feature data to fall clearly within the "clustered" field, even given this skew toward NNI > 1, accentuates the nonrandom nature of these population distributions. . NNI values for mapped Venusian volcanic feature data sets (red crosses) and the results of 50 random distribution simulations of sample sizes corresponding to the feature data set with which it is plotted (blue dots). Black circles represent mean values from the random data sets. The classifications I-III refer to the "Dome," "Coronoid," and "Volcano" families. Figure 5 shows the histograms of nearest neighbor distances for all feature types (grey) overplotted with the histograms of the random populations selected from the 50 examples to have an NNI nearest to NNI = 1 (blue). In most cases these best fit a gamma or exponential distribution, although some of the random populations are indistinguishable from a normal distribution. With the exception of shield fields, the peaks of the observed frequency distributions always occur at smaller nearest neighbor distances than those of the random distributions, again indicative of a tendency for clustering. Figures 6 and 7 show the histograms for distance to nearest rift (both red and black lines in Figure 2d), and separately for distance to nearest "postplains," or "young", rift (red lines only in Figure 2d) as proposed by Krassilnikov et al. [2012] according to the previously defined stratigraphic relationships of Basilevsky and Head [2000a], for each volcano-tectonic feature class along with fitted frequency distribution curves. All distributions fit an exponential distribution and show the affinity for and/or preferred distance from the rift axes identified in Figure 2d.   Figure 4. Plot of NNI against sample size for randomly generated point distributions. These data best fit a curve with the exponential equation y(x) = 0.953 + x −0.44 , giving a maximum R 2 = 0.98. Best fit curve selected based this R 2 value. The intercept is not quite 1 as would be expected due to the lack of high-n data sets.

Group I: Domes
This first group comprises the steep-sided and fluted domes. They show clear similarities in terms of their strong clustering, each displaying low NNIs of ∼0.72; when considered together as a single grouping, they have an NNI of ∼0.6 ( Figure 3). The exponential nearest neighbor frequency distribution and its large deviation from the normal and gamma distributions of the random populations (Figure 5a) highlight this strong tendency to occur in clusters, which is also evident on inspection of the feature density maps in Figure 8. It is likely that the features in this group share similar formation mechanisms; together, they represent a distinct subgroup in the database and warrant classification and analysis as a group feature set.
Steep-sided domes show rift distance frequency distribution curves almost identical to the random distribution curves (Figures 6aii and 7aii). Fluted domes, however, deviate more from random where "all rifts" are considered (Figures 6ai and 7ai), indicating a preference to occur near "old rifts." The tendency for fluted domes to occur near rift is apparent on inspection of Figure 8a.

Group II: Coronoids
This group comprises the features making up the coronoids: the coronae, novae, and arachnoids. Figure 3 shows that this group is marginally more uniformly distributed as a family compared with Group I while still being strongly clustered, most notably the arachnoids (also see Figures 5bii and 9b). Histograms in Figure 5b are more similar to the random distribution than the domes; however, peaks in the data show that individual features are consistently closer to their nearest neighbors in general than the randomly generated data. The Novae (Figures 5bi and 9a) tend to cluster into loosely linear groupings, as opposed to discrete clusters. The type 1 coronae appear to be randomly distributed on a small scale but are clustered across the BAT region at a large scale (Figure 9c).
The histograms in Figure 6b show another important distinction of Group II compared with Group I: whereas Group I features may be showing some tendency to deviate from the random distribution, Group II features display a much stronger deviation in the form of steepening of the distribution curve away from the random distribution. This is seen most strongly for the novae (Figure 6bi) but is also observed for the arachnoids and type 1 coronae (Figures 6bii and 6biii). Type 2 coronae are an exception and seem to occur off rift and somewhat clustered (Figure 9d), as previously noted by Stofan et al. [2001b]. The coronoids' general occurrence very close to or on rifts indicates a close link between this type of volcanic feature manifesting at the surface and the spatial association with the rifts. Figure 9 shows the feature density maps for Group II, highlighting the grouping characteristics and tendency to occur on rift.

Group III: Volcanoes
This group is more morphologically diverse, with isolated calderas and volcanoes ranging in size from small volcanoes below the imaging resolution of the radar to features of up to >1000 km in diameter.
Calderas are distributed as shown in Figure 10a. There is some clustering, with an NNI of ∼0.72 and a significant jump in the histogram at a nearest neighbor distance of <250 km showing ∼40% of all calderas occurring at least this close to another caldera (Figure 5ci). Figure 10a also shows that like the domes, most calderas lie in a near-rift setting. The tendency to occur in groups could be indicative of their occurrence over off-rift plumes, with one or two rare on-rift exceptions. The distance from rift data, in terms of both "all rifts" and "young rifts", is indistinguishable from random (two-sample Kolmogorov-Smirnov p values of 0.90 and 0.76, respectively). This indicates a formation process completely independent of rift location. Their location on both the proposed older shield plains as well as the proposed more stratigraphically recent regional and lobate plains units Distributions are shown with gamma, exponential, or normal distribution curves as determined by best fit to the data; black and blue curves correspond to the observed and random data, respectively. All data are plotted on figures of a fixed x axis interval width up to 5000 km to enable comparison between different data sets. All mapped data sets are divided into 24 bins of equal width. Where the ranges of the data sets differ between the random and observed values, the random values are binned to intervals to match those of the observed data, rather than maintaining the same number of bins overall. This, as well, is intended to enable comparison between data sets. The y axis data were fitted to plots of equal height to highlight the trends observed in the population distributions, rather than the absolute number of features. [Ivanov and Head, 2015] suggests that even assuming a directional model for Venus' evolution, these were not temporally restricted to any particular period of relative timing and continued to form throughout the proposed transition in volcano-tectonic regime. The structure of these calderas, featuring faulting and subsidence, implies a tectonic component caused by the removal of a magma source, which suggests that shallow bodies of magma may have continued to affect surface features throughout the regional plains emplacement.
The histogram for the intermediate volcanoes (≥20 km to <100 km diameter, Figures 5cii and 10b)  Distributions are shown with gamma, exponential, or normal distribution curves as determined by best fit to the data; black and blue curves correspond to the observed and random data, respectively. All data are plotted on figures of a fixed x axis interval width up to 5000 km to enable comparison between different data sets. All mapped data sets are divided into 24 bins of equal width. Where the ranges of the data sets differ between the random and observed values, the random values are binned to intervals to match those of the observed data, rather than maintaining the same number of bins overall. This, as well, is intended to enable comparison between data sets. The y axis data were fitted to plots of equal height to highlight the trends observed in the population distributions, rather than the absolute number of features.
volcanoes). The majority appear to occur close to, rather than on, rifts. Figure 6cii illustrates this with the striking peak in the data at ∼160 km, a trend echoed in the fluted dome data (Figure 6ai) suggesting that these intermediately sized features may form via some rift related, albeit off-rift, mechanism.
The shield field histograms (Figures 5ciii, 6ciii, and 7ciii) show nearest neighbor and distance from rift distributions very close to the corresponding random distributions. However, in terms of clustering, Figure 5ciii may mask some distribution characteristics observed in Figure 10c where, for example, the signals from the more uniformly distributed eastern portion of the study area and the more clustered western portion of the study area are lost when considered together. The vast majority of shield fields occur on the proposed preregional AIREY ET AL. VOLCANISM AND GLOBAL TECTONICS ON VENUS Distributions are shown with gamma, exponential, or normal distribution curves as determined by best fit to the data; black and blue curves correspond to the observed and random data, respectively. All data are plotted on figures of a fixed x axis interval width up to 5000 km to enable comparison between different data sets. All mapped data sets are divided into 24 bins of equal width. Where the ranges of the data sets differ between the random and observed values, the random values are binned to intervals to match those of the observed data, rather than maintaining the same number of bins overall. This, as well, is intended to enable comparison between data sets. The y axis data were fitted to plots of equal height to highlight the trends observed in the population distributions, rather than the absolute number of features.
plains "shield plains" stratigraphic unit as defined in the directional model of Venus geology [Ivanov and Head, 2011, in areas remaining unaffected by the regional plains emplacement along with the "old" rift features with which they commonly coincide.
The large (≥100 km diameter) volcanoes show a distribution almost identical to the randomly generated sample ( Figure 5civ); however, inspection of Figure 10d suggests that this may be due to the same "actual" distribution dichotomy that affects the shield fields, that is, both uniform and clustered subpopulations. The distance from "all rifts" data in Figure 6civ display a subtle divergence from the random sample due to a large spike in the leftmost bin with 39 large volcanoes (∼35%) occurring within ∼50 km of a rift segment. This deviation from random is even more prominent where "young" rifts are concerned (Figure 7civ) suggesting a closer association with the most recently active rift geometry according to the directional model. The large volcano classification encompasses a rather broad size range (100-1000 km diameter), so a new subclassification, very large volcanoes, is shown in Figure 10e. This subclass shows a strong tendency to occur closer to the rift features proposed to be the most recent with, of the 17 very large volcanoes, 11 (around two thirds) occur- ring within 200 km, and 13 (around three quarters) occurring within 350 km, of a "young" rift segment. It is a small data set, however.
Due to the small sample size, it is hard to derive any spatial relationships in the anemone data set (Figures 5cv,  6cv, 7cv, and 10f ), although strong clustering is evident in one particular location (south Atla Regio). The more widely distributed nature of the other anemones limits the NNI and distribution analysis, which is simply included for completeness.

Timing of Volcano-Tectonic Events
Volcano-tectonic associations have been previously used to test models of Venus' geological evolution over time [e.g., Solomon et al., 1992;Stofan et al., 1992Stofan et al., , 1995Basilevsky and Head, 1998;Guest and Stofan, 1999]. In section 3 the associations of the new database of volcanic features were plotted with both the "old" rifts AIREY ET AL. VOLCANISM AND GLOBAL TECTONICS ON VENUS and "young" rifts mapped by Krassilnikov et al. [2012]. These "old" and "young" rifts have been suggested as representative of two distinct geological time periods on Venus with the division marked by the proposed globally synchronous stratigraphic marker of the regional plains unit (combining upper and lower regional plains units, Figure 11). These ideas are fundamental to the directional model of Venus geology [Basilevsky and Head, 1998, 2000a, 2002Head, 2013, 2015]. Guest and Stofan [1999], however, argue for a nondirectional model, with processes such as plains emplacement and rifting not occurring in a fixed sequence, and instead largely randomly and interleaved over shorter periods. Another argument for a more variable resurfacing history includes the calculated mantle potential temperatures for the Venera site compositions [Shellnutt, 2016]. These are suggestive of an ambient mantle thermal regime occurring concurrent with high thermal regime regions such as Beta Regio, which is inconsistent with either end-member (catastrophic versus equilibrium) resurfacing model. In this section the implications of the current findings in terms of the timing and evolution of Venus' volcanism, tectonism, and geology are explored. Figure 11. Detail of the study area taken from the global geological map of Venus [Ivanov and Head, 2011], redrawn to emphasize the BAT study region. Units most relevant to the interpretations in this study are, in the left hand column, the interpreted postplains rift zones (dark purple) and the regional plains themselves (pale blue and midblue) and in the center column, the preplains shield fields (teal) and groove belts (analogous with "old" rifts (pale pink)).
First of all, the coronoids, with the exception of type 2 coronae, all show a strong association with rifts (Figures 6bi-6biii), especially at Parga and Hecate Chasmata, which is consistent with the subcrustal lid rejuvenation model [Ghail, 2015]. However, this association is less significant when only the "young" rifts are considered (Figure 7bi-7biii), consistent with corona formation events being broadly coincident with phases of older rifting. Although the corona formation events may not necessarily be directly genetically linked with the proposed initial rifting stage, if they are correlated in both time and space, it is likely that there is a common underlying cause of their formation. Although it remains possible, given the uncertainties inherent in the observations in this study, that coronae are associated with both "young" and "old" rifts and hence that they continued to form to significant extents throughout both rifting episodes, the suggestion here that they are Figure 12. Bar chart of coronae data from Martin et al. [2007]. Data included for coronae identified in the cited article as on rift and further classified based on formation relationship with any associated rift as postrift, postrift/rift synchronous, rift synchronous, prerift/rift synchronous, or prerift. Number of on-rift coronae identified by Martin et al. [2007] that coincide with rift mapping of Krassilnikov et al. [2012] (and reproduced in this study) upon which this figure is based number 47 on "old" rifts and 22 on "young" rifts. more strongly associated with the proposed "old" rifts is corroborated by the study of Krassilnikov et al. [2012]. They examined a subset of coronae (20% of the global database) and concluded that as few as 3% of coronae formed after the regional plains [Krassilnikov et al., 2012].
To test the consistency of coronae distribution with alternative resurfacing models of Venus, rift-relative corona formation times at Parga Chasma specified in Martin et al. [2007] may be investigated in the context of our current study. Martin et al. [2007] identified the coronae as postrift, postrift/rift synchronous, rift synchronous, prerift/ rift synchronous, or prerift (note that these classifications are independent of the "young" and "old" rift definitions of Krassilnikov et al. [2012]). Where these coronae coincide with the rifts in our study, they were included in the analysis in Figure 12. Of coronae coincident with

10.1002/2016JE005205
"young" rifts, 70% are classified as prerift, 18% as rift synchronous, and only 12% as postrift, suggesting that if the directional model is borne out, a relatively small proportion of coronae were formed after the proposed later stage of rifting. Of those occurring on the "old" rifts, these categories become more balanced at 42%, 31%, and 27%, respectively. This is consistent with the concept of on-rift corona formation preceding and continuing throughout the rifting episodes; however, there is very little evidence of purely postrift corona activity at Parga Chasma suggesting either cessation of corona formation or ongoing rifting concurrent with synchronous corona activity. The distribution of coronae therefore appears consistent with predictions of each of the different resurfacing models and does not clearly distinguish between them.
The distribution of shield fields is informative, even though the population distributions for distance to either proposed rift stage ("old" or "young") is consistent with a random distribution (Figures 6ciii and 7ciii). It is evident on inspection of Figure 10c that the more prominent shield field clusters are most strongly associated with the terrain of the "old" rifting stage (regions of preregional plains rifting, black lines), largely coincident with the shield plains unit (Figure 11). Shield fields in areas coincident with either the regional plains or rift zone units are less common, but not entirely absent (Figures 10c and 11), suggesting that the bulk of shield field-forming volcanism occurred prior to the proposed plains resurfacing in these areas.
The domes and calderas show strong clustering (Figures 5aiii and 5ci), with NNIs of 0.608 and 0.717, respectively ( Figure 3). Intermediate volcanoes show only moderate clustering (Figure 5cii), with an NNI of 0.915 (Figure 3), which could be due to the weighting on the lower value bins of the histogram from the strongly clustered part of the population being diminished by the secondary peak at ∼700 km resulting in a bimodal distribution (grey bars, Figure 5cii). The distributions with respect to the rifts of these three classes of dispersed, small to intermediately sized volcano-tectonic features do not differ greatly from random, providing no evidence for a genetic association. However, they all correlate broadly with the shield plains and groove belts regions ( Figure 11).
Finally, the large volcanoes (>100 km diameter) show a population distribution with respect to "all rifts" that is broadly consistent with the random population; there is, however, a strong spike in the leftmost bin corresponding to ∼35% of large volcanoes occurring within ∼50 km of a rift (Figure 6civ). When only "young" rifts are considered (Figure 7civ) this distribution remains strong, suggesting that it is largely this subpopulation responsible for causing the population distribution to diverge from random; otherwise, the difference from random would decrease commensurate with the reduction in the population. This is further evidence of a link between rift zones interpreted as "young" and persistent large volcano forming centers of volcanism as suggested previously [e.g., Basilevsky and Head, 2000b] and hence a clear difference in the volcanic characteristics of the rifts designated as "young" by Krassilnikov et al. [2012], which supports their classification as in some way distinct from the "old" rifts. Further evidence for enhanced volcanism in this proposed second stage of rifting is evident in postemplacement volcanism at old coronae, which most commonly occurs coincident with the shield plains units and less commonly into the postregional plains period [Krassilnikov et al., 2012], as well as the recent Venus Express data suggesting direct evidence for active volcanism at Ganis Chasma [Shalygin et al., 2015].

Implications for the Underlying Tectonics
As discussed in section 4.1, the new evidence in this study implies a period of corona formation and widespread small-scale volcanism [Ivanov and Head, 2013]. The underlying causes of these volcanic features are thought to be the interaction of subsurface mantle plumes or diapirs for the former Janes et al., 1992;Squyres et al., 1992;Stofan et al., 1992] and magma migration to the surface for the latter [Head and Wilson, 1986;Head et al., 1992;Crumpler and Aubele, 2000]. The dispersed nature of this initial phase of Venus volcanism could result from many widely dispersed small plumes impinging on the base of a globally relatively thin lithosphere in an earlier "mobile lid" regime [Moresi and Solomatov, 1998] or the initiation of subcrustal lid extension [Ghail, 2015], in either case causing buoyant uplift of the lithosphere. The crust then responds by gravitational collapse/flow, allowing magma to propagate up into the rifts. As the lithosphere thickens, the buoyant uplift will be increasingly associated with fewer, larger, plumes [Phillips and Hansen, 1994;Turcotte, 1995;Brown and Grimm, 1999]. As long as there is a flux of magma into the crust (from partial melting in the plume), then there is a mechanism for continued magma-assisted rifting, by analogy with that which is proposed to occur at the Main Ethiopian Rift on Earth [Kendall et al., 2005;Corti, 2009], with coronae and volcanism associated with rifts at the surface. , a detachment layer at the base of the crust and a CO 2 -induced "petrological" asthenosphere facilitate spreading and recycling of the mantle "lid." Regions of upwelling result in lid rejuvenation and rift formation, whereas regions of downwelling may be responsible for stress-related features such as the wrinkle ridges found to affect the regional plains occurring at localized areas of compression.
Can these different models of global resurfacing and rifting be reconciled? What processes would result in the two stages of surface activity supported by this study and the synthesis of previous work [i.e., Ivanov and Head, 2015] described in section 4.1? Ghail [2015] proposed a mechanism by which the mantle part of the stagnant lid might be able to detach from the overlying crust and be spread and recycled into the deeper mantle. This model would result in the transition from that illustrated in Figure 13a to that in Figure 13b. Figure 13a represents Venus' crust and upper mantle at a time before the episode of proposed "young" rifting, where globally dispersed centers of magma upwelling generate a correspondingly globally dispersed incidence of volcanism and minor rifts as observed in this study associated with the preregional plains rifts (c.f. Figures 10c and 11). Conductive cooling and gradual thickening of the lithosphere over time [Phillips and Hansen, 1994;Turcotte, 1995;Brown and Grimm, 1999] results in the transition to Figure 13b, where mantle upwelling is concentrated and localized at rift axes and large-scale volcanism is then associated with larger rift zones. This is illustrated in this study by the evidence provided in Figures 7civ, 10d, and 10e, and described in section 4.1.
In the absence of any major subduction-related plate driving forces, it is proposed that the predominant mechanism by which major rift formation progresses is analogous to magma-assisted rifting [Kendall et al., 2005;Corti, 2009]. Smrekar et al. [2010] do note, however, that the apparent depths of compensation beneath these chasmata do not imply large-scale upwelling, suggesting that the exact rifting mechanism may be characterized by a more complex network of enhanced magmatism. Overall, the processes described in this section may provide a satisfactory tectonic evolution model for Venus pending further study.

Synthesis of Evidence for Volcano-Tectonic Processes and Timing
The evidence presented in sections 4.1 and 4.2 may be considered broadly consistent with some observations cited in the directional model of rift evolution [Basilevsky and Head, 2000a]. Where this is the case, the findings are summarized in Figure 14, which describes the tectonic and volcanic evolution of Venus proposed by the directional model. The "supporting evidence" row has been added to highlight the new evidence in this study.
Although there is a broad consistency and despite basing the rift classification used in this study upon directional concepts, there remains a degree of nondirectional activity [Guest and Stofan, 1999], such as the rare postplains coronae identified by Krassilnikov et al. [2012], and large volcanoes persisting away from "young" rifts (e.g., Ituana "Corona", a large volcano located at 19.5 ∘ N, 154 ∘ E, Figure 10e), suggesting that the ongoing Figure 14. Summary of the tectonic history of Venus since tessera formation (oldest observed stratigraphic unit) as outlined in accordance with the directional model of Venus evolution. Modified from Ivanov and Head [2015]. The model describes three main global volcano-tectonic regimes beginning with a tectonic regime characterized by intense tectonic deformation and the formation of minor rift features (groove belts corresponding to "old rift" features in this study). This is followed by a volcanic regime characterized by sequential phases of plains emplacement. Finally, the network rifting-volcanism regime dominates along with the formation of the youngest and most intense rifting episodes. Unit names in the "rock stratigraphic units and structures" row refer to regions mapped on Figure 11: tesserae (t), densely lineated plains (pdl), ridged plains/ridge belts (pr/RB), groove belts (gb), shield plains (psh), lower and upper regional plains (rp1 and rp2), lobate plains (pl), and rift zones (rz). Lettering in the "supporting evidence" row refers to evidence cited in Figures: (a) 2; (b) 6b, 7b, and 9c; (c) 3, 5a, 5c, 8c, 10a, and 10b; (d) 3, 5c, 10c, and 11; (e) 13; (f ) 7c, 10d, and 10e.
sequence of events is more complex than either end-member model. It should also be noted that with current observational data there is also a limit in terms of to what extent we can categorically test one model against another and there are certainly other scenarios that could give rise to similar distributions of volcanic features on Venus.
Although attempts can be made to infer the relative sequence of events to a degree by studying the relationships described in this work, the absolute timescales of these processes are largely speculative. In the context of Venus' proposed geological history set out in Ivanov and Head [2015], the relative timing of events may be attributed to some broad geological time periods. A reasonable estimate of Venus' mean surface age based on the cratering history is ∼750 Ma, after McKinnon et al. [1997], and the processes occurring during the "Global Tectonic Regime" (Figure 14) are likely to have been occurring from this time. The period of time (probably of order tens of million years) in which the stagnant lid scenario might have transitioned to the proposed subcrustal spreading mechanism happened at some stage since then. This suggests that both rift-forming episodes were prevalent for order hundreds of million years, with the late-stage "old" rift scenario coinciding broadly with the "Global Volcanic Regime" (Figure 14), and the "young" rift scenario coinciding broadly with the "Network Rifting-Volcanism Regime" (Figure 14) following the transition to the proposed spreading mechanism.

A Failed Rift?
The chain of concentrated volcanic features including domes, coronae, and shield fields coincident with a region of rift axes defined as "old" in southern Llorona Planitia located in the upper left of the study area at around 10-30 ∘ N and 90-150 ∘ E (Figures 8c, 9c, 10c, and 11) is a striking feature of the study region. As previously noted in sections 4.1-4.3, coronae are strongly associated with rifting, indicating that this region may have been in the early stages of transition into a "young" rift (postregional plains). However, following the proposed regional plains emplacement, there is no evidence of further rifting or development of large-scale volcanism in this region. This suggests that this rift section failed (as described in section 4.2 and Figure 13), with the associated large-scale magmatism probably being diverted toward the more concentrated active "young" rifts as described in section 4.1 and Figure 13. This further example of evidence in support of a transition from widespread volcanism to concentrated rifting is consistent with the response to thickening lithosphere (and mobile to stagnant lid) described in Phillips and Hansen [1998].
On Earth, continental rifting is a well-studied component of plate tectonics (Figure 1a). It is possible that this process may have characteristics analogous with the processes responsible for the rifts observed on Venus. The East African rift system (EARS) has often been cited as potentially analogous to Venusian rifting, particularly the region between Beta and Phoebe Regiones [McGill et al., 1981;Campbell et al., 1984;Stofan et al., 1989;Foster and Nimmo, 1996]. This analogy is based primarily on the similarity of geomorphological characteristics in these two settings, with fault bounded rift segments connecting domed regions of enhanced volcanism. The proposed failed rift on Venus may have been the result of a plume-lithosphere interaction similar to that responsible for the observed difference between the eastern (magma-rich, active) and western (magma-poor, perhaps failing) branches of the central EARS on either side of the Tanzanian craton (centered at 3 ∘ S, 32 ∘ E, Figure 15) [Roberts et al., 2012]. Koptev et al. [2015] modeled the possible geodynamic processes responsible for this EAR rift failure as a mantle plume head impacting on the base of a significantly thickened area of lithosphere (i.e., the Tanzanian craton). The plume is deflected by the craton in their model to form magma-rich (eastern high density of recent volcanism) and magma-poor (western low density of recent volcanism) branches of the rift on either side ( Figure 15). If indeed plume-rift formation processes and interactions are in any way analogous on Earth and Venus, it seems likely that a similar event may have been responsible for the proposed failed rift on Venus. Thetis Regio (Figure 2b), the easternmost area of the highland region Aphrodite Terra, is a region of thickened lithosphere (evidenced by its high topography) that occurs in between the proposed failed section and the system of "young" chasmata running along the southern fringe of Thetis Regio. It is plausible to suggest the mechanism proposed by Koptev et al. [2015] might apply here too, though the scale differs significantly; the distance between the rifts on either side of the thickened lithosphere on Venus is around five times that on Earth (∼3500 km and ∼700 km, respectively). If this is the case, deflection or bifurcation of a plume beneath a much wider region of thickened lithosphere and the resulting surface expression implies a larger plume on Venus in order to result in the differing surface expressions on either side. Further geodynamic modeling to test the plausibility of applying this model to this region of Venus is to be encouraged.

Conclusions
With the use of new geospatial analyses performed in ArcGIS, the quantification of spatial relationships between volcanic features and rifts on Venus was achieved and applied to global volcano-tectonic concepts. The systematic study of nearest neighbor distances, distances from rifts, and distances from "young" rifts (according to the directional evolutionary model and also locally in nondirectional terms) facilitated the comparison of spatial and temporal trends of volcanic behavior.
The nearest neighbor analysis provided evidence that all the dome-type and corona-type subpopulations display a tendency for clustering behavior, having NNI values significantly lower than 1 (0.6-0.9). This is despite the propensity for truly random populations with relatively small sample sizes such as these to appear somewhat more uniformly distributed due to the edge effects. With the exception of those for shield fields and large volcanoes, all the feature population distribution histograms were more positively skewed when compared with random distributions of the same sample size, providing supporting evidence for the clustering.
Rift associations were inferred via the deviation of a feature's population distribution (as a function of distance from rift) from that of a random population. Tests were performed to see whether there was a distinction between the "old" and "young" rifts defined by Krassilnikov et al. [2012] in terms of the spatial associations of the volcanic features. It was found that dome-type features in general had no discernible relationship with rifts, corona-type features had a strong association with "all rifts" (with the exception of type 2), less so with only "young" rifts, and with the volcano class there are clear associations with "all rifts". This last class has a more complex relationship with "young" rifts; the small volcanoes (those in shield fields) have a preference to occur away from rifts; intermediate volcanoes and, most notably, large volcanoes show a preference to occur close to or, most strongly in the case of the latter, on rifts. The study found that ∼35% of large volcanoes in the study region occur within 50 km of a rift segment, and when very large volcanoes are considered, approximately two thirds occur within 200 km and approximately three quarters occur within 350 km, of a young rift segment.
When these data analyses were compared with characteristics of the directional [Basilevsky and Head, 1998] and nondirectional [Guest and Stofan, 1999] evolutionary models for Venus, a number of observations were found to be consistent with components of the directional sequence. When applying the methods independently to the two inferred relative rift ages, the differing associations discovered were indicative of a broad agreement with rift definitions as described by a directional-type model. In such models, initial widely dispersed rifting structures with associated coronae are indicative of prerift processes and dispersed shield fields dominated before the emplacement of the "regional plains" unit. In previous studies advocating the directional model this unit has been proposed to be globally synchronous. Our new data do not bear on this issue and could be consistent with local or global directionality, or indeed some other reason for different classes of rift on Venus.
New evidence in this study shows a closer association of coronae with "all rifts" than with "young rifts". We suggest that this is consistent with an early tectonic regime related to a relatively thin lithosphere allowing the dispersed, relatively small, plumes impinging on the lithosphere to generate the widespread small-scale volcanism evidenced by the shield fields. As the lithosphere cools conductively and thickens, we propose that the plume activity becomes concentrated at fewer, larger plumes and is focused along the "young" rift zones. The evidence in this study supports our suggested model, showing a stronger correlation of large volcanoes with "young rifts" than with "all rifts".
Relationships confirmed by this work are also consistent with the evolution from a stagnant lid regime to the hypothesized subcrustal spreading regime of Ghail [2015]. The transition from widespread to more localized activity could result from the concentration of plumes in response to magma-assisted rifting, a thickening lithosphere, and crustal detachment of the subcrustal lid. These processes would promote subcrustal spreading and the enhanced recent volcanism concentrated along "young" rifts. Although other processes could