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

Elena C. Reinisch; Bradley G. Henderson

doi:10.1117/1.JRS.17.044511

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

Elena C. Reinisch, Bradley G. Henderson

Author Affiliations +

Journal of Applied Remote Sensing, Vol. 17, Issue 4, 044511 (November 2023). https://doi.org/10.1117/1.JRS.17.044511

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.

1. Introduction

Subsidence signatures due to mining, especially that of copper, gold, and coal, have been a subject of literature.¹ 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.² ISL 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.² 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.³ However, ground subsidence and/or uplift can also occur due to anthropogenic activities, including geothermal energy production and mining.⁴^,⁵ Differencing 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).

Fig. 1

Overhead image of McArthur River mine with main operations outlined in red;⁶ image is from Google Earth (Landsat/Copernicus and Maxar Technologies, 2022). Spatial coverage is the same as that for the InSAR analysis.

Fig. 2

Overhead image of Four Mile mine, East section, with locations of monitoring wells marked with black and white targets;⁷ image is from Google Earth (CNES/Airbus, 2022). Spatial coverage is the same as that for the InSAR analysis.

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.⁶ Although 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⁶ (see Figs. 3 Fig. 4 Fig. 5 Fig. 6–7 and 4 Fig. 5 Fig. 6 Fig. 7 Fig. 8 Fig. 9–10, especially). Mining occurs at an average depth of 550 m.⁶ There were two major active mining zones in 2017: zone 1 (near shaft 1) and zone 4 (near shaft 3).⁶ Zone 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.

Fig. 3

General approach to deformation analysis per site in this study.

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. 5

(a) Best-fitting estimates of rates of voxel-wise volumetric change under the assumption that all deformation at McArthur River can be attributed to the ore removal process. Voxels considered to be in the deforming region (where ${\dot{V}}^{(M)} < - 2 σ_{\dot{V}}^{(M)}$ ) are outlined in white. (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). We note that the deformation fields are shown in terms of displacement rate (mm/year) along LOS of the radar, where positive values indicate uplift and negative values indicate subsidence. Map coordinates are in Easting and Northing (km) of the Universal Transverse Mercator (UTM) projection on the WGS84 ellipsoid zone 13 V.²²

Fig. 6

Histogram of best-fitting estimates of volumetric strain rates for McArthur River under a prior that attributes all observed deformation to ore removal alone. The prior probability distribution function is shown in red, and the posterior probability function is shown in blue.

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⁷ are marked with black and white targets.

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⁷ are outlined in black.

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 ${\dot{V}}^{(F)} < - 2 σ_{\dot{V}}^{(F)}$ ) 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.²²

Fig. 10

Histogram of best-fitting estimates of volumetric strain rates at Four Mile under a prior that attributes all observed deformation to ore removal alone. The prior probability distribution function is shown in red, and the posterior probability function is shown in blue.

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 specifically on Four Mile East for this study (e.g., Fig. 2), the location of which is centered on a large deposit of uranium.⁸ We consider data from 2017, which was the last year that production was predominantly focused in this area of the mine.⁹

Our 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

2.1.

SAR/InSAR Data

We use data from ESA’s S1 constellation.¹⁰ S1 consists of two satellites carrying C-band wavelength sensors ( $\sim 5.5 cm$ ); each satellite has a standard cadence of about 12 days. We use single-look 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.¹⁰ In S1 IW mode, each SAR scene consists of three subswath images (roughly $5 m \times 20 m$ ).¹⁰ In addition to the S1 data, SRTM 1-arc sec digital elevation models were obtained for topographic comparison.¹¹

We 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 $\hat{s} = [- 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 $\hat{s} = [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,¹² an open-source InSAR processing software that uses Generic Mapping Tools¹³ 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.¹⁴

Baseline plots and the final pair datasets for McArthur River and Four Mile are shown in Appendix A.

2.2.

Time-Series Analysis

To improve the SNR measured deformation signatures in our InSAR pairs, we perform time-series 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.¹⁵ This 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.

2.3.

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;¹⁶ 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

Eq. (1)

{\dot{ρ}}^{(\mod)} = Gm,

where

{\dot{ρ}}^{(\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.¹⁶

2.3.1.

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.⁶ We 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$ .¹⁷ We 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_{u 3 o 8}^{(M)} = 7.6 Mkg$ .⁶ Given an ore density of $ρ_{u 3 o 8} = 8380 kg / m^{3}$ , this equates to a volume change rate ${\dot{V}}_{ore}^{(M)}$ of

Eq. (2)

{\dot{V}}_{ore}^{(M)} = m_{u 3 o 8}^{(M)} / ρ_{u 3 o 8} = - 8.7 \times 10^{2} m^{3} / year,

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,¹⁸ and the average grade of ore for 2017 was roughly 8%.⁶ We therefore calculate the weight

m_{rock}^{(M)}

of removed host rock:

Eq. (3)

m_{rock}^{(M)} = m_{u 3 o 8}^{(M)} / 0.08 = 76.2 Mkg .

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}$ .⁶^,¹⁹ Although the subsurface is comprised of a mix of these minerals, we calculate end-member values for the volume change rate ${\dot{V}}_{rock}^{(M)}$ associated with the removal of host rock:

Eq. (4)

{\dot{V}}_{rock}^{(M)} = m_{rock}^{(M)} / [ρ_{q}, ρ_{s}] = [- 3.5, - 4.1] \times 10^{4} m^{3} / year .

From these two values of the volume change rate, we calculate a reasonable range of total volume change rates associated with mining at McArthur River:

Eq. (5)

{\dot{V}}^{(M)} = {\dot{V}}_{ore}^{(M)} + {\dot{V}}_{rock}^{(M)} = [- 3.6, - 4.2] \times 10^{4} m^{3} / year .

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_{T}^{(M)}$ is $6.9 \times 10^{8} m^{3}$ . The volume change rates in Eq. (5) are converted to volumetric strain rates $ϵ^{(M)} = {\dot{V}}^{(M)} / V_{T}^{(M)} = [5.2, 6.1] \times 10^{- 5} {year}^{- 1}$ . We find a suitable $μ - 1 σ$ prior for volumetric strain rates at McArthur River to be

Eq. (6)

{\dot{ϵ}}^{(M)} \in [- 5.6 \pm 2.8] \times 10^{- 5} {year}^{- 1} .

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.

2.3.2.

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⁷ 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.⁷ Some 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,⁷ 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 ${\dot{V}}_{bleed}^{(F)}$ associated with bleed for the given bleed percentage range:

Eq. (7)

{\dot{V}}_{bleed}^{(F)} = [0.01, 3] \times (- 1.6 \times 10^{9}) / 1000 = [- 0.16, - 4.73] \times 10^{4} m^{3} / year .

We then estimate the total rate of volume change ${\dot{V}}_{ore}^{(F)}$ associated with the removed ore. The total mass of produced ore in its final product form is $m_{u 3 o 8}^{(F)} = 1.8 \times 10^{6} kg$ with a density $ρ_{u 3 o 8} = 8380 kg / m^{3}$ .²⁰ From this and the molar masses of both the final product ore $M_{u 3 o 8} = 842 g / mole$ and mined ore $M_{u} = 238 g / mole$ , we calculate the corresponding mass of mined ore:

Eq. (8)

m_{u}^{(F)} = m_{u 3 o 8}^{(F)} / (M_{u 3 o 8} \times 10^{- 3}) \times 3 \times (M_{u} * 10^{- 3}) = 1.5 \times 10^{6} kg .

We now calculate the volume change rate ${\dot{V}}_{ore}^{(F)}$ of mined ore given its density $ρ_{u} = 1900 kg / m^{3}$ :

Eq. (9)

{\dot{V}}_{ore}^{(F)} = m_{u}^{(F)} / ρ_{u} = - 7.9 \times 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.²¹ Using 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 ${\dot{V}}_{by-prod}^{(F)} = 21.9 m^{3} / year$ . We find a range for the total volume change rate ${\dot{V}}^{(F)}$ associated with mining at Four Mile to then be

Eq. (10)

{\dot{V}}^{(F)} = {\dot{V}}_{bleed}^{(F)} + {\dot{V}}_{ore}^{(F)} + {\dot{V}}_{by-prod}^{(F)} = [- 0.24, - 4.81] \times 10^{4} m^{3} / year .

Table 1

Typical by-product elements recoverable from sulfuric acid ISL mining.21

Element	Average solution concentration (mg/L)
Rhenium	0.2 to 0.5
Molybdenum	10 to 40
Selenium	0 to 50 to 60
Vanadium	10 to 40
Scandium	0.15 to 0.6
Yttrium	0 to 15
Rare earths	10 to 40

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_{T}^{(F)} = 4.2 \times 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 ${\dot{ϵ}}^{(F)} = [- 0.46, - 9.23] \times 10^{- 4} {yr}^{- 1}$ . We establish a wide $μ \pm 1 σ$ prior for reasonable volumetric strain rates associated with mining:

Eq. (11)

{\dot{ϵ}}^{(F)} \in [- 4.7 \pm 4.4] \times 10^{- 4} {year}^{- 1} .

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.

3. Results

3.1.

McArthur River Underground Mine

3.1.1.

Time-series analysis

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.

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).

3.1.2.

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)} = {\bar{\dot{ϵ}}}^{(M)} = 5.4 \times 10^{- 6} {year}^{- 1}$ , which is significantly outside our prior range of $[- 5.6 \pm 2.8] \times 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)} = 2.6 \times 10^{- 5} {year}^{- 1}$ . Considering voxels with ${\dot{V}}^{(M)} < - 2 σ_{\dot{V}}^{(M)}$ to be significantly deforming, we calculate the volume of the deforming region to be $V = 3.1 \times 10^{7} m^{3}$ and the volume change rate of the voxels within the modeled chamber to be $\dot{V} = (- 7.7 \pm 0.1) \times 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,²³^,²⁴ where a misfit value close to one indicates that the observed data is well represented by the model. We find that $\sqrt{χ^{2}} = 2.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).

3.2.

Four Mile In Situ Leach Mine

3.2.1.

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.²⁵ We also apply a smoothing factor to weight against transient changes in topography²⁶ 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.⁸ Given that the layout at Four Mile is such that monitoring wells surround injection and production wells,⁷ 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 $\sim 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.

3.2.2.

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)} = {\bar{\dot{ϵ}}}^{(F)} = - 1.8 \times 10^{- 5} {year}^{- 1}$ , which is slightly outside our prior $μ - 1 σ$ range of $[- 4.7 \pm 4.4] \times 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 \times 10^{- 4} {year}^{- 1}$ . We find the volume of the modeled deforming region to be $V = 2.1 \times 10^{6} m^{3}$ and the volume change rate of the voxels within that region to be $\dot{V} = (- 4.8 \pm 0.01) \times 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 $\sqrt{χ^{2}} = 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).

4. 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.⁷^,⁹

To 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.

4.1.

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.⁶^,²⁷^,²⁸ Generally 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

Eq. (12)

{\dot{ϵ}}_{T} = α_{T} \dot{T},

where

α_{T}

represents the thermal volumetric expansion coefficient and

\dot{T}

represents the temperature change rate at the site.²⁹ Based on the literature,³⁰ we assume a reasonable mean estimate of

\dot{T} = - 27 K / year

and a mean estimate of

α_{T} = 5.5 \times 10^{- 6} 1 / K

. We assign a wide uncertainty of

- 27 K

for

\dot{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 \times 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.³¹ This leads to a prior of

Eq. (13)

{\dot{ϵ}}_{T}^{(M)} \in (- 1.5 \pm 1.7) \times 10^{- 4} {year}^{- 1} .

Our joint prior is defined as the sum of the two random variables characterizing the volumetric strain rate ${\dot{ϵ}}_{T}^{(M)}$ related to thermal contraction and the volumetric strain rate ${\dot{ϵ}}_{m}^{(M)}$ related to ore removal [e.g., Eq. (6)], which we find through the sum of two Gaussians:³¹

Eq. (14)

{\dot{ϵ}}_{m + T}^{(M)} = {\dot{ϵ}}_{T}^{(M)} + {\dot{ϵ}}_{m}^{(M)} \in (- 2.0 \pm 1.7) \times 10^{- 4} {year}^{- 1} .

Repeating our analysis using this prior, we find a posterior mean volumetric strain rate estimate of $μ_{m + T}^{(M)} = {\bar{\dot{ϵ}}}_{m + T}^{(M)} = - 1.2 \times 10^{- 4} {year}^{- 1}$ [Fig. 11(a)]. This value falls within the $μ + 1 σ$ prior range defined in Eq. (14), and the corresponding maximum posterior parameter uncertainty is $1.4 \times 10^{- 4} {yr}^{- 1}$ , which is smaller than our prior uncertainty (Fig. 12). We find the volume of the deforming region to be $V = 3.7 \times 10^{7} m^{3}$ with a volume change rate of $(- 1.1 \pm 0.1) \times 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 $\sqrt{χ^{2}} = 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.

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 ${\dot{V}}^{(M)} < - 2 σ_{\dot{V}}^{(M)}$ ) 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.²²

Fig. 12

Histogram 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.

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.³¹ At 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.

4.2.

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 \times 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)} = {\bar{\dot{ϵ}}}^{(F)} = - 1.4 \times 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 \times 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 $\dot{V} = (- 8.6 \pm 0.2) \times 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 $\sqrt{χ^{2}} = 1.2$ .

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 ${\dot{V}}^{(F)} < - 2 σ_{\dot{V}}^{(F)}$ ) 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.²²

Fig. 14

Histogram 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.

To determine if our wider prior provides a better fit to our previously defined prior, we run an $F$ -test for model comparison.³¹ At 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}$ , 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.

5. 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

6.1.

McArthur River S1 T107 2017 Baseline Plot and Pairs

The 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. 15 and Table 2, respectively.

Table 2

T107 Pair epochs for McArthur River 2017 analysis.

Reference (YYYY-MM-DD)	Repeat (YYYY-MM-DD)
2017-03-01	2017-03-13
2017-03-13	2017-03-25
2017-04-06	2017-04-18
2017-04-06	2017-04-30
2017-04-06	2017-05-12
2017-04-06	2017-06-05
2017-04-06	2017-06-17
2017-04-06	2017-06-29
2017-04-18	2017-04-30
2017-04-18	2017-05-12
2017-04-18	2017-06-05
2017-04-18	2017-06-17
2017-04-30	2017-05-12
2017-04-30	2017-05-24
2017-04-30	2017-06-05
2017-04-30	2017-06-17
2017-04-30	2017-06-29
2017-04-30	2017-07-11
2017-05-12	2017-05-24
2017-05-12	2017-06-05
2017-05-12	2017-06-17
2017-05-12	2017-06-29
2017-05-12	2017-07-11
2017-05-12	2017-08-04
2017-05-24	2017-06-29
2017-05-24	2017-07-11
2017-05-24	2017-08-04
2017-05-24	2017-08-16
2017-06-05	2017-06-17
2017-06-05	2017-06-29
2017-06-05	2017-07-11
2017-06-05	2017-08-16
2017-06-05	2017-08-28
2017-06-17	2017-07-11
2017-06-17	2017-08-04
2017-06-17	2017-08-16
2017-06-17	2017-08-28
2017-06-29	2017-07-11
2017-06-29	2017-09-09
2017-07-11	2017-08-16
2017-07-11	2017-08-28
2017-07-11	2017-09-09
2017-07-11	2017-10-03
2017-08-04	2017-09-09
2017-08-16	2017-09-09
2017-08-16	2017-10-03
2017-08-28	2017-09-09
2017-08-28	2017-10-03
2017-09-09	2017-10-03
2017-10-03	2017-10-15

Fig. 15

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

6.2.

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.

Table 3

T60 Pair epochs for Four Mile East 2017 analysis.

Reference (YYYY-MM-DD)	Repeat (YYYY-MM-DD)
2017-01-08	2017-01-20
2017-01-08	2017-03-09
2017-01-20	2017-02-25
2017-02-01	2017-03-09
2017-02-01	2017-04-02
2017-02-13	2017-04-02
2017-02-25	2017-03-09
2017-02-25	2017-04-02
2017-02-25	2017-04-14
2017-02-25	2017-05-08
2017-02-25	2017-05-20
2017-03-09	2017-04-02
2017-03-09	2017-04-14
2017-03-09	2017-05-08
2017-03-09	2017-05-20
2017-03-09	2017-06-01
2017-04-02	2017-04-26
2017-04-02	2017-06-01
2017-04-02	2017-06-25
2017-04-14	2017-04-26
2017-04-14	2017-06-01
2017-04-14	2017-06-25
2017-04-14	2017-07-07
2017-04-26	2017-05-08
2017-04-26	2017-05-20
2017-04-26	2017-06-01
2017-04-26	2017-07-19
2017-05-08	2017-05-20
2017-05-08	2017-06-13
2017-05-08	2017-06-25
2017-05-20	2017-06-01
2017-05-20	2017-06-25
2017-05-20	2017-07-19
2017-06-01	2017-07-07
2017-06-01	2017-08-24
2017-06-13	2017-08-24
2017-06-13	2017-09-05
2017-06-25	2017-07-19
2017-06-25	2017-09-17
2017-07-07	2017-07-19
2017-07-07	2017-07-31
2017-07-07	2017-08-12
2017-07-07	2017-08-24
2017-07-07	2017-09-05
2017-07-19	2017-08-24
2017-07-19	2017-09-05
2017-07-31	2017-08-24
2017-07-31	2017-09-05
2017-07-31	2017-10-11
2017-08-12	2017-08-24
2017-08-12	2017-09-05
2017-08-12	2017-10-11
2017-09-05	2017-10-23
2017-09-05	2017-11-28
2017-09-17	2017-09-29
2017-09-17	2017-11-04
2017-09-17	2017-11-28
2017-09-17	2017-12-10
2017-09-29	2017-12-22
2017-10-11	2017-12-22
2017-10-23	2017-11-04
2017-11-04	2017-11-28
2017-11-04	2017-12-10
201711-16	2017-12-10

7. 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 InSAR-observed 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,¹⁶^,³² 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.³³ Both stations are part of the MAGNET network with publicly accessible time series of relative position analyzed using standard procedures.³⁴^,³⁵ Their relative proximity ( $\sim 3 km$ separation) allows deformation around both stations to be measured by the same interferogram (Fig. 16).

Fig. 16

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

We conduct our validation similar to previous works,¹⁶^,³² using data from Sentinel-1’s track T144 ( $\hat{s} = [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 $\hat{s}$ . 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.¹⁷

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.

Pair	InSAR ρ˙ (mm/year)	GPS ρ˙ (mm/yearr)
2016-03-26 to 2016-06-06	$2.5 \pm 0.8$	$3.5 \pm 8.2$
2016-05-13 to 2016-06-06	$4.8 \pm 0.3$	$3.2 \pm 8.1$
2016-06-06 to 2016-06-30	$2.8 \pm 0.1$	$1.7 \pm 9.0$

Code 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³⁷ as well as through ESA’s Copernicus Data Space Ecosystem.³⁸ The digital elevation models used in this study may be freely obtained via GMTSAR resources.³⁹ GPS data used for validation are also freely available from the Department of Energy’s Geothermal Data Repository.³⁴ InSAR processing software GMTSAR¹²^,³⁶ was used in the generation and time-series 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¹³ 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).

Acknowledgments

We thank Timothy Johnson of Pacific Northwest National Laboratory for his project leadership and two anonymous reviewers for their helpful comments. We additionally thank the European Space Agency and the Alaska Satellite Facility for open access to Sentinel-1 data, the Nevada Geodetic Laboratory for open access to GPS time series, and the GMTSAR and GMT development teams for making their InSAR processing and geospatial software suites freely available. The authors have no relevant financial interests in the manuscript and no other potential conflicts of interest to disclose.

References

1.

M. Eneva, “Detection of surface deformation at mining and geothermal sites using Satellite Radar Interferometry (InSAR),” in 44th US Rock Mech. Symp. and 5th US-Canada Rock Mech. Symp., (2010). Google Scholar

2.

M. Abzalov, Applied Mining Geology, 12 Springer( (2016). Google Scholar

3.

D. Massonnet and K. L. Feigl, “Radar interferometry and its application to changes in the Earth’s surface,” Rev. Geophys., 36 (4), 441 –500 https://doi.org/10.1029/97RG03139 REGEEP 8755-1209 (1998). Google Scholar

4.

G. Falorni, J. Morgan and M. Eneva, “Advanced InSAR techniques for geothermal exploration and production,” GRC Trans., 35 1661 –1666 (2011). Google Scholar

5.

Z. Yang et al., “Use of SAR/InSAR in mining deformation monitoring, parameter inversion, and forward predictions: a review,” IEEE Geosci. Remote Sens. Mag., 8 (1), 71 –90 https://doi.org/10.1109/MGRS.2019.2954824 (2020). Google Scholar

6.

L. Bray, G. M. Murdock and A. D. Renaud, McArthur river operation Northern Saskatchewan, Canada National Instrument 43-101 Technical Report (2018). Google Scholar

7.

C. Bartels and H. R. P. Ltd, Four Mile uranium mine program for environment protection and rehabilitation, and radioactive waste management plan (2018). Google Scholar

8.

S. B. Hore and S. M. Hill, “Palaeoredox fronts: setting and associated alteration exposed within a key section for understanding uranium mineralisation at the Four Mile West deposit,” MESA J., 55 34 –39 (2009). Google Scholar

9.

G. Birch and H. R. P. Ltd, 2019 mining compliance report—Beverley, Beverley North & Four Mile mineral leases (2020). Google Scholar

10.

D. Geudtner et al., “Sentinel-1 system capabilities and applications,” in IEEE Geosci. and Remote Sens. Symp., 1457 –1460 (2014). https://doi.org/10.1109/IGARSS.2014.6946711 Google Scholar

11.

T. G. Farr et al., “The shuttle radar topography mission,” Rev. Geophys., 45 (2), RG2004 https://doi.org/10.1029/2005RG000183 REGEEP 8755-1209 (2007). Google Scholar

12.

D. Sandwell et al., “GMTSAR: an InSAR processing system based on generic mapping tools,” (2011). Google Scholar

13.

P. Wessel et al., “Generic mapping tools: improved version released,” EOS Trans. Am. Geophys. Union, 94 (45), 409 –410 https://doi.org/10.1002/2013EO450001 EOSTAJ 0096-3941 (2013). Google Scholar

14.

C. W. Chen and H. A. Zebker, “Network approaches to two-dimensional phase unwrapping: intractability and two new algorithms,” J. Opt. Soc. Am. A, 17 (3), 401 –414 https://doi.org/10.1364/JOSAA.17.000401 (2000). Google Scholar

15.

P. Berardino et al., “A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms,” IEEE Trans. Geosci. Remote Sens., 40 2375 –2383 https://doi.org/10.1109/TGRS.2002.803792 IGRSD2 0196-2892 (2002). Google Scholar

16.

E. C. Reinisch, M. Cardiff and K. L. Feigl, “Characterizing volumetric strain at Brady Hot Springs, Nevada, USA using geodetic data, numerical models and prior information,” Geophys. J. Int., 215 (2), 1501 –1513 https://doi.org/10.1093/gji/ggy347 GJINEA 0956-540X (2018). Google Scholar

17.

Z. Sadeghi et al., “Benchmarking and inter-comparison of Sentinel-1 InSAR velocities and time series,” Remote Sens. Environ., 256 112306 https://doi.org/10.1016/j.rse.2021.112306 (2021). Google Scholar

18.

B. Jamieson, “Mining the high grade McArthur River uranium deposit,” (2002). Google Scholar

19.

D. Bronkhorst et al., “McArthur river operation northern Saskatchewan, Canada National Instrument 43-101 Technical Report,” (2012). Google Scholar

20.

B. C. Reed, The History and Science of the Manhattan Project, Springer( (2014). Google Scholar

21.

V. Benes et al., “Manual of acid in situ leach uranium mining technology,” (2001). Google Scholar

22.

J. P. Snyder, Map Projections–A Working Manual, 1395 US Government Printing Office( (1987). Google Scholar

23.

G. Strang and K. Borre, Linear Algebra, Geodesy, and GPS, SIAM( (1997). Google Scholar

24.

R. C. Aster, B. Borchers and C. H. Thurber, Parameter Estimation and Inverse Problems, Elsevier( (2018). Google Scholar

25.

X. Xu et al., “Tectonic and anthropogenic deformation at the Cerro Prieto geothermal step-over revealed by Sentinel-1A InSAR,” IEEE Trans. Geosci. Remote Sens., 55 (9), 5284 –5292 https://doi.org/10.1109/TGRS.2017.2704593 IGRSD2 0196-2892 (2017). Google Scholar

26.

X. Tong and D. Schmidt, “Active movement of the Cascade landslide complex in Washington from a coherence-based InSAR time series method,” Remote Sens. Environ., 186 405 –415 https://doi.org/10.1016/j.rse.2016.09.008 (2016). Google Scholar

27.

G. Newman et al., “Artificial ground freezing: an environmental best practice at Cameco’s uranium mining operations in Northern Saskatchewan, Canada,” in 11th Int. Mine Water Assoc. Congr.–Mine Water–Managing the Challenges, 113 –118 (2011). Google Scholar

28.

M. Semin, I. Golovatyi and A. Pugin, “Analysis of temperature anomalies during thermal monitoring of frozen wall formation,” Fluids, 6 (8), 297 https://doi.org/10.3390/fluids6080297 (2021). Google Scholar

29.

H. F. Wang, Theory of Linear Poroelasticity with Applications to Geomechanics and Hydrogeology, 2 Princeton University Press( (2000). Google Scholar

30.

T. Smith, D. Milne and Z. Szczepanik, “Interpretation of stressmeter data in variable temperature conditions,” in 1st Canada-US Rock Mech. Symp., (2007). Google Scholar

31.

D. Wackerly, W. Mendenhall and R. L. Scheaffer, Mathematical Statistics with Applications, Cengage Learning( (2014). Google Scholar

32.

E. C. Reinisch et al., “Time-series analysis of volume change at Brady Hot Springs, Nevada, USA, using geodetic data from 2003–2018,” J. Geophys. Res. Solid Earth, 125 (9), e2019JB017816 https://doi.org/10.1029/2019JB017816 (2020). Google Scholar

33.

K. L. Feigl and L. M. Parker, “PoroTomo final technical report: poroelastic tomography by adjoint inverse modeling of data from seismology, geodesy, and hydrology,” Madison, WI, United States (2019). Google Scholar

34.

C. Kreemer, “BRAD BRDY and BRD1 GPS station RINEX files 01-19-2018 [data set],” (2018). Google Scholar

35.

G. Blewitt et al., “Terrestrial reference frame NA12 for crustal deformation studies in North America,” J. Geodyn., 72 11 –24 https://doi.org/10.1016/j.jog.2013.08.004 JOGEE7 0264-3707 (2013). Google Scholar

36.

D. Sandwell et al., “Open radar interferometry software for mapping surface deformation,” EOS Trans. Am. Geophys. Union, 92 (28), 234 –234 https://doi.org/10.1029/2011EO280002 EOSTAJ 0096-3941 (2011). Google Scholar

37.

Alaska Satellite Facility, “ASF data search vertex,” https://search.asf.alaska.edu (2023). Google Scholar

38.

European Space Agency, “Copernicus data space ecosystem,” https://dataspace.copernicus.eu (2023). Google Scholar

39.

GMTSAR, “Generate DEM files for use with GMTSAR,” https://topex.ucsd.edu/gmtsar/demgen/ (2010). Google Scholar

Biography

Elena C. Reinisch is a staff scientist in the Space Remote Sensing and Data Science Group (ISR-6) at Los Alamos National Laboratory. She received her PhD in geophysics from the University of Wisconsin–Madison in 2019, having developed methods to characterize InSAR-observed subsidence at geothermal fields. Her research interests are in the combination of remote sensing, data science, and computational geodesy, and she has expertise in developing algorithms for SAR analysis of Earth observation studies.

Bradley G. Henderson is a staff scientist in ISR-6 at Los Alamos National Laboratory. He received his PhD in geophysics from the University of Colorado in 1995, investigating thermal-infrared emission spectra of airless planetary bodies. His main research interest is the physics and techniques of remote sensing. His experience includes algorithm development and data analysis for hyper- and multi-spectral remote sensing, radiative transfer modeling of light in particulate solids, polarization, and analysis of SAR imagery.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 International License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Citation Download Citation

Elena C. Reinisch and Bradley G. Henderson "Spatio-temporal analysis and volumetric characterization of interferometric synthetic aperture radar-observed deformation signatures related to underground and in situ leach mining," Journal of Applied Remote Sensing 17(4), 044511 (14 November 2023). https://doi.org/10.1117/1.JRS.17.044511

Received: 21 June 2023; Accepted: 16 October 2023; Published: 14 November 2023

Access the abstract

JOURNAL ARTICLE
23 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

KEYWORDS

Deformation

Mining

Voxels

Interferometric synthetic aperture radar

Global Positioning System

Data modeling

Uranium