Radar model and conditions for the detection of the interior of lava tubes
We consider a SAR sensor operating from orbit. We assume that after focusing the system range and azimuth resolutions are better than the cavity width \(w\). The radar received power \({P}_{R}\) is directly proportional to the radar backscattering coefficient \({\sigma }^{0}\). The value of \({P}_{R}\) is given by24:
$${P}_{R}\left(\theta \right)=\frac{K{\sigma }^{0}(\theta )}{{R}^{4}}.$$
(1)
The system constant \(K\) in the equation depends on factors such as the transmitted power, the antenna gain, and other relevant parameters. The variable \(R\) denotes the slant range. The term clutter refers to all unwanted radar echoes that have the potential to obscure the desired signal. In the SAR image, the cavity interior (which is the signal of interest) near a collapse point can be detected if the radar echo power originating from the cavity is larger than the echo power arising from the surface (i.e., clutter) for the same image pixel. This is considering the cavity interior and surface scattering points at the same slant range that results in the superimposition of two radar echoes. Therefore, the cavity can be detected only when the radar cross section of the cavity interior \({\sigma }_{c}^{0}\) is greater than that of the surface \({\sigma }_{s}^{0}\) (i.e., \({\sigma }_{c}^{0}\left({\theta }_{c}\right) > {\sigma }_{s}^{0}({\theta }_{s})\)) for any given pixel. Here, \({\theta }_{c}\) and \({\theta }_{s}\) represent the local incidence angles inside the cavity and over the surface, respectively, for each pixel in the image.
Radar model for the determination of skylight and lava tubes geometric parameters
The lava tube interior is illuminated if the radar look angle \({\theta }_{L}\) satisfies the following condition21,24 (see Supplementary Fig. 5):
$${\theta }_{L}\le {{{\rm{atan}}}}\left(\frac{{D}_{1}}{{h}_{t}}\right)$$
(2)
where \({D}_{1}\) and \({h}_{t}\) are the pit diameter (range direction) and the cavity roof thickness, respectively. This implies that \({\theta }_{L}\) should be sufficiently low to allow the radar signal to propagate inside the lava tube, rather than being solely reflected by the cavity roof exposed by the skylight. The collapse depth \(h={h}_{t}+{h}_{c}\) is estimated as (see Supplementary Fig. 5):
$$h={h}_{t}+{h}_{c}={\Delta R}_{g,1}\sin {\theta }_{I}\cos {\theta }_{L}$$
(3)
where \({\Delta R}_{g,1}\) is the ground range interval of the radar shadow cast from the pit’s edge, \({h}_{c}\) is the cavity height measured from the top of the rock pile and \({\theta }_{I}\) is the radar incidence angle. The radar signal propagation \({P}_{1}\) inside the lava tube is equal to24:
$${P}_{1}=\left({\Delta R}_{g,2}-{\Delta R}_{g,3}\right)\sin {\theta }_{I}\sin {\theta }_{L}$$
(4)
where \({\Delta R}_{g,2}\) is the ground range extent of the radar anomaly originating from the lava tube interior measured from the relevant pit edge and \({\Delta R}_{g,3}\) is the ground range interval measured between the reflection from the pit’s edge and the first reflection from the roof (see Supplementary Fig. 5). The radar signal propagation inside the lava tube can be also measured with reference to the pit’s edge. In this case, the measurement is denoted as \({P}_{2}\) and it is estimated to be:
$${P}_{2}\le {\Delta R}_{g,2}\sin {\theta }_{I}\sin {\theta }_{L}.$$
(5)
The maximum depth \({H}_{t}\) based on radar signal propagation inside the cavity is equal to:
$${H}_{t}={h}_{c}+{h}_{{pile}}+{h}_{t}\ge {\Delta R}_{g,2}\sin {\theta }_{I}\cos {\theta }_{L}.$$
(6)
The roof thickness \({h}_{t}\) could be estimated in case of vertical roof wall as:
$${h}_{t}=\frac{{\Delta R}_{g,3}\sin {\theta }_{I}}{\cos {\theta }_{L}}.$$
(7)
However, the unknown inclination of the roof wall makes it impossible to reliably estimate its thickness through equations. Indeed, when the roof wall is not vertical, the previous equation could yield an overestimation of the real value. It is possible to obtain an approximation of the inclined roof thickness as:
$${h}_{t}={\Delta R}_{g,3}\sin {\theta }_{I}\cos {\theta }_{L}.$$
(8)
Nonetheless, in this case, the obtained value may represent an underestimate of the real one.
Therefore, due to the uncertainty about the actual geometry of the cavity roof, the thickness \({h}_{t}\) has been estimated numerically, starting from the knowledge of the slant range interval measurement \({\Delta R}_{g,3}\sin {\theta }_{I}\) as follows:
$${h}_{t}={{\arg }}{\min }_{|{y}_{p}|}\left\{|\Delta R\left({x}_{p},{y}_{p}\right)-{\Delta R}_{g,3}\sin {\theta }_{I}|\right\}$$
(9)
where the parametric slant interval measurement \(\Delta R\left({x}_{p},{y}_{p}\right)\) as function of the roof cartesian reflection point \(({x}_{p},{y}_{p})\) is equal to:
$$\Delta R\left({x}_{p},{y}_{p}\right)= \, {R}_{1}-{R}_{0}=\sqrt{{\left({h}_{s}-{y}_{p}\right)}^{2}+{\left({h}_{s}\tan \left({\theta }_{L}\right)-{x}_{p}\right)}^{2}}\\ -\sqrt{{\left({h}_{s}\right)}^{2}+{\left(-{h}_{s}\tan \left({\theta }_{L}\right)\right)}^{2}}$$
(10)
where \({h}_{s}\) is the spacecraft height.
The cavity height \({H}_{i}\) within the lava tube is equal to:
$${H}_{i}=({\Delta R}_{g,2}-{\Delta R}_{g,3})\sin {\theta }_{I}\cos {\theta }_{L}.$$
(11)
The cavity width \(w\) and the pit major and minor axes, denoted as \({D}_{x}\) and \({D}_{y}\), respectively, are directly measured on the image. The lava tube floor slope is estimated as:
$${\theta }_{s}={{{\rm{atan}}}}\left[\frac{{H}_{i}-{h}_{c}}{{P}_{1}}\right].$$
(12)
To take into account the Magellan pixel resolution, all the geometric measurements of the conduit have been quantized as follows:
$$\underline{y}={{{\rm{r}}}}{{{\rm{o}}}}{{{\rm{u}}}}{{{\rm{n}}}}{{{\rm{d}}}}\left(\frac{y}{75}\right)\cdot 75$$
(13)
where y is a generic cavity parameter (e.g., \(h\)). Accordingly, the measurement uncertainty is equal to ±75 m.
The collapse linear volume \(V\) is equal to2:
$$V=\frac{\pi }{2}h\,{D}_{y}.$$
(14)
The collapse asymmetry ratio \({{{\rm{AR}}}}\) is defined as2:
$${{{\rm{AR}}}}=\frac{{D}_{y}}{h}.$$
(15)
Analysis of the uncertainties in the determination of skylight and lava tubes geometric parameters from Magellan data
A radar measurement taken at a slant range resolution equal to \({R}_{s}={R}_{r}\sin ({\theta }_{i})\), where \({R}_{r}\) is the ground range resolution, has a root-mean-square error \(\delta R\) equal to21:
$${\delta R=R}_{s}/\sqrt{2{{{\rm{SNR}}}}}.$$
(16)
Given that, for Magellan’s system, \({R}_{s}=88\) m and \({{{\rm{SNR}}}}=8\) dB, a value of 24.77 m is obtained for \(\delta R\). Thus, we have that the uncertainty due to SNR is lower than the one given by the pixel resolution, which is equal to p = 75 m (\(\delta R < p\)). Knowing the values of the incidence angle, \({\theta }_{i}\), and look angle, \({\theta }_{l}\), corresponding to the analyzed SAR measurement, we can bound the errors for the different estimated quantities. The pit examined in this paper, namely pit A (see Fig. 3), was imaged by Magellan’s SAR system with an incidence angle of 42.4° and a look angle of 40°. Consequently, the following uncertainties can be evaluated for the estimated geometric parameters (see Methods and Supplementary Fig. 5). The estimated collapse depth uncertainty is evaluated as \(p\sin {\theta }_{i}\cos {\theta }_{l}=38.74{{{\rm{m}}}}\). The radar signal propagation uncertainty is equal to \(p\sin {\theta }_{i}s{{{\rm{in}}}}{\theta }_{l}=32.75{{{\rm{m}}}}\). The lava tube roof thickness (worst-case overestimation model) uncertainty is equal to \(p\sin {\theta }_{i}/\cos {\theta }_{i}=66.02{{{\rm{m}}}}\). The cave height uncertainty is computed as \(p\sin {\theta }_{i}\cos {\theta }_{l}=38.74{{{\rm{m}}}}\).
Magellan radar system specifications and image data format
Magellan is a SAR system operating in S-band (2.385 GHz). The transmitted peak power is 325 W and the swath width is 25 km (variable). The Full Resolution Radar Mosaicked images (FMAP) used in this study are full-resolution (75 m/pixel) global mosaics produced by the U.S. Geological Survey from Magellan Full-resolution Basic Image Data Record (F-BIDR) data. FMAP products are reprojected from the radar coordinates (delay-doppler) onto the surface of Venus by using information about the spacecraft location based on a priori predictions and earth-based tracking. For this study, the left-look FMAP mosaic has been employed. The radar resolution ranges from 120 m × 120 m (range × azimuth) to 280 m × 120 m, depending on the spacecraft altitude.
Sentinel 1 radar system specifications and image data format
The SAR images used in this paper were acquired by Sentinel 1. Sentinel 1 is a SAR system operating in C-band (5.405 GHz). The Sentinel 1 Ground Range Detected (GRD) products used in this work encompass focused SAR data that have undergone detection, multi-looking, and projection onto the ground range utilizing an Earth ellipsoid model. The ellipsoid projection of the GRD products is adjusted based on the terrain height. Ground range coordinates refer to the slant range coordinates that have been projected onto the Earth’s ellipsoid. The pixel values represent the detected magnitude. The resultant product exhibits an approximately square spatial resolution and square pixel spacing of 10 m × 10 m, with reduced speckle owing to the multi-look processing. The spatial resolution is about 20.5 m × 22.5 m in range and azimuth, respectively.
Capella radar system specifications and image data format
The VHR SAR spotlight image of Lanzarote of size 5 km × 5 km used in this paper was acquired with Capella Space X-band radar systems. The radar central frequency and bandwidth are 9.65 GHz and 500 MHz, respectively. Data has been acquired in HH polarization with a pulse repetition frequency of 10 KHz. After the geocoding and orthorectification, the ground range resolution can vary between 0.5 and 0.8 m depending on the look angle. The azimuth resolution is 0.5 m without sidelobe filtering. The standard look angle range for Capella Space X-Band data is between 25° and 40°. Given the small size (5 km × 5 km) of the considered image, the look angle does not appreciably vary across the imaged scene. No post-processing has been applied to the data. We exploited the GEO (Geocoded Terrain Corrected) data products, in which pixel values contain the radiometrically calibrated intensity in linear scale. The image is multi-looked nine times in the azimuth direction to enhance their radiometric resolution. GEO data products are geocoded and terrain-corrected using a Digital Elevation Model to improve the geolocation accuracy beyond what is achievable with only considering the ellipsoid.