Decision-Grade Risk and Cost Mapping for Illegal Gold Mining at Crucitas, Costa Rica: Prioritisation, Phased Remediation Portfolios, and Uncertainty-Aware Policy Ranking

1. Introduction

Artisanal and illegal gold extraction in ecologically sensitive tropical landscapes presents a persistent governance challenge because environmental damages, public-health risks, and fiscal liabilities can continue to accumulate even where industrial mining is formally prohibited (Intergovernmental Forum on Mining, Minerals, Metals and Sustainable Development [IGF], 2019; Hattingh et al., 2021). In such settings, policy debate often polarizes around normative positions (prohibition versus regularisation) while lacking quantified evidence on three decision-critical questions: where impacts are concentrated, what remediation would cost over time, and how alternative legislative pathways shift net public liabilities once enforcement, residual risk, and potential revenue or offset mechanisms are treated explicitly rather than implicitly.

Artisanal and informal gold mining has expanded rapidly in tropical forest regions; in the Peruvian Amazon, rising gold prices have been linked to mining-driven deforestation and to increased national mercury imports used in artisanal extraction (Swenson et al., 2011). In remote frontiers, multi-sensor satellite archives enable multi-decadal mapping of artisanal mining expansion and can be paired with satellite-derived water-siltation metrics (e.g., total suspended solids, TSS) to support consistent monitoring (Lobo et al., 2016). Site-scale remote-sensing studies also note that limited information on the magnitude of damage and pollution can hinder remediation and mine-site formalization (Adamek et al., 2021). Because mercury and other ASGM contaminants occur across multiple media and exposure pathways, standardized guidance recommends integrated site characterization and sampling designs that explicitly link contaminant profiles to exposure and health-impact analyses (World Bank, 2021).

The Crucitas mining landscape (Cutris, northern Costa Rica) exemplifies this tension in a transboundary, connectivity-relevant setting near the San Juan River and embedded within the San Juan–La Selva corridor context, where land-use decisions can propagate beyond the immediate footprint through hydrologic and ecological pathways (Monge & Chassot, 2009). Public technical commentary from the period of the halted industrial project emphasized that the proposed open-pit configuration would have been situated in a high-rainfall environment with shallow groundwater and plausible downstream exposure routes for sediment and process-related contaminants (Astorga, 2009). Contemporaneous reporting documented how quickly irreversible landscape change could occur under enabling decisions, including rapid forest clearing and the prospect of much larger authorized disturbance envelopes (Biamonte, 2011). These attributes motivate an assessment approach that is spatially explicit, hydrologically attentive, and capable of representing how liabilities scale under alternative access and extraction dynamics.

Independent field evidence further indicates that illegal alluvial practices in the Crucitas setting are associated with contamination mechanisms directly relevant to hydrologic prioritisation and remediation design. Field sampling and synthesis have reported mercury contamination concerns and acid-drainage processes linked to illegal mining activity, underscoring that the governance problem is operationally environmental and can translate into persistent, water-facing liabilities (Gazel et al., 2019).

Current policy debate in Crucitas is not well represented by a binary “mine vs. no mine” framing. Public reporting describes multiple concurrent pathways spanning (i) concession-based industrial regularisation with structured state participation, (ii) regulated artisanal/cooperative arrangements, and (iii) non-extractive “development pole/geopark” designs that seek to finance recovery and livelihoods through institutional finance mechanisms and remediation-linked gold recovery rather than opening a new extraction pathway (Martínez, 2024). Coverage further suggests that legal eligibility under the executive pathway could extend across the full district of Cutris, implying that policy design may expand exposure and tail-risk considerations well beyond the historic concession geometry (Pomareda García, 2025). In parallel, investigative reporting indicates that illegal-mining supply chains may displace processing impacts beyond the extraction footprint through off-site informal processing nodes, complicating “site-only” remediation accounting unless explicitly parameterised in sensitivity analysis (Esquivel Solano & Rivera, 2024). Related synthesis has flagged inconsistencies in gold-export statistics and traceability narratives, reinforcing that enforcement effectiveness, leakage, and provenance controls are not background assumptions but central determinants of pathway performance (Delfino, 2025).

Despite rapid advances in Earth observation and environmental risk mapping, two gaps limit decision relevance in contested mining landscapes. First, remote-sensing products are frequently used to demonstrate disturbance or land-cover change, but are seldom translated into decision-grade prioritisation surfaces that link spatial signals to actionable intervention packages and phased implementation (Maus et al., 2022). Second, economic assessments often report aggregate damages or remediation costs without a transparent mechanism for coupling spatial prioritisation to present-value (PV) financial planning and for comparing policy pathways on a common fiscal basis—particularly when long-horizon closure and water-treatment liabilities dominate cost profiles and are sensitive to discounting, updating, and tail-risk treatment (Chambers, 2024; Bocoum et al., 2021). These gaps constrain practical decision-making, including: (i) estimation of remediation liability under current landscape conditions, (ii) identification of priority intervention locations and sequences, and (iii) evaluation of how specific legislative pathways shift expected net PV liability once enforcement costs, residual risk, and potential revenue or offset terms are incorporated. Prior work on Crucitas has largely emphasized environmental evidence and governance/enforcement dynamics (Monge & Chassot, 2009; Astorga, 2009; Gazel et al., 2019), while many Earth-observation applications for artisanal/illegal gold mining remain oriented toward detection/monitoring and disturbance mapping rather than translating spatial signals into phased, budgetable intervention portfolios (Moomen et al., 2022); in parallel, mine-closure and financial-assurance literature develops PV-consistent approaches for long-horizon liabilities (notably water-treatment and post-closure obligations) but typically without explicit coupling to EO-derived prioritisation in contested illegal-mining settings (Chambers, 2024). This study contributes a transferable, decision-grade framework that integrates Earth-observation signals with field-anchored screening to deliver spatial prioritisation, phase-timed remediation portfolios, and uncertainty-aware PV liability and policy ranking in an auditable, reproducible workflow.

This study presents an integrated, multi-model framework for the Crucitas landscape that converts remote-sensing change signals into phased financial liabilities and policy-pathway comparisons under uncertainty. Three research questions guide the analysis: (1) Which grid cells exhibit the strongest multi-domain signals consistent with intervention prioritisation in land, water, and hydrology dimensions? (2) What are the PV-consistent gross and net public liabilities of a phase-timed remediation portfolio derived from those prioritisation surfaces? (3) How do alternative policy pathways shift the distribution of net PV liabilities and exceedance risk under unit-cost uncertainty?

The framework comprises five linked models. Model 1 constructs baseline-to-recent change proxies in terrestrial disturbance indicators and riparian water proxies using consistent windows (baseline January 2019–December 2020; recent February 2024–January 2025). Model 2 integrates these signals with terrain and hydrological predictors to compute component indices (land, water, hydrology) and a combined 0–100 prioritisation surface at grid-cell resolution. Model 3 provides field-anchored screening by evaluating whether remotely sensed and terrain-derived predictors recover observed variability in mercury-related indicators under leave-one-out cross-validation, thereby bounding interpretation of remote-sensing–only screening. Model 4 translates the prioritisation surface into a phase-timed remediation portfolio with explicit intervention packages and unit-cost envelopes, computing PV gross and net liabilities and decomposing costs by phase and package; policy pathways are implemented as model-ready levers that modify costs, revenues/offsets, timing, and residual provisioning. Model 5 propagates key uncertainties through Monte Carlo simulation to quantify exceedance risk and evaluate how policy pathways shift the full distribution of net PV liabilities. Collectively, the framework provides an end-to-end translation from environmental signals to fiscal decision variables, enabling transparent comparison of legislative pathways on a common PV and probabilistic-risk basis. An overview of the integrated framework is provided in Figure 1.

2. Methodology

This section describes the geospatial data acquisition and Google Earth Engine–based export pipeline used to generate ring-compatible predictor stacks for downstream analysis.

2.1. Study Area Definition and Ring-Based Analysis Zones

Spatial inference is anchored at the Crucitas reference location (centroid of the project feature) and evaluated using concentric annular zones (“rings”) to support distance-structured comparisons (Gorelick et al., 2017). Rings are defined as a donut-shaped region between an inner and outer radius, ensuring that feature extraction remains comparable across distance intervals. To ensure transparency and reproducibility of the spatial partitioning used for downstream modeling, ring geometry is formalized in Equation (1):

$${\Omega }_{r_1:r_2} = \mathit{Buffer}(x_0,r_2) \setminus \mathit{Buffer}(x_0,r_1)$$

where:

Ω_r1:r2 = annular analysis zone between inner radius r₁ and outer radius r₂

x₀ = Crucitas reference point (centroid)

Buffer(x₀,r) = circular buffer of radius r around x₀

r₁,r₂ = ring radii (e.g., 1–5 km, 5–10 km, 10–20 km)

\setminus = set-difference operator (outer minus inner)

2.2. Context Covariates: Topography and Rasterized Hydrographic/Settlement Layers

A multiband context stack is generated to provide consistent spatial covariates for interpretation and downstream modeling, including elevation-derived terrain metrics and rasterized vector context layers. Elevation is sourced from NASADEM, from which slope and hillshade are derived (NASA Jet Propulsion Laboratory, 2020). Hydrographic network lines, settlement points, and administrative boundaries are rasterized into binary presence layers to provide structured contextual predictors at a consistent footprint for ring-based feature extraction (Gorelick et al., 2017).

2.3. Sentinel-2 Indices and Monthly Compositing (TIF23 Disturbance Stack; TIF22 Water-Proxy Stack)

Sentinel-2 Surface Reflectance imagery is filtered to the export footprint and screened using scene classification to reduce cloud and shadow contamination (Drusch et al., 2012; European Union/ESA/Copernicus, n.d.). Disturbance-sensitive indices are computed from surface reflectance, including NDVI, NDMI, and BSI (Equations (2)–4). Monthly composites are constructed using pixel-wise median statistics, producing a fixed-schema, time-indexed multiband stack for disturbance characterization.

$$\mathit{NDVI}={\displaystyle \frac{{\rho }_{\mathit{NIR}}-{\rho }_{\mathit{RED}}}{{\rho }_{\mathit{NIR}}+{\rho }_{\mathit{RED}}}}$$

where:

ρ_NIR is near-infrared surface reflectance (Sentinel-2 B8).

ρ_RED is red surface reflectance (Sentinel-2 B4).

$$\mathit{NDMI}={\displaystyle \frac{{\rho }_{\mathit{NIR}}-{\rho }_{\mathit{SWIR}1}}{{\rho }_{\mathit{NIR}}+{\rho }_{\mathit{SWIR}1}}}$$

where:

ρ_SWIR1 is shortwave infrared 1 surface reflectance (Sentinel-2 B11).

ρ_NIR is near-infrared surface reflectance (Sentinel-2 B8).

$$\mathit{BSI}= {\displaystyle \frac{({\rho }_{\mathit{SWIR}1}+{\rho }_{\mathit{RED}})-({\rho }_{\mathit{NIR}}+{\rho }_{\mathit{BLUE}})}{({\rho }_{\mathit{SWIR}1}+{\rho }_{\mathit{RED}})+({\rho }_{\mathit{NIR}}+{\rho }_{\mathit{BLUE}})}}$$

where:

ρ_BLUE is blue surface reflectance (Sentinel-2 B2).

ρ_RED,ρ_NIR,ρ_SWIR1 are as defined above.

For water-proxy monitoring, a riparian corridor mask is used to restrict analysis to areas proximate to streams. A turbidity/suspended-sediment proxy is computed using NDTI (Equation (5)) (European Union/ESA/Copernicus, n.d.). Stable surface water is identified using the Global Surface Water occurrence product (Pekel et al., 2016) to constrain interpretation to persistent water bodies.

$$\mathit{NDTI}={\displaystyle \frac{{\rho }_{\mathit{RED}}-{\rho }_{\mathit{GREEN}}}{{\rho }_{\mathit{RED}}+{\rho }_{\mathit{GREEN}}}}$$

where:

ρ_GREEN is green surface reflectance (Sentinel-2 B3).

ρ_RED is red surface reflectance (Sentinel-2 B4).

NDTI is treated as an optical proxy for relative turbidity/suspended-sediment conditions and is not a direct contaminant concentration metric.

2.4. Hydrology and Erosion-Driver Layers (Ring-Compatible Footprint)

Hydrologic and terrain-based predictors are derived from MERIT Hydro and NASADEM, including upstream-area proxies, distance-to-stream, and related terrain metrics (Yamazaki et al., 2019; NASA Jet Propulsion Laboratory, 2020). All layers are exported over a buffer that fully covers the maximum ring extent to prevent outer-ring truncation and to ensure consistent sampling support across Ωr₁:r₂.

2.5. Integrated Impact-Driver Stack (TIF24) and Baseline-to-Recent Change Metrics

A consolidated impact-driver GeoTIFF integrates multi-sensor indicators relevant to land disturbance and environmental change, including Sentinel-2 baseline versus recent composites (Drusch et al., 2012; European Union/ESA/Copernicus, n.d.), Sentinel-1 VV/VH backscatter changes (Torres et al., 2012; European Union/ESA/Copernicus, n.d.), precipitation sums from CHIRPS (Funk et al., 2015), Dynamic World bare-ground indicators (Brown et al., 2022), and forest-loss indicators from the Hansen global forest change lineage (Hansen et al., 2013). Change features are operationalized as signed baseline-to-recent differences (Equation (6)), ensuring consistent directionality across indices and facilitating cross-layer interpretation.

$$\Delta V=V_{\mathrm{recent}}-V_{\mathrm{base}}$$

where:

V_base is the baseline-window statistic (e.g., median over 2019–2020).

V_recent is the recent-window statistic (e.g., median over 2024–2025).

ΔV is the signed change feature used for screening (e.g., ΔNDVI, ΔBSI, ΔVV).

2.6. Export Conventions and Reproducibility Controls

All outputs are exported as GeoTIFFs using consistent coordinate reference settings and deterministic band naming conventions to enable automated feature extraction (Gorelick et al., 2017). Monthly stacks enforce a stable schema across time by preserving masked values in periods with insufficient observations, preventing column drift during downstream ingestion. Vector-derived context products are rasterized at analysis scale to minimize inconsistencies between gridded predictors and ring extraction. A complete inventory of geospatial inputs, Earth Engine products, temporal coverage, spatial resolution, and exported file schemas is provided in Supplementary Table S1. Following raster export, predictors are operationalized as AOI-level feature matrices via zonal reduction to support consistent downstream inference.

2.7. Zonal Feature Extraction over Analysis Zones

To generate ring-compatible predictors for downstream analysis, raster stacks are summarized over each analysis zone Ω (e.g., administrative unit, buffer, or annulus) using zonal reduction. For border-adjacent contexts, analysis zones are spatially constrained to the national boundary prior to reduction to ensure that extracted values reflect only the intended jurisdictional footprint. Zonal features are computed by applying a reducer R_j to the distribution of valid pixel values within Ω, formalized as Equation (7):

$$s_{k,j}(\Omega ,t)={\mathrm{R}}_j\left(\left\{V_k(x,t) | x\in \Omega , \mathrm{M}(x,t)=1\right\}\right)$$

where:

s_k,j(Ω,t) is the summary statistic j for predictor band k evaluated over analysis zone Ω at time t (or for a static composite when t is not time-indexed).

V_k(x,t) is the raster value of the predictor band k at pixel location x and time t.

Ω is the analysis zone used for aggregation (e.g., Ωr₁:r₂ as defined in Equation (1)).

M(x,t)∈{0,1} is the validity mask (1 for retained pixels after QA/masking; 0 otherwise).

R_j(⋅) is the reduction operator defining statistic j (e.g., mean, median, percentiles, standard deviation, min–max, sum, count), applied to the set of valid pixels within Ω.

This formulation supports both time-collapsed composites (e.g., baseline vs. recent windows and their deltas) and time-indexed monthly composites, while maintaining a stable feature schema under missing-observation months via explicit masking conventions described in Section 2.6.

2.8. Pixel-Wise Temporal Feature Construction from Monthly Predictor Stacks

Model 1 compresses monthly predictor stacks into a fixed set of pixel-wise temporal descriptors used for downstream mapping and modelling. For each predictor band k and pixel x, descriptors include distributional moments over the available series, a linear trend term estimated over valid observations, windowed means for baseline and recent periods, and a baseline-to-recent change term. All descriptors are computed using validity masking to ensure consistency under missing observations (masked months) and optional AOI masking. The full mathematical specification of the temporal feature vector and the trend estimator is provided in the Supplementary material (equation-compendium table), and is referenced later in the Methodology to avoid repeated cross-referencing.

2.9. Model 2: Impact Scoring via Robust Normalization and Percentile Aggregation

Model 2 converts Model 1 predictors into spatially comparable impact scores by: (i) computing grid-cell zonal summaries, (ii) transforming predictor deviations using robust dispersion statistics, (iii) mapping deviations to empirical percentiles to obtain unitless comparability across predictors, and (iv) aggregating predictor percentiles into grouped indices (and a combined index) on a 0–100 scale.

Grid-cell predictor summaries are computed via zonal reduction (Section 2.7; Equation (7)) using the mean reducer. For each predictor band k and grid cell g, robust standardized deviation is computed relative to the distribution across grid cells using the median and median absolute deviation (MAD) Equation (8):

$$z_k(g)={\displaystyle \frac{\left|{\bar{V}}_k\left(g)-{\mathit{median}}_{g’}\right\{{\bar{V}}_k(g’\left)\right\}\right|}{1.4826\hspace{0.17em}{\mathit{MAD}}_{g’}\left\{{\bar{V}}_k\right(g’\left)\right\}}}$$

where:

${\bar{V}}_k\left(g\right)$ = grid-cell zonal mean for predictor k (computed using Equation (7) with mean reducer)

G_k = set of grid cells with finite ${\bar{V}}_k$(⋅) for predictor k

MAD{${\bar{V}}_k$}=median(∣${\bar{V}}_k$−median(${\bar{V}}_k$)∣) over g′∈G_k

1.4826 scales MAD to be consistent with the standard deviation under normality

Robust deviations are then mapped to a unit interval using the empirical percentile (empirical CDF) across grid cells Equation (9):

$$p_k(g)={\displaystyle \frac{1}{\left|G_k\right|}}\sum _{g’\in G_k}I\left(z_k\left(g’){ \le z}_k\right(g)\right)$$

where:

p_k(g)∈[0,1] = percentile rank of grid cell g for predictor k

F_k(⋅) = empirical CDF of z_k over G_k

Grouped indices are computed as the mean percentile rank across predictors assigned to group GGG, scaled to 0–100 Equation (10):

$$I_G(g)=100\cdot {\displaystyle \frac{1}{\left|K_G\right|}}\sum _{k\in K_G}{p}_k\left(g\right)$$

where:

I_G(g)∈[0,100] = group-level index for group G at grid cell g

K_G = set of predictors assigned to group G (e.g., land, water, hydrology)

Finally, the combined index is computed as the mean across the set of groups S used in the analysis Equation (11):

$$I_{\mathrm{COMB}}(g)={\displaystyle \frac{1}{\left|S\right|}}{\sum }_{G\in S}{I}_G\left(g\right)$$

where:

I_COMB(g)∈[0,100] = combined impact-prioritization score

S = set of included groups (e.g., {LAND,WATER,HYDRO}

2.10. Scientific Validation and Triangulation Using Field Sediment Observations

To evaluate whether mapped impact patterns are directionally consistent with independent field evidence, sediment sampling points (Hg and co-measured metals, when available) are used for triangulation. This step does not infer causality; it tests whether field concentrations exhibit spatial gradients that are qualitatively compatible with model-derived scores.

Distance from each sampling location to the mine reference point is computed in a metric projection Equation (12):

$$d(x)={\displaystyle \frac{1}{1000}}\hspace{0.17em}{\Vert x-x_m\Vert }_2$$

where:

d(x) = distance from sampling location x to mine reference point x_m (km)

∥⋅∥₂ = Euclidean norm in projected coordinates (m); division by 1000 converts to km

Sampling points are then assigned to annular distance classes (e.g., 0–5, 5–10, 10–20 km) for near–far comparisons. Ring-stratified summaries of Hg are computed to support ring-wise reporting, and directional agreement is assessed by relating ring-level Hg summaries to model-derived ring/AOI scores (association tests reported as correlation metrics).

2.11. Model 3: Field-Based Triangulation Using Multi-Scale Raster Sampling

Model 3 calibrates the remote-sensing predictor space against independent field observations by extracting multi-scale neighborhood summaries from candidate predictor rasters at each sampling location and fitting a Random Forest regressor under leave-one-out cross-validation (LOOCV). Preprocessing (median imputation and robust scaling) is performed fold-internally. Dimensionality is controlled via coverage screening and correlation-based prefiltering (TOPK).

2.11.1. Multi-Scale Feature Extraction at Field Locations

For each sampling location x_i and predictor raster V_k, neighborhood-aggregated means are computed over circular windows of radius r (in the raster’s native CRS/grid) Equation (13):

$${\mu }_{k,r}(x_i)={\displaystyle \frac{1}{N_{i,k,r}}}\sum _{u\in W_r\left(x_i\right)}{m}_k\left(u)\hspace{0.17em}{V}_k\right(u)$$

where:

μ_k,r(x_i) = neighborhood mean of predictor k around site i using radius r

W_r(x_i) = set of raster pixels u within radius r of x_i (in the raster’s native grid/CRS)

m_k(u)∈{0,1} = validity mask (1 if finite and not nodata; 0 otherwise)

N_i,k,r=∑_u∈Wr(xi)m_k(u) = number of valid pixels in the window

V_k(u) = raster value for predictor k at pixel u

Local variability is summarized using the neighborhood standard deviation Equation (14):

$${\sigma }_{k,r}(x_i)=\sqrt{{\displaystyle \frac{1}{N_{i,k,r}}}\sum _{u\in W_r\left(x_i\right)}{m}_k\left(u)\hspace{0.17em}(V_k\right(u)-{\mu }_{k,r}\left(x_i\right){)}^2}$$

where:

σ_k,r(x_i): neighborhood standard deviation of predictor k around site i at radius r

μ_k,r(x_i): neighborhood mean (Equation (13))

W_r(x_i): pixels within radius r of x_i

m_k(u): validity mask (1 valid, 0 nodata/invalid)

N_i,k,r: number of valid pixels in the window

V_k(u): raster value at pixel u

To represent local gradients, a contrast feature is computed between an inner and outer radius (r₁<r₂) Equation (15):

$$C_{k,r_1,r_2}\left(x_i)={\mu }_{k,r_1}\right(x_i)-{\mu }_{k,r_2}(x_i), r_1<r_2$$

where:

C_k,r1,r2(x_i): contrast feature for predictor k at site i between radii r₁ and r₂

μ_k,r1(x_i), μ_k,r2(x_i) neighborhood means at inner/outer radii (Equation (13))

r₁<r₂: inner and outer window radii (same raster grid/CRS)

2.11.2. Predictor Screening and Dimensionality Control

Candidate rasters are retained only if they meet a minimum extraction-coverage requirement across field locations Equation (16):

$$I_k=I\left(\sum _{i=1}^{n}I ({\mu }_{k,r_0}(x_i)\in R)\ge n_{\mathit{min}}\right)$$

where:

I_k∈{0,1} indicates whether predictor k passes the coverage screen

μ_k,r0(xi): extracted point-support value for predictor k at site i (Equation (13) with radius r₀)

r₀ = point-support (pixel-scale) extraction radius

n = number of field sites with valid target measurements

n_min = minimum required number of valid extractions across sites

I(⋅): indicator function (1 if condition holds, 0 otherwise)

Given the high dimensionality of predictors, a fast filter selects the top K predictors by absolute correlation with the calibration target y (e.g., log₁₀(Hg) and/or a composite contamination index), computed over available paired observations Equation (17):

$$K^{\star }=\arg \underset{K, |K|=K_0}{\max }\sum _{m\in M}\left|\mathit{corr}\left(X_p,y\right)\right|$$

where:

P = set of screened predictors after Equation (16)

K⋆ = selected subset of predictors

K₀: target subset size (number of predictors retained)

X_p = feature vector for predictor p across sampled points (from Equations (13)–15)

y: calibration target (e.g., log₁₀(Hg) and/or ContamIndex)

corr(⋅,⋅) = Pearson correlation computed on pairwise-complete observations

2.11.3. Cross-Validated Calibration and Significance Screening

Out-of-sample performance is estimated using LOOCV. For each point i, the prediction is generated by a model trained on all other points, with preprocessing (imputation and robust scaling) performed within each training fold Equation (18):

$${\widehat{y}}_i=f^{(-i)}(x_i), i=1,\dots ,n$$

where:

${\widehat{y}}_i$: LOOCV prediction for site i

x_i: feature vector at site i (constructed from Equations (13)–15)

f⁽⁻ⁱ⁾(⋅): model fit using all sites except i (with fold-internal preprocessing)

n: number of field sites used in calibration

To test whether observed predictive skill exceeds chance structure, a permutation test is applied by shuffling targets and recomputing LOOCV R₂ across B permutations Equation (19):

$$p={\displaystyle \frac{1+{\sum }_{b=1}^{B}I\left(R_b^2{\backslash \mathit{geR}}_{\mathit{obs}}^2\right)}{B+1}}$$

where:

p: permutation-test p-value for predictive skill

R_obs²: observed LOOCV R² using the true targets

R_b²: LOOCV R² from permutation b (targets shuffled)

B: number of permutations

I(⋅): indicator function (1 if condition holds, 0 otherwise)

Model 3 extracts multi-scale predictor features at field locations, produces LOOCV predictions for the calibration target(s), reports model-based importance summaries for interpretability, and saves the preprocessing-plus-regressor artifacts to enable reproducible reruns under fixed seeds and documented settings. Candidate predictors comprise disturbance and hydro-geomorphic rasters from the main export pipeline, complemented by a fixed core pack of static covariates (soil SOC/pH/clay, surface-water occurrence, MERIT Hydro proxies, slope, and WorldCover tree/wetland masks) exported as single-band GeoTIFFs under deterministic naming conventions to support automated feature extraction. Implementation details that are not essential for the main narrative—such as exact feature definitions, extraction-coverage screening and correlation-based prefiltering for dimensionality control, and the LOOCV/permutation-testing formulas with symbol definitions—are consolidated in Supplementary Table S3.

2.12. Model 4: Phase-Timed Financial Costing and Present-Value Screening

Model 4 converts spatial impact scores into a (PV) liability profile at grid-cell resolution by: (i) assigning each cell to an implementation phase using ranked score cut-offs; (ii) mapping dominant impact groups to intervention packages and associated bill-of-quantities measures; (iii) discounting one-time capital expenditures (CAPEX) to the start of the assigned phase; (iv) discounting recurring operating expenditures (OPEX) as phase-activated streams; and (v) netting eligible credits on a PV-consistent basis (credits may exceed costs, but net PV liabilities are bounded at zero).

Cells are assigned to implementation phases using score-ranked cut-offs. Let π_i denote the normalized descending rank of the total score s_i (higher score ⇒ smaller π_i), and let (τ1,τ2,τ3)=(0.30,0.60,1.00) denote cumulative phase thresholds. The phase assignment is Equation (20):

$$\mathrm{p}\mathrm{h}\mathrm{a}\mathrm{s}\mathrm{e}(i)=\left\{\begin{array}{ll}1& if{\pi }_i\le {\tau }_1,\\ 2& if{\tau }_1<{\pi }_i\le {\tau }_2,\\ 3& if{\tau }_2<{\pi }_i\le {\tau }_3,\end{array}\right.({\tau }_1,{\tau }_2,{\tau }_3)=(\mathrm{0.10,0.30,1.00})$$

where:

phase(i) = assigned implementation phase for cell i

π_i = normalized descending rank of s_i (e.g., Nπi=rank↓(s_i)/N, with N evaluated cells)

τ_p =cumulative phase thresholds (top 30%, top 60%, all cells)

A common discount factor is used throughout Model 4 Equation (21):

$$\mathit{DF}(t)={\displaystyle \frac{1}{(1+r{)}^t}}$$

where:

DF(t) = discount factor at t years from the reference year

r = real annual discount rate

One-time intervention costs (CAPEX) are phase-timed to the start of the assigned phase and discounted accordingly. Let t_i^start denote the start year (in years from the reference year) of the phase assigned to cell i, and let M denote the set of one-time measures in the bill of quantities Equation (22):

$${\mathit{CAPEX}}_i^{\mathit{PV}}=\mathit{DF}\left(t_i^{\mathit{start}}\right)\sum _{m\in M}{q}_{i,m}\hspace{0.17em}{c}_m$$

where:

CAPEX_i^PV = PV of one-time costs for cell i

q_i,m = quantity of measure m applied in cell i (e.g., derived from cell area/length and package rules)

c_m = unit cost of measure m

Total PV cost and net PV cost (after credits) are computed as follows Equation (23):

$${\mathit{TC}}_i^{\mathit{PV}}={\mathit{CAPEX}}_i^{\mathit{PV}}+{\mathit{OPEX}}_i^{\mathit{PV}}, {\mathit{NC}}_i^{\mathit{PV}}=\mathit{max} \left(0,{\mathit{TC}}_i^{\mathit{PV}}-{\mathit{CR}}_i^{\mathit{PV}}\right)$$

where:

TC_i^PV = total PV cost for cell i

NC_i^PV = net PV cost for cell i after credits (bounded below by zero via max(0,·))

OPEX_i^PV = PV of recurring costs (phase-activated; e.g., monitoring and water-related operations where applicable)

CR_i^PV = PV of eligible credits applied to cell i

To keep the main manuscript concise while preserving reproducibility, the explicit stream formulations for (i) phase-activated OPEX streams and (ii) PV-consistent finite credit streams are provided in Supplementary Table S3.

Policy pathway overlays are implemented as additive PV adjustments applied to the portfolio-level sum ∑_iNC_i^PV, with revenue and offset terms treated as negative costs. For Exp. 24.717, the royalty term is parameterized as a minimum 5% take on annual gross sales (“ventas anuales”) and discounted over the PV horizon (see Supplementary Table S7 for the pathway lever definitions and uncertainty structure).

2.13. Model 5: Sensitivity Analysis, Monte Carlo Uncertainty, Exceedance Risk, and 2060 Projections

Model 5 quantifies uncertainty in the phase-timed present-value (PV) liability estimates by propagating uncertainty in key unit-cost drivers through (i) one-way sensitivity screening and (ii) Monte Carlo simulation. Uncertainty is summarized using distributional percentiles of the portfolio net PV, and communicated to decision-makers using exceedance probabilities for policy-relevant budget thresholds. In addition, a reporting-only indexation projection expresses 2026-valued summary metrics in 2060 terms under alternative annual escalation assumptions.

To preserve internal consistency with Model 4, all PV streams in Model 5 use the same discount rate and the same phase activation logic (phase start years and phase durations). For a constant annual cost rate aaa that begins at year t₀ (years from the reference year) and persists for Y years, the PV is Equation (24):

$$\mathit{PV}(a;\hspace{0.17em}{t}_0,Y)=a\sum _{k=1}^Y{\displaystyle \frac{1}{(1+r{)}^{\hspace{0.17em}{t}_0+k}}}$$

where:

a= annual cost rate (USD·y⁻¹)

$t_0$ = start year of the active phase (years from the reference year)

Y= number of active years in the phase

r = real annual discount rate

k = year index within the active period (1,…,Y)

Uncertain unit costs are represented using a triangular distribution parameterized by low, central (mode), and high values. For unit-cost item m and Monte Carlo draw b Equation (25):

$$c_m^{\left(b\right)}\sim \mathit{Triangular}\left(c_m^L,\hspace{0.17em}{c}_m^C,\hspace{0.17em}{c}_m^H\right), b=1,\dots ,B$$

where:

cm(b) = sampled unit cost for item m in draw b

c_m^L,c_m^C,c_m^H = low, central (mode), and high parameters

B = number of Monte Carlo draws

For each draw b, uncertainty is propagated by sampling unit-cost ratios and applying draw-specific multipliers to the Model 4 present-value cost components at the grid-cell level; PV-consistent credits are then subtracted and non-negativity is enforced per cell Equation (26):

$$T{C}_i^{\left(b\right)} ={\lambda }_{\mathit{CAPEX}}^{\left(b\right)}\hspace{0.17em}\mathit{CAPE}{X}_i^{\mathit{PV},C}+{\lambda }_{\mathit{OPEX}}^{\left(b\right)}\hspace{0.17em}\mathit{OPE}{X}_i^{\mathit{PV},C}, N{C}^{\left(b\right)} =\sum _{i=1}^N\left(0, T{C}_i^{\left(b\right)}-C{R}_i^{\mathit{PV}}\right) .$$

where:

$N{C}^{\left(b\right)}$= simulated total net present-value (PV) cost (USD) in Monte Carlo draw b.

$T{C}_i^{\left(b\right)}$= simulated gross PV cost (before credits) for grid cell i in draw b.

i=1,…,N indexes grid cells

c(b) = vector of sampled unit costs in draw b

$\mathit{CAPE}{X}_i^{\mathit{PV},C}$,$\mathit{OPE}{X}_i^{\mathit{PV},C}$ are the central (Model 4) PV CAPEX and PV OPEX for cell i, computed under Model 4 timing (OPEX discounted via Equation (24)).

CR_i^PV = PV credits for cell i (held fixed when credits are already phase-timed in Model 4)

${\lambda }_{\mathit{CAPEX}}^{\left(b\right)}\hspace{0.17em}$,)${\lambda }_{\mathit{OPEX}}^{\left(b\right)}$ = draw-specific unit-cost multipliers (ratios relative to central unit-cost values) applied to CAPEX and OPEX PV components.

max(0,⋅)=enforces non-negative net liabilities at the grid-cell level (no negative net PV costs).

Exceedance risk for a policy-relevant budget threshold T is estimated empirically from Monte Carlo samples Equation (27):

$$\widehat{P}\left(\mathit{NC}>T\right)={\displaystyle \frac{1}{B}}\sum _{b=1}^{B}I\left({\mathit{NC}}^{\left(b\right)}>T\right)$$

where:

$\widehat{P}\left(\mathit{NC}>T\right)$ = estimated probability that portfolio net PV exceeds T

T= threshold (USD, PV terms)

I(⋅) = indicator function (1 if true; 0 otherwise)

For long-horizon reporting only (not for PV accounting), 2026-valued summary metrics (e.g., mean, P10, P50, P90 of {NC(b)}) are expressed in 2060 terms using an annual indexation rate γ Equation (28):

$${\mathit{NC}}_{2060}(\gamma )={\mathrm{NC}}_{2026}\hspace{0.17em}(1+\gamma {)}^{2060-2026)}$$

where:

NC₂₀₂₆ = selected summary metric expressed in 2026 USD terms

NC₂₀₆₀(γ) = corresponding indexation projection to 2060

γ = assumed annual indexation rate (dimensionless)

Implementation details not essential to the main narrative—such as the full pixel-wise temporal descriptor definitions, the explicit empirical-percentile formulation, and the phase-activated PV stream definitions for recurring costs and credits—are consolidated in Supplementary Table S3 to preserve reproducibility while keeping the main Methods decision-focused.

Model 4 and Model 5 are parameterised using external quantitative anchors and policy levers compiled in Supplementary Tables S4–S7. These tables provide (i) Costa Rica–specific credit benchmarks (e.g., PSA schedules) used to construct PV-consistent credit streams, (ii) monitoring tariff anchors used to parameterise recurring OPEX intensities, (iii) mine-water treatment and water-recovery cost envelopes used to define scenario bounds for water-related OPEX (where applied), (iv) inflation/indexation conventions used for consistent long-horizon reporting to 2060, and (v) policy pathway levers—including the Exp. 24.717 minimum royalty take (≥5% of annual gross sales), enforcement costs, implementation lags, and residual tail-risk terms—used to implement delta-versus-baseline comparisons and define uncertainty ranges propagated in Model 5.

2.14. Limitations

The proposed workflow should be interpreted as a reproducible, spatially explicit screening and costing framework rather than a definitive contaminant quantification system. First, several predictors are remote-sensing proxies (e.g., vegetation condition, bare soil, SAR backscatter, optical turbidity indices) and hydrology-derived surrogates rather than direct measurements of Hg or dissolved metals; therefore, mapped outputs represent relative disturbance and exposure signals contingent on proxy validity and local process coupling. Second, masking and compositing choices (cloud/shadow filtering, monthly medians, and retention of fully masked periods to preserve band schema stability) can yield uneven temporal support and residual observation bias, while baseline–recent differencing may conflate anthropogenic change with interannual hydroclimatic variability. Third, zonal aggregation across rings, buffers and grid cells introduces scale dependence and modifiable areal unit effects; although national clipping reduces cross-border pixel leakage, truncated outer geometries and boundary-adjacent cells can still generate edge effects that complicate direct comparability across rings. Fourth, hydrology- and riparian-restricted analyses depend on global hydrography products and user-defined thresholds (e.g., upstream-area cutoffs and corridor widths) and may not capture local drainage alterations, fine channels, engineered diversions or seasonal connectivity in small headwaters. Fifth, the scoring and prioritisation logic is inherently relative—robust normalisation, percentile ranking and phase assignment are distribution-dependent and do not represent absolute hazard or risk metrics; transfer to other domains or time windows would require recalibration of thresholds, package rules and weighting assumptions. Sixth, calibration and scientific validation are constrained by sample density, representativeness, temporal alignment and spatial autocorrelation; conventional LOOCV and random splits can overstate generalisation when samples are clustered unless spatial blocking and independent hold-out areas are available. Seventh, the financial and policy modules are assumption-driven: unit costs, package definitions, phase durations, credit eligibility (including PSA-related offsets), discounting and horizons are simplified representations that may diverge from procurement realities, regulatory requirements and implementation constraints, and several policy levers (e.g., illegal inflow reduction effectiveness, leakage, enforcement intensity, royalty basis, tail-risk severity/probability, implementation lag) are scenario-conditional and should be treated as sensitivity parameters rather than point estimates. Eighth, the uncertainty analysis propagates parameterised cost uncertainty only (via unit-cost bounds and the selected Model 5 propagation mode), and does not fully quantify uncertainty from proxy fidelity, hydrography error, masking artifacts, AOI design, behavioural responses, or model-structure alternatives; potential correlations among cost items are also not explicitly modelled. Finally, any 2060 reporting is an indexation-based translation of 2026-valued cost summaries under assumed annual escalation rates and should not be interpreted as a structural forecast of technology evolution, institutional change, land-use trajectories or future contaminant loads.

3. Results

3.1. Model 1 — Baseline Conditions and Recent-Change Signals

Model 1 characterises baseline-to-recent change signals using pixel-wise Δ surfaces defined as recent minus baseline (baseline: January 2019–December 2020; recent: February 2024–January 2025). Across terrestrial disturbance proxies, both Cutris and the 5-km buffer exhibit modest net “greening” signals. In Cutris, mean changes indicate small increases in vegetation indices (ΔNDVI = +0.008; ΔNDMI = +0.032) alongside a decrease in bare-soil exposure (ΔBSI = −0.024). The 5-km buffer shows a stronger NDVI increase (ΔNDVI = +0.067) with a smaller NDMI increase (ΔNDMI = +0.024) and a comparable decline in BSI (ΔBSI = −0.031), consistent with a more pronounced net shift in vegetation-related signals at the buffer scale.

Riparian water proxies display stronger spatial differentiation. Cutris riparian areas show a decrease in NDTI (ΔNDTI = −0.015) and an increase in RedNIR (ΔRedNIR = +0.037), whereas the riparian 5-km buffer shows decreases in both NDTI (ΔNDTI = −0.031) and RedNIR (ΔRedNIR = −0.056). These contrasting directions indicate heterogeneous riparian proxy responses between the Cutris domain and its surrounding buffer. Figure 2 summarises the spatial Δ surfaces and the area-level mean contrasts used as upstream evidence for subsequent impact scoring.

Together, these baseline-to-recent Δ proxies provide the upstream spatial evidence used to parameterise Model 2 impact scoring and to define subsequent phase-timed costing

3.2. Model 2 — Impact Scoring and Prioritization Surfaces

Model 2 constructs 0–100 prioritisation indices at 1-km grid resolution across the Crucitas landscape (n = 1,324 cells), producing component scores for land, water, and hydrology sensitivity and a combined prioritisation index (Figure 3). A total of 923 cells contained sufficient valid input data to be scored; the remaining 401 cells lacked valid pixels across the predictor stack and are therefore best interpreted as no-data (rather than low sensitivity) in spatial visualisations. For downstream portfolio costing (Model 4), missing component values are set to zero prior to phase ranking so that fiscal totals remain defined over the full 1,324-cell grid; these cells therefore populate the lowest-priority tail. Across scored cells the combined index is centred near the mid-scale (median 48.9; mean 50.0) with a pronounced upper tail (95th percentile 76.9; maximum 92.6), indicating a limited set of high-priority hotspots rather than a uniform gradient. Component indices show comparable dispersion (land median 50.8; water median 46.5; hydrology median 49.3), and land and water sensitivity co-vary at the grid-cell scale (Spearman ρ ≈ 0.58 where both indices are defined), consistent with spatial co-occurrence of terrestrial and aquatic vulnerability. The mapped surfaces in Figure 3b–e provide the spatial basis for downstream portfolio design and costing in Models 4–5.

A broadly mid-range prioritisation field is observed with a limited set of high-score hotspots, alongside moderate co-occurrence of land and water sensitivity across scored grid cells. Cells with missing predictor coverage (n = 401) are treated as no-data in mapping and should not be interpreted as low-priority areas; for portfolio costing they are retained with scores set to zero to avoid dropping area from fiscal totals.

3.3. Model 3 — Field-Anchored Screening Performance for Mercury and a Composite Contamination Index

Model 3 assesses whether remotely sensed and terrain-derived predictors recover observed variability in (i) sediment mercury expressed as log₁₀(Hg) and (ii) a composite contamination index (ContamIndex) using leave-one-out cross-validation (LOOCV) on the available field dataset (n=13). Observed Hg is strongly right-skewed (39–8,800 ppb; median 330 ppb), motivating analysis on the log₁₀ scale. Predictive skill is limited under this small-sample setting: LOOCV R² is 0.070979 for log₁₀(Hg) and 0.096187 for ContamIndex, with corresponding errors of MAE = 0.402611 and 0.296410, and RMSE = 0.533339 and 0.363161, respectively. Permutation testing (80 permutations) indicates that apparent skill does not exceed chance structure at conventional thresholds (p=0.074074 for log⁡₁₀(Hg); p=0.111111 for ContamIndex), consistent with small-sample constraints. Figure 4 summarises sampling coverage, LOOCV predicted–observed relationships, and the leading predictors identified by the permutation-importance analysis.

Across runs, importance rankings most consistently emphasised soil chemistry (pH), terrain/relief, and spatial context variables as contributors to the weak but non-zero signal.

3.4. Model 4 — Deterministic Financial Assurance Costs by Phase and Remediation Package (Baseline Portfolio)

Model 4 converts the prioritisation surface into a phased remediation portfolio across 1,324 cost grid cells, yielding a gross PV liability of US$548.0 million (10-year PV horizon; 5% discount rate). PSA-related credits total US$167.3 million PV; applying the cell-level non-negativity floor for net liabilities (Equation (23)) yields a portfolio net PV of US$408.0 million (−25.6% relative to gross; without the floor, simple gross-minus-credits would be US$380.8 million). Gross costs are dominated by remediation CAPEX (US$477.5 million; 87.1%), while monitoring contributes a material share (US$69.1 million; 12.6%); the baseline water-response module is comparatively small (US$1.36 million; 0.25%).

Liabilities are front-loaded in gross terms, with Phase 1 = US$265.0 million (48.4%), Phase 2 = US$196.6 million (35.9%), and Phase 3 = US$86.5 million (15.8%). After PSA credits (with the Equation (23) floor applied per cell), net PV is concentrated in Phase 1 (US$233.3 million) and Phase 2 (US$159.6 million), while Phase 3 remains credit-heavy but contributes a residual US$15.1 million net PV because credits exceed gross Phase 3 costs in aggregate (raw Phase 3 net = −US$12.1 million) and negative cell-level nets are floored to zero. By intervention package, Riparian and water protection is the largest net contributor (US$188.1 million; 46.1%), followed by Runoff stabilization (US$136.5 million; 33.5%) and Land revegetation and erosion control (US$83.3 million; 20.4%), consistent with PSA credits being concentrated in revegetation measures. Figure 5 synthesises the phase and package decomposition, policy-pathway overlays, and the unit-cost ranges used to parameterise the bill of quantities.

3.5. Model 4 — Policy-Pathway Effects (Deterministic Overlays)

Model 4 evaluates legislative pathways as deterministic overlays on the baseline net PV portfolio (US$408.0 million), adding expected enforcement OPEX and tail-risk provisioning and subtracting applicable revenue/offset terms. Policy-adjusted values use the same 10-year PV horizon and 5% discount rate, so differences across pathways reflect incremental policy terms rather than changes in the underlying remediation bill of quantities (Figure 5c–d).

Under the central parameterisation, the auction-based concession pathway (Exp. 24.717; mining permitted) yields a policy-adjusted net PV of US$336.1 million (Δ = −US$71.9 million vs baseline). The reduction is driven by modeled royalty revenues (PV US$93.8 million) and avoided illegal-extraction costs (PV US$7.77 million), which exceed incremental tail-risk (PV US$25.9 million) and enforcement (PV US$3.70 million) costs. In this run, the royalty is implemented as a minimum 5% take applied to an annual gross-sales stream; where operational sales forecasts are not specified, the central estimate uses a documented gross-sales proxy (Supplementary Table S7) to preserve internal consistency.

The recovery-focused pathway without industrial mining (Exp. 24.675) increases the policy-adjusted net PV to US$418.6 million (Δ = +US$10.6 million), reflecting enforcement (PV US$3.89 million) and tail risk (PV US$11.7 million) partially offset by avoided illegal-extraction costs (PV US$4.90 million). The baseline legal-constraint pathway (Law No. 8904) yields US$421.7 million (Δ = +US$13.7 million), driven by enforcement (PV US$4.08 million) and tail risk (PV US$12.2 million) with a smaller offset from avoided illegal-extraction costs (PV US$2.57 million). In the deterministic central case, the fiscal ranking is governed by the balance between enforceability and residual tail risks versus the presence of a material revenue stream (Figure 5d), while the baseline remediation portfolio is held constant across pathways.

3.6. Model 5 — 2060 Projections and Policy Ranking Summary

Model 5 propagates uncertainty by applying triangular stochastic multipliers to the exported Model 4 PV cost components (CAPEX, monitoring OPEX, and water OPEX), holding the PV credit stream fixed; simulated cell-level net costs are floored at zero and then aggregated to a portfolio total. The same draw-level treatment is applied to policy overlay components where present (enforcement and tail-risk costs versus royalty and other offset terms). For the baseline remediation portfolio, the net PV distribution is right-tailed with P10–P90 = US$385.4–519.1 million and P50 = US$450.1 million (mean US$451.0 million; simulated max US$648.1 million). Exceedance probabilities indicate that outcomes above US$400 million are common (P = 0.8357), and outcomes above US$500 million occur in a material share of draws (P = 0.1786). To report 2060 projections, portfolio outcomes are scaled from the PV base year to 2060 using (1+γ)³⁴, giving baseline P50₂₀₆₀ = US$0.883 billion (γ=2%) and US$1.708 billion (γ=4%) (from US$450.1 million at γ=0%).

Policy-pathway uncertainty is evaluated by adding a stochastic policy delta to each baseline draw. Under Exp. 24.717, the policy-adjusted net PV is P10–P90 = US$312.4–457.2 million with P50 = US$382.8 million (mean US$383.9 million), a median reduction of ~US$67.3 million versus the baseline median; P(>US$400 million) = 0.3833 and P(>US$500 million) = 0.0179, and the policy reduces net PV in 100% of draws. Under Exp. 24.675, the policy-adjusted net PV increases to P50 = US$460.8 million (P10–P90 US$396.2–529.8 million), with P(>US$400 million) = 0.8846 and P(>US$500 million) = 0.2325. Under Law No. 8904, the policy-adjusted net PV is P50 = US$463.9 million (P10–P90 US$399.0–533.4 million), with P(>US$400 million) = 0.8956 and P(>US$500 million) = 0.2497. Under the same 2060 scaling, P50₂₀₆₀ is US$0.748B / 1.447B (Exp. 24.717), US$0.902B / 1.746B (Exp. 24.675), and US$0.908B / 1.757B (Law 8904) for γ=2%/4%, respectively.

Figure 6 summarises portfolio-level uncertainty through the cumulative distribution of net PV cost, the corresponding exceedance curve, a pseudo–one-way sensitivity ranking of unit-cost parameters, and a projection of net cost to 2060 under alternative indexation rates

Across scenarios, Exp. 24.717 consistently lowers exceedance risk relative to baseline in all draws, whereas the non-mining pathways shift the distribution upward, increasing exceedance probabilities at common liability thresholds.

4. Discussion

This study’s primary contribution is an integrated decision pipeline that converts observable environmental signals into (i) spatially explicit prioritisation surfaces, (ii) a phase-timed remediation portfolio, and (iii) present-value (PV) fiscal comparisons of alternative policy pathways under uncertainty. The key implication is that Crucitas is not well represented by either a purely descriptive disturbance narrative or an aggregate “total damages” figure. Liabilities concentrate spatially and temporally, and the fiscal performance of legislative pathways is governed by whether added enforcement and residual-risk obligations are credibly offset by enforceable revenues or other verified offsets.

Interpreting Models 1–3 together supports a clear distinction between robust spatial triage and predictive certainty about contaminant magnitudes. Remote-sensing proxies and multi-criteria indices provide a defensible basis for prioritising where interventions and verification should occur first. By contrast, the field-anchored screening (Model 3) provides structured but limited evidence under small-sample constraints and heterogeneous exposure pathways. This evidentiary architecture is methodologically appropriate for contested landscapes: it prevents over-claiming while still producing operationally actionable outputs. The decision interpretation is therefore “prioritise and validate,” not “predict and declare.” Consistent with this role, Model 3 yields only weak LOOCV skill (R²≈0.07–0.10) and does not meet conventional permutation-based significance thresholds, so outputs are interpreted as triangulation evidence rather than predictive mapping.

Model 4’s portfolio logic is as important as the PV magnitude itself. Converting a prioritisation surface into intervention packages and phases makes the liability eligible for budgeting, sequencing, and assurance design. Portfolio decomposition clarifies which components dominate early capital exposure versus long-horizon operating and monitoring burdens and provides a practical accounting rule for policy use: credits can reduce gross totals, but net liabilities should be bounded at zero at the grid-cell level (and thus cannot be negative in aggregate) to avoid planning artefacts. This portfolio framing also enables transparent comparison of legislative pathways. The policy question is not only what a bill permits, but what it implies for enforceable funding streams, implementation lags, monitoring obligations, and residual provisioning.

The results position Crucitas less as a one-off cleanup invoice and more as a closure-governance and financial-assurance problem under persistent illegality: costs become decision-relevant only when paired with credible enforceability, transparent cost baselines, and explicit residual-risk provisioning. International guidance converges on the same logic—closure planning should be integrated and iterative, linked to robust cost estimation and long-term governance, with financial assurance designed to prevent liabilities being externalized to the state and to enable third-party execution if obligations are not met (World Bank, 2021; Hattingh et al., 2021; ICMM, 2019). In contested or illegal settings where an operator cannot be reliably bonded, the evaluated policy pathways operationalize a public-sector analogue: ring-fenced funding and enforceable offsets, paired with monitored effectiveness and uncertainty bounds, so that “policy choice” is assessed in fiscally comparable terms rather than rhetorically. This aligns with integrated closure good practice emphasizing early planning, progressive implementation, and governance for relinquishment (ICMM, 2019; World Bank, 2021). For illegal landscapes specifically, this study adds a required extension: enforcement performance is treated as a first-order uncertainty term that must be priced, stress-tested, and carried into residual risk.

Within this accounting logic, pathway ranking reflects a balance between added obligations (enforcement OPEX, governance overhead, residual tail-risk provisioning) and offsets (royalty-like revenues, avoided illegal-extraction damages where plausibly achieved, and other verified financing). Formalisation is often justified on the basis of potential revenues and enforceable repair obligations, contingent on strong oversight and credible institutional controls (Redacción La República, 2023). Fiscal outcomes also depend on concession design and operator quality; weak screening can reduce compliance performance even when headline royalty terms appear attractive (Martínez, 2025). Timing is a first-order determinant of net PV: long-lead legal, technical, and financing steps can prolong interim enforcement and remediation burdens before offsets materialise (Arrieta, 2024).

Model 5 adds decision value by shifting interpretation from deterministic rankings to conditional risk management: under what parameter conditions does a pathway reduce liability, and what is the probability of exceeding policy-relevant thresholds? This perspective is particularly relevant in Crucitas because right-tail drivers—enforcement effectiveness, residual risk severity, and the reliability of any offset stream—are precisely the parameters most exposed to governance slippage. Reporting that illegal extraction activity may be expanding to additional sites underscores why enforcement should be treated as spatially distributed and adaptively dynamic rather than as a one-time control action (Hidalgo, 2025). Operational interventions can suppress capacity in the near term, but they also motivate treating effectiveness and persistence as uncertain inputs rather than fixed constants (Colindres Lagos, 2025).

A second boundary condition concerns leakage and provenance. If material can be moved, processed elsewhere, and exported under weak traceability, “site-only” remediation portfolios can understate total liabilities or misattribute where costs accrue. Investigative reporting linking illegal supply chains to laundering routes and export flows provides a rationale for treating traceability failure and off-site processing as explicit sensitivity cases, because these mechanisms can shift tail outcomes even when median costs appear stable (Bolaños & Ubaque Calixto, 2025). Relatedly, civil-society position documents emphasizing district-scale footprint expansion and long-lived liabilities are best interpreted as stakeholder-aligned evidence that the debate concerns long-duration risk and governance credibility, rather than as empirical magnitude estimates (Observatorio de Bienes Comunes UCR, 2025). This alignment supports the exceedance-risk framing while preserving a clear distinction between scenario assumptions and measured quantities.

A decision-grade interpretation must also recognize that fiscal-liability accounting operates within a socio-environmental conflict system where legitimacy, inclusion, and institutional capacity condition what is implementable. Comparative scholarship on Costa Rican gold-mining conflicts shows that outcomes are shaped by collective action, legal contestation, and the state’s enforcement and regulatory capacity—constraints that can dominate feasibility even when a pathway appears fiscally favourable on paper (Blanco Obando, 2021). Contested environmental decisions also distribute participation and exposure asymmetrically, highlighting why fiscal rankings should be presented as conditional and complemented by transparent benefit and risk allocation criteria (Rodríguez-Soto, 2025). Mercury risk reduction is similarly mediated through public expectations and multi-level governance roles, reinforcing that enforcement performance, monitoring capacity, and risk communication are policy-dependent levers rather than background assumptions (Campos-Guevara et al., 2022).

The framework relies on remote-sensing proxies of disturbance and hydrologic context rather than direct contaminant concentration surfaces, and Model 3 is constrained by small-sample field data that support triangulation but not predictive mapping. Uncertainty propagation is implemented primarily through unit-cost envelopes; enforcement effectiveness, behavioural adaptation, leakage, and proxy fidelity are treated as scenario assumptions rather than fully stochastic drivers. Spatial boundaries also matter: if processing and disposal are displaced beyond the mapped footprint, system-wide liabilities can exceed “site-only” portfolio totals. Long-horizon indexation and scenario projections should therefore be interpreted as stress tests for fiscal exposure rather than forecasts.

A bounded, decision-oriented conclusion follows. Across all legislative pathways, Crucitas faces substantial remediation obligations requiring credible financing and sequenced implementation. Any fiscal advantage from a mining-allowed pathway is contingent on enforceable and transparent offsets, competent implementation, and credible residual-risk provisioning. Conversely, any non-extractive pathway relying on remediation-linked offsets or institutional finance requires comparable transparency about feasibility, timing, and governance controls. The most immediate value of the framework is not to “settle” the debate, but to clarify which assumptions drive rankings, how uncertainty concentrates in the right tail, and which empirical evidence—expanded field validation, traceability and provenance assessment, and cost engineering—would most reduce decision uncertainty.

5. Conclusions

Crucitas faces substantial remediation obligations across legislative pathways, requiring credible financing, transparent sequencing, and explicit treatment of residual risk. The integrated framework translates remote-sensing change signals into a reproducible prioritisation surface, converts that surface into a phase-timed remediation portfolio with PV-consistent accounting, and evaluates policy pathways on a common fiscal and exceedance-risk basis. Field-anchored screening (Model 3) provides triangulation evidence rather than predictive certainty about contaminant magnitudes under the current small-sample constraint, reinforcing a “prioritise and validate” interpretation. Portfolio results indicate a large gross PV liability with partial offset from PSA-related credits, while pathway rankings are driven by enforceability and timing of offsets and by the magnitude of enforcement and tail-risk provisions rather than by pathway labels alone. Under the central parameterisation, Exp. 24.717 reduces net PV relative to baseline through modeled revenues; non-mining pathways increase net PV through enforcement and residual-risk terms. Under uncertainty, Exp. 24.717 reduces the probability of exceeding policy-relevant thresholds across draws, whereas non-mining pathways shift the liability distribution upward.

The results should be interpreted as conditional accounting, not a definitive forecast. Decision relevance is maximised when load-bearing assumptions are made testable: enforcement effectiveness and persistence, traceability and leakage, credibility of revenue/offset collection, and tail-risk severity. Priority evidence to reduce decision uncertainty includes expanded field sampling designed for spatial representativeness, supply-chain and provenance assessment, and site-specific engineering cost studies capable of narrowing unit-cost envelopes. With these additions, the same pipeline can support periodic re-ranking of intervention priorities and policy-pathway performance as conditions and governance constraints evolve. If sustained enforcement cannot be guaranteed, non-extractive pathways require dedicated, verifiable financing; if an extractive pathway is pursued, it should be conditioned on ring-fenced funds, enforceable guarantees, and independent audit from the outset