ABSTRACT Title of Thesis: DEFORMATION AND FAILURE OF LUNAR LAVA TUBES Edward A. Williams, Master of Science, 2021 Thesis Directed By: Professor Laurent G. J. Montesi, Department of Geology Lava tubes are long void spaces left by lava flows. Terrestrial tubes have been extensively studied. Evidence indicates lunar tubes also exist, potentially ideal environments for lunar exploration, habitation, and study. This thesis presents models of elastic deformation around lunar lava tubes to determine the expected stresses, failure regions, and surface features associated with tubes and their dimensions. Failure on internal surfaces is extensive in the modeled tubes, leaving only small usable floor regions (maximum 500 m) near tube edges. Cracks and debris may be problematic for utilizing tubes. Shape influences failure extents; narrower, semicircular-floored tubes suffer less floor failure. Linear surface cracks and bulges are expected, at distances from the tube that depend on tube dimensions, and might be used to locate tubes and determine their dimensions without entering. Such cracking seems to be present near the Southern end of Rima Mairan, suggesting the rille continues into a tube. DEFORMATION AND FAILURE OF LUNAR LAVA TUBES by Edward A. Williams Thesis submitted to the Faculty of the Graduate School of the University of Maryland, College Park, in partial fulfillment of the requirements for the degree of Master of Science 2021 Advisory Committee: Professor Laurent G. J. Montesi, Chair Professor Wenlu Zhu Professor Nicholas C. Schmerr Professor Megan E. Newcombe ? Copyright by Edward A. Williams 2021 Acknowledgements I would like to thank my advisor, Dr. Laurent Montesi, for his guidance and support. I would also like to thank Ernie Bell for his comments and insights on terrestrial tubes and the failure observed inside them. This work was supported by the NASA SSERVI GEODES grant 80NSSC19M0216. ii Table of Contents Acknowledgements ....................................................................................................... ii Table of Contents ......................................................................................................... iii List of Tables ............................................................................................................... iv List of Figures ............................................................................................................... v List of Abbreviations ................................................................................................. xiv 1. Introduction ........................................................................................................... 1 1.1. Lava Tubes on Earth ..................................................................................... 1 1.2. Lunar Observations and Evidence of Lunar Lava Tubes.............................. 7 1.2.1. Sinuous Rilles ....................................................................................... 7 1.2.2. Holes ................................................................................................... 12 1.2.3. Gravity ................................................................................................ 16 1.2.4. Ground Penetrating Radar................................................................... 18 1.3. Uses of Lunar Lava Tubes .......................................................................... 20 1.4. Prior Modeling Work .................................................................................. 24 1.4.1. General Comparisons .......................................................................... 24 1.4.2. Blair et al. (2017) ................................................................................ 25 1.4.3. Theinat et al. (2020) ............................................................................ 27 1.4.4. Modiriasari et al. (2018) ..................................................................... 29 2. Methods............................................................................................................... 32 2.1. Modeling Process ........................................................................................ 32 2.1.1. Main Model Suites .............................................................................. 32 2.1.2. Model Suite Variants .......................................................................... 37 2.2. Model Analysis ........................................................................................... 40 2.3. Lunar Surface Observations ........................................................................ 49 3. Results ................................................................................................................. 52 3.1. Regions of Failure on Tube Internal Surfaces ............................................ 52 3.1.1. Failure Regions ................................................................................... 52 3.1.2. Usable Floor Extents ........................................................................... 56 3.2. Surface Features .......................................................................................... 59 3.3. Effects of Tube Shapes ............................................................................... 64 3.4. Lunar Observation Results .......................................................................... 72 4. Discussion ........................................................................................................... 79 4.1. Failure on Tube Internal Surfaces and Usable Floor Extents ..................... 79 4.2. Surface Features .......................................................................................... 82 4.3. Effects of Different Tube Shapes ................................................................ 85 4.4. Lunar Surface Observations ........................................................................ 86 4.5. Testing Robustness of Conclusions with Different Model Setups ............. 89 4.5.1. Block Size ........................................................................................... 89 4.5.2. Mesh .................................................................................................... 91 4.5.3. Grid ..................................................................................................... 94 5. Conclusions ......................................................................................................... 97 Bibliography ............................................................................................................... 98 This Table of Contents is automatically generated by MS Word, linked to the Heading formats used within the Chapter text. iii List of Tables Table 1: Width and roof thickness values used for the modeled tubes. Each combination of dimensions was used. Table 2: Variables and constant parameters used in tensile and compressive failure criteria. iv List of Figures Figure 1: The interior of N?huku, a lava tube in Hawaii that is open for visitors to walk through. Image credit: NPS Photo/Janice Wei (J. Wei). National Park Service photograph, available at https://www.nps.gov/havo/learn/news/20180202_nakuku_lighting.htm (accessed May 5 2021). Figure 2: Chain of holes into a lava tube, Lava Beds National Monument, CA. Coordinates at the center of the image are 41.75?N, 121.50?W. Satellite imagery obtained from Google Maps Figure 3: Cross sections of ??inahou Ranch lava tube at increasing distances from the vent (from a to d) over 2.62 km of the tube, showing where the flowing lava has downcut through a paleosol layer into a pre-existing flow (Figure 8 from Greeley et al., 1998) Figure 4: Section of Cassone lava tube roof where lava layers plastered onto the roof have collapsed, revealing linings with mm-thick layers. The ruler shown is 25 cm long. The arrows indicate ?pull-apart? stalactites caused by the detachment of a lava layer (Figure 7 from Calvari and Pinkerton, 1999). Figure 5: Sinuous rilles in Northern (a) and Southern (b) Oceanus Procellarum, mapped by Hurwitz et al. (2013). The background maps are Lunar Reconnaissance Orbiter Camera (LROC) Wide Angle Camera (WAC) global mosaics; topography contours on these are displayed in increments of 500 meters. Features identified as sinuous rilles are shown in black. Features identified as collapsed lava tubes are v indicated with arrows and labels. The authors note that the examples shown are not a comprehensive survey of collapsed lava tubes (Figure 26 from Hurwitz et al., 2013). Figure 6: The Rima Marius lunar sinuous rille with a discontinuity labeled in the figure that may be a lava tube segment. Image credit: Lunar Reconnaissance Orbiter Camera (LROC) Narrow Angle Camera (NAC), Image M135507533R Figure 7: Examples of terrestrial (a.) and lunar (b.) holes connected to lava tubes. The example in a. is found in the Corona volcano lava tube system (centered on 29.166? N, 13.457? W) (Satellite image from Google Earth). The example in b. is the Marius Hills Hole (MHH) in the Marius Hills region of the Moon (Modified from Figure 8 from Robinson et al., 2012. LROC NAC image M155607349R). Figure 8: Lunar Reconnaissance Orbiter Camera (LROC) Narrow Angle Camera (NAC) images of the Marius Hills Hole (MHH). (Modified from Figure 8 from Robinson et al., 2012). A-D show images of the hole with various conditions. A: 0.5? emission angle, 25? incidence angle (LROC NAC Image M122584310L). B: 9? emission angle, 13? incidence angle (LROC NAC Image M155607349R). C: 45? emission angle, 34? incidence angle (LROC NAC Image M137929856R). D: 30? emission angle, 13? incidence angle (LROC NAC Image M155614137R). Scale for D not given in Robinson et al., 2012. Figure 9: Map of the correlation between the free-air and Bouguer gravity anomalies for the region around the MHH, using a cross-correlation technique (modified from Figure 2 in Chappaz et al., 2016). Cold colors designate mass surplus; hot colors designate mass deficits. This combination of anomalies highlights mass deficits buried beneath the surface. The MHH is shown as a red circle; the long buried mass vi deficit near the MHH is the candidate lava tube (Chappaz et al., 2016). See details in Chappaz et al., 2016 for explanations of data analysis and units. Figure 10: Locations where LRS data showed signatures indicative of void spaces, overlaid on a Bouguer gravity map, made using a cross-correlation method to highlight where the Bouguer anomaly appears consistent with the expected anomaly for a long linear mass deficit such as a tube (gravity map used by Kaku et al., 2017 was developed by Chappaz et al., 2017). Cold colors on the gravity map designate mass surplus; hot colors designate mass deficits. The colors of the circular points show ????, the difference in echo power between the first and second echo peaks; void spaces are most likely at points with low ???? values (black, purple, or red). Most strong void space echo features are present on the Bouguer mass deficit, indicating a tube; the other strong echo features may indicate individual voids surrounded by mass surpluses. The MHH is next to point T1. Surface features, including the rille, are also shown. (Figure 5 from Kaku et al., 2017) Figure 11: Schematic of one possible design of a lunar base built inside a 100 meter scale lunar lava tube (Figure 3 from Daga et al., 1992). Figure 12: Blair et al. (2017) stability results for tubes of different widths and roof thicknesses in lithostatic (upper panel) and Poisson (lower panel) stress states (Figure 3 from Blair et al., 2017). The bold line in the Poisson stress state data divides two types of failure through the tube roof: compressive below the line, and tensile above it (Blair et al., 2017). Figure 13: Stability results for tubes of different widths and roof thicknesses (Figure 13 from Theinat et al., 2020). vii Figure 14: Modiriasari et al. (2018) results for tubes of different roof thicknesses and aspect ratios. Overburden yield (%) is the percentage of the thickness of each tube roof undergoing failure above the tube center. Tubes are considered stable (0% overburden yield), quasi-stable (some failure, but less than 50%), or unstable (50% or more) (Figure 1 from Modiriasari et al., 2018). Figure 15: (a.) Diagram of tube model geometry, showing the roof thickness, the tube half-width, and the direction of gravity. The model is symmetric across the dotted line. Only a portion of the model near the tube is shown in a.; b. shows the entire block for this example. c. shows the mesh used over the entire block; a significantly finer mesh is used near the tube itself than for other regions of the block Figure 16: Representative geometries of the tube shapes considered. This selection includes tubes with all three top-half width-to-height ratios and a subset of the floor curvatures used, including the deepest-dipping and the flat floors. Figure 17: Example unprocessed results from one COMSOL? model. The plotted stress is labeled as the First Principal Stress, and is positive; due to COMSOL?s? sign and stress naming conventions, this is the tensile stress, considered in this work as the Third Principal Stress and negative. Appropriate sign conversions are applied when reading data from COMSOL?. Stresses are only plotted for the upper region of the modeled block, to avoid plotting large-magnitude stresses far from the tube (such as at the bottom of the block, caused by the entire block?s weight) and making the colorbar less appropriate for the stresses close to the tube. Figure 18: Example of the grid used to interpolate data from a model. The grid in a. is the original, uniform grid. The grid in b. has been reshaped to fit the outline of the viii tube (green); reshaped grids like this are used for interpolation. The plots are zoomed in on the region of the block containing the tube. Figure 19: Model output for one modeled lunar tube, with a width of 6.5 km and a roof thickness of 650 m. The filled colored contours show the Hoek-Brown ratio described in Equation 5; compressive failure is present where this ratio is at least one, shown in white on the color scale. Shaded red designates where the rock is in tension. Circles plotted on the tube outline show tensile (red) and compressive (blue) failure. Figure 20: Region of a sinuous rille including the MHH viewed using LROC Quickmap with different lighting conditions. Location: 14.1?N, 56.8?W. 12.04 m/pixel. a.: Large solar incidence angle, East illumination only. b.: Large solar incidence angle, West illumination only. c.: Small solar incidence angle (0 ? 40 degrees). d.: Medium solar incidence angle (40 ? 55 degrees). e.: Large solar incidence angle (55 ? 80 degrees). Images were from the LROC NAC. Figure 21: Fractional extent of each tube floor that undergoes tensile failure vs. tube width. Colors and marker sizes denote different roof thicknesses used. The inset shows a simplified visualization of where this failure takes place. Figure 22: Fractional extent of each tube roof that undergoes tensile failure vs. tube width. Colors and marker sizes denote different roof thicknesses used. The inset shows a simplified visualization of where this failure takes place. Figure 23: Usable floor extents vs. tube widths for lunar (a.) and terrestrial (b.) tubes. Colors and line widths denote different roof thicknesses used. The inset shows a simplified visualization of where the failure types take place and where usable floor (green) is found. The green line represents the center of the tube. ix Figure 24: Sketch showing failure regions and distances to surface features. Figure 25: Horizontal distance from the center of the tube to the maximum tensile stress point on the surface, D_s, vs. tube width. Colors and marker sizes denote different roof thicknesses used. The inset shows a simplified visualization of the tube geometry and this distance on the surface. Figure 26: Horizontal distance from the center of the tube to the peak of the elastic bulge on the surface, D_b, vs. tube width. Colors and marker sizes denote different roof thicknesses used. The inset shows a simplified visualization of the tube geometry and this distance on the surface. Figure 27: Ratio of D_b over D_s vs. tube width (km) for modeled lunar (a.) and terrestrial (b.) tubes. The average ratio value is 1.73 ? 0.10 for lunar and 1.80 ? 0.14 for terrestrial tubes. Figure 28: Comparison of results of tubes with the same top-half width-to-height ratio (6:2; original half-ellipse) but different floor curvatures: (a.) 33.3%, (b.) 66.6% (symmetric), and (c.) 100% (semicircle). The models with more curving floors have less tensile floor failure; the 100% floor curvature tube has no floor tensile failure at all. Tube width: 5.5 km; tube roof thickness: 950 m. The filled colored contours show the Hoek-Brown ratio described in Equation 5; compressive failure is present where this ratio is at least one, shown in white on the color scale. Shaded red designates where the rock is in tension. The tube shape is outlined with a black line. Circles plotted on the tube outline show tensile (red) and compressive (blue) failure. Horizontal axis: horizontal distance (km); vertical axis: depth (km). x Figure 29: Comparison of results of tubes with the same floor curvature (66.6%) but different top-half width-to-height ratios: (a.) 6:1 (Flatter half-ellipse), (b.) 6:2 (Original half-ellipse), and (c.) 6:3 (Semicircle). Both models with a half-elliptical top half have tensile floor failure, although the 6:2 tube has less failure than the 6:1 tube; in the 6:3 tube, floor failure is prevented completely. Figure 30: Model output for a modeled lunar tube, with a width of 5.5 km and a roof thickness of 950 m. The width-to-height ratio of the top half of the tube is 6:1. The degree of floor curvature is 100%. The filled colored contours show the Hoek-Brown ratio described in Equation 5; compressive failure is present where this ratio is at least one, shown in white on the color scale. Shaded red designates where the rock is in tension. The tube shape is outlined with a black line. Circles plotted on the tube outline show tensile (red) and compressive (blue) failure. Figure 31: Model output for a modeled lunar tube, with a width of 5.5 km and a roof thickness of 950 m. The width-to-height ratio of the top half of the tube is 2:1 (half- circular). The degree of floor curvature is 33.3%. The filled colored contours show the Hoek-Brown ratio described in Equation 5; compressive failure is present where this ratio is at least one, shown in white on the color scale. Shaded red designates where the rock is in tension. The tube shape is outlined with a black line. Circles plotted on the tube outline show tensile (red) and compressive (blue) failure. Figure 32: Tube shape effects on maximum floor tensile stress. The x-axis shows floor curvature depth. The dashed black line shows the tensile strength of the material; any points below this line have failure. Colors designate different top-half aspect ratios used: Half-Circular (6:3 top-half width-to-height ratio) in red, Half- xi Elliptical (6:2 ratio; tube shape used for original models) in blue, and Flatter Half- Elliptical (6:1 ratio) in green. ?W? and ?RT? designate tube Width and Roof Thickness respectively. Diagrams show color-coded visualizations of some of the tube shapes used. Figure 33: Examples of features noted on or near lunar sinuous rilles. For each figure two versions are shown; one with a graticule and one without to avoid obscuring features. All panels are from regions of Rima Sharp. a. and b. show a region where apparent linear features follow part of the rille path on one side. c. and d. show possible linear cracking patterns inside the rille trough and on the rille edge. Imagery from LROC. Figure 34: Southern region of Rima Mairan. The rille ends in two aligned depressions following the path of the rille, connected by two lines of linear cracks. These could be indications of a partially collapsed lava tube continuing the path of the rille. Imagery from LROC. Figure 35: Imagery and data for the region of the lunar surface including the MHH and part of the surrounding rille. Figure a. is the LROC imagery of the area, including graticule and scale; the boundaries of the region shown and the scale are the same in b., c., d., and e. The insets in a. and c. are zoomed in on the MHH itself; both insets cover the same region with the same scale. The outline of the rille is shown in dashed yellow in a., and the MHH is circled in solid yellow; the same outlines are shown in the other figures, in either black or white depending on the image. b. shows SLDEM2015 slope data. Clear patterns follow the rille path, especially along the western section, and the MHH itself, but no clear indications of linear cracking (such xii as cracking causing broken blocks to change position and affect the slope) are visible to either side of the rille. c. shows SELENE Multiband Imager FeO abundance data. Patterns in the data can be seen on the left edge following the path of a deeper part of the rille. A bright spot is visible in or near the MHH. d. shows LRO Diviner silicate minerology data. A pattern of roughness in the data follows the entire rille, distinct from the smooth North-South linear patterns seen outside the rille. This roughness includes the shallower section in which the MHH is found, which was not followed by the patterns in some other datasets. e. shows Diviner rock abundance data, the fraction of rocks as opposed to fine soil present. This data is based on comparison of daytime and nighttime temperatures. No clear pattern is seen following the rille. Figure 36: Comparison of model results for (a.) original model setup (20 km block dimensions) and (b.) model setup using a larger block (140 km block dimensions). The roof thickness used was 0.65 km, and the tube width used was 3.5 km. Interpretations of contours, colors, and red and blue circles are the same as in Figure 19 and similar figures. Figure 37: Comparison of model results for (a.) original model setup and (b.) model setup using a finer COMSOL Multiphysics? mesh. The roof thickness used was 0.65 km, and the tube width used was 3.5 km. Interpretations of contours, colors, and red and blue circles are the same as in Figure 19 and similar figures. Figure 38: Comparison of model results for (a.) original model setup and (b.) model setup using a finer MATLAB? sampling grid. The roof thickness used was 0.65 km, and the tube width used was 3.5 km. Interpretations of contours, colors, and red and blue circles are the same as in Figure 19 and similar figures. xiii List of Abbreviations Geological Strength Index (GSI) Gravity Recovery and Interior Laboratory (GRAIL) Ground Penetrating Radar (GPR) Lava Beds National Monument (LBNM) LIght Detection And Ranging (LIDAR) Lunar Orbiter Laser Altimeter (LOLA) Lunar Radar Sounder (LRS) Lunar Reconnaissance Orbiter (LRO) Lunar Reconnaissance Orbiter Camera (LROC) Marius Hills Hole (MHH) Narrow Angle Camera (NAC) Selenological and Engineering Explorer (SELENE) Wide Angle Camera (WAC) xiv 1. Introduction 1.1. Lava Tubes on Earth A lava tube is a long, quasi-linear subsurface void space inside a lava flow (Huff and Owen, 2015). A lava tube forms when a lava flow develops a solid crust of rock on its sides and roof, by one of multiple formation mechanisms, and the interior lava subsequently flows out and drains the space. If the crust is strong enough to support its own weight, this process leaves a long, empty tube in the path of the flow. Lava tubes are present in basaltic flows in many locations on the Earth, including the lava flow systems of Hawaiian volcanoes (L?veill? and Datta, 2010) (Figure 1) and the Lava Beds National Monument (LBNM) in California (Waters et al., 1990). These tubes typically have widths from a few meters to tens of meters. Some terrestrial tubes can reach widths of approximately 30 meters (Greeley, 1971) and lengths of 65 kilometers (Pipan and Culver, 2019). Terrestrial lava tubes often have holes, where part of the roof is missing and the tube interior is accessible from the surface. These holes are referred to as skylights if they are places where the roof has never fully formed, or collapses if they are caused by the roof partially collapsing. As the formation mechanisms of holes into lunar lava tubes are not known, ?holes? is used here to refer to any opening through a lava tube roof that does not cover the entire tube width, without implying a specific formation mechanism. Some tubes on Earth have a distinct chain of holes, visible on the surface, along their lengths (Paris et al., 2019), including some LBNM tubes as shown in Figure 2. Tubes can form in ?a?? or p?hoehoe flows (Calvari and Pinkerton, 1999). Complex tube systems can 1 form where multiple tubes exist, such as the LBNM in California (Waters et al., 1990) and Mt. Etna in Sicily (Calvari and Pinkerton, 1999). Multiple tubes can merge in some cases to form compound tubes. Tubes can also split, partially collapse, become blocked with lava where they did not drain, form lavafalls where the lava flowed over a sharp drop, and change in width and height at different points (Waters et al., 1990; Calvari and Pinkerton, 1999). Figure 1: The interior of N?huku, a lava tube in Hawaii that is open for visitors to walk through. Image credit: NPS Photo/Janice Wei (J. Wei). National Park Service photograph, available at https://www.nps.gov/havo/learn/news/20180202_nakuku_lighting.htm (accessed May 5 2021). 2 Figure 2: Chain of holes into a lava tube, Lava Beds National Monument, CA. Coordinates at the center of the image are 41.75?N, 121.50?W. Satellite imagery obtained from Google Maps Four distinct formation mechanisms for lava tubes are frequently discussed (Valerio et al., 2008; Coombs and Hawke, 1992). Three of these mechanisms involve channelized flows, where the flow of lava is confined to a relatively narrow flow 3 path, such as by levees on either side. One mechanism occurs when a channelized lava flow develops regions of solidified crust on its surface, spreading from the edges of the flow to eventually meet at the center. A similar but distinct mechanism happens when plates of rafted solidified crust carried on the flow collide and merge to form a roof covering the entire flow and reaching the edges. A channelized flow bounded by levees can also form a roof if the levees grow gradually by lava overflows and spattering, until they meet above the flow (Valerio et al., 2008; Coombs and Hawke, 1992). The fourth formation mechanism does not depend on a channelized flow or levees; this mechanism takes place when an inflated flow on the surface, with solidified crust covering its top surface and sides, progresses through the continual extension of lobes. Lobes forming at the end of an inflated, crusted-over flow allows lava to continue to flow through it while keeping the crust largely intact. If the extending lobes also inflate and form solidified crusts, and then form lobes themselves to extend further, this process can form a long tube (Wentworth and MacDonald, 1953). Inflation by lava can expand tubes after they have been created. Tubes can also be enlarged by the flowing lava downcutting into the surface through erosion. Depending on characteristics of the lava such as the velocity, temperature, and flow type (laminar or turbulent) and the properties of the surface material being eroded, thermal or mechanical erosion may play a dominant role. Either erosion type can cut downwards into the surface, deepening and possibly widening the tube. Evidence of downcutting can be seen in some tubes, where layers of the materials that the flow cut through can be seen (Figure 3) (Greeley et al., 1998). 4 Figure 3: Cross sections of ??inahou Ranch lava tube at increasing distances from the vent (from a to d) over 2.62 km of the tube, showing where the flowing lava has downcut through a paleosol layer into a pre-existing flow (Figure 8 from Greeley et al., 1998) Lava tubes can allow lava to flow much farther from the source vent than it would otherwise be able to (Calvari and Pinkerton, 1998); the crust of a tube decreases heat loss to the surface significantly, insulating the interior lava and 5 keeping it liquid long after it would otherwise have cooled and solidified (Huff and Owen, 2015; Sakimoto and Zuber, 1998). This property can have significant effects on the flow systems of terrestrial volcanoes and the hazard assessments of lava flows, allowing them to reach places they would otherwise be unable to (Calvari and Pinkerton, 1998). Tubes can often have other lava flow through them after the initial flow that formed them. These later flows can leave distinct layers of lava plastered to the tube walls, of varying thicknesses based on the characteristics of the flows (Calvari and Pinkerton, 1999). Plastered layers of lava can fail separately from the rest of the tube. Failure of these onion-skin layers can create distinct tube features and debris (Waters et al., 1990). This can also lead to areas of the tube roof where a collapse has left a cross-section of the distinct lava layers visible (Figure 4) (Calvari and Pinkerton, 1999). 6 Figure 4: Section of Cassone lava tube roof where lava layers plastered onto the roof have collapsed, revealing linings with mm-thick layers. The ruler shown is 25 cm long. The arrows indicate ?pull-apart? stalactites caused by the detachment of a lava layer (Figure 7 from Calvari and Pinkerton, 1999). 1.2. Lunar Observations and Evidence of Lunar Lava Tubes 1.2.1. Sinuous Rilles Sinuous rilles are winding troughs seen in many locations on the lunar mare. These features are found on the flood basalt plains of the lunar maria, and are believed to have been created by ancient lava flows (Hurwitz et al., 2013; Chappaz et al., 2016). These features can have widths up to several kilometers and stretch for hundreds of kilometers (Hurwitz et al., 2013). One such feature, Hadley Rille, was investigated by Apollo 15 (Spudis et al., 1988). As the lavas that formed the lunar 7 maria were very low-viscosity basalts, their paths generally followed the slope of the pre-existing surface topography (Hurwitz et al., 2013), similar to rivers on Earth; this results in the sinuous shapes of these rilles. As lava tubes can form by the same mechanism as sinuous rilles, with the addition of the formation of a crust, it is believed that lunar rilles and tubes are related and can both exist on one lava flow path. Some rilles may lead into underground extensions in the form of tubes (Greeley, 1971; Hurwitz et al., 2013). Several sinuous rilles on the lunar surface show signs of possibly leading into lava tubes. These include locations where a rille appears to have a roofed-over segment, such as near one end of Rima Marius. There, the rille path has a discontinuity (Figure 6). Discontinuities such as this may be caused by a formerly larger tube that has undergone roof collapse over all but a small section, or a rille that only formed a roof and became a tube over a short stretch. The interpretation of this roofed section of Rima Marius as an intact tube segment is strengthened by the analysis of Roberts and Gregg (2019), who conclude by studying nearby smaller, cross-cut rilles that Rima Marius likely formed as a lava tube. Partially roofed rille segments are also present in other locations on the Moon. Hurwitz et al. (2013) surveyed lunar sinuous rilles and located completely or partially collapsed lava tubes in Northern and Southern Oceanus Procellarum and Mare Fecunditatis (Figure 5). They noted that multiple sinuous rilles show signs of being connected to tubes. The authors also noted that, as the focus of their study was on rilles and not tubes, the tubes that they listed are not a comprehensive survey. 8 Figure 5: Sinuous rilles in Northern (a) and Southern (b) Oceanus Procellarum, mapped by Hurwitz et al. (2013). The background maps are Lunar Reconnaissance Orbiter Camera (LROC) Wide Angle Camera (WAC) global mosaics; topography contours on these are displayed in increments of 500 meters. Features identified as sinuous rilles are shown in black. Features identified as collapsed lava tubes are indicated with arrows and labels. The authors note that the examples shown are not a comprehensive survey of collapsed lava tubes (Figure 26 from Hurwitz et al., 2013). In some cases a clear rille discontinuity is not present, but the rille path shows other signs of being connected to a tube. Rima Mairan terminates in two large aligned depressions connected by cracking following the path of the rille, which may be indications of partial collapse of a lava tube that extends from the rille; this rille and its relationship to the model results is discussed in more detail in later sections. An unnamed rille in the Marius Hills region of the Moon is also believed to be related to a lava tube, although it does not show obvious discontinuities. The path of this rille 9 contains the Marius Hills Hole (MHH), which is considered likely to be a hole into a concealed lava tube (Haruyama et. al., 2009; Robinson et al., 2012; Chappaz et al., 2016; Kaku et al., 2017). 10 Figure 6: The Rima Marius lunar sinuous rille with a discontinuity labeled in the figure that may be a lava tube segment. Image credit: Lunar Reconnaissance Orbiter Camera (LROC) Narrow Angle Camera (NAC), Image M135507533R 11 1.2.2. Holes Lava tubes on Earth often have holes through their roofs, which connect the tube interior to the surface environment. These holes have been observed in many locations in terrestrial lava tubes, such as in Hawaiian tube systems (Kauahikaua et al., 1998), at the Lava Beds National Monument in California (Waters et al., 1990), and the Mt. Etna tube system (Calvari and Pinkerton, 1998). Chains of holes occur in some places on Earth, following the path of a tube (Sauro et al., 2020); these can in some cases be used to trace the path of an unexplored tube (Torrese et al., 2021). Holes in a tube roof (Figure 7) can be formed either when a flow does not completely roof over into a tube, leaving a small uncovered area, or when a tube with a complete roof suffers collapse of a small roof section. These areas are not completely uncovered, as a sinuous rille is; the tube roof remains largely intact, with the hole only taking up a portion of its width. 12 Figure 7: Examples of terrestrial (a.) and lunar (b.) holes connected to lava tubes. The example in a. is found in the Corona volcano lava tube system (centered on 29.166? N, 13.457? W) (Satellite image from Google Earth). The example in b. is the Marius Hills Hole (MHH) in the Marius Hills region of the Moon (Modified from Figure 8 from Robinson et al., 2012. LROC NAC image M155607349R). Holes have been observed on the surfaces of the Moon and Mars that are believed to lead into lava tubes. Some of these holes are individual, while others form chains similar to those seen on Earth (Sauro et al., 2020). The sinuous paths that some of these hole chains follow do not appear to be consistent with any other formation mechanism, such as a tectonic origin, besides lava tubes. Additionally, some of these features have a visible vent at the end of the chain, further indicating a volcanic origin (Sauro et al., 2020). In some cases, suspected holes into tubes are associated with sinuous rilles (Hurwitz et al., 2013). The MHH in the Marius Hills volcanic region of the Moon is a particularly well-studied example of a hole on the path of a rille that is believed to lead into a tube (Figure 8). Observations of the MHH using high- 13 resolution imagery from the SELENE (Selenological and Engineering Explorer) spacecraft led researchers to suspect that it was an opening into a larger void space beneath the surface (Haruyama et. al., 2009). Later angled observations using the Lunar Reconnaissance Orbiter Camera (LROC) Narrow Angle Camera (NAC) on the Lunar Reconnaissance Orbiter (LRO) confirmed that the hole leads into a wider underground void space (Robinson et al., 2012). This discovery supports the theory that the MHH is a hole into a lunar lava tube. The dimensions of the hole have been determined through imagery observations; the hole is approximately 53 m wide and 45 m deep (Robinson et al., 2012). The width of the void space below the hole is at least 12 m beyond the width of the hole (Robinson et al., 2012), but only a lower bound estimate of the width can be determined by seeing how much floor space is visible through the hole. As imagery from SELENE and LROC suggested that the MHH could be a hole into a lunar lava tube, the hole and the region around it were later studied using several other methods to investigate this candidate tube. 14 Figure 8: Lunar Reconnaissance Orbiter Camera (LROC) Narrow Angle Camera (NAC) images of the Marius Hills Hole (MHH). (Modified from Figure 8 from Robinson et al., 2012). A-D show images of the hole with various conditions. A: 0.5? emission angle, 25? incidence angle (LROC NAC Image M122584310L). B: 9? emission angle, 13? incidence angle (LROC NAC Image M155607349R). C: 45? emission angle, 34? incidence angle (LROC NAC Image M137929856R). D: 30? emission angle, 13? incidence angle (LROC NAC Image M155614137R). Scale for D not given in Robinson et al., 2012. 15 1.2.3. Gravity Gravity data from the GRAIL (Gravity Recovery and Interior Laboratory) mission also indicates the existence of lava tubes beneath the lunar surface. GRAIL consisted of two satellites, Ebb and Flow, which orbited the Moon and measured its gravity field to an extremely high resolution. The distance between the two satellites was measured using a laser as they orbited the Moon; as the spacecraft passed near a gravitational anomaly such as a mass concentration or mass deficit, the distances between the spacecraft changed as one reacted to the anomaly before the other. Analysis of these precisely measured distances allowed the gravity field of the Moon to be calculated. The data can reveal lunar features that are invisible from the surface. Large void spaces beneath the surface are distinguishable in the gravity data as mass deficits (Chappaz et al, 2016). GRAIL data for multiple areas on the Moon, including the area around the MHH, was studied to search for signs of lava tubes. The data around the MHH revealed a long mass deficit that was consistent with a long void space. The authors of this analysis also created lava tube models to determine the gravity signature that a lava tube would produce; the observed gravity deficit was consistent with the expected signature of a lava tube. This feature was interpreted as a strong lava tube candidate (Figure 9) (Chappaz et al., 2016). 16 Figure 9: Map of the correlation between the free-air and Bouguer gravity anomalies for the region around the MHH, using a cross-correlation technique (modified from Figure 2 in Chappaz et al., 2016). Cold colors designate mass surplus; hot colors designate mass deficits. This combination of anomalies highlights mass deficits buried beneath the surface. The MHH is shown as a red circle; the long buried mass deficit near the MHH is the candidate lava tube (Chappaz et al., 2016). See details in Chappaz et al., 2016 for explanations of data analysis and units. Similar gravity deficits indicating candidate tubes were found in other regions of the Moon. Some such candidates were found near sinuous rilles such as Rima Sharp, while others were in unmarked areas. The authors of the GRAIL data analysis concluded that multiple strong lunar lava tube candidates were detected. Although 17 this method detected multiple tube candidates, exact tube dimensions were difficult to calculate. This method also could only detect relatively large tubes; tubes with widths less than approximately 1 km were undetectable (Chappaz et al., 2016). 1.2.4. Ground Penetrating Radar Ground Penetrating Radar (GPR) data of the lunar subsurface also indicates the presence of lunar lava tubes. The Lunar Radar Sounder (LRS) instrument on the SELENE spacecraft used GPR to investigate the lunar subsurface to depths up to several km. Data from the LRS was studied for the area around the MHH by Kaku et al. (2017), to determine if a lava tube was present in this area. The data showed echo features indicative of underground void spaces in multiple locations in the target area. Most void-space-indicating echo features were associated with gravity deficits seen in GRAIL data; some were associated with the rille path and its possible extension into a tube. A strong void space echo signature was seen near the MHH itself (Figure 10). The authors concluded that a lunar lava tube is present in this location (Kaku et al., 2017). 18 Figure 10: Locations where LRS data showed signatures indicative of void spaces, overlaid on a Bouguer gravity map, made using a cross-correlation method to highlight where the Bouguer anomaly appears consistent with the expected anomaly for a long linear mass deficit such as a tube (gravity map used by Kaku et al., 2017 was developed by Chappaz et al., 2017). Cold colors on the gravity map designate mass surplus; hot colors designate mass deficits. The colors of the circular points show ????, the difference in echo power between the first and second echo peaks; void spaces are most likely at points with low ???? values (black, purple, or red). Most strong void space echo features are present on the Bouguer mass deficit, indicating a tube; the other strong echo features may indicate individual voids surrounded by mass surpluses. The MHH is next to point T1. Surface features, including the rille, are also shown. (Figure 5 from Kaku et al., 2017) 19 The authors could not determine the dimensions of the tube from this analysis. Two distinct interpretations of the radar echo data were discussed (Kaku et al., 2017); although both indicate the existence of a tube, the tube dimensions change based on which interpretation is used, leaving these values uncertain. 1.3. Uses of Lunar Lava Tubes Lunar lava tubes are potentially very useful locations for future lunar exploration and habitation, possibly including the construction of long-term lunar bases (Figure 11), and could also present unique scientific opportunities. The surface of the Moon is constantly bombarded by micrometeorites and radiation, which pose significant risks to unprotected astronauts and equipment on the surface (H?rz, 1985). The environment inside a lava tube is expected to be protected from these surface hazards, as the rock of the roof provides effective shielding (H?rz, 1985; De Angelis et al., 2002). The lunar surface also undergoes extreme, potentially hazardous temperature changes between day and night of approximately 280?C. In contrast, the temperature inside a lunar lava tube is expected to be relatively constant, at approximately -20?C (H?rz, 1985). The shielded environment inside a lunar lava tube would protect anything inside from surface hazards and eliminate the need for manmade shielding built into constructed structures. A lunar habitat constructed on the surface would require heavy shielding to be protected from these hazards, which would increase the weight of the structure, increase the cost and complexity of construction, and limit what materials could be used; inside a tube these issues are not present (H?rz, 1985). 20 Lunar dust on the surface could also interfere with mission equipment and operations, and may be difficult to keep out of a lunar habitat. As most of the interior of a tube is not connected to the surface, a tube could keep large volumes of dust out of the environment around a habitat. Additionally, the walls of a lava tube could be used for additional structural support for large habitats if required (H?rz, 1985). Figure 11: Schematic of one possible design of a lunar base built inside a 100 meter scale lunar lava tube (Figure 3 from Daga et al., 1992). The protected environment inside a lunar lava tube could provide scientific as well as practical benefits, by providing unique opportunities to study a preserved lunar environment. Basalts inside a lava tube would be protected from processes that effect the surface such as solar wind, micrometeorites, and larger impacts, and 21 therefore would not have been altered over time as basalts exposed on the surface would have been. These preserved rocks would allow the volcanic behavior of the Moon at the time of tube formation to be studied in previously impossible ways (Haruyama et al., 2012). In this way, new insights could be gained into the sources of the basalts that formed the tube, how they were erupted, and the overall geological history of the Moon. Lunar lava tubes may also protect volatiles, such as water ice deposits, that could not survive on the exposed lunar surface (Haruyama et al., 2012). A stable temperature consistently below 0?C, and protection from sunlight and surface hazards, could allow such deposits to survive for extremely long periods of time; the same principle has been observed on Earth, where lava tubes such as Skull Cave in the LBNM can preserve water ice for many years despite seasonal temperature changes on the surface (Kern and Thomas, 2014). It is theorized that the same effect could also preserve water ice in Martian lava tubes (Williams et al., 2010). Exploration inside a lunar tube could discover such volatile deposits; these would be of great use for studying how volatiles have evolved on the lunar surface, and the history of the Moon overall. Additionally, tubes may have nearby pyroclastic deposits and other materials that could be used as resources for the long-term maintenance of a lunar base (Coombs and Hawke, 1992). Lunar lava tubes present ideal environments for future lunar exploration and habitation and could offer unique and highly valuable scientific opportunities; however, the usefulness and usability of these tubes could be reduced by failure on 22 their internal surfaces. Failure in these tubes is unlikely to take place in the future, as the tubes have similar ages to the surrounding mare basalts of several billion years (Hiesinger et al., 2010). However, significant roughness and debris fields on tube floors, remaining from ancient collapses, could hinder mission activities and construction. Cracks on the tube floor and debris from roof collapses could interrupt areas of flat floor, making traversal and construction more difficult. Debris and cracks could also have jagged edges that would pose risks to space suits and inflatable structures such as envisioned by Daga et al. (1992). These hazards should be taken into consideration and, if necessary, remediation strategies identified, when planning mission activities inside tubes; it would therefore be of great use for mission planning to be able to determine how much failure has occurred inside a tube and how much usable floor is present before entering it. Being able to determine the dimensions of a tube without entering it would also be extremely useful for planning uses of the tube and what could be constructed inside it; although evidence shows the existence of lunar lava tubes, it is not currently possible to determine their exact dimensions from these methods. In this project, I will describe where internal failure is likely in lunar lava tubes with various dimensions, report how much of the tube floors are expected to be uninterrupted and easily usable, and introduce a method for determining tube dimensions from the surface features they cause. 23 1.4. Prior Modeling Work 1.4.1. General Comparisons Prior works by other researchers have also modeled lunar lava tubes to determine their stability; however, these projects focused on determining if each tube was stable or unstable overall, and not on where failure occurs inside the tubes or what surface features would be produced. These include the works of Blair et al. (2017), Theinat et al. (2020), and Modiriasari et al. (2018). The setup of this research was informed by the assumptions and the results of these prior works. The setup of the modeling work of Blair et al. (2017), Theinat et al. (2020), and Modiriasari et al. (2018) shares aspects with the model setup used in this research. Tubes were simulated using two-dimensional finite-element models and plane strain assumptions, representing the entire tube with one slice and greatly reducing computational costs. The models used in all cases cover only half of the tube, using an axis of symmetry to represent the entire tube, further decreasing computational costs. The modeled tubes were set in larger block of material, with parameters chosen to simulate a lunar rock mass. The mesh covering the entire block was made finer near the tube itself in all cases, although the specific mesh setup was different between projects. The boundary conditions imposed are similar for all the projects: the research described here and in Theinat et al.?s (2020) and Modiriasari et al.?s (2018) work allowed the vertical boundary of the block to move vertically but not horizontally, while Blair et al. (2017) did not allow it to move in either direction. The bottom surface of the block in all projects was not allowed to move in any direction. 24 The same lava tube shape was used in Blair et al.?s (2017) and Theinat et al.?s (2020) work: tubes were idealized as half-ellipses with a 3:1 width to height ratio. This shape is almost identical to the shape used in the main model suites of this research; the tube shape used in this research was informed by these prior projects. The tubes used by Blair et al. (2017) and Theinat et al. (2020) had sharp corners, as opposed to the rounded corners used here; the stress effects of sharp corner are significant to this research but not for these prior works, as Blair et al. (2017) and Theinat et al. (2020) studied the overall stability of the tubes and not regions of internal failure. Blair et al. (2017) and Theinat et al. (2020) both varied the same tube dimensions that were varied in this research: width and roof thickness. However, these dimensions were varied over different ranges, and the authors of these works also changed other parameters that were not varied in this project. Tube shape, which was varied in this research, was not changed in these works. Modiriasari et al. (2018) kept tube width constant, but varied tube shape and roof thickness. 1.4.2. Blair et al. (2017) Blair et al. (2017) modeled lunar lava tubes using the Abaqus? software suite. The authors used tube widths ranging from 250 meters to 10 kilometers, and roof thicknesses from 1 meter to 500 meters. The authors also used both Poisson and lithostatic stress states to determine their results on tube stability; the research introduced here uses only a Poisson stress state. Blair et al. (2017) also introduced 25 tectonic strains into their models by forcing the far vertical edge of the block to move by set distances, putting the block in compression or tension. Blair et al. (2017) studied the behavior of the lava tube roofs using a Mohr- Coulomb plastic failure envelope to determine tube overall stability. Wherever the rock exceeds this envelope and plastic strains are present, the material is considered to have failed. The authors divided their modeled tubes into three stability categories, based on the percentage of the roof thickness through the center of the roof that exceeded the failure envelope in each case: stable (0%), quasi-stable (less than 50%), and unstable (50% or more). The authors found that kilometer-wide tubes could be stable with a roof thickness as small as 2 meters, and that stable tubes with widths up to 5 kilometers are possible. Quasi-stable tubes up to 6.75 kilometers wide were also found to be possible (Blair et al., 2017). The full set of stability results is shown in Figure 12. The results found by Blair et al. (2017) showed that tubes with very large widths could be stable on the Moon, and could exist with relatively thin roofs. This result supported conclusions from GRAIL gravity data that lunar lava tubes with widths of several kilometers may be present (Chappaz et al., 2016), by indicating that such tubes could feasibly exist. Comparison of Blair et al.?s (2017) stability results with patterns of failure seen in the research described here are given in Section 4.1. 26 Figure 12: Blair et al. (2017) stability results for tubes of different widths and roof thicknesses in lithostatic (upper panel) and Poisson (lower panel) stress states (Figure 3 from Blair et al., 2017). The bold line in the Poisson stress state data divides two types of failure through the tube roof: compressive below the line, and tensile above it (Blair et al., 2017). 1.4.3. Theinat et al. (2020) Theinat et al. (2020) also used the Abaqus? software suite. The authors used tube widths ranging from 300 meters to 4 kilometers, and roof thicknesses from 1 27 meter to 2 kilometers. In addition to varying tube dimensions, the authors also varied the tensile strength of the material used. To determine tube stability the authors studied the values of convergence, or the radial displacement difference between two diametrically opposed points chosen on the perimeter of the opening, for each tube (Theinat et al., 2020). The authors found that a tube 300 meters wide could be stable with a 1-meter-thick roof, a kilometer-wide tube could be stable with a minimum roof thickness of 100 meters, and the largest modeled tube, 4 kilometers wide, was stable with a minimum roof thickness of 1 kilometer (Theinat et al., 2020). The full set of results is given in Figure 13. Theinat et al. (2020) overall found tubes to be less stable than Blair et al. (2017). Both show a general trend of wider tubes requiring thicker roofs to remain stable, but Theinat et al. (2020) found larger minimum roof thicknesses required for tube stability. The authors also used roof thickness values above 500 meters, which Blair et al. (2017) did not, and found all tubes with these thicknesses to be stable (Theinat et al., 2020). 28 Figure 13: Stability results for tubes of different widths and roof thicknesses (Figure 13 from Theinat et al., 2020). 1.4.4. Modiriasari et al. (2018) Modiriasari et al. (2018) also used the Abaqus? software suite. The authors used a constant tube width of 4 kilometers, and roof thicknesses of 20, 50, 100, and 200 meters. The authors used three different tube shapes; one was a half-ellipse with a 3:1 width-to-height ratio, as used by Blair et al. (2017) and Theinat et al. (2020), and modified for this research. The other two shapes used had width-to-height ratios of 2:1 (semicircular) and 4:3 (half-ellipse, for which the vertical axis of the full ellipse would be longer than the horizontal). 29 The authors used a Mohr-Coulomb plastic failure envelope to determine where failure is present, as Blair et al. (2017) did, and used the same categorizations for stable, quasi-stable, and unstable tubes as Blair et al. (2017). The authors found that 4-kilometer-wide, half-elliptical tubes with a 3:1 width-to-height ratio were quasi-stable with roofs 50 to 200 meters thick, but unstable with a 20-meter-thick roof. These specific results agreed with the results found by Blair et al. (2017), and agreed with the trend of stability increasing with increasing roof thickness. Modiriasari et al. (2018) found that decreasing the tube width-to-height ratio decreased the amount of failure through the thickness of the tube roof. Of the tubes that they modeled, the percentage of roof thickness that failed was largest for the 3:1 ratio tubes and smallest for the 4:3 ratio tubes, for every roof thickness used (Figure 14). Comparison of these results to the results of the study presented here concerning the effects of tube shapes are discussed in Section 4.3. 30 Figure 14: Modiriasari et al. (2018) results for tubes of different roof thicknesses and aspect ratios. Overburden yield (%) is the percentage of the thickness of each tube roof undergoing failure above the tube center. Tubes are considered stable (0% overburden yield), quasi-stable (some failure, but less than 50%), or unstable (50% or more) (Figure 1 from Modiriasari et al., 2018). 31 2. Methods 2.1. Modeling Process 2.1.1. Main Model Suites I have produced and analyzed models simulating lava tubes with lunar conditions to determine the deformation, stress, and failure that they create. The tube models were created and run using COMSOL Multiphysics?, a finite-element simulation software, and its LiveLink? for MATLAB? feature. The Livelink? for MATLAB? feature allows the software to be manipulated by MATLAB? code instead of manually, making it possible to automate the creation, running, and analysis of models. This allows model suites to be created much more reliably and efficiently than by hand, and allows entire new sets of models to be introduced and run easily. The lava tube models are two-dimensional, using plane strain assumptions, and use an idealized shape of a cross-section of a lava tube. Each model simulates half of a tube, and is symmetric along the left edge. As such, each model effectively functions as a full tube that extends infinitely in both directions out of the modeling plane, without needing to model the entire three-dimensional structure (Figure 10). The momentum conservation equations are solved assuming steady-state in the region surrounding the tube. The rheology is elastic with shear modulus G = 12.163 GPa (derived from equations in Marinos and Hoek, 2000) and Poisson ratio ? = 0.3 (Theinat et. al, 2020; Schultz, 1995). The model is loaded simply by gravity, assuming a lunar gravitational acceleration of -1.662 m.s-2. Each model structure 32 consists of a 20-km-by-20-km block of material, with a void space cut out of it to form the tube. The material making up each block is made to simulate a lunar rock mass, instead of a completely intact lunar basalt. A rock mass is partially fractured, such as by cooling, tectonic strains, or nearby impacts, leading to degrees of blockiness or fragmentation. The model block features a reflective condition on the left side to serve as an axis of symmetry, a fixed constraint at the bottom of the block, and a surface on the right side, away from the tube, which is allowed to move vertically but not horizontally; these conditions are used to make the block behave as if it is surrounded on the sides and below by more lunar material, as a tube would be in reality. The remaining surfaces, consisting of the top surface and the interior tube surfaces, are free, and the model is allowed to deform. The shape of the tube used for the main model suites was idealized as a half- ellipse with rounded corners (Figure 15). The aspect ratio used for these suites of half-elliptical tubes was a 3:1 width-to-height ratio, giving a full ellipse of the same aspect ratio a 3:2 width-to-height ratio. The half-ellipse shape and ratio follow previous modeling works on lunar lava tubes (Blair et al., 2017; Theinat et al., 2020). The interior of the tube is left empty and the outline of the tube is allowed to deform. The difference in density between the empty tube interior and the surrounding material generates a complex stress field in which the tube tends to rise relative to the rest of the rock mass. Rounded corners were added after the use of sharp corners in initial models caused unrealistically high stress to occur near the tube edges; artificially sharp corners were also not consistent with observations of real tubes. After introducing rounded corners, the anomalous high-stress concentrations at the 33 corners were reduced, and the stresses throughout the rest of the tube were not greatly affected. 34 Figure 15: (a.) Diagram of tube model geometry, showing the roof thickness, the tube half-width, and the direction of gravity. The model is symmetric across the dotted line. Only a portion of the model near the tube is shown in a.; b. shows the entire block for this example. c. shows the mesh used over the entire block; a significantly finer mesh is used near the tube itself than for other regions of the block 35 Tube models with systematically varying dimensions were run to determine the effects of tube dimensions on the deformation, stress, and failure produced. The tube full width and the tube roof thickness, or burial depth, were varied within the limits shown in Table 1. Roof thickness represents a combination of the thickness of the original layer that formed the tube roof and the thicknesses of additional layers of mare emplaced over the tube after it formed. While the tube shape is kept unchanged, these two dimensions completely describe the tube geometry. A parametric sweep was used to construct and run multiple models at a time with different dimensions, keeping the shape of the tube otherwise unchanged. Multiple different combinations of dimension ranges were tested to observe the results before a set of dimensions was chosen to focus on. The final set of dimensions chosen for the main model suites were tube widths of 0.5 through 6.5 km, with increments of 1 km, and roof thicknesses of 50 m through 950 m, with increments of 300 m. These dimensions were chosen based on the widths of lunar features such as sinuous rilles, the prior modeling results of Blair et al. (2017) and Theinat et al. (2020) concerning stable tube dimensions, and my experimentation with different ranges of dimensions. 36 Table 1: Width and roof thickness values used for the modeled tubes. Each combination of dimensions was used. 2.1.2. Model Suite Variants The focus of this modeling was on lunar lava tubes; however, some models with different starting parameters to simulate terrestrial conditions were also used, for comparison with the modeled lunar tubes and real terrestrial tubes. This was done by increasing the gravitational acceleration to -9.8 m.s-2, using smaller and shallower tube models to match observations of terrestrial tubes, reducing the size of the surrounding block to follow the smaller tube sizes, and adjusting the mesh used to account for smaller tubes and thinner roofs without compromising accuracy. The tube shape was also changed for some models to determine its effect on stress and failure. The shape used for most models was a half-ellipse with a 3:1 width-to-height ratio and rounded corners and a flat floor; this shape was used as the basis for the different variations of tube shape used. For some models, this shape was 37 kept unchanged except for the width-to-height ratio, or aspect ratio; this ratio was decreased to 2:1 and increased to 6:1, resulting in semicircular tubes and half- elliptical tubes with a flatter shape, respectively. In other models, the original width-to-height ratio was kept unchanged, but a curving floor was used instead of a flat floor. Curving floors were implemented by making the tube floor concave, with the geometry of a half-ellipse; a curving floor dips to its maximum depth at the tube center. Degrees of floor curvature were designated and categorized by their depth at the tube center as a percentage of the tube half-width: a tube with a 0% floor curvature has a completely flat floor, as in the original models, and a tube with 100% floor curvature has a floor that dips deeply enough to form a semicircle. These shape alterations were chosen to investigate results seen in initial test models with full-circle and full-ellipse shapes; the initial full-ellipse models had a 3:2 full-width to full-height ratio. After models with these shape variations were tested, they were combined to make a set of models with each of the three chosen aspect ratios over a range of floor curvature values between 0% and 100%; as an example, using a tube with a 2:1 width-to-height ratio and applying a floor curvature of 100% results in a completely circular tube. The values of the width-to-height ratios were defined for tubes with flat floors, and do not take into account the added tube height from the floor curving downwards. For tubes with curving floors, this ratio is effectively the width-to-height ratio of the top half of the tube alone, without floor curvature; as such, this value is sometimes referred to as the top-half width-to-height ratio. 38 A subset of the different tube shapes run is shown for illustration in Figure 16. Each different tube shape model was run with large and small widths and roof thicknesses, to determine the effects of tube dimensions on differently shaped tubes. Figure 16: Representative geometries of the tube shapes considered. This selection includes tubes with all three top-half width-to-height ratios and a subset of the floor curvatures used, including the deepest-dipping and the flat floors. All tubes shown here have roof thicknesses of 950 m and widths of 5.5 km. 39 2.2. Model Analysis In addition to being used to create the model suites, MATLAB? code is also used to export the data of interest from the models and to analyze it. A MATLAB? master script connects separate functions and scripts in the modeling process and controls a set of parameters and model aspects that can be changed between runs. After a suite of models have been created and run successfully, the master script applies data analysis code to the models. Running models in COMSOL? allows stresses across the rock to be found and displayed; the analysis code allows data values at points of interest to be determined and allows the results to be analyzed and interpreted. Figure 17 shows the variation of tensile stress obtained for a tube of width 5.5 km and a roof thickness of 950 m. This is an example of the raw stress data produced by COMSOL Multiphysics?, before analysis code is applied. 40 Figure 17: Example unprocessed results from one COMSOL? model. The plotted stress is labeled as the First Principal Stress, and is positive; due to COMSOL?s? sign and stress naming conventions, this is the tensile stress, considered in this work as the Third Principal Stress and negative. Appropriate sign conversions are applied when reading data from COMSOL?. Stresses are only plotted for the upper region of the modeled block, to avoid plotting large-magnitude stresses far from the tube (such as at the bottom of the block, caused by the entire block?s weight) and making the colorbar less appropriate for the stresses close to the tube. The analysis scripts, as well as the scripts for defining model geometries, are changed for different tube shapes. The analysis scripts first extract and then analyze the desired data from the completed models. The data points of particular interest are the values of the principal stresses, in the bulk of the rock and especially at its edges where it meets the tube or the surface above. For each model a uniform grid is created 41 to cover the block, and reshaped to meet the outline of the tube (Figure 18); the reshaping of this grid is one aspect of the analysis that must be adjusted for different tube shapes. Stress values across the bulk rock are interpolated at the points of this reshaped grid. This reshaping is done to ensure that the tube edges themselves are sampled instead of nearby points, so that the stress and failure of tube internal surfaces are determined correctly. These stress values are plotted for visual interpretation. Figure 18: Example of the grid used to interpolate data from a model. The grid in a. is the original, uniform grid. The grid in b. has been reshaped to fit the outline of the tube (green); reshaped grids like this are used for interpolation. The plots are zoomed in on the region of the block containing the tube. Once the relevant data has been exported, other code compares the stress values to two failure criteria, tensile and compressive, to determine if failure occurs 42 and where. These failure criteria require values for strength parameters of the rock mass to be calculated. A Geological Strength Index (GSI) is used to describe the strength of a rock mass. The GSI is a classification created to describe the behavior of rock masses and allow their strength and deformation properties to be estimated. GSI values are chosen based on the structure of the rock mass (how intact or disintegrated it is) and its surface conditions (how rough or weathered the surfaces are, and whether coatings or fillings are present) (Marinos and Hoek, 2000). GSI values range from 100, for rock that is intact or has widely spaced discontinuities that have a rough surface with no evidence of weathering, to 0, for a rock mass that has no blockiness but highly weathered surfaces with coatings or fillings of soft clay (Marinos and Hoek, 2000). The GSI value chosen to best represent a lunar rock mass was 70, based on examination of terrestrial tubes (Theinat et al., 2018), assuming zero aqueous weathering on the Moon, and a ?blocky? rock mass structure (Blair et al., 2017). For some models, GSI was varied to 60 and 80 to represent various degrees of damage that could be present in the mare and to test the effects of reasonable GSI variations on the model results. Material parameter values for a completely intact lunar basalt were also chosen, and combined with the GSI to calculate material parameters of the basaltic rock mass. The intact rock parameters required for these calculations were ???, the compressive strength, and ??, a dimensionless strength parameter analogous to the rock?s frictional strength. Using high values of ?? results in steeper Mohr envelopes and higher instantaneous friction angles than for lower ?? values for low effective normal stresses (Eberhardt, 2012). For intact rock with an ?? value much 43 ? greater than 1, ?? is approximately equal to the ratio ? (Hoek and Brown, 1980). |??| The value of ??? was chosen to be 266 MPa (Schultz, 1995). The value of ?? was chosen to be 25, based on similar analyses (Marinos and Hoek, 2000; Blair et. al., 2017). ?? values of 20 and 22 were also used in some models, to test the effects of different values chosen in analyses (Schultz, 1995; Blair et al., 2017). The GSI value and these intact rock parameter values were used to calculate two dimensionless parameters for the rock mass, ? and ??, using Equations 1 and 2: GSI?100 ? = exp ( ) (1) 9 GSI?100 ?? = ?? exp ( ) (2) 28 ?? is a version of the ?? parameter adapted for a rock mass using GSI; more fractured and damaged rock masses have lower ?? values (Eberhardt, 2012). ? represents how fractured the rock mass is. ? = 1 for a completely intact rock, and ? = 0 for a very fractured material with a tensile strength of zero. The cohesion of the rock mass is related to this parameter (Eberhardt, 2012). The calculated ? and ?? values used for the rock mass for the central suite of models were ? = 0.0357 and ?? = 8.5630 (see Table 2). A value of ? for the intact rock did not need to be chosen, as ? = 1 for all intact rock. The density used for the material was ?? = 3100 kg. m?3 (Blair et al., 2017; Kiefer et al., 2012). The tensile failure criterion is the Von Mises criterion ?3 = ??, where ?3 is the third principal stress (most tensile stress) and ?? is the tensile strength of the rock 44 mass (Table 2). ?? was calculated from the material parameters for the rock mass and the intact rock by Equation 3. ? ? = ??? (?? ? ?? 2 ? + 4?) (3) 2 The tensile strength value calculated was 1.11 MPa; see Table 2 for a summary of calculated and chosen strength parameters. Tensile failure is considered to be present if ?3 meets or exceeds Equation 3. For compressive failure, the Hoek-Brown failure criterion is used. This criterion is given by Equation 4, ?1 = ?3 + ???????3 + ? ? 2 ?? (4) where ?1 is the first principal stress. The Hoek-Brown failure criterion is an empirically derived relationship describing the strengths of intact rocks or rock masses (Eberhardt, 2012; Hoek and Brown, 2019). Unlike the Mohr-Coulomb failure criterion, the Hoek-Brown criterion describes a non-linear relationship of rock shear strength on normal stress. This nonlinear behavior was one reason why the Hoek- Brown failure criterion was chosen for this analysis, as this allows the criterion to better fit observed data for small normal stresses than the Mohr-Coulomb criterion (Eberhardt, 2012). Small normal stresses are expected in some regions of the models, such as near the surface of the material where the stress created by surrounding and overlying rock is expected to be small. This criterion was also chosen because of its applicability to rock masses. The environment around lunar lava tubes is expected to be a fractured rock mass, as is true for terrestrial tubes, instead of completely intact 45 rock (Theinat et al., 2020). The Hoek-Brown failure criterion can be applied easily to rock masses using the Geological Strength Index (GSI). The GSI was developed to use observations of a rock mass, and strength parameters of the intact rock, to calculate appropriate values of the parameters used in the criterion (Eberhardt, 2012; Hoek and Brown, 2019). As GSI was developed to be applied to this criterion (Hoek and Brown, 2019), the criterion can be applied to a rock mass in a straightforward way that has been specifically designed for the purpose and tested by extensive use by many researchers (Eberhardt, 2012). 46 Variable or Parameter Value Source First principal stress -- -- Third principal stress -- -- Tensile strength 1.11 MPa Calculated value (Equation 3) Compressive strength 266 MPa Schultz, 1995 (intact rock) Rock strength parameter 0.0357 Calculated value (Equation 1) Rock strength parameter 8.5630 Calculated value (Equation 2) Table 2: Variables and constant parameters used in tensile and compressive failure criteria After failure criterion calculations have been completed, the interpolated stresses inside the rock are plotted with overlaid regions of calculated failure for visual interpretation. The code displays the stresses inside the rock using one of several plotting options, with regions of failure overlayed in red (tensile failure) and blue (compressive failure) circles, as in Figure 19. Regions where the rock is under tension are also displayed as semitransparent red shading. The value plotted in 47 contours in the bulk rock in this figure is a ratio (Equation 5) derived from the Hoek- Brown criterion. ?1?? ? 3HB = (5) ???? 2 ???3+? ??? This value shows how close the rock is to Hoek-Brown failure at any point, and indicates failure when it is greater than or equal to one (shown in white in Figure 19). Figure 19: Model output for one modeled lunar tube, with a width of 6.5 km and a roof thickness of 650 m. The filled colored contours show the Hoek-Brown ratio described in Equation 5; compressive failure is present where this ratio is at least one, shown in white on the color scale. Shaded red designates where the rock is in tension. Circles plotted on the tube outline show tensile (red) and compressive (blue) failure. When running through a parametric sweep of tubes of varying dimensions, in addition to producing plots for user analysis of each tube model individually the code also calculates multiple diagnostic variables for each model run and compiles them 48 for the entire model suite. This list of diagnostic variables summarizes the behavior of each tube. These variables include the maximum tensile stresses on the tube roof and floor, the percentage of the tube floor that undergoes tensile failure, the location of maximum tensile stress on the surface, the location of maximum upwards deformation on the surface, and other quantities. These variables are collected for each tube and stored along with the tube dimensions, allowing the overall behavior of tubes of different dimensions to be easily viewed, compared, and analyzed. These variables are later plotted to observe how variables describing tube behavior are affected by the tube dimensions and shapes used. The analysis of multiple diagnostic variables with respect to tube dimensions and tube shape is described in later sections. 2.3. Lunar Surface Observations In addition to the modeling work described in Sections 1 and 2, observations were also conducted of the lunar surface, studying data from several lunar missions, to determine if the surface features predicted by the models could be located in actual examples. This study was conducted using the LROC Quickmap online application, which allows access to many datasets of the lunar surface from LROC and other sources; these include other instruments on the LRO, GRAIL, SELENE, and Lunar Prospector, among others. Individual datasets include, among others, WAC (Wide Angle Camera) and NAC (Narrow Angle Camera) images from LROC, compositional data from the Multiband Imager on SELENE, and surface roughness data from LOLA (Lunar Orbiter Laser Altimeter) on LRO. 49 Multiple sinuous rilles and their surrounding areas were examined in LROC images. LROC imagery taken with several different solar incidence angles was referred to, to better visualize the lunar surface under a range of lighting conditions (examples shown in Figure 20). The models predict linear patterns of surface cracking (described more in Sections 3.2. and 4.2.); particular attention was paid to features such as this in the surface observations. Figure 20: Region of a sinuous rille including the MHH viewed using LROC Quickmap with different lighting conditions. Location: 14.1?N, 56.8?W. 12.04 m/pixel. a.: Large solar incidence angle, East illumination only. b.: Large solar incidence angle, West illumination only. c.: Small solar incidence angle (0 ? 40 degrees). d.: Medium solar incidence angle (40 ? 55 degrees). e.: Large solar incidence angle (55 ? 80 degrees). Images were from the LROC NAC. 50 A more in-depth examination was conducted of the area around the MHH; this is an area of particular interest as it contains a sinuous rille, a likely hole into a tube, and a strong candidate lava tube (see Section 1.2.). The MHH region was investigated in visual imagery in more detail than other sinuous rilles, and was studied using multiple other available datasets from different instruments. This in-depth study was intended to find features that may not be visible in LROC images, such as if cracking near a tube was only visible in LOLA data showing the roughness of the surface, or if fractured regions had different thermal behaviors from intact regions in a way that could be observed in Diviner temperature data. 51 3. Results 3.1. Regions of Failure on Tube Internal Surfaces 3.1.1. Failure Regions The code calculates where tensile and compressive failures occur on the interior surfaces of each tube using the failure criterion and material parameters discussed in Section 2. The models predict extensive regions of failure present inside almost all of the tubes simulated. To allow the behavior of these internal failure regions to be effectively compared and studied, they were categorized into three types: tensile failure on the tube floor, tensile failure on the tube roof, and compressive failure on the tube wall. The diagnostic variables calculated and saved for each tube include the fractional extents of each of these failure types; fractional extent is defined as the horizontal extent of the failure as a fraction of the tube width. Individual tube model results were also examined manually. The most common type of failure found in the models is tensile failure of the tube floor, which is present in all but one modeled tube (Figure 21). Tensile failure occurs centered on the center of the tube floor and spreading out towards the sides, leaving only relatively small unaffected regions near the tube walls. This tensile failure occurs as the tube floor bows upward due to the lack of pressure from the empty tube compared to the rock outside of it. As the floor flexes upwards and takes on a convex shape, the rock is stretched and put under tension. The largest flexure, and strongest tensile stresses, develops at the center of the tube floor. This tensile stress eventually leads to cracking, in a type of failure sometimes referred to as 52 extrados cracks (e.g. Borghi et al., 2015). In all cases where it is present, this failure covers at least 65% of the tube floor. It covers over 80% of the floor in most cases and reaches up to almost 95% of the floor in some, generally wider and deeper, cases. This type of failure is generally the most extensive of the three types. Figure 21: Fractional extent of each tube floor that undergoes tensile failure vs. tube width. Colors and marker sizes denote different roof thicknesses used. The inset shows a simplified visualization of where this failure takes place. 53 The second-most commonly found type of failure is tensile failure of the tube roof (Figure 22). Similarly to tensile floor failure, tensile roof failure occurs centered at the peak of the tube roof and spreads outwards towards the edges. This failure is caused by the tube roof bowing inwards into the void space under its own weight, as it is not supported from below as the surrounding rock is. This flexing stretches the bottom surface of the roof, and may cause tensile failure if the tensile stresses are large enough (intrados cracks) (Borghi et al., 2015). As with tensile floor failure the flexure, and the magnitude of the tensile stress caused, is largest at the tube center. This type of failure is generally less extensive than tensile floor failure, with a maximum fractional extent of under 40% of the tube width. Tensile roof failure occurs mostly in tubes that are relatively wide and deep, leaving the majority of the shallowest and narrowest tubes unaffected. 54 Figure 22: Fractional extent of each tube roof that undergoes tensile failure vs. tube width. Colors and marker sizes denote different roof thicknesses used. The inset shows a simplified visualization of where this failure takes place. Compressive wall failure occurs in the most limited subset of the modeled tubes; it is present only in large and deep tubes (widths at least 4.5 km; roof thicknesses at least 350 m). This type of failure initiates at the tube corner and spreads upwards along the wall, increasing in extent with increasing width and roof thickness. The horizontal extent of this failure is less overall than that of either of the other failure types, even for the largest and deepest tubes. This type of failure, like tensile 55 roof failure, may create debris that would drop onto the floor below, causing significant debris fields. 3.1.2. Usable Floor Extents From the three failure types described above, the usable floor extent for each modeled tube is determined (Figure 23). Usable floor is the region of the floor that has is expected to be smooth and uninterrupted by cracks or debris, and easily usable for traversal, exploration, and construction. The usable floor is defined as any region of the floor that has not undergone tensile failure and is not directly underneath any region of tube roof that is covered by tensile or compressive failure. The floor undergoing tensile failure is considered unusable because this failure is expected to produce significant cracking; the floor underneath any type of roof failure is considered unusable to account for debris fields created by fallen fragments of rock. Both fallen debris and cracked surfaces would be obstacles for movement and construction, and could have jagged edges that would pose risks to space suits and pressurized structures. 56 Figure 23: Usable floor extents vs. tube widths for lunar (a.) and terrestrial (b.) tubes. Colors and line widths denote different roof thicknesses used. The inset shows a simplified visualization of where the failure types take place and where usable floor (green) is found. The green line represents the center of the tube. In most lunar tubes, tensile failure on the floor leaves only a small unaffected region near the tube edge; in some cases, compressive failure spreading from the tube edge renders some of this small region unusable. The combination of these factors leaves relatively small extents of usable floor near the tube edges (Figure 23 a). The maximum usable floor extent found in all the lunar tubes modeled is only 470 m. This extent could be reduced further by other sources of debris and failure that are not captured in these models, such as the failure of distinct lava layers pasted onto the wall, if they are present. The widths of usable floor regions are generally greater in shallower tubes than in tubes with thicker roofs, but do not depend in a predictable 57 way on tube widths. The increased widths of usable floor regions in shallower tubes are interpreted as being due to the decreased weight of rock surrounding the tube causing smaller stresses and less deformation on the tube floor, therefore leading to less extrados floor cracking. The finding that large regions of debris and cracking are likely to be present would be highly relevant for planning missions to explore and utilize these tubes. In contrast to lunar cases, no failure of any kind occurs on the internal surfaces of any modeled terrestrial tube. As a result, the entirety of the floor of each terrestrial tube is considered usable (Figure 23 b). Observations from inside real terrestrial tubes show that significant debris and failure are often present, in contrast to these modeled results (Waters et al., 1990). Examples of failure in terrestrial tubes can be found in the lava tube systems of Lava Beds National Monument (LBNM) in California, among other locations; Waters et al. (1990) provide detailed maps of some of the tube systems in LBNM, showing many regions of debris. However, this contrast between observed and predicted failure could be explained by several factors that are not captured in the current models, such as thin lava layers, deposited after the tube was formed, failing independently from the main rock. Failure of such veneer layers is observed in multiple terrestrial locations. Examples are documented by Calvari and Pinkerton (1999) in Mt. Etna tubes and by Waters et al. (1990) in LBNM tubes, such as in Arch Cave. Arch Cave has layers of pasted lava veneer that have partially or completely peeled off along a significant portion of its length, resulting in curling benches and onion-skin-like debris (Waters et al., 1990). Other unaccounted for factors that could introduce failure, including weathering and seismic 58 shaking, may be significantly less important on the Moon as compared to the Earth; these factors may lead to failure in terrestrial tubes but not be expected to cause failure in lunar tubes. Another such factor is subcritical crack growth, which can allow materials to fail at stress conditions significantly below those usually required for failure. Subcritical cracking is believed to be a prevalent and significant factor contributing to rock failure on Earth (Atkinson, 1984). As this cracking is facilitated by the presence of water or humidity (Atkinson, 1984; Eppes and Keanini, 2017), it may not be a significant factor capable of triggering failure in the extremely dry environment of the Moon. 3.2. Surface Features The modeled tubes predict several surface features caused by the tubes on the lunar surface above them. Such features would be potentially observable around real lunar tubes. Some quantities describing these surface features were found to depend approximately linearly on the widths of the tubes that produce them. These quantities include the distance from the center of the tube to the peak of the elastic surface bulge and to the point of maximum surface tensile stress (shown in Figure 24). The bulges and regions of tensile stress are caused by the tube roof and surrounding surface acting similarly to an elastic plate (Turcotte and Schubert, 2014); as the roof flexes downwards into the void space under its own weight, regions of the surface to either side are forced to flex upwards, stretching the surface and putting it under tensile stress. Where this stress exceeds the tensile strength of the rock mass, tensile failure is expected to occur. This process results in two types of potentially observable surface features: an elastic bulge and a line of tensile cracking. Both are expected to occur 59 symmetrically on either side of the tube, running parallel to the tube path. Tubes are therefore predicted to produce distinctive and possibly informative patterns of surface features, consisting of a long depression over the center of the tube and a line of surface cracking and a surface bulge running parallel to the tube at predictable distances away from it. As these predictable distances were found to depend on tube width by linear relationships, these features if observed could be used with these relationships to deduce approximate values of tubes? widths from the surface. Figure 24: Sketch showing failure regions and distances to surface features. The distances from the tube center to these surface features are defined as ??, the distance to the location of maximum tensile stress on the surface, and ??, the 60 distance to the peak of the elastic surface bulge. The relationships between these distances and tube widths are shown in Figures 25 (??) and 26 (??); both quantities increase linearly as tube width increases. These quantities also increase linearly with tube roof thickness; the relationships of ?? and ?? to both width and roof thickness are given in Section 4.2. Figure 25: Horizontal distance from the center of the tube to the maximum tensile stress point on the surface, ??, vs. tube width. Colors and marker sizes denote different roof thicknesses used. The inset shows a simplified visualization of the tube geometry and this distance on the surface. 61 Figure 26: Horizontal distance from the center of the tube to the peak of the elastic bulge on the surface, ??, vs. tube width. Colors and marker sizes denote different roof thicknesses used. The inset shows a simplified visualization of the tube geometry and this distance on the surface. The ratio of these two distances, ??/??, was also calculated for both terrestrial and lunar cases. This ratio was found to be constant in each case: 1.73 ? 0.10 for lunar and 1.80 ? 0.14 for terrestrial tubes (Figure 27). Of the two surface 62 features discussed, the peak tensile stress locations are more likely to be observable, and hold the best promise for constraining tube width. Surface bulges are expected to be shallow, and may be lost in the roughness and other features of volcanic fields, while surface tensile stress maxima are expected to produce significant cracking in almost all cases. However, it is also possible that surface tensile cracking could be covered by subsequent lava flows on the surface and rendered undetectable, if later flows occurred near the tube. If both kinds of features can be observed, however, the constant ??/?? ratio could be used to deduce the location of an otherwise undetectable tube. Figure 27: Ratio of ?? over ?? vs. tube width (km) for modeled lunar (a.) and terrestrial (b.) tubes. The average ratio value is 1.73 ? 0.10 for lunar and 1.80 ? 0.14 for terrestrial tubes. 63 3.3. Effects of Tube Shapes For the suite of models run with different tube shapes, the maximum tensile stresses on the tube floors were examined to determine what aspects of tube shape encouraged or discouraged tensile floor failure. This type of failure in particular was studied because it was the most extensive type found in the main model suites and was generally the most important factor controlling usable floor extents. The models show that both the width-to-height ratios and floor curvatures of the tubes affect their floor stresses significantly. The maximum tensile stress on the floor decreases with increasing floor curvature. Using a deeper curvature for a tube?s floor therefore decreases the extent of tensile floor failure (Figure 28). The smallest stress values are found in tubes with deeply-dipping floors, including those that form complete semicircles. Some of these tubes with large floor curvatures have stresses below the tensile strength of the rock, and therefore do not experience any tensile floor failure. The largest floor stress values are found in tubes with completely flat floors, such as those shown in Section 2.1.1. 64 Figure 28: Comparison of results of tubes with the same top-half width-to-height ratio (6:2; original half-ellipse) but different floor curvatures: (a.) 33.3%, (b.) 66.6% (symmetric), and (c.) 100% (semicircle). The models with more curving floors have less tensile floor failure; the 100% floor curvature tube has no floor tensile failure at all. Tube width: 5.5 km; tube roof thickness: 950 m. The filled colored contours show the Hoek-Brown ratio described in Equation 5; compressive failure is present where this ratio is at least one, shown in white on the color scale. Shaded red designates where the rock is in tension. The tube shape is outlined with a black line. Circles plotted on the tube outline show tensile (red) and compressive (blue) failure. Horizontal axis: horizontal distance (km); vertical axis: depth (km). Tubes with a larger width-to-height ratio for their top halves have generally greater magnitudes of the maximum tensile floor stress. This relationship is true for all ratios used; tubes with a top-half width-to-height ratio of 6:3 (semicircle) generally have lower floor tensile stresses than tubes with a 6:2 ratio (half-ellipse, as used in main model suites), which have generally lower stresses than 6:1 ratio (flatter half- 65 ellipse) tubes. This ratio can affect floor stress enough to control whether floor failure occurs or not in some cases. Even in cases where floor failure occurs, changing this ratio can control the extent of the tube floor that undergoes failure (Figure 29). Figure 29: Comparison of results of tubes with the same floor curvature (66.6%) but different top-half width-to-height ratios: (a.) 6:1 (Flatter half-ellipse), (b.) 6:2 (Original half-ellipse), and (c.) 6:3 (Semicircle). Both models with a half-elliptical top half have tensile floor failure, although the 6:2 tube has less failure than the 6:1 tube; in the 6:3 tube, floor failure is prevented completely. Floor curvature has a greater effect on floor stress than width-to-height ratio, for the ranges used. A flatter half-elliptical tube with 33.3% floor curvature (top and bottom halves are symmetrical) has significant tensile floor failure, while a purely circular tube does not. Changing the floor curvature of the flatter half-elliptical model to match the circular successfully prevents floor failure (Figure 30). In contrast, 66 changing the aspect ratio of the model to match the circular, while leaving floor curvature unchanged, still results in failure (Figure 31). 67 Figure 30: Model output for a modeled lunar tube, with a width of 5.5 km and a roof thickness of 950 m. The width-to-height ratio of the top half of the tube is 6:1. The degree of floor curvature is 100%. The filled colored contours show the Hoek-Brown ratio described in Equation 5; compressive failure is present where this ratio is at least one, shown in white on the color scale. Shaded red designates where the rock is in tension. The tube shape is outlined with a black line. Circles plotted on the tube outline show tensile (red) and compressive (blue) failure. 68 Figure 31: Model output for a modeled lunar tube, with a width of 5.5 km and a roof thickness of 950 m. The width-to-height ratio of the top half of the tube is 2:1 (half- circular). The degree of floor curvature is 33.3%. The filled colored contours show the Hoek-Brown ratio described in Equation 5; compressive failure is present where this ratio is at least one, shown in white on the color scale. Shaded red designates where the rock is in tension. The tube shape is outlined with a black line. Circles plotted on the tube outline show tensile (red) and compressive (blue) failure. 69 The maximum tensile floor stresses for the full suite of models run with different tube shapes are shown in Figure 32. The dashed black line indicates the minimum tensile stress required for tensile failure. The effects of both the width-to- height ratios and floor curvatures of the tubes on floor stresses are shown; these effects can be seen to apply over almost all of the modeled tubes. The trend of floor tensile stress increasing with width-to-height ratio is true for the significant majority of models run, and only does not apply to the extremes of the floor curvature ranges used, where the floor is very deep-dipping with very low stress values or is completely flat. The dimensions used for the first group of models with varying floor shapes were relatively large, near the upper limits of the ranges of tube width and roof thickness used in the main model suites. Keeping roof thickness unchanged but using a significantly smaller width, or vice versa, resulted in decreased floor stress magnitudes in the majority of cases. The observed trends in tube shape?s effects on stress also apply to the tubes with different dimensions, with the only exceptions being at the extreme ends of the floor curvature range. 70 Figure 32: Tube shape effects on maximum floor tensile stress. The x-axis shows floor curvature depth. The dashed black line shows the tensile strength of the material; any points below this line have failure. Colors designate different top-half aspect ratios used: Half-Circular (6:3 top-half width-to-height ratio) in red, Half- Elliptical (6:2 ratio; tube shape used for original models) in blue, and Flatter Half- Elliptical (6:1 ratio) in green. ?W? and ?RT? designate tube Width and Roof Thickness respectively. Diagrams show color-coded visualizations of some of the tube shapes used. 71 3.4. Lunar Observation Results Multiple lunar sinuous rilles were examined in LROC imagery; overall, clear patterns of cracking running parallel to sinuous rilles or candidate tube sections could not be found. Some sinuous rilles were observed to have nearby linear patterns of features that may be cracks or bulges, but these did not consistently follow the path of the rille and were not symmetric across the rille, and could have been a result of other processes. Apparent linear cracking patterns were also observed in some cases inside the troughs of rilles, close to the upper edges; however, these were not considered to be possible tube-related cracks, as they were not a distance away on flatter ground and may have been caused by material collapsing and sliding down the rille slopes. 72 Figure 33: Examples of features noted on or near lunar sinuous rilles. For each figure two versions are shown; one with a graticule and one without to avoid obscuring features. All panels are from regions of Rima Sharp. a. and b. show a region where apparent linear features follow part of the rille path on one side. c. and d. show possible linear cracking patterns inside the rille trough and on the rille edge. Imagery from LROC. Although cracking was not observed in most cases, one case showed possible strong cracking similar to the predicted surface features. The southernmost section of Rima Mairan in Oceanus Procellarum, shown in Figure 34, ends in two large, aligned depressions, the arrangement of which follows the path of the rille. These features 73 may be collapses into a concealed lava tube that the rille extends into underground. Visible between these depressions are linear features that resemble cracks, following the path that the proposed lava tube would follow. These features may be the kind of surface tensile failure caused by a tube that the models predict. 74 Figure 34: Southern region of Rima Mairan. The rille ends in two aligned depressions following the path of the rille, connected by two lines of linear cracks. These could be indications of a partially collapsed lava tube continuing the path of the rille. Imagery from LROC. Patterns of surface bulges running parallel to tubes are also predicted by the models; no such patterns were clearly observed next to tube candidates or rilles in the LROC Quickmap analysis. However, it is possible that the combination of bulges to 75 either side of the tube and depression above the tube create a surface expression similar to a shallow sinuous rille above the tube, making these features harder to distinguish as tube surface expressions. The MHH is found in a depression that appears to be part of a sinuous rille; it is possible that this feature is instead the surface expression of a tube connected to the rille. However, the heights and depths of the predicted surface features may not be sufficient for this to be true, depending on the depth of the observed depression and the dimensions of the candidate tube. The region around the MHH was examined in multiple datasets in addition to LROC imagery (Figure 35). SELENE Multiband Imager data, showing different surface compositional abundances, showed patterns following some deeper parts of the rille, but not the shallower section associated with the candidate tube. This data also showed a possible FeO abundance (wt%) bright spot in the MHH itself. LRO Diviner silicate minerology data showed a possible pattern of roughness following the entire rille, including the shallower section in which the MHH is found. LOLA data showing slope at 100 m, roughness at 100 m, and roughness at 25 m were studied for the area; some data showed patterns following the path of the rille, but no patterns similar to the predicted surface features were seen. Data from Digital Elevation Maps made using LOLA and SELENE data (SLDEM2015) showing slope, azimuth, and topography all showed patterns following the rille, but also did not show the expected cracking. 76 77 Figure 35 (previous page): Imagery and data for the region of the lunar surface including the MHH and part of the surrounding rille. Figure a. is the LROC imagery of the area, including graticule and scale; the boundaries of the region shown and the scale are the same in b., c., d., and e. The insets in a. and c. are zoomed in on the MHH itself; both insets cover the same region with the same scale. The outline of the rille is shown in dashed yellow in a., and the MHH is circled in solid yellow; the same outlines are shown in the other figures, in either black or white depending on the image. b. shows SLDEM2015 slope data. Clear patterns follow the rille path, especially along the western section, and the MHH itself, but no clear indications of linear cracking (such as cracking causing broken blocks to change position and affect the slope) are visible to either side of the rille. c. shows SELENE Multiband Imager FeO abundance data. Patterns in the data can be seen on the left edge following the path of a deeper part of the rille. A bright spot is visible in or near the MHH. d. shows LRO Diviner silicate minerology data. A pattern of roughness in the data follows the entire rille, distinct from the smooth North-South linear patterns seen outside the rille. This roughness includes the shallower section in which the MHH is found, which was not followed by the patterns in some other datasets. e. shows Diviner rock abundance data, the fraction of rocks as opposed to fine soil present. This data is based on comparison of daytime and nighttime temperatures. No clear pattern is seen following the rille. 78 4. Discussion 4.1. Failure on Tube Internal Surfaces and Usable Floor Extents The model results concerning internal failure on tube internal surfaces are potentially applicable to planning future lunar exploration and use of lunar lava tubes. In the modeled tubes, a large portion of the tube floors are expected to be covered by cracks or debris from internal failure, leaving a small region of usable floor near the tube edges. This conclusion could be highly relevant to any future missions to explore or build inside a lava tube, to inform mission planning and predict the limits on what size of structure could be constructed. The upper limit of usable floor found here, 500 meters, could affect general planning of what kind of structures could be built in a tube. Some schematics of possible bases inside lunar tubes include large inflatable structures, covering most of the tube floor (Daga et al., 1992). As most of the floor is expected to be covered in large cracks with potentially sharp edges, the use of large inflatable structures may be impossible, or may need to be significantly adapted to accommodate these hazards. Such plans could be changed to cover only the usable regions and not the damaged floor, or to include specifically reinforced material to prevent tearing or puncturing. Knowing that the usable floor in a tube is expected to be near the tube edges, and the floor in the center of the tube is likely to be highly cracked, could also inform plans of how to enter a tube through a hole in its roof. If an exploration vehicle or manned mission entered a tube through a relatively narrow hole that covers only a small portion of the tube roof near the center, they would likely find themselves in a 79 heavily fractured area that would need to be crossed to reach usable floor. Planning for this circumstance could be important for such missions. The stress and failure patterns found in this work generally show agreement with the model examples given in prior works. Although Blair et al. (2017) only analyzed stresses and computed failure through the tube roof above the tube apex, they gave some example figures showing horizontal stresses through the rock; these seem to agree with results found in this research. Theinat et al. (2020) similarly focused on whether or not tubes were stable overall, but did give an example diagram showing where tensile stress was present in two of their models, and where tensile failure would occur depending on the tensile strength of the rock in these tubes. They observed patches of tensile failure at the surface, tube roof apex, and tube floor, which agree with the results presented here, and at the tube wall where it meets the tube floor, which does not agree. This may be related to the shape of the tube corners, which were rounded in this research but not in Theinat et al.?s (2020). The authors did not analyze in detail the extents of these failure regions. The tensile roof failure found in these modeled tubes agrees with a stability pattern seen by Blair et al. (2017). The authors found an overall trend of wider tubes requiring thicker roofs to remain stable, but found an exception to this trend in some tubes with an initial Poisson stress state, as used here. Some wide tubes with this stress state were stable with small roof thicknesses but became only quasi-stable with roof thicknesses of 200 meters or more due to tensile failure at the apex of the tube roofs. Examples included the 3.5-kilometer-wide tubes, which had tensile roof apex failure in relatively deep cases (200 meter roof thickness), no roof failure in some 80 shallower tubes (50 meter roof thickness), and compressive surface failure in very shallow tubes (20 meter roof thickness). This behavior agrees with the results of this work, where tensile roof failure is present spreading outwards from the roof apex in most tubes, but does not affect most shallow (50 meter roof thickness) tubes. This behavior was not shown in Theinat et al.?s (2020) results, possibly due to the different criterion for failure used. Blair et al. (2017) also note that failure starting from the apex of the roof is observed in terrestrial tubes. A small number of models run for this work with thin roofs also had compressive surface failure above the center of the tube, similar to very shallow tubes used by Blair et al. (2017). Model stability results found here could be used to inform decisions of which tubes to explore. Shallower modeled tubes have more usable floor than tubes with thicker roofs. These tubes may therefore be prioritized as more valuable tubes to use. Shallower tubes would also have thinner rock layers protecting everything inside from surface hazards. However, a radiation safety analysis of lunar lava tubes predicted that a roof thickness on the order of one to two meters would provide acceptable shielding and a roof thickness of six meters would effectively block all hazardous radiation, assuming an additional regolith layer thickness of five meters (De Angelis et al., 2002). Tube roof thicknesses of a few tens of meters are expected to also shield effectively against meteorite impacts and provide a near-constant temperature environment (H?rz, 1985). Even the shallowest of the lunar lava tubes modeled here have roof thicknesses significantly thicker than these minimum values, and could therefore provide effective shielding to their contents. Roof thickness may be observable through imaging of a tube hole at an angle, as for the MHH (Robinson 81 et al., 2012). If a tube?s roof thickness could be determined from the surface, from angled imaging, ground-penetrating radar, or a different source, shallower tubes could be located and this conclusion could be used to prioritize using them. The conclusions stated above are generally applicable to lunar tubes. Additionally, if the dimensions of a specific target tube could be found, these models could be adapted to use these dimensions to predict more precisely the internal conditions of that specific tube. 4.2. Surface Features The surface feature relationships observed in these models may be applied to the study and exploration of real tubes on the Moon. The predicted relationships between surface feature locations and tube widths provide a new way to determine tube widths from the surface, while the constant ratio of these distances may provide a method to locate an otherwise unobservable tube. If a lava tube is located under the lunar surface, possibly by observations of a hole chain, the path of a partly covered sinuous rille, or a topographic dip, and is accompanied by the surface features that these models predict, the relationships found here could be used to deduce the tube?s approximate width. The distance from the tube center to the peak of the elastic surface bulge, ??, depends on tube width and roof thickness by ?? = (0.71 ? 0.0055)? + (1.77 ? 0.033)?? + (?0.057 ? 0.027) (6) 82 where ?is the tube width in km and ?? is the tube roof thickness in km. The distance from the tube center to the location of maximum tensile stress, ??, depends on tube width and roof thickness by ?? = (0.43 ? 0.023)? + (1.09 ? 0.14)?? + (?0.11 ? 0.12) (7) where ?is the tube width in km and ?? is the tube roof thickness in km. Uncertainties in both equations are 95% confidence bounds, given by the MATLAB? ?fit? function. As failure is expected to occur at the maximum tensile stress location, this feature may be observable as an area of cracking. If observations of the area around a tube were made by satellite, it is likely that holes or other indications that a tube is present could be observed, but the width and depth of the tube could not be determined. High-resolution data or observations from rovers or astronauts may be able to detect cracking patterns and possibly a topographic bulge near the tube. If either cracks or a bulge running parallel to a tube were observed on the lunar surface, the relationships described above could be used to determine the width of the tube. This method only requires imagery of the area to be used, and does not necessitate astronauts placing and using equipment on the surface; this is an advantage over other methods such as GPR scans conducted on the surface (e.g. Esmaeili et al., 2020). The proposed method might also be usable even if the tube was too deep for such other geophysical methods to be able to image it. The dependence of the surface distances on roof thickness as well as width could be used to help determine tube characteristics. If the roof thickness of a tube could also be found from the surface, such as by angled imagery, this could be combined with the observed surface features 83 to decrease the uncertainty of the result and determine the tube width more precisely. Being able to determine tube width from the surface could be of great use. Without such methods, tube dimensions could only be found by descending into a tube through a hole, which could be a difficult and resource-intensive process. Surface features may be observed without the tube itself being visible; in this case, the patterns dictating surface feature distances could be applied in a different way. If both a topographic bulge and a line of cracking were found running parallel to a suspected tube, but no holes or other direct indicators of the tube?s location were present, the ratio of ??/?? could be used to deduce the tube?s location. This ratio was found to be constant regardless of tube dimensions, with values of 1.73 ? 0.10 for lunar tubes and 1.80 ? 0.14 for terrestrial ones. The distance between the two surface features could be used with this ratio to calculate the distance from either surface feature to the tube center, and therefore locate the concealed tube itself. The deduced distances could then also be applied to the relationships discussed earlier to find the width of the concealed tube. A potential caveat with this method is that the elastic surface bulge may be very difficult or impossible to detect. The amplitude of the bulge is expected to be very shallow, on the scale of approximately a meter, and could be lost amid the other, unrelated surface features present and the roughness of the volcanic field. However, on a flat and undisturbed area, this bulge may still be visible. Even if the elastic bulge is not observable, the location of the line of tensile cracking alongside the tube can be used by itself to draw conclusions about the tube width. 84 4.3. Effects of Different Tube Shapes The effects of tube shape on floor stability described in Section 3.3. indicate that the maximum tensile stress on a tube?s floor depends on multiple aspects of the tube?s geometry. This stress generally decreases with increasing depth of floor curvature, decreasing width-to-height ratio, and decreasing width and depth. The patterns in tube shape?s effect on floor stress could be a relevant consideration when choosing which tubes to explore and utilize for construction. Tube shape could be factored into choices of which tube to use if it could be determined from the surface, possibly by methods such as ground-penetrating radar, angled imaging at a tube hole, or a LIDAR (LIght Detection And Ranging) scanner being lowered into the tube. Choosing a tube with a curving floor could help to ensure failure is not present, but could also pose difficulties for tube use that would not be present with a flat-floored tube; for example, a sloping floor could be more difficult to build or walk on than a flat surface. Additionally, fallen debris could slide or roll from its original position on a curving floor, possibly making the locations of debris fields more difficult to predict and causing additional floor damage. Both the benefits and drawbacks of tubes with curving floors should be considered when selecting a target tube. Choosing tubes with higher roofs compared to their widths could also help avoid encountering floor failure, and may be a valuable consideration to prioritize; this aspect of tube shape could be determined from the surface if tube height could be determined from images of a hole and tube width could be determined from a surface feature relationship. 85 Tube shapes were also altered in Modiriasari et al. (2018)?s models, to use different half-ellipse width-to-height ratios. Floor curvatures were not added. The authors calculated the overburden yield, or percent of the thickness of a tube?s roof above the tube center that undergoes any kind of failure, for each tube. The authors found that decreasing the width-to-height ratio of the tubes decreased the overburden yield in each case, and caused the tubes to be more stable. Although Modiriasari et al. (2018) studied failure through the thickness of the center of the roof and this work studies stress and failure on the tube floor, the trend of stability increasing with taller, narrower tubes is the same. Some tube shapes were shared, while others were not; only Modiriasari et al. (2018) used half-elliptical tubes with a 4:3 width-to-height ratio, and only this work used such tubes with a 6:1 ratio. The results from these models indicate that this trend continues for these tube shapes. These trends may be present because a wider floor or roof relative to the tube can flex upwards or downwards respectively more than for a narrower tube, allowing larger stresses to develop. 4.4. Lunar Surface Observations The tube surface features predicted by the models were not clearly observed in most sinuous rilles studied. An exception to this is Rima Mairan; features that appear to be linear surface cracking are present between the two depressions at the endpoint of the rille, where a partially collapsed tube could be present (Figure 34). These cracks would be consistent with the predicted cracking for a tube with a width of 86 approximately 200 m. It is possible that the rille transitions into a tube in this location, and these cracks are the tube?s surface expression. This would support the model predictions regarding what surface features should develop, and would allow the approximate width of this candidate tube to be determined. In addition to the surface cracking, the existence of a tube in this location would also explain the aligned orientations and shapes of the two depressions. This rille was the only one found to have clear linear cracking over a tube candidate section, and the only one to have studied that had such large aligned depressions. It is possible that the changes in stress state caused by these collapses made the surface cracks expand and become more visible, or that the existence of cracking caused these collapses to occur. The study of GRAIL gravity data by Chappaz et al. (2016) revealed a mass deficit indicative of a strong candidate tube in this area, south of the end of Rima Mairan. This deficit is not in the same exact location or orientation as the aligned depressions, and is significantly wider than the 200 m predicted tube width. However, this could be an indication of another, larger tube in the area, suggesting that this area has the right geophysical characteristics to promote lava tube formation. This would strengthen the interpretation that Rima Mairan extends into a tube. The predicted 200 m wide tube is not shown in Chappaz et al.?s (2016) results; however, since this study cannot detect features less than 1 km in width, this feature would not be expected to be detected if it was present. Clear patterns of surface bulges running parallel to candidate tubes were not observed; however, as the predicted surface bulges are in many cases very shallow, they could be present but not observable due to being obscured by other nearby 87 topographic features. Additionally, it is possible that surface features that appeared similar to sinuous rilles were instead the surface dips and bulges created by large underlying tubes. These surface depressions would not have the sharp edges observed in many rilles, but could resemble rilles in other aspects. The depression depths predicted by the models are on the scale of tens of meters in many cases, and the depths of sinuous rilles can range from 4.8 m to 534 m, with a median depth of 49 m, based on a survey of over 200 rilles by Hurwitz et al. (2013). Therefore, although some depressions are too large to be created by any of the modeled tubes, some rille- like features could potentially be caused by surface depressions above tubes. Although bulges and cracks were not clearly observed in many locations, these features could still be present. The buildup of regolith and degradation of surface features over time could have rendered existing features impossible to observe using remote measurements. Observations from the lunar surface by future manned or rover missions may be able to detect these features, allowing the relationships found to be applied. Additionally, the present study has only covered a subset of all lunar sinuous rilles, and has not used all available datasets on each one. An in-depth study using datasets beyond the LROC imagery was only conducted on the MHH region. Using other data sets on further rilles may reveal undiscovered features, such as if LOLA roughness data revealed evidence of cracks partially covered by regolith that would otherwise be invisible. 88 4.5. Testing Robustness of Conclusions with Different Model Setups Several aspects of the model setup were varied on a representative model to determine if the results were significantly affected by them or if they were robust regardless of these changes. These aspects of model setup were the size of the block containing the tubes, the fineness of the mesh used in COMSOL Multiphysics?, and the fineness of the grid used to extract data from the models and analyze it in MATLAB?. One representative model was chosen on which to test the effects of these changes. The tube dimensions chosen for this representative model were a roof thickness of 0.65 km and a width of 3.5 km; this model showed patterns of internal failure similar to those seen in the majority of models, and was in approximately the center of the ranges of widths and roof thicknesses used. This model was also chosen to allow unusual overlapping points in the plot of fractional extents of floor tensile failure (Figure 21) to be studied. There are multiple instances in Figure 21 where points are found at the same locations. This model and its counterpart with the same width but a 0.95 km roof thickness are an example of such overlapping points. When testing the effects of a finer mesh, both this representative case and its deeper counterpart were modeled, allowing this unusually precise agreement to be investigated (see Section 4.5.2). 4.5.1. Block Size In the main model suites, a block with width and depth of 20 km was used in all cases. To determine if stress effects from the edges of the block were affecting the stress fields around the tubes, test models were run with much larger block 89 dimensions. Each block is a square, with equal width and depth. Prior research modeling lunar lava tubes, such as Blair et al. (2017) and Theinat et al. (2020), used block dimensions of 20 times the widths of the tubes, to keep block edges far enough away from the tubes to not influence the results. Based on these choices, the larger block size for this research was chosen as 140 km, which is over 20 times the largest tube width used. In the models with the larger blocks, the mesh setup parameters were adjusted to ensure that the mesh around the tube was kept as fine as in the original models, and not automatically changed to be coarser when applied to the larger block. Apart from the changes to block dimensions and mesh parameter adjustments, the models were kept unchanged from the originals. Using a large-block model as opposed to the original models caused little change to the results. The patterns of stress and failure shown in the model output was very similar to the original results (Figure 36). The extent of usable floor found in the representative tube changed very little; the change was less than one meter, which was significantly smaller than what was considered the precision of the models. This change was therefore considered negligible. 90 Figure 36: Comparison of model results for (a.) original model setup (20 km block dimensions) and (b.) model setup using a larger block (140 km block dimensions). The roof thickness used was 0.65 km, and the tube width used was 3.5 km. Interpretations of contours, colors, and red and blue circles are the same as in Figure 19 and similar figures. 4.5.2. Mesh All models have a fine mesh covering the top section of the modeled block, which contains the tube and surrounding area. To define this fine mesh in the main model suites, the finest of the available element size settings predefined by COMSOL 91 Multiphysics? was used; this predefined setting is referred to in the software as ?Extremely Fine?. A finer mesh was used in test models to determine if this had a significant influence on the results. As no finer predefined mesh setting was provided by COMSOL Multiphysics?, a custom mesh was created. This custom mesh used the same element size parameters as the ?Extremely Fine? mesh with the exception of the maximum allowed element size, which was decreased from 200 m to 20 m. Changing from the ?Extremely Fine? mesh to the finer custom mesh caused slight changes, but no significant effects on the results. The stress and failure plots for each tube tested are very similar for the original and finer mesh (Figure 37); using the finer mesh increases or decreases the number of compressive failure points near the tube edges slightly, but does not change the overall sizes or distributions of failure regions. The usable floor extent for the representative tube changes by less than one meter; similarly to the results of testing different block sizes, this change is considered negligible. Using a finer mesh caused the floor tensile failure percentages for the representative tube and its deeper counterpart, which had previously been identical, to separate. The two data points had previously occupied the same point in the tensile floor failure fractional extent plot (Figure 21). It was hypothesized that these tubes had exactly the same fractional floor failure extents because of mesh effects, and a finer mesh would show a slight difference between these points. Models of these tubes with the custom finer mesh gave results that were similar but not identical to one another, supporting this interpretation. The new percentages were close to the 92 original values; the points were separated but significant changes to the results were not introduced. Figure 37: Comparison of model results for (a.) original model setup and (b.) model setup using a finer COMSOL Multiphysics? mesh. The roof thickness used was 0.65 km, and the tube width used was 3.5 km. Interpretations of contours, colors, and red and blue circles are the same as in Figure 19 and similar figures. 93 4.5.3. Grid The grid used by MATLAB? to sample points across the models was constructed using set parameters for the numbers of horizontal and vertical increments. These parameters are defined in the code and can be altered for different model sets. Test models were run using five times the numbers of horizontal and vertical increments for these grids than were used in the main model suites, to determine if this change influenced the results significantly. Using a finer grid caused more significant effects than changing the block size or mesh parameters. The extent of compressive failure at the tube corner increased when a finer grid was used (Figure 38). This may be a result of the finer grid sampling points near the tube edge that were under enough stress to fail but were not sampled in previous models. However, although using a finer grid gave a larger extent of compressive failure, the change to this failure region is still relatively small compared to the overall patterns seen and does not change the trends that were observed, such as relatively little usable floor being present in each tube and these usable floor regions being found near tube edges. The extent of tensile floor failure was not affected significantly by the grid change, with a difference between test and original models of approximately one meter. 94 Figure 38: Comparison of model results for (a.) original model setup and (b.) model setup using a finer MATLAB? sampling grid. The roof thickness used was 0.65 km, and the tube width used was 3.5 km. Interpretations of contours, colors, and red and blue circles are the same as in Figure 19 and similar figures. Changes in block size or COMSOL Multiphysics? mesh caused very small effects on the failure patterns inside the modeled representative tube. Using a different MATLAB? sampling grid caused larger effects on failure, but did not change the overall patterns observed. Therefore, patterns of failure seen inside the 95 tubes are considered to be robust overall with regards to the changes in model setup applied. 96 5. Conclusions We found that lunar lava tubes are expected to have extensive failure on their internal surfaces (tensile floor failure, tensile roof failure, and compressive wall failure). This failure leaves only relatively small regions of contiguous usable floor near the tube walls, with a maximum extent of 470 m. This could be an important consideration for planning future missions to explore and utilize lunar tubes. Tensile floor failure is the most extensive type, and is affected by the shape of tube used. Tubes with smaller width-to-height ratios (more circular) and more curving floors generally have less tensile floor stress than other tubes. Tubes are also predicted to cause distinctive patterns of surface features, specifically a line of tensile cracking and a topographic bulge running parallel to the tube at predictable distances away from it. These distances depend linearly on tube width, and the ratio between them is constant; if observed, these features could therefore be used to deduce tube width from the surface or locate a concealed tube. An examination of locations where lava tubes may be present in various datasets and imagery types showed apparent linear cracking, of the type predicted by the models, on the southern end of Rima Mairan; this location also has two aligned possible collapse features, and seems likely to be a lava tube extension of the rille. These features were not clearly present in other locations studied, but may be visible in different data sets, different locations, or by future manned or rover surface observations. 97 Bibliography Angelis, D. G., Wilson, J. W., Clowdsley, M. S., Nealy, J. E., Humes, D. H., & Clem, J. M. (2002). Lunar lava tube radiation safety analysis. Journal of radiation research, 43(Suppl), S41-S45. https://doi.org/10.1269/jrr.43.S41 Atkinson, B. K. (1984). Subcritical crack growth in geological materials. Journal of Geophysical Research: Solid Earth, 89(B6), 4077-4114. https://doi.org/10.1029/JB089iB06p04077 Blair, D. M., Chappaz, L., Sood, R., Milbury, C., Bobet, A., Melosh, H. J., Howell, K. C., & Freed, A. M. (2017). The structural stability of lunar lava tubes. Icarus, 282, 47-55. https://doi.org/10.1016/j.icarus.2016.10.008 Borghi, A., Renard, P., Fournier, L., & Negro, F. (2015). Stochastic fracture generation accounting for the stratification orientation in a folded environment based on an implicit geological model. Engineering Geology, 187, 135-142. https://doi.org/10.1016/j.enggeo.2014.12.019 Calvari, S., & Pinkerton, H. (1998). Formation of lava tubes and extensive flow field during the 1991?1993 eruption of Mount Etna. Journal of Geophysical Research: Solid Earth, 103(B11), 27291-27301. https://doi.org/10.1029/97JB03388 Calvari, S., & Pinkerton, H. (1999). Lava tube morphology on Etna and evidence for lava flow emplacement mechanisms. Journal of Volcanology and Geothermal Research, 90(3-4), 263-280. https://doi.org/10.1016/S0377-0273(99)00024-4 Chappaz, L., Sood, R., Melosh, H. J., Howell, K. C., Blair, D. M., Milbury, C., & Zuber, M. T. (2017). Evidence of large empty lava tubes on the Moon using GRAIL gravity. Geophysical Research Letters, 44(1), 105-112. https://doi.org/10.1002/2016GL071588 Coombs, C. R., & Hawke, B. R. (1992, September). A search for intact lava tubes on the Moon: Possible lunar base habitats. In NASA CONFERENCE PUBLICATION (pp. 219-219). NASA. Daga, A. W., Daga, M. A., & Wendel, W. R. (1992, September). Evolving concepts of lunar architecture: The potential of subselene development. In NASA CONFERENCE PUBLICATION (pp. 281-281). NASA. Eberhardt, E. (2012). The hoek?brown failure criterion. In The ISRM Suggested Methods for Rock Characterization, Testing and Monitoring: 2007-2014 (pp. 233- 240). Springer, Cham. DOI: 10.1007/978-3-319-07713-0_20 Eppes, M. C., & Keanini, R. (2017). Mechanical weathering and rock erosion by climate?dependent subcritical cracking. Reviews of Geophysics, 55(2), 470-508. https://doi.org/10.1002/2017RG000557 98 Esmaeili, S., Kruse, S., Jazayeri, S., Whelley, P., Bell, E., Richardson, J., Garry, W. B., & Young, K. (2020). Resolution of Lava Tubes with ground penetrating radar: The TubeX project. Journal of Geophysical Research: Planets, 125(5), e2019JE006138. https://doi.org/10.1029/2019JE006138 Greeley, R. (1971). Lava tubes and channels in the lunar Marius Hills. The Moon, 3(3), 289-314. Greeley, R., Fagents, S. A., Harris, R. S., Kadel, S. D., Williams, D. A., & Guest, J. E. (1998). Erosion by flowing lava: Field evidence. Journal of Geophysical Research: Solid Earth, 103(B11), 27325-27345. https://doi.org/10.1029/97JB03543 Haruyama, J., Hioki, K., Shirao, M., Morota, T., Hiesinger, H., van der Bogert, C. H., ... & Pieters, C. M. (2009). Possible lunar lava tube skylight observed by SELENE cameras. Geophysical Research Letters, 36(21). https://doi.org/10.1029/2009GL040635 Hiesinger, H., Head III, J. W., Wolf, U., Jaumann, R., & Neukum, G. (2010). Ages and stratigraphy of lunar mare basalts in Mare Frigoris and other nearside maria based on crater size?frequency distribution measurements. Journal of Geophysical Research: Planets, 115(E3). https://doi.org/10.1029/2009JE003380 Hoek, E., & Brown, E. T. (2019). The Hoek?Brown failure criterion and GSI?2018 edition. Journal of Rock Mechanics and Geotechnical Engineering, 11(3), 445-463. https://doi.org/10.1016/j.jrmge.2018.08.001 Horz, F. (1985). Lava tubes-potential shelters for habitats. In Lunar bases and space activities of the 21st century (pp. 405-412). 1985lbsa.conf..405H Huff, W.D., & Owen, L.A. (2015). Volcanic Landforms and Hazards, Reference Module in Earth Systems and Environmental Sciences, Elsevier, ISBN 9780124095489, https://doi.org/10.1016/B978-0-12-409548-9.09512-9 Hurwitz, D. M., Head, J. W., & Hiesinger, H. (2013). Lunar sinuous rilles: Distribution, characteristics, and implications for their origin. Planetary and Space Science, 79, 1-38. https://doi.org/10.1016/j.pss.2012.10.019 Kaku, T., Haruyama, J., Miyake, W., Kumamoto, A., Ishiyama, K., Nishibori, T., ... & Howell, K. C. (2017). Detection of intact lava tubes at Marius Hills on the Moon by SELENE (Kaguya) Lunar Radar Sounder. Geophysical Research Letters, 44(20), 10-155. https://doi.org/10.1002/2017GL074998 Kern, Z., & Thomas, S. (2014). Ice level changes from seasonal to decadal time- scales observed in lava tubes, Lava Beds National Monument, NE California, USA. Geografia Fisica e Dinamica Quaternaria, 37(2), 151-162. DOI 10.4461/GFDQ.2014.37.14 99 Kiefer, W. S., Macke, R. J., Britt, D. T., Irving, A. J., & Consolmagno, G. J. (2012). The density and porosity of lunar rocks. Geophysical Research Letters, 39(7). https://doi.org/10.1029/2012GL051319 L?veill?, R. J., & Datta, S. (2010). Lava tubes and basaltic caves as astrobiological targets on Earth and Mars: a review. Planetary and Space Science, 58(4), 592-598. https://doi.org/10.1016/j.pss.2009.06.004 Marinos, P., & Hoek, E. (2000, November). GSI: a geologically friendly tool for rock mass strength estimation. In ISRM international symposium. International Society for Rock Mechanics and Rock Engineering. Modiriasari, A., Theinat, A. K., Bobet, A., Melosh, H. J., Dyke, S. J., Ramirez, J., Maghareh, A., & Gomez, D. (2018, March). Size and Structural Stability Assessment of Lunar Lava Tubes. In Lunar and Planetary Science Conference (No. 2083, p. 2803). https://ui.adsabs.harvard.edu/#abs/2018LPI....49.2803M/abstract Paris, A. J., Davies, E. T., Tognetti, L., & Zahniser, C. (2019). PROSPECTIVE LAVA TUBES AT HELLAS PLANITIA: LEVERAGING VOLCANIC FEATURES ON MARS TO PROVIDE CREWED MISSIONS PROTECTION FROM RADIATION. Washington Academy of Sciences. Journal of the Washington Academy of Sciences, 105(3), 13-36. 10.13140/RG.2.2.23401.03689 Pipan, T. & Culver, D. C. (2019). Chapter 107 - Shallow subterranean habitats. White, W. B., Culver, D. C., & Pipan, T. (3rd Ed.). Encyclopedia of Caves (Third Edition, Pages 896-908). Academic Press. ISBN 9780128141243. https://doi.org/10.1016/B978-0-12-814124-3.00107-2 Roberts, C. E., & Gregg, T. K. (2019). Rima Marius, the Moon: Formation of lunar sinuous rilles by constructional and erosional processes. Icarus, 317, 682-688. https://doi.org/10.1016/j.icarus.2018.02.033 Robinson, M. S., Ashley, J. W., Boyd, A. K., Wagner, R. V., Speyerer, E. J., Hawke, B. R., Hiesinger, H., & Van Der Bogert, C. H. (2012). Confirmation of sublunarean voids and thin layering in mare deposits. Planetary and Space Science, 69(1), 18-27. https://doi.org/10.1016/j.pss.2012.05.008 Sakimoto, S. E. H., & Zuber, M. T. (1998). Flow and convective cooling in lava tubes. Journal of Geophysical Research: Solid Earth, 103(B11), 27465-27487. https://doi.org/10.1029/97JB03108 Sauro, F., Pozzobon, R., Massironi, M., De Berardinis, P., Santagata, T., & De Waele, J. (2020). Lava tubes on Earth, Moon and Mars: A review on their size and morphology revealed by comparative planetology. Earth-Science Reviews, 103288. https://doi.org/10.1016/j.earscirev.2020.103288 100 Schultz, R. (1995). Limits on strength and deformation properties of jointed basaltic rock masses. Rock Mechanics and Rock Engineering, 28(1), 1-15. https://doi.org/10.1007/BF01024770 Theinat, A. K., Modiriasari, A., Bobet, A., Melosh, J., Dyke, S., Ramirez, J., Maghareh, A., & Gomez, D. (2018). Geometry and structural stability of lunar lava tubes. In 2018 AIAA SPACE and Astronautics Forum and Exposition (p. 5185). https://doi.org/10.2514/6.2018-5185 Theinat, A. K., Modiriasari, A., Bobet, A., Melosh, H. J., Dyke, S. J., Ramirez, J., Maghareh, A., & Gomez, D. (2020). Lunar lava tubes: Morphology to structural stability. Icarus, 338, 113442. https://doi.org/10.1016/j.icarus.2019.113442 Torrese, P., Pozzobon, R., Rossi, A. P., Unnithan, V., Sauro, F., Borrmann, D., Lauterbach, H., & Santagata, T. (2021). Detection, imaging and analysis of lava tubes for planetary analogue studies using electric methods (ERT). Icarus, 357, 114244. https://doi.org/10.1016/j.icarus.2020.114244 Turcotte, D., & Schubert, G. (2014). Geodynamics (3rd ed.). Cambridge: Cambridge University Press. Valerio, A., Tallarico, A., & Dragoni, M. (2008). Mechanisms of formation of lava tubes. Journal of Geophysical Research: Solid Earth, 113(B8). https://doi.org/10.1029/2007JB005435 Waters, A. C., Donnelly-Nolan, J. M., & Rogers, B. W. (1990). Selected caves and lava-tube systems in and near Lava Beds National Monument, California. U.S. Geological Survey Bulletin 1673 Wei, J., Day Hike- N?huku (Thurston Lava Tube). National Park Service. Date accessed: May 5, 2021. https://www.nps.gov/havo/planyourvisit/day-hike- nahuku.htm Wentworth, C. K., & Macdonald, G. A. (1953). Structures and forms of basaltic rocks in Hawaii (No. 994). US Govt. Print. Off.,. https://doi.org/10.3133/b994 Williams, K. E., McKay, C. P., Toon, O. B., & Head, J. W. (2010). Do ice caves exist on Mars?. Icarus, 209(2), 358-368. https://doi.org/10.1016/j.icarus.2010.03.039 101