Spatio-temporal analysis and volumetric characterization of interferometric synthetic aperture radar-observed deformation signatures related to underground and in situ leach mining

Abstract. The effect of uranium mining on ground deformation is a relatively unexplored area, especially in terms of surface subsidence related to subsurface ore removal. We use interferometric synthetic aperture radar and spatiotemporal techniques to characterize subsidence signals at the McArthur River underground mine in Canada and the Four Mile in situ leach mine in Australia. We enhance the signal-to-noise ratio of our datasets via time-series techniques and compare results from active periods with results during inactivity to establish a baseline for mining-related signals. We then relate observed surface subsidence to subsurface volumetric strain rates via a voxel parameterization and Bayesian, geostatistical inversion. We use priors on our volumetric strain rates to identify whether these rates are best attributed to ore removal or if additional factors are contributing to subsidence at these sites. We find that the subsidence at McArthur River is best explained by a combination of ore removal and thermal contraction resulting from ground freezing practices. Ore removal via solution extraction alone explains the subsidence at Four Mile, although the localized subsidence pattern and resulting strain rates suggest an intricate combination of sinks and sources in the field, possibly from injection and production well locations and the subsequent flow of solution.


Introduction
Subsidence signatures due to mining, especially that of copper, gold, and coal, have been a subject of literature. 1 However, there is room for expanding on such signatures related to uranium mining, especially in terms of monitoring.In particular, the effects of underground or in situ leach (ISL) mining are relatively unexplored.As opposed to open-pit mining, in which materials are removed directly from ground excavation, underground mining operations happen beneath the surface such that any surface deformation would be due to the physical removal of subsurface material. 2ISL mining is a particular type of underground mining that removes material by pumping a solution underground to dissolve mineral-concentrated ore, which is then pumped to the surface and treated. 2 This results in a minimally invasive method for mining.Satellite interferometric synthetic aperture radar (InSAR) is a powerful geodetic remote sensing technique that makes use of the phase difference between two synthetic aperture radar (SAR) images to determine whether surface deformation has occurred at a site.Such deformation can be due to natural hazards, such as earthquakes, landslides, or volcanic unrest. 3However, ground subsidence and/or uplift can also occur due to anthropogenic activities, including geothermal energy production and mining. 4,5Differencing the phase between two SAR images taken over the same location but at different times measures the change in radar line-of-sight (LOS) distance between the satellite sensor and the ground pixel between passes of the satellite, indicating whether surface deformation has occurred.The resolution of such displacement measurement is dependent on the half-wavelength of the satellite, allowing for centimeter-scale resolution in displacement along LOS of the radar.
Given the availability of repeat, global-coverage SAR data via the European Space Agency (ESA)'s Sentinel-1 (S1) constellation, we aim to determine whether C-band, Terrain Observation with Progressive Scans SAR (TOPSAR) data can discern subsidence at underground and ISL uranium mining sites and whether such detected deformation is reasonably attributed to ore removal activities.We focus on two mining sites: the McArthur River mine in Canada for underground mining (Fig. 1) and the Four Mile mine in Australia for ISL mining (Fig. 2).
Owned and operated by Cameco, the McArthur River uranium mine is located in Northern Saskatchewan, Canada, and houses the largest deposit of high-grade uranium in world.In conjunction with the nearby Key Lake mine, McArthur River produced roughly 250.6 million pounds of U 3 O 8 between 2000 and 2013. 6Although upkeep at the mine is still ongoing, operations at McArthur River were suspended in July 2018 due to low uranium prices.Thus we focus on the 2017 period for measuring signals of interest.Details of the layout of operations at McArthur River can be found in publicly available operations records 6 (see Figs. 3-7 and  4-10, especially).Mining occurs at an average depth of 550 m. 6 There were two major active mining zones in 2017: zone 1 (near shaft 1) and zone 4 (near shaft 3). 6Zone 2, which is near zone 1, was in development and had maintenance activities, such as ground freezing, occurring during this time.In these active zones, ore is mined and then mixed into a slurry of surrounding rock and water for transportation to the surface via the shafts.
The Four Mile uranium mine is located in the Frome Basin of South Australia, about 5 to 10 km northwest of the Beverley uranium mine.Currently owned and operated by Quasar, Four Mile commenced operations on June 28, 2015, and is still currently operational.The mine uses the ISL process to extract uranium; during ISL, a weakly acidic solution is pumped into the subsurface to dissolve uranium-containing ore.The dissolved ore solution is then extracted and treated for uranium, with the cleaned fluid being reinjected for continued use.We focus Fig. 3 General approach to deformation analysis per site in this study.
Fig. 2 Overhead image of Four Mile mine, East section, with locations of monitoring wells marked with black and white targets; 7 image is from Google Earth (CNES/Airbus, 2022).Spatial coverage is the same as that for the InSAR analysis.specifically on Four Mile East for this study (e.g., Fig. 2), the location of which is centered on a large deposit of uranium. 8We consider data from 2017, which was the last year that production was predominantly focused in this area of the mine. 9ur general approach is shown in Fig. 3.We first examine each candidate site to determine if surface deformation is detectable via InSAR.Given the extreme environmental conditions at both sites and the likelihood of any observable signal being small in magnitude due to the subsurface   nature of the activity, we enhance the signal-to-noise ratio (SNR) of our InSAR dataset by first performing time-series analysis; such analysis reduces noise in individual pairs by considering average deformation trends across all pairs over the timespan of the dataset.We compare our results with baseline time-series analyses at times/locations where no operations were taking place to confirm the correlation between signal presence and mining operations.We then perform spatial deformation analysis to further identify whether the observed signals are related to ore removal.Using a voxel parameterization and Bayesian geostatistical inversion, we test scenarios explaining our observed subsidence in terms of subsurface volume change.We analyze whether such subsidence is due primarily to ore removal or if there are other factors contributing to deformation at each site.
2 Materials and Methods

SAR/InSAR Data
We use data from ESA's S1 constellation. 10 S1 consists of two satellites carrying C-band wavelength sensors (∼5.5 cm); each satellite has a standard cadence of about 12 days.We use singlelook complex SAR data from the interferometric wide (IW) swath mode, which is a TOPSAR scanning mode that uses a sweeping azimuthal antenna pattern to achieve larger image swath areas at the cost of azimuthal resolution. 10In S1 IW mode, each SAR scene consists of three subswath images (roughly 5 m × 20 m). 10 In addition to the S1 data, SRTM 1-arc sec digital elevation models were obtained for topographic comparison. 11e use a bounding region of ½57.75°; 57.77° in latitude and ½−105.08°;−105.03° in longitude for our McArthur River analysis, which spatially correlates with ascending track T107 (frame 1163).To avoid processing data taken during suboptimal conditions (e.g., snow cover), we only consider data from March to October of 2017.This corresponds to 20 scenes in our dataset.Due to the relatively small size of our bounding box, we only utilize the second subswath of each scene.The unit vector ŝ ¼ ½−0.63; −0.11; 0.77 for this track describes the track's pointing vector from the ground to the satellite sensor.
We perform a similar selection of data at Four Mile using a bounding box of ½−30.1645°;−30.1202° in latitude and ½139.4734°;139.5467° in longitude.We consider data during the entire year of 2017.We identify 29 suitable scenes from descending track T60 (frame 4197).We utilize only the second subswath of these scenes as it completely encompasses our region of interest.The unit vector for this track is ŝ ¼ ½0.58; −0.16; 0.80.
For each site, we form a dataset of interferometric pairs using different combinations of scenes.To maximize coherence in our dataset, we impose a maximum orbital separation (the distance between the physical location of the satellite at the first epoch versus its location at the second epoch) of 120 m and a maximum temporal separation of 90 days.Pair formation and further down-selecting of pairs that contained poor coherence or atmospheric contamination resulted in a final time-series dataset of 50 pairs and 18 scenes spanning from March 1st to October 15th for McArthur River.We similarly arrive at a final time-series dataset of 64 pairs for Four Mile after down-selecting pairs during postprocessing and removing pairs that had poor coherence throughout.We process our interferometric pairs using GMTSAR, 12 an open-source InSAR processing software that uses Generic Mapping Tools 13 for visualization.We include ionospheric correction via the split-spectrum method, and pairs are converted from wrapped phase to unwrapped phase using the Snaphu method. 14aseline plots and the final pair datasets for McArthur River and Four Mile are shown in Appendix A.

Time-Series Analysis
To improve the SNR measured deformation signatures in our InSAR pairs, we perform timeseries analysis on each dataset prior to applying deformation modeling.This helps to eliminate transient noise in individual pairs while enhancing persistent deformation signals throughout the dataset.
We use GMTSAR's implementation of small-baseline subsets (SBAS) analysis to estimate the average deformation velocity from our pairs.SBAS performs pixel-wise time-series analysis on a stack of interferograms using weighted least squares to solve for average velocity per pixel. 15his process converts a stack of pair-wise displacement estimates to estimates of cumulative displacement at each scene (i.e., epoch) in the dataset.From this, an average velocity map is derived for the time period of the dataset.SBAS relies on maximizing the possible coherence among a network of interferograms by selecting pairs with small orbital and temporal baselines.We accordingly select our pairs using orbital and temporal separation benchmarks of 120 m and 90 days, respectively, as described above.

Spatial Deformation Modeling
Using the average deformation rate field derived from time-series analysis for each site, we perform spatial deformation analysis on these fields to characterize the subsurface volume change related to the observed surface deformation signature at McArthur River and Four Mile.We apply Bayesian geostatistical methods first introduced for geothermal activity; 16 we summarize this approach below.
This method parameterizes a deformation field as a gridded layer of voxels, where each voxel is representative of a volumetric sink (or source) with an associated volumetric strain rate.Linear, geostatistical inversion is used to solve for the best-fitting estimates of volumetric strain rate for each voxel, which is then converted to estimates of the volume change rate.The corresponding set of matrix equations to solve by linear inversion is where ρðmodÞ represents a vector of modeled LOS displacement rates and G represents the design matrix relating the displacement within each voxel to each observed value of the LOS displacement rate.We interpret m, a vector of volumetric strain rates, by assigning prior models for reasonable estimates of voxel-wise volumetric strain rates based on possible causes of volume change at the site (e.g., ore removal).A prior uncertainty for the volumetric strain rates is defined as part of these priors, which is then used to define a model covariance matrix Q using characteristic distance scales l x and l y , which describe spatial correlation lengths of the observed deformation signals.We solve for the best-fitting estimates of voxel-wise volumetric strain rates by minimizing the negative logarithm of the posterior probability distribution. 16

Parameterization and prior for McArthur River
We parameterize the 2017 deformation rate map from McArthur River as a gridded layer of ð100 mÞ 3 cubic voxels centered at 550 m depth according to the depth of mining operations. 6e use an exponential geostatistical model covariance matrix with characteristic distance scales of l y ¼ 1700 m along strike of 45°N and l x ¼ 500 m across strike of 45°N to account for the geometry of the mining and corresponding deformation signal.Sentinel-1 has been shown to have measurement uncertainty of roughly 1.1 mm∕year. 17We use this value to assign an inherent data uncertainty to our measurements after verification with our own GPS validation analysis (Appendix B).We assign priors to our volumetric strain rates.We consider mining (e.g., ore removal) as the sole contributor of deformation at the site.To assign a reasonable prior on volumetric strain rates associated with ore removal, we first estimate the volume change that we would expect due to mining operations alone.For McArthur River, the production of fully processed ore for 2017 is listed at roughly m ðMÞ u3o8 ¼ 7.6 Mkg. 6Given an ore density of ρ u3o8 ¼ 8380 kg∕m 3 , this equates to a volume change rate VðMÞ ore of E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 4 ; 1 7 6 where the negative value for the volume change rate is used to represent the removal process.However, we must also take into consideration the amount of host rock removed during the mining process.The slurry at McArthur River is cited to have roughly 50% solid material by weight, 18 and the average grade of ore for 2017 was roughly 8%. 6We therefore calculate the weight m ðMÞ rock of removed host rock: We know that the surrounding host rock of McArthur River is primarily composed of sandstone with density ρ s ¼ 2323 kg∕m 3 and quartzite with density ρ q ¼ 2711 kg∕m 3 . 6,19Although the subsurface is comprised of a mix of these minerals, we calculate end-member values for the volume change rate VðMÞ rock associated with the removal of host rock: We convert this range to a volumetric strain rate prior by dividing by the total volume of our voxel grid.With 690 voxels at ð100 mÞ 3 each, the total volume V ðMÞ T is 6.9 × 10 8 m 3 .
The volume change rates in Eq. ( 5) are converted to volumetric strain rates ϵ ðMÞ ¼ VðMÞ ∕V ðMÞ T ¼ ½5.2; 6.1 × 10 −5 year −1 .We find a suitable μ − 1σ prior for volumetric strain rates at McArthur River to be E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 7 ; 5 3 7 Although the uncertainty in our prior will be used as a weighting factor for our model covariance matrix, the mean of our prior will not be used when finding our best-fitting estimates.It will instead be used as a comparison with our posterior mean found through Bayesian, geostatistical inversion.

Parameterization and prior for Four Mile
We apply the above approach to the 2017 average deformation rate map from the Four Mile ISL mine.We use a grid depth of 120 m in accordance with mining procedures at the site 7 and a voxel size of ð50 mÞ 3 .Due to the localized, circular deformation patterns at Four Mile, we do not impose strict characteristic length scales for the geostatistical model covariance matrix.Instead, we assign characteristic length scales l x ¼ l y ¼ 10 m, recognizing that this favors independence among the voxel parameters.
ISL mining at Four Mile involves dissolving the ore in a solution (lixiviant) that is pumped into the ground.The sulfuric acid-based lixiviant is designed to leach the ore from the subsurface such that it can be recovered via a solution.The nature of ISL mining signifies that all of the solution pumped into the subsurface is recovered in a closed flow circuit, with the exception of a small amount of bleed such that a groundwater gradient exists in the system. 7Some additional elements may be leached as well, resulting in by-products, although the concentrations of such elements in the produced solution are typically small.Given the bleed percentage, the volume of ore produced, and estimates of additional leached mineral concentrations, we can calculate the volume change rate expected at the site due to mining processes.
According to records, 7 the average production flow rate at Four Mile for 2017 was 1.6 GL∕year, with production mostly attributed to the Four Mile East and Four Mile Northeast sections and at a bleed rate of 0.01% to 3% relative to said flow.
From this, we find the total volume change rate VðFÞ bleed associated with bleed for the given bleed percentage range: We then estimate the total rate of volume change VðFÞ ore associated with the removed ore.The total mass of produced ore in its final product form is m ðFÞ u3o8 ¼ 1.8 × 10 6 kg with a density ρ u3o8 ¼ 8380 kg∕m 3 . 20From this and the molar masses of both the final product ore M u3o8 ¼ 842 g∕mole and mined ore M u ¼ 238 g∕mole, we calculate the corresponding mass of mined ore: We now calculate the volume change rate VðFÞ ore of mined ore given its density ρ u ¼ 1900 kg∕m 3 : E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 4 ; 5 0 9

VðFÞ
ore ¼ m ðFÞ u ∕ρ u ¼ −7.9 × 10 2 m 3 ∕year: Finally, we estimate the volume change rate associated with by-product recovery.There are several commonly associated by-products of sulfuric acid ISL technology, which are shown in Table 1. 21Using the concentrations and densities listed for each element as well as the total production flow listed above, we find the total volume change rate associated with by-product recovery to be VðFÞ by-prod ¼ 21.9 m 3 ∕year.We find a range for the total volume change rate _ V ðFÞ associated with mining at Four Mile to then be E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 1 1 4 ; 4 0 5 VðFÞ by-prod ¼ ½−0.24; −4.81 × 10 4 m 3 ∕year: We note that this is technically a measure of reasonable volume change rates for all of Four Mile, whereas we are only considering Four Mile East.Given that production at Four Mile West was limited during 2017, we consider a range of 50% to 100% of the rates in Eq. ( 10) to be likely attributable to Four Mile East specifically.Now that we have a range of reasonable volume change rates, we convert this range to volumetric strain rates.There are 416 voxels in our parameterization, each with a volume of ð50 mÞ 3 .This equates to a total volume V ðFÞ T ¼ 4.2 × 10 6 m 3 for the deforming region.Thus the range of volume change rates in Eq. ( 10) is converted to volumetric strain rates as εðFÞ ¼ ½−0.46; −9.23 × 10 −4 yr −1 .We establish a wide μ AE 1σ prior for reasonable volumetric strain rates associated with mining: As mentioned above, only the uncertainty of this prior will be used during inversion; the mean will be used only as a comparison after our inversion is completed.To determine if the signal that we are measuring at McArthur River can be reasonably attributed to mining activities, we perform another time-series analysis during a time when the mine was no longer operating (2019).Figure 4(b) shows the resulting displacement rate measured during this period of no operations.Although some faint variations in background signal are apparent, most likely due to atmospheric effects and/or the wetland-like conditions at McArthur River, we see a distinct lack of deformation signal compared with Fig. 4(a).

Spatial deformation modeling and volume change analysis
We apply the voxel parameterization described above and use Bayesian, geostatistical inversion to estimate the optimal volumetric strain rates for each voxel (Fig. 5).We find the mean of the posterior volumetric strain rates attributed to ore removal to be μ ðMÞ ¼ εðMÞ ¼ 5.4 × 10 −6 year −1 , which is significantly outside our prior range of ½−5.6 AE 2.8 × 10 −5 year −1 for reasonable values of volumetric strain rates.This suggests that our prior, which interprets the deformation signal in terms of ore removal alone, does not accurately represent the subsurface conditions at McArthur River.This is also reflected when comparing the resulting posterior distribution of best-fitting strain rate estimates with the prior distribution for volumetric strain rates that we considered reasonable (Fig. 6).The largest individual parameter uncertainty is σ ðFÞ to be significantly deforming, we calculate the volume of the deforming region to be V ¼ 3.1 × 10 7 m 3 and the volume change rate of the voxels within the modeled chamber to be V ¼ ð−7.7 AE 0.1Þ × 10 4 m 3 .We calculate the dimensionless misfit of the model to our data using the square root of the reduced χ 2 -test statistic, 23,24 where a misfit value close to one indicates that the observed data is well represented by the model.We find that ffiffiffiffi ffi

Time-series analysis
We perform a similar analysis at the Four Mile ISL mine as was done for McArthur River.We correct for ionospheric effects during processing, although gradient effects still prevail in some pairs; thus we include additional atmospheric correction by common point stacking. 25We also apply a smoothing factor to weight against transient changes in topography 26 due to the arid conditions at Four Mile.Because we already down-selected to maximize the quality of our dataset, we do not weight temporally by coherence when applying SBAS.The results from the SBAS analysis showing average velocity at Four Mile over 2017 compared with the site location are shown in Fig. 7.We can see very clear subsidence centered on the mining site, with a majority of the deformation occurring between monitoring well locations and in correlation with the location of ore deposits. 8Given that the layout at Four Mile is such that monitoring wells surround injection and production wells, 7 we can infer that this subsidence is occurring near and/or around locations of injection and production activities.
To more accurately interpret the signal captured at Four Mile, we compare our results with results taken over an inactive well field belonging to Beverley Mine, a sister site located ∼8 km SE that ceased main operations in 2014 and was inactive in 2017.We choose to use this sister site at the same time period as a baseline as opposed to using a period of inactivity at Four Mile (as was done for our McArthur River comparison) due to Four Mile's continued operations and the lack of S1 data prior to the mine's opening.Using the same SAR dataset as for the period of operations at Four Mile, we form corresponding pairs over a section of Beverley mine.We then perform SBAS analysis using the same optimized pair subset as those from the analysis at Four Mile.The results are shown in Fig. 8 using the same color scale as Fig. 7.We see no apparent deformation at Beverley, especially where the production wells are located.

Spatial deformation modeling and volume change analysis
Applying the voxel parameterization described in Sec. 2, we find the best-fitting estimates of volumetric strain rates (and accordingly volume change rates) for each voxel in our gridded layer [Fig.9(a)].We find the mean of the posterior volumetric strain rates μ ðFÞ ¼ εðFÞ ¼ −1.8 × 10 −5 year −1 , which is slightly outside our prior μ − 1σ range of ½−4.7 AE 4.4 × 10 −4 year −1 for reasonable values of volumetric strain rates (Fig. 10).The largest individual uncertainty associated with our parameter estimates is σ ðFÞ ϵ ¼ 4.1 × 10 −4 year −1 .We find the volume of the modeled deforming region to be V ¼ 2.1 × 10 6 m 3 and the volume change rate of the voxels within that region to be V ¼ ð−4.8 AE 0.01Þ × 10 3 m 3 ∕year, which is within the range that we considered reasonable for volume change rates due to ISL mining at Four Mile.We find the misfit of the model to be ffiffiffiffi ffi χ 2 p ¼ 1.5.The corresponding deformation field modeled from the best-fitting estimates of volumetric strain rates and observed and residual deformation fields are shown in Figs.9(b)-9(d).

Discussion
Our time-series results indicate clear subsidence at both mining locations.Furthermore, comparisons with baseline time periods or locations when operations were halted suggest that the observed deformation at both sites may be attributed to the presence of mining activities.The signal location at the McArthur River underground mine is aligned along the surface with the known location of subsurface activity.The majority of the signal is oblong and along strike of the mining chamber beneath, which is consistent with the idea that this observed surface signal is related to subsurface activity at the site.The observed deformation at the Four Mile ISL mine is localized and centered within areas where we know monitoring wells are present.This suggests that the observed deformation is correlated to injection and production wells at the site, as we know that the layout of Four Mile is such that injection and production wells are surrounded by monitoring wells. 7,9o further analyze whether volume change associated with the observed amount of surface deformation could be reasonably attributed to ore removal, we additionally perform spatial deformation modeling using the results from our time-series analyses.We discuss our results in detail for each site below.

McArthur River
We find in our spatial analysis that our results under the prior assumption that observed deformation can be attributed to ore removal alone are unreasonable and unrealistic.The posterior mean estimate of volumetric strain rates lies significantly outside our prior range, and the uncertainty in the posterior estimates is larger than the uncertainty in our prior.Furthermore, the estimates correspond to an inflation scenario, which is in direct conflict with the subsidence event that we observed taking place at the site.This indicates that our prior assumption that the deformation at McArthur River can be solely described by ore removal is incorrect; there are additional factors contributing to the observed subsidence at the site, and we need to redefine a prior to incorporate such additional factors.After careful reconsideration of mining operations at McArthur River, we examine the possibility of ground freezing contributing to subsidence at the site via thermal contraction of the rock matrix.Ground freezing is done as a safety precaution to protect against harmful materials removed during the mining process; a more in-depth discussion as to how this procedure works is discussed in the literature. 6,27,28Generally speaking, ground freezing is accomplished by pumping cooled brine into an underground system (e.g., around mining zone locations) at temperatures between −25°C (248 K) to −35°C (238 K).The brine freezes the surrounding rock matrix into which it is pumped, essentially creating a protective barrier around harmful materials.Due to the temperature exchange that this operation causes, there is the possibility that the affected rock matrix may experience thermal contraction, in turn causing additional subsidence.
We examine the possibility that thermal contraction is contributing to volume change at McArthur River by repeating our spatial analysis with a prior updated to take into account additional effects from thermal contraction due to ground freezing.Strain rate due to thermal contraction of the rock matrix is characterized as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 1 1 7 ; 7 1 2 εT ¼ α T Ṫ; (12)   where α T represents the thermal volumetric expansion coefficient and _ T represents the temperature change rate at the site. 29Based on the literature, 30 we assume a reasonable mean estimate of Ṫ ¼ −27 K∕year and a mean estimate of α T ¼ 5.5 × 10 −6 1∕K.We assign a wide uncertainty of −27 K for _ T to allow for the possibility of either little to no temperature change of the surrounding rock matrix or even greater change.We assign an uncertainty of 3 × 10 −6 1∕K for α T .We find our corresponding prior for reasonable volumetric strain rates associated with thermal contraction using the product of two Gaussians. 31This leads to a prior of E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 1 1 7 ; 6 0 0 Our joint prior is defined as the sum of the two random variables characterizing the volumetric strain rate εðMÞ T related to thermal contraction and the volumetric strain rate εðMÞ m related to ore removal [e.g., Eq. ( 6)], which we find through the sum of two Gaussians: 31 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 4 ; 1 1 7 ; 5 3 3 Repeating our analysis using this prior, we find a posterior mean volumetric strain rate estimate of μ 11(a)].This value falls within the μ þ 1σ prior range defined in Eq. ( 14), and the corresponding maximum posterior parameter uncertainty is 1.4 × 10 −4 yr −1 , which is smaller than our prior uncertainty (Fig. 12).We find  22 the volume of the deforming region to be V ¼ 3.7 × 10 7 m 3 with a volume change rate of ð−1.1 AE 0.1Þ × 10 5 m 3 ∕year.This values lies within the range considered reasonable when converting the prior in Eq. ( 14) to the volume change rate.The misfit of the model is ffiffiffiffi ffi χ 2 p ¼ 1.8, which is improved from our original run.Our resulting modeled and residual deformation fields are shown in Figs.11(c) and 11(d), respectively.
To determine if considering thermal contraction due to ground freezing as an additional contributor to subsidence significantly improves the model fit, we run an F-test for model comparison. 31At an alpha level of 0.05, we test the null hypothesis that our joint prior does not provide a significantly better fit to the data than our first prior against the alternate hypothesis that our joint prior does provide a significantly better fit.We find F calc ¼ 1.39, which is greater than F α¼0.05 ¼ 1.02.Thus we reject the null hypothesis and conclude at 95% confidence that incorporating thermal contraction into our prior provides a significantly better posterior fit to the data.Thus we infer that thermal contraction due to ground freezing is also contributing a detectable level of subsidence to the deformation observed at McArthur River.
Although some residual signal exists along the edges of the deforming region in Fig. 11(d), the model does a suitable job of fitting most of the observed deformation.The best-fitting parameterization [Fig.11(a)] highlights three major areas of significantly deforming voxels.We note that these areas are similar in location to the mining zones that we know to be active as well as shaft locations.

Four Mile
The spatial analysis at McArthur River resulted in a posterior average volumetric strain rate that was slightly smaller than expected in magnitude.This could be due in part to Four Mile comprising three segments (Four Mile West, Four Mile East, and Four Mile Northeast), and although we only imaged Four Mile East, the data that we had to inform our priors came from Four Mile as a whole.Closer examination of the best-fitting volumetric strain rates in Fig. 9(a) suggest that there are localized areas of volumetric strain rate expansion, which may also contribute to the smaller average.The presence of such voxels further suggest the presence of such localized subsurface expansion, which may be attributed to injection at the site.
Although knowledge of well locations at the site is necessary for true confirmation, we examine the potential for some of the voxels to be experiencing positive volumetric strain rates.We rerun our analysis using a prior standard deviation of 1 × 10 −3 year −1 , which allows more flexibility in the range of best-fitting estimates, including considering positive volumetric strain rates associated with local volumetric expansion (e.g., injection) to be reasonable.With this prior, we find the mean of our posterior volumetric strain rates to be μ ðFÞ ¼ εðFÞ ¼ −1.4 × 10 −5 year −1 , which is less (in magnitude) than our previous estimate [Fig.13(a)].However, this estimate now falls within the μ − 1σ range that we use to define reasonable volumetric strain rates.Our posterior uncertainty is now σ ðFÞ ϵ ¼ 8.8 × 10 −4 year −1 , which is an order of magnitude smaller than our prior and suggests a more informed posterior fit (Fig. 14).The corresponding estimate of the volume change rate of the deforming region is V ¼ ð−8.6 AE 0.2Þ × 10 3 m 3 ∕year, which is still within the range that we considered reasonable for volume change rates due to ISL mining at Four Mile.We find the misfit of the model to be ffiffiffiffi ffi χ 2 p ¼ 1.2.To determine if our wider prior provides a better fit to our previously defined prior, we run an F-test for model comparison. 31At an alpha level of 0.05, we test the null hypothesis that our broader prior does not provide a significantly better fit to the data than our first prior against the alternate hypothesis that our broader prior does provide a significantly better fit.We find F calc ¼ 1.30, which is greater than F α¼0.05 ¼ 1.02.Thus we reject the null hypothesis and conclude at 95% confidence that broadening our prior provides a significantly better posterior fit to the data.The corresponding modeled and residual deformation fields are shown in Figs.13(b)-13(d).
Although we found that broadening the prior for volumetric strain rates provided a better fit to the data than our original prior at Four Mile, the spatial analyses resulting from both priors are unable to fully resolve the individualized patterns of deformation observed at the site [e.g., Fig. 13(d)].We attempted to address this by decreasing our voxel size to ð25 mÞ 3 ,  22 but such residual signatures still persisted in our results.The deformation signature may be better parameterized by a non-cubic voxel field, although such work is outside the scope of this study.
Both Four Mile spatial analyses found that regions of significantly deforming voxels seemed to be grouped in areas where wells may be located.However, further analysis with well locations is needed before drawing this conclusion.

Conclusions
We have identified subsidence signatures at both an underground and an ISL mining site using Sentinel-1 data.Using SBAS time-series techniques for signal-to-noise improvement in our dataset, we have detected subsidence that spatially correlates with locations of mining operations at each site, and comparison with baseline times/locations when no operations were taking place confirmed that our observed deformation signals also correlate temporally with mining operations.By employing Bayesian, geostatistical inversion in a voxel parameterization framework, we found the best-fitting estimates of the subsurface volume change rate corresponding to the observed surface deformation.Through established priors, we tested whether it is reasonable to attribute such observed subsidence entirely to ore removal at each site.We found that the subsidence observed at McArthur River is best explained through a combination of ore removal and thermal contraction resulting from ground freezing practices at the site.The subsidence at Four Mile ISL mine may be reasonably explained by ore and solution removal alone, although the localized pattern of subsidence and resulting best-fitting volumetric strain rate patterns suggest a more intricate pattern of sinks and sources potentially due to injection and production wells and the flow of solution between them.Further inference may be drawn with more knowledge of well layouts at the site.
6 Appendix A: Pair Datasets

Four Mile East S1 T60 2017 Baseline Plot and Pairs
The baseline plot for the epochs considered in our 2017 Four Mile analysis as well as the list of pairs in our final time-series analysis are shown in Fig. 16 and Table 3, respectively.

Appendix B: InSAR Measurement Validation with GPS
To verify the uncertainty in our InSAR-derived measurements, we perform a validation check with global positioning survey (GPS) measurements.There are no available GPS measurements at either site (McArthur River or Four Mile) during our time period of interest, so a direct    comparison with our dataset is not possible.Instead, we perform a validation analysis at a proxy site that is both actively deforming due to subsurface processes and that has GPS data available to measure such deformation to determine the uncertainty that we have in Sentinel-1 InSARobserved deformation due to underground processes.We choose Brady Hot Springs geothermal field in Nevada, a site that has been used in previous studies to validate InSAR measurements within an actively deforming region due to subsurface processes, 16,32 for such analysis.This site is an ideal choice due to the existence of a continuous GPS station (BRDY) in a non-deforming region as well as a continuous GPS station (BRD1) located in the middle of the subsiding region, courtesy of a 2016 field study led by the University of Wisconsin in support of a US Department of Energy's Geothermal Technologies Office project. 33Both stations are part of the MAGNET network with publicly accessible time series of relative position analyzed using standard procedures. 34,35Their relative proximity (∼3 km separation) allows deformation around both stations to be measured by the same interferogram (Fig. 16).We conduct our validation similar to previous works, 16,32 using data from Sentinel-1's track T144 (ŝ ¼ ½0.65; −0.11; 0.75) to compute our pairs.To address the relative nature of displacement measured by InSAR, we treat BRDY as a baseline for deformation and difference InSAR-measured range change around BRD1 from range change measured around BRDY to estimate InSAR-observed range change within the deforming region relative to BRDY.To find comparable measurements from GPS, we first convert GPS displacements to range (i.e., displacement along LOS of the SAR sensor) via the unit pointing vector ŝ.We then difference the GPS measurements at each station over time according to each epoch of the InSAR pair to find measurements of displacement over the span of the InSAR pair.Finally, we find the GPS-derived range changes from BRD1 with respect to those from BRDY to similarly estimate GPS-observed range change from the deforming region only.We consider three InSAR pairs spanning a period of normal operations at the site (i.e., when deformation should be occuring): 2016-03-26 to 2016-06-06, 2016-05-13 to 2016-06-06, and 2016-06-06 to 2016-06-30.Our results are shown in Table 4 with uncertainties taken from standard deviation measurements and uncertainty error propagation for our InSAR and GPS estimates, respectively.After comparison, we find a mean difference of 1.2 mm∕year between InSAR-observed measurements of surface deformation and GPS-observed measurements of surface deformation due to subsurface processes, which is in good agreement with the quantified Sentinel-1 uncertainty from a benchmark study. 17de and Data Availability SAR data used in this study are made publicly available by the ESA under the Sentinel-1 constellation.It may be accessed free-of-charge from the Alaska Satellite Facility using the Distributed Active Archive Center 37 as well as through ESA's Copernicus Data Space Ecosystem. 38The digital elevation models used in this study may be freely obtained via GMTSAR resources. 39GPS data used for validation are also freely available from the Department of Energy's Geothermal Data Repository. 34InSAR processing software GMTSAR 12,36 was used in the generation and timeseries analysis of the Sentinel-1 InSAR pairs in this work; it is an open-source software package that is available on GitHub (https://github.com/gmtsar).Google Earth and Generic Mapping Tools 13 were used for visualization.The spatial deformation modeling scripts used for our volume change analysis were based on those available from GitHub's General Inversion of Phase Technique repository (Brady branch, https://github.com/feigl/gipht/tree/Brady).

Fig. 4 Fig. 5 Fig. 6
Fig. 4 (a) Average displacement rate along radar LOS at the McArthur River mine during normal operations in 2017 found from SBAS analysis of 50 pairs.Results are overlaid on Google Earth image (Landsat/Copernicus and Maxar Technologies, 2022).Positive displacement indicates uplift, and negative displacement indicates subsidence.The location of the main operations is outlined in red.(b) Average displacement rate at the McArthur River mine during period of inactivity in 2019 found from SBAS analysis of 70 pairs.

Fig. 8
Fig. 8 Average displacement rate along radar LOS at the inactive Beverley mine during 2017.Results are overlaid on Google Earth image (CNES/Airbus, 2022).Positive displacement indicates uplift, and negative displacement indicates subsidence.Production well locations 7 are outlined in black.

Fig. 7
Fig. 7 Average displacement rate along radar LOS at the Four Mile mine during 2017.Results are overlaid on Google Earth image (CNES/Airbus, 2022).Positive displacement indicates uplift, and negative displacement indicates subsidence.Monitoring well locations from records 7 are marked with black and white targets.

Fig. 9 Fig. 10
Fig.9(a) Best-fitting estimates of rates of voxel-wise volumetric change under the assumption that all deformation at Four Mile can be attributed to the ore removal process.Voxels considered to be in the deforming region (where _V ðF Þ < −2σ ðF Þ _ V )are outlined in double white lines.(b)-(d) Deformation fields corresponding to the voxel parameterization and Bayesian, geostatistical analysis done at Four Mile: (b) observed, (c) modeled, and (d) residual (observed minus modeled).The deformation fields are shown in terms of displacement rate (mm/year) along line of sight of the radar, where positive values indicate uplift and negative values indicate subsidence.Map coordinates are in Easting and Northing (km) of the UTM projection on the WGS84 ellipsoid zone 54 J.22

Figure 4 (
Figure 4(a) shows the results from the SBAS analysis of average LOS velocity at McArthur River over 2017.Comparing these results with the site location in Google Earth (e.g., Fig.1), we can see very clear subsidence centered on the mining site, where most of the operations are taking place underground.

5 .
The corresponding best-fitting estimates of the volume change rate are shown in Fig. 5(a).The modeled deformation field resulting from such estimates of the volume change rate, as well as the observed deformation field and residual deformation field, are shown in Figs.5(b)-5(d).

Fig. 11 (
Fig. 11 (a) Best-fitting estimates of rates of voxel-wise volumetric change under the assumption that deformation at McArthur River can be attributed to both thermal contraction and the ore removal process.Voxels considered to be in the deforming region (where _ V ðMÞ < −2σ ðMÞ _ V ) are outlined in double white lines.(b)-(d) Deformation fields corresponding to the voxel parameterization and Bayesian, geostatistical analysis done at McArthur River: (b) observed, (c) modeled, and (d) residual (observed minus modeled).The deformation fields are shown in terms of LOS displacement rate (mm/year), where positive values indicate uplift and negative values indicate subsidence.Map coordinates are in Easting and Northing (km) of the UTM projection on the WGS84 ellipsoid zone 13 V.22

Fig. 12
Fig.12Histogram of best-fitting estimates of volumetric strain rates for McArthur River under a prior attributing observed deformation to both ore removal and thermal contraction of the rock matrix.The prior probability distribution function is shown in red, and the posterior probability function is shown in blue.

Fig. 13 (
Fig. 13 (a) Best-fitting estimates of rates of voxel-wise volumetric change at Four Mile using a wide prior.Voxels considered to be in the deforming region (where V ðF Þ < −2σ ðF Þ V ) are outlined in white.(b)-(d) Deformation fields corresponding to the voxel parameterization and Bayesian, geostatistical analysis done at Four Mile: (b) observed, (c) modeled, and (d) residual (observed minus modeled).The deformation fields here are shown in terms of the LOS displacement rate (mm/year), where positive values indicate uplift and negative values indicate subsidence.Map coordinates are in Easting and Northing (km) of the UTM projection on the WGS84 ellipsoid zone 54 J.22

6. 1
McArthur River S1 T107 2017 Baseline Plot and PairsThe baseline plot for the epochs considered in our 2017 McArthur River analysis as well as the list of pairs in our final time-series analysis are shown in Fig.15and Table2, respectively.

Fig. 14
Fig.14Histogram of best-fitting estimates of volumetric strain rates at Four Mile under a prior that allows for localized contraction and expansion.The prior probability distribution function is shown in red, and the posterior probability function is shown in blue.

Fig. 15
Fig. 15 Baseline plot connecting epochs in our McArthur River 2017 dataset via temporal and orbital separation.

Fig. 16
Fig. 16 Baseline plot connecting epochs in our Four Mile 2017 dataset via orbital and temporal separation.

Table 1
21pical by-product elements recoverable from sulfuric acid ISL mining.21

Table 2
T107 Pair epochs for McArthur River 2017 analysis.

Table 3
T60 Pair epochs for Four Mile East 2017 analysis.

Table 4
Results from InSAR-GPS validation at Brady Hot Springs, Nevada, using InSAR data from Sentinel-1 Track T144 and GPS data from the MAGNET network's BRDY and BRD1 continuous GPS stations.