Mining–climate interactions reshape soil–water storage in a coastal sandy basin (Taolagnaro, SE Madagascar), 1993–2022

1 Mining–climate interactions reshape soil–water storage in a coastal sandy basin (Taolagnaro, SE Madagascar), 1993–2022 https://doi.org/10.5281/zenodo.17516357 Zo R. Rasoanaivo1 1 Doctoral Host Team: Geo-resources and Environment, Doctoral School: Earth and Evolutionary Sciences, University of 5 Antananarivo, Antananarivo, 101, Madagascar – ORCID: 0009-0003-0725-3764 – email to: [email protected] Abstract. Large-scale mineral-sand extraction can alter infiltration and soil–water storage in coastal dune systems. Basin-scale changes in soil–water storage (ΔSW) were evaluated in Taolagnaro, SE Madagascar, over 1993–2022. The framework integrates long-term meteorological trends (1941–2023), field measurements of hydrophysical properties, satellite-derived 10 spatial indicators, and Soil and Water Assessment Tool (SWAT+) water-balance diagnostics within a before–after/control– impact (BACI) design (1993–2007 vs. 2008–2022; mining vs. reference zones). Field data indicate a 66 % decline in saturated hydraulic conductivity from 184 ± 42 to 61 ± 15 cm h⁻¹ after mining disturbance (extraction/restoration), implying reduced infiltration. Regional warming (+0.018 °C yr⁻¹; ≈ +0.27 °C over 2009–2023) and an increase in potential evapotranspiration estimated with the Thornthwaite method (≈ +2 mm yr⁻¹) have raised atmospheric demand. Basin-wide ΔSW was near 15 equilibrium before mining (slope ≈ +0.29 mm yr⁻¹) but shifted to persistent drying after 2008 (−1.49 mm yr⁻¹). Net recharge mapped for 1993–2007 (~4,759 ha) is nearly absent in 2008–2022, while drying expands to ~8,586 ha, including ~2,132 ha within the mining zone. Under a BACI framework with counterfactual extrapolation, the mining-specific effect centres on a median offset of about −0.5 mm yr⁻¹ (area-mean over the mining zone); the mean is sensitive to interannual variability associated with El Niño–Southern Oscillation (ENSO) and the Indian Ocean Dipole (IOD). Bootstrap and threshold-sensitivity 20 analyses consistently support the post-2008 drying. Reduced infiltration capacity, together with increased evaporative demand, has shifted the basin towards persistent drying since 2008, with implications for wetland resilience, seasonal baseflow, and restoration strategies that enhance infiltration. 1 Introduction Coastal mineral-sand extraction modifies topography, soil structure and vegetation cover, with potential impacts on infiltration 25 and storage processes. In sandy dune systems, soil–water storage (ΔSW) is sensitive to structural soil changes and regional climate variability. This study investigates Taolagnaro (Fort-Dauphin), SE Madagascar, where large-scale ilmenite mining launched in 2009. We aim to: (i) quantify changes in infiltration capacity (Ksat) and soil hydrophysics; (ii) characterize regional climate trends and variability; (iii) diagnose shifts in basin-scale ΔSW before and after mining; and (iv) attribute part of the observed drying to mining using a before–after/control–impact (BACI) design formalized as a BACI counterfactual model. 30 2 2 Study area and datasets 2.1 Study area The study area is a low-relief coastal sandy basin near Taolagnaro, including a delineated mining zone and an adjacent reference zone with comparable biophysical conditions but outside mining disturbance (Figure 1). The landscape comprises stabilized dunes, interdunal wetlands, and short flashy drainage networks. 35 Figure 1: Aerial view of the Mandromondromotra watershed and the Mandena extraction site The mining footprint was manually digitized from Google Earth imagery dated 2023. The watershed boundary was delineated from the ALOS World 3D—30 m (AW3D30, JAXA) DEM accessed via OpenTopography and processed hydrologically in QGIS. CRS: WGS 84 (EPSG:4326). Base imagery: © Google 2024; imagery © Maxar Technologies and Airbus (2024). 40 Map data © 2024 Google; SIO, NOAA, U.S. Navy, NGA, GEBCO. Imagery accessed via Google Earth on 2 January 2024. 3 2.2 Climate data Station records at Taolagnaro airport (DGM, 2023; Weather&Climate, 2024) provide temperature and rainfall series for 1941– 2023. Potential evapotranspiration (PET) diagnostics use the Thornthwaite method applied to station temperatures (Thornthwaite, 1948). Annual PET summaries indicate a linear trend of ≈ +2 mm yr⁻¹ (≈ +30 mm over 2009–2023). Sensitivity 45 analyses contrast simpler formulations when this clarifies the result. 2.3 Soil and field measurements Field campaigns measure hydrophysical properties with double-ring infiltrometers and laboratory analyses (texture, porosity) at replicated plots before and after extraction/restoration. Protocols specify ring diameters, infiltration stages, initial moisture, and outlier criteria; the sampling plan balances mining and reference zones. 50 2.4 Remote sensing and spatial layers DEM (USGS, 2023) supports hydrological preprocessing and mapping, and land-cover plus ancillary rasters guide stratification and SWAT+ parameterization (Arnold et al., 2012). The workflow harmonizes all rasters to a common projection and resolution with standardized masking and nodata handling. 2.5 Hydrological model and diagnostics 55 SWAT+ defines watershed boundaries, HRUs, and annual water-balance components. The analysis summarizes annual components and extracts ΔSW from model diagnostics for temporal and spatial assessments. 3 Methods 3.1 Experimental design: phases and zones The design splits 1993–2007 (pre-mining) and 2008–2022 (post-mining) and contrasts a mining zone with a reference zone. 60 Analyses aggregate annual ΔSW by zone and phase. 3.2 Climate trend and variability analysis Analyses estimate linear temperature trends with confidence intervals; non-parametric diagnostics include Mann–Kendall and Sen’s slope together with Pettitt change-point tests. PET trends derive from temperature series, and sensitivity analyses examine alternative PET formulations. 65 4 3.3 Soil hydrophysics: measurement and synthesis Field campaigns measure Ksat at replicated sites before and after mining. Reporting includes mean ± SD, median (IQR), and Hedges g with 95% CIs. Non-parametric tests address departures from distributional assumptions. 3.4 SWAT+ configuration and verification The SWAT+ setup follows standard procedures. Simulations undergo water-balance closure checks against acceptable 70 residuals. Outputs emphasize annual components and ΔSW diagnostics, and the documentation records parameter choices for sandy coastal systems. 3.5 ΔSW definition and spatial classification The analysis defines ΔSW as the annual change in soil-water storage (mm yr⁻¹) from SWAT+ diagnostics at basin/zone scale. Each pixel–year receives a label—net recharge or net drying—by applying thresholds to ΔSW. The workflow tallies area 75 (hectares) per class by phase and zone. Uncertainty reflects threshold choice and interannual variability. 3.6 Attribution framework (BACI / counterfactual extrapolation) The attribution framework estimates the mining-specific effect under a BACI design with counterfactual extrapolation for the mining zone. The workflow constructs a ΔSW series expected without mining from pre-mining relationships and contemporaneous conditions in the reference zone. Annual deviations (observed minus counterfactual) quantify the effect. For 80 2008–2022, the median offset equals −0.5 mm ha⁻¹ yr⁻¹ and the mean equals +0.2 mm ha⁻¹ yr⁻¹, influenced by 2011 (+24 mm). A pre-mining ANOVA assesses phase comparability. 3.7 Uncertainty and sensitivity Uncertainties stem from: (i) climate data coverage at Taolagnaro airport; (ii) remote-sensing spatial and seasonal constraints (Landsat 30 m; May scenes in 2013, 2018, 2023); (iii) soil and SWAT+ parameter uncertainty (SoilGrids and legacy maps; 85 semi-distributed HRUs; lack of discharge observations for calibration); (iv) ΔSW classification thresholds and temporal aggregation; and (v) interannual variability driven by ENSO/IOD. The report presents central estimates with indicative ranges. The analysis does not produce formal 95% bootstrap confidence intervals. Sensitivity checks recompute class areas under alternative thresholds, test ΔSW trends and zone contrasts with inclusion/exclusion of each ENSO/IOD-linked outlier set and assess the effect of alternative land-cover years and soil inputs on SWAT+ outputs. 90 5 4 Results 4.1 Climate trends Temperature records indicate a significant warming trend of ≈ +0.018 °C yr⁻¹ since 1941 (Sen’s slope; 95 % CI), amounting to ≈ +1.5 °C over the full record and ≈ +0.27 °C over 2009–2023. PET increases by ≈ +2 mm yr⁻¹. Rainfall exhibits strong interannual variability consistent with ENSO/IOD phases. Figure 2 summarises these results. 95 Figure 2: Temperature trends. 4.2 Soil hydrophysics Ksat decreases from 184 ± 42 to 61 ± 15 cm h⁻¹ (−66%; Table 1). We report sample sizes (n), Hedges g, and 95% CIs; nonparametric tests confirm the shift. The magnitude accords with compaction/structure alteration after extraction and restoration. 100 Table 1: Ksat distributions prevs post-mining Plot Ksat before (cm/h) Ksat after (cm/h) Reduction (%) 1 210 ± 12 65 ± 8 69% 2 185 ± 15 55 ± 7 70% 3 240 ± 10 80 ± 9 67% 4 130 ± 14 45 ± 6 65% 5 155 ± 11 60 ± 7 61% Average 184 ± 42 61 ± 15 66% 6 4.3 ΔSW dynamics at basin scale ΔSW time series show a near-equilibrium pre-mining slope of ≈ +0.29 mm yr⁻¹ and a post-mining decline of ≈ −1.49 mm yr⁻¹. 105 The analysis reports OLS slopes with 95% confidence intervals and standard diagnostics. Figure 3 displays these fits, and the evidence indicates a structural break around 2008 consistent with the study design. Figure 3: ΔSW time series with OLS slopes (pre/post) and 95% CIs. 4.4 Spatial patterns of recharge and drying 110 From 1993–2007, net recharge covers ≈ 4,759 ha overall (≈ 1,139 ha in the mining zone), while drying spans ≈ 3,827 ha (≈ 993 ha in the mining zone). During 2008–2022, mapped recharge drops to ≈ 0 ha, whereas drying expands to ≈ 8,586 ha, including ≈ 2,132 ha within the mining zone. Figure 4 summarises these shifts. Sensitivity checks provide indicative ranges; the analysis does not produce formal 95% bootstrap confidence intervals. 115 Figure 4: Spatial distribution of cumulative ΔSW by Hydrological Response Unit (HRU) for the periods 1993–2007 (pre-mining) and 2008–2022 (post-mining). 7 4.5 Zone contrasts and attribution Pre‑mining ANOVA indicates no significant ΔSW difference between zones, which supports parallel pre‑phase behaviours. Over 2008–2022, the counterfactual analysis for the mining zone yields a median ΔSW offset of −0.5 mm ha⁻¹ yr⁻¹ and a mean 120 of +0.2 mm ha⁻¹ yr⁻¹, with 2011 exerting strong positive influence (+24 mm). Year-to-year rainfall anomalies coincide with major ENSO/IOD transitions: 1998 and 2013 are very wet years, while 2003, 2008, and 2018 are dry years; the 2008 deficit stands out despite the prevailing La Niña (Table 2). Table 2: Correlation between ENSO/IOD events and break point years 125 Break Point Year (detected) Annual Rainfall (mm) Corresponding ENSO/IOD Events Rainfall Impact relative to Interannual Mean (1585.2 mm) 1998 2059.8 End of El Niño 1997-98 (intense) & Start of La Niña 1998-00 (prolonged) & IOD Positive 1997-98 Surplus (+474.6 mm). Very wet year, potentially marking the rapid transition from a dry El Niño to a wet La Niña. 2003 1435.8 El Niño 2002-03 (moderate) Deficit (-149.4 mm). Drier year, which is consistent with the influence of an El Niño on southern Madagascar. 2008 1213.2 La Niña 2007-08 (strong) & IOD Mixed/Negative 2007-08 Very deficient (-372.0 mm). This is counter-intuitive for a La Niña (which is often associated with more abundant rainfall). This suggests other predominant local or regional factors, or a complexity of the impact of these phenomena for Taolagnaro. 2013 2203 End of La Niña 2010-12 (multiple) Very surplus (+617.8 mm). This break point marks the end of a La Niña period. It could reflect persistent wet La Niña conditions before a regime change, or a year of heavy rainfall just after the end of its influence. 2018 1434.2 End of La Niña 2017-18 & Start of El Niño 2018-19 (weak) Deficit (-151.0 mm). Coincides with the transition from a La Niña to an El Niño, often synonymous with drier conditions. 4.6 Uncertainty and sensitivity Despite uncertainties in PET formulation and model structure, the sign and magnitude of post-2008 drying are robust to reasonable choices. Excluding 2011 does not change the qualitative attribution; ranges widen but remain negative. 5 Discussion 130 5.1 Process interpretation Reduced Ksat limits soil intake during high-intensity rainfall, amplifying quick flow and curbing recharge. Warming elevates PET, further suppressing ΔSW. The near disappearance of mapped net recharge after 2008 is consistent with this mechanism. 8 5.2 Comparison with related studies Studies on mined and restored sandy systems report Ksat declines comparable to those measured here. In warm coastal 135 climates, higher PET shifts water-partitioning toward evaporation. Related work likewise reports ΔSW responses of similar magnitude in disturbed dune fields and sandy catchments undergoing vegetation and structural change. 5.3 Alternative explanations and confounders The BACI counterfactual design mitigates three confounders: land-cover change outside the mining footprint, groundwater abstraction, and large-scale climate modes. 140 5.4 Implications for water management and restoration This study recommends restoration targeting infiltration (soil loosening, organic amendments, deep-rooting vegetation) and systematic monitoring (Ksat recovery, ΔSW maps). 5.5 Reproducibility and transparency All processing steps are scripted. A reproducible repository with versioned code, configuration files, and representative inputs 145 will enable reproduction. 6 Conclusions 1) Field Ksat decreased by ~66% after mining began. 2) Regional warming (~+0.018 °C yr⁻¹) and PET rise (~+2 mm yr⁻¹) increased atmospheric demand. 3) Basin-scale ΔSW shifted from near-equilibrium (+0.29 mm yr⁻¹) to structural drying after 2008 (−1.49 mm yr⁻¹). 150 4) Net recharge area collapsed to near-zero, while drying expanded (~8586 ha; ~2132 ha within the mining zone). 5) A BACI counterfactual analysis indicates a mining-specific ΔSW deficit around −0.5 mm ha⁻¹ yr⁻¹; results are robust in sensitivity tests. Data availability Derived datasets (ΔSW rasters and CSV summaries), SWAT+ configuration files, and analysis scripts are archived at Zenodo: 155 - https://doi.org/10.5281/zenodo.16965306, - https://doi.org/10.5281/zenodo.16414860, - https://doi.org/10.5281/zenodo.16088264. Third-party raw data (Taolagnaro airport station records from DGM; USGS/Landsat; WAPOR) are not redistributed; scripts list scene IDs and retrieval steps. Correspondence and data requests: razorivo[email protected]. 160 9 Author contributions Z.R.R.: conceptualisation, methodology, investigation, formal analysis, visualisation, writing – original draft, and writing – review & editing. Competing interests The author declares no conflict of interest. 165 Acknowledgements The author thanks the Doctoral Host Team (Geo-resources and Environment, University of Antananarivo) and local agencies for facilitating field access. Any remaining errors are the author’s own. Financial support This research received no specific grant from funding agencies in the public, commercial, or not-for-profit sectors. 170