A Global Survey of Lithospheric Flexure at Steep-Sided Domical Volcanoes on Venus Reveals Intermediate Elastic Thicknesses

Topographic ﬂexure in response to vertical loads reveals key lithospheric properties, including elastic thickness and the heat ﬂow from the interior. Flexural stresses may also control volcano morphology. One previous study predicted that steep-sided domes on Venus usually form where the elastic thickness is ˜15-40 km. We surveyed ﬂexural signatures around steep-sided domes and conﬁrmed this hypothesis. We determined elastic thickness from topographic proﬁles with a curve-ﬁtting algorithm and a plate bending model in Cartesian and axisymmetric geometry. We used a yield stress envelope to convert elastic thickness and plate curvature into mechanical thickness and surface heat ﬂow. The average elastic thickness for domes not near coronae is ˜30 km, corresponding to a heat ﬂow of ˜60 mW/m 2 . Coronae on Venus are typically associated with elastic thicknesses of < 10-15 km. Domes near coronae yielded elastic thicknesses in this range, and higher heat ﬂows than domes not near coronae.


Introduction
Properties of the lithosphere set the upper boundary condition for the evolution of planetary interiors.Interior dynamics in turn govern surface conditions over geologic time via volcanism, tectonics, and atmospheric outgassing (e.g., Foley & Driscoll, 2016;Smrekar et al., 2018).Venus's lithosphere is currently poorly understood, but we can use surface observations to learn about the interior and evolution of Earth's "evil twin".In particular, the elastic thickness of the lithosphere might control the morphology of volcanoes.McGovern et al. (2013) used two different magma ascent criteria to explore how coronae, steep-sided domes, and large conical edifices would form.Steep-sided domes, also called "pancake domes," appear nearly circular and flat from above (e.g., Gleason et al., 2010).This peculiar morphology is usually attributed to siliceous (dacitic and/or rhyolitic) volcanism (e.g., Head et al., 1991;Pavri et al., 1992).However, other studies suggest that the radar properties of these domes are more similar to basaltic or andesitic flows (e.g., Ford, 1994;Plaut et al., 2004;Stofan et al., 2000).In any case, steep-sided domes may most likely form at regions with elastic thicknesses of ~10-40 km as shown in Figure 1 (McGovern et al., 2013).However, this hypothesis has not yet been confirmed via measurements of the actual elastic thickness.
Coronae are oval-shaped features with annular fractures and complex topography (e.g., Smrekar & Stofan, 1997;Stofan et al., 1992).According to McGovern et al. (2013), volcanic features resembling some coronae were most likely to form at regions with the lowest elastic thicknesses: ~15 km or less.Previous studies derived thin elastic thicknesses at coronae showing signs of lithospheric flexure.O'Rourke & Smrekar (2018) used newly available elevation models derived from stereo images (Herrick et al., 2012) to analyze topographic profiles of flexural signatures around coronae.That study used plate bending models of both broken and continuous (unbroken) elastic plates to estimate elastic thickness and heat flow.They determined the heat flow at all of these locations based on a conversion from elastic to mechanical thicknesses using an assumed rheology for the lithosphere.Overall, they found elastic thicknesses of ~5 to 15 km at most coronae, corresponding to heat flows >95 mW/m 2 .This result was somewhat unexpected because most models for the thermal evolution of Venus predict that the average surface heat flow is ~40 mW/m 2 today with spikes to ~60 mW/m 2 possible in the last billion years (e.g., Armann & Tackley, 2012;Gillmann & Tackley, 2014;Noack et al., 2012;O'Rourke & Korenaga, 2015;Weller & Kiefer, 2020).In contrast, the average surface heat flow for Earth is ~86 mW/m 2 now (Jaupart et al., 2007).
Large conical volcanoes are predicted to form mostly at regions with elastic thicknesses of 40 km or greater.We did not explore large conical volcanoes in this study because the topographic signatures of flexure necessary for our analysis are not visible at those locations.McGovern & Solomon (1997) explained that flexural moats around these large volcanoes were likely filled in by lava flows from the large volcanoes themselves.Even using the stereo-derived topography, we were unable to find flexural signatures that would permit the analyses used for coronae previously and for steep-sided domes in this study.Russell & Johnson (2019) previously studied lithospheric flexure at one steep-sided dome.Specifically, they obtained topographic profiles radiating from a flexural moat around Narina Tholi.Their pioneering analyses revealed a very thin lithosphere at this location, with an elastic thickness of ~2-3 km.This result was somewhat expected given that Narina Tholi is located directly on the rim of Aramaiti Corona, which have been found to be associated with <15 km-thick lithosphere (O'Rourke & Smrekar 2018;Johnson & Sandwell 1994).They concluded that Narina Tholi was emplaced during the late stages of the formation of Aramaiti Corona when the heat flow from the interior to the surface was particularly elevated.Recently, Gülcher et al. (2020) argued that plume-lithosphere interaction at Aramaiti Corona could be actively occurring today.Generally speaking, the geologic context of individual volcanic features is critical to interpreting measurements of elastic thickness from lithospheric flexure.In this study, we conducted a global survey of lithospheric flexure at steep-sided domes to assemble a database of elastic thicknesses and other properties for this class of features.

Topographic Profiles
We accessed a database of volcanic features on Venus (Pavri et al., 1992) to find the latitude and longitude of steep-sided domes.We used JMARS to locate these domes using Magellan Synthetic Aperture Radar (SAR) data.Only the domes we were able to definitively identify in both the SAR imagery and topography (75 out of ~150 candidates) were used in this study.The remaining candidates were disregarded if we were unable to locate them in the imagery, if there was a gap in the topography data, or if the topography data was not of a high enough quality to distinguish the dome.We drew 8 topographic profiles at increments of ~45° around 75 steep-sided domes.The profile extending north (0°) is labelled profile 1, and the labelling numbers increase clockwise through profile 8 (Figure 2a).Wherever possible, the topographic profiles were taken using stereo-derived topography, which covers about 20% of the surface at horizontal and vertical resolutions of ~1 km and ~100 m, respectively (Herrick et al., 2012).Topography from areas without stereo coverage was extracted from the global topographic data record (GTDR) collected from the Magellan mission.The horizontal resolution is more than an order-of-magnitude worse than the stereo dataset at ~10-20 km (Ford & Pettengill, 1992).The profiles were drawn starting at approximately the edge of the dome and extended out ~800 km.We imported the topographic profiles into MATLAB and trimmed them to begin at the lowest point in the moat around the dome.The elevation values were normalized to 0 m at the start by subtracting the first elevation value.
Only the profiles that showed clear evidence of flexural signatures were chosen for further analysis.These signatures appear as a moat surrounding the dome, followed by a rise in elevation leading to a forebulge, before the topography eventually flattens out (Figure 2b).Some profiles encountered interfering topography such as craters or other mountains and were either discarded or truncated.We found flexural signatures in 29 profiles from 14 different domes.These domes were chosen for further analysis (Figure 2c).We noted the type of terrain in which each dome was located: tesserae or plains.Additionally, we determined which domes were in close proximity to other features such as coronae.Fotla Corona lies south of dome 75, interacting with only profiles 3, 4, 5, and 6.We thus divide dome 75 into 75a (profiles 1, 2, and 8) and 75b (profiles 3-6) to indicate which profiles intersect the corona and which do not.

Flexure and Elastic Thickness
We used standard models of lithospheric flexure to infer the elastic thickness from topography.Steep-sided domes are ~1 km tall (e.g., Gleason et al., 2010;Pavri et al., 1992) and thus produce vertical loads of order ~1 MPa, which are not sufficient to break the lithosphere by themselves.In general, we assume that the lithosphere is continuous (unbroken) and that in-plane (horizontal) forces are negligible.A differential equation describes the resulting deformation under a vertical load (Turcotte & Schubert, 2002): where D is the flexural rigidity, w is the displacement from the pre-flexure elevation, q is the vertical load per unit area, and r is the radial distance along track.Table S1 defines the constants used in our analysis such as the density contrast across the lithosphere between the atmosphere and upper mantle (Dr) and gravitational acceleration (g).We can solve the governing equation in both Cartesian and axisymmetric (cylindrical) geometry.The axisymmetric version is preferred because the steep-sided domes are relatively small with diameters comparable to the elastic thickness.In other words, the horizontal curvature of the load is important to assessing the flexural response of the plate.In contrast, switching the assumed geometry for flexural analyses of much larger features such as coronae changes the derived elastic thicknesses by ~20% or less (e.g., Johnson & Sandwell, 1994;O'Rourke & Smrekar, 2018).
We fit analytic solutions to Equation 1 to the topographic profiles.In axisymmetric geometry, q(r) = q0 is a constant for r less than the diameter of the dome (R) and zero elsewhere.The solution for r ≥ R, outside of the dome, is (Johnson & Sandwell, 1994): Here a is the flexural parameter and s accounts for any regional slope.This equation uses the Kelvin functions ber, ker, bei, and kei with apostrophes indicating first-derivatives.The Cartesian model treats the dome as an end load so q(r) = 0 Pa everywhere along the profile.The standard Cartesian solution is (Turcotte & Schubert, 2002): where w0 is the maximum displacement (i.e., at the start of the profile in the flexural trough).We used the fminsearch curve-fitting algorithm in MATLAB (Lagarias et al., 1998) to find the best-fit parameters for each model.We seed the curve-fitting algorithm with initial guesses: the distance along track to the forebulge for a; the maximum normalized elevation divided by the total length of the profile for s; the difference between the maximum and minimum elevations for w0 in the Cartesian model, and 1 MPa for q0 in the axisymmetric model.
The elastic thickness is extracted from the best-fit flexural parameters derived in both geometries.Specifically, the elastic thickness (he) is (Turcotte & Schubert 2002): where E is Young's modulus (100 GPa) and n is Poisson's ratio (0.25).This equation is derived from the standard definitions of the flexural rigidity, D = (Ehe 3 )/[12(1 -n 2 )], and the flexural parameter, a = [4D/(Drg)] 1/4 .The curve-fitting algorithm returns formal uncertainties for all parameters.We use a Monte Carlo method to propagate the uncertainties on a to estimate the "1sigma" uncertainties on the best-fit elastic thickness.O'Rourke & Smrekar (2018) used a sophisticated Markov chain Monte Carlo to fit flexural models to topography, but also demonstrated that simpler curve-fitting algorithms yield equivalent results.

Mechanical Thickness and Heat Flow
We used a yield stress envelope to calculate the mechanical thickness, surface heat flow, and thermal gradient at each dome.In an idealized elastic plate, the horizontal stress varies linearly between the maximum extensional stress at the top of the plate to the maximum compressional stress at its lower boundary.The real, mechanical plate is thicker than the elastic lithosphere because brittle and ductile failure occur at the top and bottom of the plate, respectively.The total bending moments for the idealized elastic and mechanical plates should be equal, allowing us to convert elastic thickness estimates to mechanical thickness estimates (e.g., Brown & Grimm, 1996;McNutt, 1984).We used the maximum curvature (d 2 w/dr 2 ) of the elastic plate in these calculations.We determined that the axisymmetric geometry was the most accurate and chose to move forward analyzing the mechanical thicknesses for only this model.
Because we know that the mechanical plate begins at the surface and ends where the temperature is high enough to cause ductile flow, we can use the mechanical thickness to determine the thermal gradient.The equation used to obtain the thermal gradient is: TS = 740 K is the surface temperature, and hm is the mechanical thickness.TC is the critical temperature above which ductile flow begins, which depends on the rheology of the material.Table 1 reports the results calculated using a dry olivine rheology (TC ~ 1013 K).Table S2 lists values calculated using two rheological models for dry diabase (Mackwell et al., 1998).We then calculate the surface heat flow: with kc = 4 W/m/K being the thermal conductivity of the crust (Turcotte & Schubert, 2002).

Agreement with Gravity Data
We determined the level of agreement between our derived elastic thicknesses and those in the global map derived from admittance, which is the ratio of gravity to topography (Figure 14 in Anderson & Smrekar, 2006).By searching within an oval the size of the resolution of the gravity, we counted the number of pixels for which the values match those from our study.Domes where >50%, 0-50%, and 0% of the pixels agree within ±20 km are considered to have "excellent", "good", and "no" agreement, respectively.Two domes (21 and 51) are partially within areas where elastic thickness could not be retrieved from admittance data, which might artificially lower the level of agreement.

Table 1
List of all domes with at least one topographic profile amenable to a flexural interpretation.From left to right, the first four columns show the index number, diameter, longitude, and latitude of each dome.Elastic thicknesses (he) were calculated using plate-bending models with axisymmetric and Cartesian geometry.Mechanical thicknesses (he) and surface heat flows (FS) were derived using a yield stress envelope for dry olivine.Our estimates of elastic thickness were compared with the map in Anderson & Smrekar (2006)

Results
Table 1 collects the results for all profiles in this study.Out of the 14 domes that showed evidence of flexure, 6 were located within the swaths of stereo data.In order to not overrepresent domes for which multiple profiles showed flexural signatures, we calculated the average values and uncertainties for all profiles around each dome using a Monte Carlo technique.The values in Table 1 for domes 2, 5, 7, 75a, and 75b are the averages of multiple profiles around those domes that were amenable to flexural interpretations.Table S3 lists the parameters derived from each individual profile.When considering all profiles, we found an average elastic thickness of ~16 and 15 km for a continuous plate using the Cartesian and axisymmetric models, respectively.
There was a significant difference in derived elastic thicknesses based on whether or not the domes were located in close proximity to coronae, which have been found mostly in regions with thin elastic thicknesses (Figure 3a).Profiles around domes 2, 7, 51, 64, and 75b were near coronae (within ~800 km, the length of the topographic profiles).Most of these profiles did have elastic thicknesses below the average expected for steep-sided domes.The average elastic thickness for profiles that are near coronae ranged from ~2-16 km using Cartesian geometry, and ~3-14 km using axisymmetric geometry.Profiles around domes 1, 5, 13,16, 20, 21, 22, 26, 32, and 75a were not located near coronae.These profiles yielded average elastic plate thicknesses ranging from ~8-41 km using the Cartesian model and ~9-38 km using the axisymmetric model.Our derived mechanical thicknesses were ~20% larger than the elastic thicknesses, while mechanical thicknesses can be a factor of ~2 larger than elastic thicknesses derived at coronae (O'Rourke & Smrekar 2018).This difference is due to the fact that domes are associated with smaller loads and thicker lithosphere than coronae, thus leading to a smaller plate curvature.
We reproduced the results around Narina Tholi from Russell & Johnson (2019).This dome is directly on the margin of Aramaiti Corona and thus has a low elastic thickness.In our analysis, Narina Tholi is referred to as dome 2. We trimmed our profiles to extend ~100 km from the flexural moat to match that of Russel & Johnson (2019).These topographic profiles also show signs of lithospheric flexure.The elastic thickness at the location of dome 2 is low: ~3 km using Cartesian geometry and ~4 km using axisymmetric geometry.This low estimate of elastic thickness for the dome agrees with the values of ~2-3 km from Russell & Johnson (2019).
We also found a significant variance in surface heat flow dependent on whether the profiles were near coronae (Figure 3b).Those profiles that were close to coronae yielded surface heat flows ranging from ~104-279 mW/m 2 , (~197 mW/m 2 on average).Analyses of flexural signatures from topographic profiles that did not encounter coronae yielded heat flows of ~29-109 mW/m 2 .The average heat flow of ~57 mW/m 2 for these domes is close to that predicted for Venus globally in some thermal evolution models (e.g., Gillmann & Tackley, 2014).
We determined the level of agreement between our derived values and those derived from gravity data.Out of 14 domes, 12 showed at least some agreement between our elastic thickness values and the gravity-derived values.Of these, 7 showed excellent agreement (>50%), another 5 domes showed good agreement (0-50%), and 2 showed no agreement.There was no bias in the levels of agreement between domes that are near coronae or not.

Discussion & Conclusions
We verified the hypothesis that steep-sided domes typically form at regions of elastic thickness between ~10 and 40 km.The distribution of steep-sided domes can thus tell us about the elastic thickness of the lithosphere on Venus (McGovern et al., 2013).Flexural signatures are not always visible as moats may have been filled in by later lava flows or other interfering events.Additionally, it might not be possible to collect reliable topographic profiles around domes that are close to other interfering topography such as craters or mountains.However, one may use the presence of steep-sided domes to infer that the elastic thickness should be between ~10 and 40 km at a given location.Lithospheric thickness can then act as a probe into the interior of the planet.Thus, the distribution of steep-sided domes on Venus is also a tool that future studies could use to constrain the regional and global heat flow out of the deep interior.
Domes located near coronae were associated with thin elastic lithosphere.Our results thus provide further evidence for thin elastic thickness at coronae.Our findings for elastic thickness were consistent with the regional elastic thicknesses derived from gravity data (Anderson & Smrekar, 2006) regardless of the proximity of domes to coronae.Profiles that are located near coronae had significantly lower mechanical thickness and thus higher average surface heat flows than those that were not near coronae.
While volcanic features are likely to form where the elastic thickness is within the hypothesized range, they appear to can form on thinner lithosphere as well.We did not consider factors such as tectonic history, geologic history, and rheology that could affect volcano morphology (e.g.Watts & Zhong, 2000).Where domes occur near coronae, it is generally impossible to determine relative ages, making the role of any localized heating associated with corona formation difficult to assess.
Further spacecraft exploration of Venus is urgently needed to improve our understanding of the Venusian interior.Stereo-derived topography proved to be an important dataset in the search for flexural signatures.Stereo data is only available for about 20% of Venus' surface, yet >40% of the flexural signatures around steep-sided domes were found using stereo data.The disproportionate number of signatures identified using stereo-derived topography suggests that flexural signatures do exist around other steep-sided domes but could not be seen in the lowerresolution Global Topographic Data Record.This study also supported the idea that flexural signatures around steep-sided domes provide reliable estimates of lithospheric thickness.Improved data from new missions would reveal different features and allow us to conduct more precise studies.Future missions to Venus to gather high-resolution radar images, topography, and gravity data on a global scale are necessary to understand the interior and evolution of Earth's nearest neighbor.

Introduction
This Supporting Information includes tables that complement the main text, along with datasets that could be used to reproduce our results.All topographic profiles were extracted from the Magellan global topographic data record (GTDR) and, where available, stereo-derived digital elevation maps using JMARS.

Table S1
Values of constants used in curve-fitting functions and calculations, which were taken from Johnson & Sandwell (1994) except the rheological laws for dry olivine and diabase are from McNutt (1984) and Mackwell et al. (1998)
Data Set S2.Raw data (elevation versus distance along track) that was extracted from each topographic profile that was amenable to a flexural interpretation.These data were fit to the equations listed in the main text using the curve-fitting algorithm that is included with Matlab.

Figure S1.
Our derived estimates for elastic thickness generally agree with those derived from admittance data (Anderson & Smrekar, 2006).Ellipses representing domes for which we derived elastic thickness overlay a map of global topography.The ellipses are color and pattern coded based on the agreement between our derived elastic thicknesses and those from a global map constructed from admittance (the ratio of gravity to topography).The size of each ellipse represents the resolution of gravity data at that point.Green triangles, solid yellow, and red striped ovals indicate domes with excellent (>50% within ±20 km), good (>0-50%), and no (0%) agreement, respectively.

Figure 1 .
Figure 1.Steep-sided domes are predicted to form where the elastic thickness is thicker than at coronae but thinner than at large conical volcanoes (McGovern et al, 2013).These images were taken in JMARS using SAR imagery from NASA Magellan illuminated from the left.Black regions indicate missing data.The corona featured is Aramaiti Corona (82°E, 25.5°S), located next to Dome 2 in this paper.The dome featured is Dome 1 in this paper (151°E, 3°S).The large conical volcano featured is Tsityostinako Mons (215°E, 15.5°S).

Figure 2 .
Figure 2. Profiles that show signs of flexure are chosen for further analysis with the curve-fitting algorithm.(a) Profiles drawn around dome 16 (66 °E, 37.5 °N) in JMARS.The profile that shows signs of flexure (profile 6) is highlighted in red, which is the only profile that does not cover heavily fractured or tectonized terrain.(b) All topographic profiles, including those that are not amenable to a flexural interpretation (grey).The profile that was later analyzed is again in red.The profiles begin at the lowest elevation point in the flexural moat surrounding the dome.(c) Cartesian and axisymmetric curve fits to the shape of profile 6 around dome 16.

Figure 3 :
Figure 3: Derived properties at all domes analyzed in our study.Error bars indicate the average standard deviation ("1-sigma") calculated from all profiles around the dome.Domes near coronae are represented as pink triangles, while blue circles show domes that are not near coronae.(a) Elastic thicknesses at each dome calculated using the axisymmetric model.Color shading indicates the predicted ranges for domes (blue) and coronae (grey) from Figure 1.(b) Surface heat flows calculated using elastic thicknesses and curvatures from the axisymmetric model and a yield stress envelope for dry olivine.
using the scale explained in the main text.Quoted uncertainties represent one standard deviation.Results from multiple profiles were averaged using Monte Carlo error propagation.