Introduction

The storage of industrial wastes is a worldwide environmental problem. Among the most prolific industries generating solid wastes is the phosphate fertilizer industry. This generates large amounts of an acid residue called phosphogypsum (PG) during the production of phosphoric acid from phosphatic rocks. With a ratio of 5 t of PG per ton of phosphoric acid (USEPA 2002), the estimated world PG production is 100 to 280 Mt per year (Parreira et al. 2003; Yang et al. 2009). Of this, only 15% is recycled, mostly in the elaboration of building and road materials, fertilizers and soil remediation compounds (Zhang et al. 2020 and references therein). The remaining 85% is stored in large retaining stockpiles, either under dry or humid conditions (Rutherford et al. 1994).

The major environmental problem concerning these stacks is their high contamination potential because PG accumulates naturally occurring toxic and radioactive elements in its mineral structure. While this environmental hazard represents a serious concern for the scientific community (Papanicolaou et al. 2009), little attention is paid, from the public domain, to the structural stability of the stacks, despite such stability exerts direct control on the safety storage of toxic PG. This disparity is well exemplified in the Huelva PG stack, SW Spain.

The Huelva PG stack is located in an estuarine environment, very close to the Atlantic Ocean, and less than 0.5 km from the city of Huelva. Reaching heights of up to 25 m, this wet embankment contains 100 Mt of wastes occupying an area of about 1000 ha. The stack received residues of up to five nearby factories and was in operation since 1967 until 2010, when it was concealed after a Spanish National Court decision.

In the last 20 years, but particularly since the closure of the disposal activity, several research studies revealed that the PG stack has an important contamination potential for the estuary, the Atlantic Ocean and the citizen of Huelva. In that sense, two main research lines have been developed. One deals with the geochemical characterization of the edge outflows discharged from the PG stack directly to the estuary. Studies like that of Papaslioti et al. (2018a, b) and Perez-López et al. (2010, 2011, 2015, 2016) have characterized leachates with pH values below 2 and anomalously elevated concentrations of P, S, F, NH4+, Fe, Zn, U, Cr, Cu and Cd. They also found that the water that percolates the PG stack and form edge outflows is essentially estuarine, most likely tidal-induced.

The other research line focuses on the potential radioactivity of stored PG. Considered a naturally occurring radioactive material (NORM) by the European Union (European Union 2013; International Atomic Energy Agency 2003), these wastes incorporate high concentrations of U-series radionuclides, including the highly radiotoxic 226Ra, 210Po and 210Pb. The concentration of these natural radionuclides in the Huelva PG stack and their communication pathways with surrounding environments, including atmosphere, water, neighbouring surface sediments and basement, have been assessed by several researches (Bolivar et al. 1998, 2002, 2009; Guerrero et al. 2019, 2020, 2021a, b; Hierro et al. 2013; Mas et al. 2006). These studies concluded that their atmospheric impact is negligible, the surrounding and underlying sediments are contaminated only in the close proximity of the PG, and that the water is by far the major vector for spreading radioactive contamination into the Huelva estuarine environment.

Surprisingly, there is an absence of public research evaluating the structural stability of the PG stack and monitoring the deformation suffered by the unconsolidated salt-marsh substratum due to PG embankment overpressure. Following specific recommendations reported by the Regional Government (Junta de Andalucía 2009), the major concessionary company have developed some studied in this regard, but this information is unavailable for the public. One of the scarce researches published on this issue is that of Valverde-Palacios et al. (2011). After examining the mechanical properties of mud samples from six core drills, they obtained a positive safety factor under non-seismic and seismic conditions for the stability of the PG, and calculated a maximum period of 3.2 years to reach 90% of expected basement consolidation.

In order to fill the gap of knowledge on this issue, and provide new insights on the ground deformation induced by the Huelva PG stack, the present study includes a remote sensing analysis of the area based on SAR images.

Analysing Sentinel-1 images with a specific differential SAR interferometry (DInSAR) algorithm called Parallel Small Baseline Subset (P-SBAS) and combining ascending and descending orbital datasets of this satellite, we have been able to retrieve vertical and horizontal displacements accounted during 4.7 years in the entire area covering the PG stack. The use of DInSAR techniques for monitoring surface deformation is widely documented, especially in areas affected by earthquakes, volcanoes, landslides and subsidence (i.e., Aslan et al. 2018; Biggs et al. 2021; Brunori et al. 2015; Burgman et al. 2000, Dang et al. 2014; Fiaschi et al. (2017); Zhang et al. 2020). This microwave remote sensing technique computes the phase difference between SAR images acquired at different times in order to measure, with centimetric to millimetric accuracy, the surface deformation affecting wide areas (up to 400 km2). The development of different multitemporal DInSAR processing techniques (Berardino et al. 2002; Ferreti et al. 2001) has favoured the analysis of fairly larger SAR datasets and the retrieval of long and accurate deformation time-series (Burgmann et al. 2000; Canuti et al. 2004; Massonnet and Feigl 1998; Metternicht et al. 2005).

For the first time, this study generates displacement maps and time series deformation curves for the entire area covering a PG stack. This information, essential for a comprehensive evaluation of the stack stability, provides critical data for accurately assessing the suitability of any proposed remediation plan and for the prediction and prevention of potential structural and/or environmental hazards.

Study area

The estuary of Huelva is conformed by the confluence of the Tinto and Odiel rivers. These extend to the Atlantic Ocean through major estuarine channels bordered by relatively narrow tidal flats and surrounded by wide areas of inter-tidal salt marshes (Fig. 1). The Huelva PG stack was placed at the right margin of the Tinto estuarine channel, settled directly over salt-marsh deposits, with no impermeable and protective cap in between. In detail, the deposits immediately below the PG basement consist of an interlaying of silty clays and sandy silts crosscutted by an intricate network of sandy channels. The thickness of this salt-marsh sequence increases to the SE, from 7.5 to 8 m thick at the NW margin of the PG stack to 36 m in the vicinity of the Tinto channel (Valverde-Palacios et al. 2011). These unconsolidated soft sediments laid directly over the Neogene massive marly formation that crops out to the NW of the city of Huelva (Fig. 2). With more than 100 m thick, this unit conforms the general basement of the estuary. Hydrodynamically, the Estuary of Huelva is located in a mixed-energy, tide-dominated coast, in the sense of Hayes (1979), characterized by mean tidal ranges between 3.06 and 1.70 m, (averaging 2.2 m) and waves with significant height of 0.7 m and about 40-m wavelength (Borrego 1992; Carro et al. 1998). The internal sector of the estuary is isolated from the direct action of waves and receives the fresh water inflow from the Tinto and Odiel rivers. With an average total water flux of 81.75 Hm3/year, both rivers are characterized by strong seasonal and inter-annual fluctuations.

Fig. 1
figure 1

Image of the Iberian Peninsula showing the location of the study area at the upper left corner. Rectangles labelled with track numbers indicate radar data coverages used in this study. The larger image shows the four zones comprising the Huelva PG stack and the two cross sections (A–A′ and B–B′) in Fig. 7. At the bottom, westward aerial view of zone 2

Full size image
Fig. 2
figure 2

Geological sketch of the study area showing the location of the Huelva PG stack including an idealized geological section (A and B) immediately below the stack

Full size image

The disposal activity in the Huelva PG stack commenced in 1967 and was definitely closured at the end of 2010, after a decision of the Spanish National Court. During this time, the stack received wastes of up to five acid phosphoric factories located in the close industrial area of the Harbour of Huelva. Until 1998, the wastes, unprocessed and mixed with seawater, were accumulated over several disposal located at the salt-marshes, reaching an average high of 10 m. Then, the industrial factories were forced to accommodate new environmental restrictions, that encompassed (1) the development of a closed circuit with freshwater to pump PG and (2) the PG accumulation into pyramidal piles above previously deposited wastes. For the closed circuit, two decanting water pounds (one of them in reserve) and a bordering channel to recover lixiviates were also required.

The PG stack is divided in four zones (Fig. 1). Zone 1 received urban and mining wastes in addition to PG until 1988. Then, it was restored covering the residues with 25–30 cm of compacted soil and vegetation. This restoration plan ended in 1992. Total amount of waste stocked was estimated in 12 Mt, occupying an area of 400 ha and a maximum high of 2–3 m.

Zone 2 was the only zone that remained active after the environmental regulation adopted in 1998. It received a total of 25 Mt of PG, in an area of 300 ha, up to 25 m high. The stack consists of a succession of truncated pyramidal piles surrounded by a peripheral channel for lixiviates. At the top, a freshwater pond was used for the closed pumping system. Since zone 2 cancelation, in 2010, no restoration measures have been adopted. There is, however, a restoration plan for this zone elaborated by the concessionary company that entails the PG cover with liner, protective bare soil and vegetated top soil. Such restoration plan, presented in 2014, still pending final government approval.

Zone 3 contains 15 Mt PG, occupying 200 ha with a height of 5–8 m. This zone received PG until 1998. Since then, it accommodated the security water pond demanded for the closed pumping system. The lack of restoration measures other than the peripheral channel for leachates built around the water pound has resulted in the acceleration of weathering processes and in the proliferation of contaminated leaching outflows that reach salt-marshes and estuary channels. The environmental impact of this problem has been exposed and evaluated in several studied, including those of Pérez-López et al. (2011, 2016). The environmental plan mentioned above also entails the liner and soil cover of zone 3.

Zone 4, with 30 Mt, 150 ha and a height of 8–10 m, stored PG until 1987. Somewhat later, in 1993, it commenced to be covered with a multilayer made of building wastes, industrial wastes and a vegetation cover. The restoration plan elapsed in time for more than 10 years. Although nowadays zone 4 is considered restored, concerns about the appropriateness of the cover has forced the concessionary factory to elaborate a new restoration plan specific to this area.

Methodology

Automatic PSBAS-InSAR procedure

In order to analyse the surface deformation of the Huelva PG stack time after its closure, two sets of Sentinel-1 SAR images have been employed, one including 137 images from an ascending track, covering from 10 Oct 2016 to 22 Jun 2021, and another with 142 images from a descending track, spanning from 4 Oct 2016 to 28 Jun 2021. Table 1 includes the ancillary details of both set of images while Fig. 1A shows their respective outline. Figure 3 illustrates the workflow of the automatic PSBAS-InSAR procedure and the subsequent non-automatic 2D decomposition.

Table 1 Main characteristics of the Sentinel-1 datasets used for the analysis
Full size table
Fig. 3
figure 3

Workflow of the automatic PSBAS-InSAR (orange blocks) and 2D decomposition (blue block) procedures followed in this study

Full size image

Satellite data were processed using the service CRN-IREA P-SBAS, hosted by the web-based Geohazard Exploration Platform (https://geohazards-tep.eu/#!) of the European Space Agency (ESA). This web tool performs unsupervised multi-temporal InSAR analyses of Sentinel-1 IWS SLC images applying a specific DInSAR algorithm called P-SBAS. This algorithm is a parallel implementation of the SBAS approach (Casu et al. 2014; De Luca et al. 2015), a technique that strategically combines differential interferograms to reduce the spatial decorrelation, thus enabling analyses with independent SAR datasets and disparate spatial scales (Berardino et al. 2002). P-SBAS is thus largely used in multi-temporal DinSAR analysis (Bonano et al. 2012; Fiaschi et al. 2017; Lanari et al. 2004, 2010). Working on a distributed computing system environment from the ESA called Grid Processing on Demand (G-POD), and benefiting from the free availability of Sentinel-1 images, CRN-IREA P-SBAS considerably decreases the time and cost required for the generation of a deformation time series analysis (Galve et al. 2017).

The CRN-IREA P-SBAS web-based interface asks for the ingestion of all the input data required for the DinSAR analysis, which include an area of interest (AOI), a set of Sentinel-1 RAW images, which can be selected from a spatial and temporal framework available within the platform, the coordinates of a predictably stable point within the AOI, and the preferred polarization type, temporal coherence threshold and digital elevation model (DEM) for topography phase correction and geocoding.

For each of the two jobs comprising this study, i.e., the ascending and descending datasets, the application initiated the processing chain workflow unpacking the orbital information associated to the selected S1 IWS SLC images. Then, the master image, the one used as reference geometry for the registration of the whole SAR data set, is identified together with the potential interferometric SAR data pairs that will be processed in the forthcoming interferometric generation step. In this early stage, the AOI is taken in consideration in order to select all the bursts within the SAR dataset involved in the analysis. Considering the orbital characteristic of Sentinel-1, with a stringent orbit and very reduced revising times, it is possible to generate interferograms with very small baselines, helpful to minimize decorrelation effects (Berardino et al. 2002; Zebker and Villasenor 1992). In consequence, no perpendicular baseline constrain was required for the selection of interferometric SAR data pairs. The selected DEM, that for this analysis was the 90-m model (1 arcsec) derived from the NASA Shuttle Radar Topography Mission (SRTM), was converted into radar coordinates (range and azimuth) following the range-Doppler equation (Curlander and McDonough 1991). Then, the so-called range azimuth matrices were generated. These include, for every pixel of each SAR image involved in the analysis, the target-to-sensor range and azimuth distances with respect the orbital position of the master image (Sansosti et al. 2006). At that time, each SLC SAR image was precisely corregistrated with respect the master image according to the spectral diversity methods by Scheiber and Moreira (1999). Based on the phase information of the SAR signal, this technique effectively determines the relative misregistration between interferometric pairs on a pixel basis.

In the next step, the inferferograms of the data pairs previously established were computed. In the present analysis, 394 and 409 interferograms were used for the ascending and descending dataset, respectively. This was followed by a spatial complex average (multi-look) operation (Franceschetti and Lanari 1999) destined to mitigate the decorrelation noise affecting the interferograms and reduce the amount of data to process (Zebker and Villasenor 1992). In both ascending and descending dataset, the number of looks employed in azimuth and range were 5 and 20, respectively. The chain workflow was followed by removal of the phase residual ramps from the generated interferometric pairs in order to reduce possible residual miss-registration errors. The differential interferograms are composed by mean of a mosaicking operation aligning the corrected multi-look interferograms with their corresponding spatial coherence maps of adjacent bursts. A noise reduction operation is then executed applying the interferogram noise filtering technique proposed by Goldstein and Werner (1998) and Baran et al. (2003). This widely used low-pass filtering algorithm applies a variable bandpass sensitive to local phase noise and fringe rate. In areas highly correlated, the phase is smoothed, whereas in areas poorly correlated, it remains broad-band.

The next step involves the phase unwrapping (PhU). This was executed applying the extended minimum cost flow (EMCF) algorithm (Pepe and Lanari 2006), a technique that estimates the unwrapped phase gradients computing the MCF algorithm (Costantini 1998) in the temporal/perpendicular baseline, and applying again the same algorithm to these estimates in the azimuth/range plane. The EMCF was applied to all those pixels with a temporal coherence value higher than the threshold fixed at the input stage, which in this case was 0.85. This high value was selected in order to ensure reliable results, despite the final density of PS points was somewhat diminished. All the unwrapped interferograms are then grouped to form the final differential interferogram sequence. This is decomposed using the single value decomposition (SVD) method (Berardino et al. 2002) in order to generate the displacement time series and the mean deformation velocity map. At the same time, the temporal coherent factor of each pixel (Pepe and Lanari 2006), which quantitatively measure the reliability of the data generated, is also retrieved. Finally, the quality of the time series obtained is enhanced by applying a filter used for detecting and removing potential atmospheric artefacts (Berardino et al. 2002).

Decomposition into 2D displacement rates

The displacement of each pixel, as calculated with the P-SBAS processing chain, quantifies the movement along the satellite line of sight (LOS). To retrieve the vertical and horizontal components, this was decomposed according to the following equation:

$${d}_{los}= {d}_{v} \mathrm{cos}\theta -\left({d}_{e} \mathrm{cos}\alpha - {d}_{n} \mathrm{sin}\alpha \right) \mathrm{sin}\theta$$
(1)

where \({d}_{los}\), \({d}_{v}\), \({d}_{e}\) and \({d}_{n}\) are the displacement along LOS, vertical, east–west and north–south directions, respectively, \(\alpha\) is the satellite heading, and \(\theta\) is the incident angle.

Given the near polar orbit characterizing the Sentinel-1 sensor, LOS measurement acquired is highly insensitive to movements in the north–south direction. For this reason, the component \({d}_{n}\) was neglected from Eq. 1. The two remaining variables,\({d}_{v}\) and \({d}_{e}\), can be calculated if two or more independent measurements are accounted for each pixel. This has been achieved by combining ascending and descending measurements into the following matrix equation (Motagh et al. 2017):

$$\left(\begin{array}{c}{d}_{los}^{a}\\ {d}_{los}^{d}\end{array}\right)= \left(\begin{array}{cc}{\mathrm{cos}\theta }^{a}& -{\mathrm{cos}\alpha }^{a}{\mathrm{sin}\theta }^{a}\\ {\mathrm{cos}\theta }^{d}& -{\mathrm{cos}\alpha }^{d}{\mathrm{sin}\theta }^{d}\end{array}\right) \left(\begin{array}{c}{d}_{v}\\ {d}_{e}\end{array}\right)$$
(2)

where superscripts a and d stand, respectively, for ascending and descending measurements and geometries.

As pixels provided by ascending and descending datasets were geographically coincident (with sub-millimetric accuracy), the north–south and east–west decomposition was performed in the vector domain, thus providing a more straightforward procedure and avoiding the uncertainties introduced by raster interpolators (Foumelis 2016).

Results

Average LOS velocity maps calculated from ascending and descending InSAR time-series analysis of the city of Huelva and nearby PG stack, spanning the Oct 2016 to Jun 2021 time interval, are, respectively, illustrated in Fig. 4a, b. In both figures, negative values (yellow to red) indicate movements away from the satellite whereas positive values (cyan to blue) correspond to movement towards the satellite. Green colour represents stable zones, with movements within the interval − 0.25 to 0.25 cm/year.

Fig. 4
figure 4

Ascending (a) and descending (b) average LOS velocity maps resulted from the InSAR time-series analysis of Sentinel-1 images. Positive values (blue colours) indicate ground displacements towards the satellite while negative values (red colours) indicate ground displacements against the satellite. (c) Average vertical velocity map retrieved from the decomposition of ascending and descending LOS velocity maps. Positive values (blue colours) indicate upward displacements while negative values (red colours) indicate downward displacements. (d) Average horizontal (east–west) velocity map retrieved from the decomposition of ascending and descending LOS velocity maps. Blue colours indicate westward displacements while red colours indicate eastward displacements

Full size image

In general, the city of Huelva and most of the PG stack have a high density of persistent scattered (PS) points, with up to 16 PS points/km2. Other areas in the image are characterized by very sparse density of PS points due to dense vegetation coverage, water-table fluctuation and/or possible temporal decorrelation. Concerning the PG stack, zones 2 and 3 show high density of PS points, contrary to what occurs in zones 1 and 4, already rehabilitated and vegetated. However, there is an area of about 2 km2 in the NE margin of zone 4 covered by PS points regularly spaced (Fig. 4).

Deformation showed in both figures is very similar, suggesting dominance of vertical motion and absence of significant atmospheric effects. The only PS points with discrepant values among both figures are located along the E and W taluses of zones 2 and 3. Assuming scarce atmospheric distortions, such discrepancy seems to be related with the occurrence of horizontal movements. In order to quantitatively evaluate the quality and consistency of the InSAR time-series analysis, an inter-comparison of ascending and descending LOS displacement rates has been performed in the “Discussion” section.

In terms of velocity rates, both maps show absence of deformation in the city of Huelva and strong negative displacements in the PG stack. Here, deformation is not only intense but mostly heterogeneous, with velocity rates above 2–3 cm/year in zone 3 and NW margin of zone 4, and above 5 cm/year in zone 2 for most PS points.

Decomposition of ascending and descending LOS velocity maps into vertical and horizontal (E-W) displacement rates is represented in Fig. 4c, d. Vertical mean velocity map in Fig. 4c corroborates the existence of strong vertical negative displacement rates in the PG stack. Most values are above 5 cm/year in zone 2 and 2 cm/year in zone 3 and NW margin of zone 4, although highly deformed areas in zones 2 and 3 reach displacement rates of 16 cm/year and 3.5 cm/year, respectively. Apart from this enormous vertical movements, the horizontal (E-W) mean velocity map (Fig. 4d) indicates eastward (in blue) and westward (in red) displacement rates > 2.5 cm/year in the E and W margins of zones 2 and 3, respectively. The northeast margin of zone 4 also experiences eastward movements of this magnitude.

Figure 5 depicts the displacement time-series attributed to six selected points located in zones 2 and 3. Left graphics show the vertical motion of points P1 to P3, whereas right graphics contain the east–west deformation measured in points located near the NE (P4) and SW (P4 an P5) margins of the PG stack.

Fig. 5
figure 5

Deformation time-series retrieved from the selected points indicated at the top image. P1 to P3, vertical deformation. P4 to P6, east–west deformation

Full size image

Vertical displacement velocities are disparate in the three selected cases, ranging from 14.6 cm/year in P1 to 2.8 cm/year in P3. This difference is also illustrated in the vertical velocity map of Fig. 4c. The displacement magnitudes finally achieved are also different, with 66.8 cm in P1, 30.6 cm in P2 and 12.2 cm in P3. Additionally, the three selected points depict an almost linear trend, indicating constant rates of vertical displacement during the entire time interval analysed. This linearity has been detected virtually in all pixels from zones 2 and 3 that depict vertical motion, regardless their displacement velocities.

Regarding the east–west movements, the magnitude in the selected points oscillates between 5.3 cm of westward displacement in P4 and 14 cm of eastward displacement in P5. In general, the east–west displacement rates and magnitudes are smaller than those measured vertically, whereas the dispersion of data in the time-series is usually higher. This phenomenon, well illustrated in P4 and P6, is probably an artifact caused by the proximity between the average displacement rates calculated for these points (1.2 and 1.3 cm/year, respectively) and the accuracy of the east–west displacement time series (close to 1 cm), suggested by Pepe et al. (2016). Despite this artifact, it is possible to distinct a trend roughly linear in P4 and P6 indicative of constant westward and eastward deformation, respectively. In P5, the trend is somewhat parabolic, indicating a vague deceleration.

The spatial and temporal variations suffered by the Huelva PG stack jointly with the potential factors responsible of the reported vertical and horizontal displacements are discussed in the following section.

Discussion

Self-consistency between ascending and descending LOS measurement

In the absence of available temporal GNNS data from the study area, an appropriate way to quantitatively assess the consistency of the InSAR mean velocity measurements is the inter-comparison between ascending and descending mean displacement rates (Aslan et al. 2018, 2019; Xu et al. 2016). When pixels from both datasets fall in different locations, due to uncertainty in SAR geolocations, resolution of results is normally down-sampled and measurements are spatially mismatched in order to be compared. However, in this case, PS points obtained in the two orbits fall at the exact same locations, with millimetric accuracy. In consequence, the resampling step was not required. Figure 6 shows the comparison between ascending and descending datasets of deforming areas detected in the PG stacks and immediate surroundings. The value of the root mean square error (RMSE) was 1.34 cm/year whereas the coefficient of determination (R2) was 0.87. Both statistics reveal good agreement between ascending and descending datasets, lending consistence and confidence to the InSAR measurement.

Data pairs poorly correlated, like those laying away from the fitting correlation line in Fig. 6a, are commonly explained by insufficient acquisitions, different temporal coverage, InSAR processing errors, horizontal displacements, seasonal biases, long temporal gaps or different looking geometries (Ferretti et al. 2001; Ge et al. 2014; Xu et al. 2014). In this case, the major factor responsible of the poorly correlated pairs is the occurrence of horizontal displacements, (Fig. 4d), because the larger the east–west horizontal displacement, the higher the difference between ascending and descending LOS values. When data showing horizontal motion (i.e., those with eastward displacement rates above 0.1 cm/year and below − 0.1 cm/year) are removed from the inter-comparison, the correlation fits much better, reaching RMSE and R2 values of 0.03 cm/year and 0.99, respectively (Fig. 6b). The very high level of fitness of datasets exclusively dominated by vertical motion represents another signal of reliability of the InSAR analysis.

Conditioning factors of vertical displacement

For a better evaluation of the vertical deformation suffered by the PG stack (Fig. 4c), different profiles perpendicular to the major PG stack axis (NE-SW) have been performed (Fig. 7). These cross-sectional views cover different areas of zones 2 and 3, confronting the topographic profiles of the stack with the vertical displacement accumulated along each profile during the time interval considered in the analysis. These profiles, as well as the vertical velocity map in Fig. 4c and the vertical deformation time-series in Fig. 5, indicate that the highest values of vertical displacement occur in zone 2. Here, the accumulated vertical deformation in 4.7 years reaches 73 cm, for almost 15 cm in zone 3. Apart from the extraordinarily high value of the vertical movement, it is also important to highlight that deformation perfectly correlates with the amount of PG accumulated. In both profiles, the largest vertical deformation coincides with the maximum PG thickness, 21 m thick in zone 2 and approximately 9 m thick in zone 3. Peripherally, where the PG pile is thinner, the deformation decreases. The NE margin of zone 4 shows signal of vertical deformation that decrease peripherally, as occur in zones 2 and 3. The rest of zone 4 is virtually devoid of PS points because this zone was vegetated during a restoration plan finished in 2003.

Fig. 6
figure 6

a Inter-comparation between ascending and descending LOS displacement rates retrieved in Huelva PG stack and surroundings. b Same inter-comparison using PS points with no horizontal (E-W) motion (i.e., with eastward and westward displacements rates between 0.1 and − 0.1 cm/year).

Full size image
Fig. 7
figure 7

Cumulative subsidence versus PG stack thickness for section A–A′ in zone 2 and B–B′ in zone 3. Location of stack profiles is in Fig. 1

Full size image
Fig. 8
figure 8

a Time consolidation curve of a silty clay marsh sample from the San Francisco Bay (reproduced from Root 1950), with indication of its three evolutive stages. b Cumulative subsidence at logarithmic time, for PS points P1 to P3 in Fig. 5

Full size image

The vertical deformation detected in the Huelva PG stack can be represented as the sum of the self-weight consolidation settlement of the PG pile plus the subsidence of the unconsolidated salt-marsh substrate. In a numerical simulation of a wet PG embankment, Karoui et al. (2020) predict a maximum long-term settlement of 8 cm during 10 years of consolidation for a PG pile 10 m high. Projecting these results to the Huelva PG stack, the maximum self-weight consolidation settlement for the same time interval in zones 2 and 3 would be around 20 and 7 cm, respectively. Nevertheless, the maximum vertical movements detected in the present time-series analysis far exceed the deduced consolidation settlement in both zones (Fig. 7).

Also, the time interval considered here (4.7 years) is considerably shorter than the self-weight consolidation time predicted by Karoui et al. (2020). Moreover, in the case of zone 3, the time elapsed from the closure of the disposal activity to the oldest date considered in this study (19 years) almost doubles the predicted self-weight consolidation time deduced by these authors. All this indicate that the effect of self-weight consolidation must be minimal in the Huelva PG stack. In consequence, the measured vertical deformation is largely caused by substrate consolidation (i.e., subsidence).

In order to explore the timing of subsidence, the vertical displacement time-series of points P1 to P3 in Fig. 5 has been compared with the curve of a consolidation test performed on a silty clay marsh sample from the San Francisco Bay, USA (Root 1950). For the sake of completeness, the comparison is made in a logarithmic scale (Fig. 8). The consolidation curve of the San Francisco Bay sample fits well with the typical consolidation curves of unconsolidated and undisturbed fine-grained sediments (Das and Sobhan 2014). These curves contain three stages. The upper one defines the initial elastic compression and is related to the increasing pressure at time of initial loading. The second stage (primary consolidation) is caused by interstitial pore water pressure dissipation. Finally, the third stage (second consolidation) entails plastic adjustment of soil fabrics at constant effective stress, after excess pore water pressure dissipation. The curve in this second consolidation stage occurs at a diminishing rate, decreasing the consolidation logarithmically and becoming asymptotic to the time horizontal axis (Das and Sobhan 2014).

Although the displacement time-series and the consolidation test in Fig. 8 have been performed at different time scales, it is possible to categorise the consolidation stage currently in progress below the PG stack when comparing the shape of their curves. In that sense, the convex shape of the three vertical displacement time-series indicates that subsidence below the Huelva PG stack is still caused by primary consolidation. The absence of any signal of deceleration along the period used in the analysis strongly suggests that the primary consolidation is not nearing completion and that subsidence, especially in zone 2, is not close to end.

At the other hand, Fig. 7 also shows that subsidence under PG stacks is differential, and that is directly controlled by PG load. While differential subsidence increases the risk of cracks, sinkholes and collapses at PG stack, continued high rates of subsidence will cause the permanent sink of high volumes of PG under the water table. Both processes will increase the inflow of high-risk pollutants, including metals, metaloids and radionuclides (Bolivar et al. 1996; Pérez-López et al. 2007) from the PG stack to the Huelva estuary.

The excessive long-term consolidation can be related with the nature, essentially lutitic, of the salt-marsh substrate (Das and Sobhan 2014; Knappett and Craig 2012), but the high rates of subsidence are hardly explained only by lithology. In a recent multibeam echosound analysis performed at the right margin of the Tinto estuarine channel, just at the eastern border of the PG stack (Fig. 1), Carro et al. (2021) proved that overpresured sediments underneath the PG stack are being injected to the channel right margin, causing uplift. The release of sediments from below the PG stack to be injected into the channel floor could explain the prolonged and anomalous high subsidence in zones 2 and 3.

Conditioning factors of horizontal E-W displacement

The horizontal E-W displacement detected in the PG stack mostly includes eastward and westward movements at rates higher than 1.5 cm/year. The most striking feature in zones 2, and in minor extend in zone 3, is the symmetry of displacements with respect the main PG stack axis. According to this displacement pattern, the horizontal movements should be mostly orthogonal to the main axis, and Fig. 4d would represent the E-W component of the horizontal displacement vector.

Horizontal movements commonly accompany subsidence in areas affected by mining, consolidation of underlaying compressive layers and withdrawal of groundwater and hydrocarbons (Lee and Shen 1969). Several features are generally present in these cases, including the occurrence of differential subsidence, the absence of horizontal displacement both at the maximum subsidence points and far beyond subsiding areas and direction of horizontal movement towards maximum subsiding zones (Lee and Strauss 1969).

The horizontal movements described in the Huelva PG stack conform with this pattern although directions of horizontal displacement are outwards the main subsiding areas. As indicated in Fig. 4d, displacements to the east and west of the main axis in zones 2 and 3 are, respectively, eastward and westward. In a combined analysis of subsidence and slope stability of an area above underground mining, Diao et al. (2019) reveal that subsidence-induced cracks are able to trigger slow slippage in a surface soil layer located on a slope. Such mechanism entails horizontal movements against the main subsiding area. This mean that cracks generated by differential subsidence can generate talus displacements in unconsolidated layers, and this necessarily involves reduction of slope stability. Considering that the marsh sediments below the PG stack are unconsolidated and are strongly affected by differential subsidence, this is the most plausible explanation for the horizontal displacements detected in present study.

The NE margin of zone 4 contains PS points reporting eastward displacement rates of the same magnitude as in zones 2 and 3. In this case, the reason seems to be related with an extreme weather event. According to the horizontal displacement time-series recorded in one of this PS points (P7 in Fig. 5), no movement was detected initially, but in March 2018, the eastward displacement commenced at a constant rate of 1.75 cm/year (Fig. 9). Throughout the end of February and the beginning of March 2018, the entire Iberian Peninsula suffered the effects of Storm Emma. The first of March, the day before the equinoctial spring tide, heavy rain and wind speed of 95 km/h coincided in the study area with high tide (Agencia Española de Meteorología 2018; Puertos del Estado 2018; Table 2), generating a storm surge that caused coastal flood damages all along the Gulf of Cádiz (Agencia Española de Meteorología 2018). This weather phenomenon also triggered the destabilization of the eastern margin of zone 4, which is constantly displacing eastward since then.

Table 2 Meteorological and tidal stations used for the analysis indicating the proximity to the study area
Full size table
Fig. 9
figure 9

a Temporal evolution of the eastward displacement of a selected PS point to the NE of zone 3 (P7 in Fig. 5). b Monthly highest tides measured at the Port of Huelva during the time span selected for the present analysis. Highest tide recorded the day before Storm Emma is highlighted in red. c Monthly accumulated precipitation in the city of Huelva during the time-span selected for the present analysis. Rain accumulated during Storm Emma is highlighted in red

Full size image

Conclusions

The InSAR time-series analysis of Sentinel-1 products in the Huelva PG stack, between Oct 2017 and Jun 2021, has retrieved strong vertical displacements in zones 2 and 3. This deformation has been unambiguously attributed to subsidence, with rates > 5 cm/year in zone 2 and 2 cm/year in zone 3. The subsidence detected is strongly differential, with values of up to 16 cm/year and 3.5 cm/year in maximum subsiding areas of zones 2 and 3, respectively. In both zones, subsidence is still active 11 and 23 years after cancelation of PG storage. During the time-span covered by the InSAR time-series analysis, subsidence shows no signal of deceleration, indicating that settlement is still at its primary consolidation stage and is far from being finalised. Subsidence is also accompanied by horizontal movements linked to talus destabilization, most prominently manifested in zone 2. High and long-lasting subsidence and constant talus movements will reduce the PG stack stability, increasing the risk of collapse and contamination by underground leaching. The vulnerability of the PG stack to adverse weather conditions has been also demonstrated in this study.

The employ of time-series analysis based on SAR imagery proves to be an effective tool for comprehensively monitoring the stability and ground deformation of large waste stockpiles.

Those accumulating toxic compounds should benefit from this method in a more regular basis.

In environmental fragile areas, these studies should be even required by the administrations.