Quantile Delta Mapping Downscaling of CMIP6 Daily Precipitation Under SSP1-2.6 to SSP5-8.5: Model Structure Dominates Projection Uncertainty at a Tropical High-Altitude Station

1. Introduction

Global freshwater availability is increasingly affected by anthropogenic climate change, with precipitation regimes undergoing shifts in intensity, frequency, and seasonality that challenge both engineering design and ecosystem functioning [1,2]. The Andean region of South America hosts some of the world’s most climate-sensitive headwater areas, supplying drinking water and irrigation to millions while simultaneously experiencing accelerating glacier retreat, altered El Niño–Southern Oscillation (ENSO) dynamics, and intensifying extreme rainfall events [3,4]. Across Latin America, hydrological extremes represent among the costliest natural hazards, with documented annual economic losses exceeding USD 20 billion during 2000–2019 according to EM-DAT global disaster records [5], disproportionately affecting mountainous areas where runoff responds non-linearly to precipitation inputs [6]. Understanding how precipitation will evolve under different greenhouse-gas forcing scenarios is therefore a societal imperative with direct implications for infrastructure sizing, flood-risk mapping, and agricultural planning [7].

Southern highland Ecuador occupies a particularly sensitive hydroclimatic niche, straddling the Intertropical Convergence Zone and subject to overlapping atmospheric teleconnections from Pacific and Atlantic sea-surface temperature anomalies [8,9]. The regional precipitation regime emerges from the convergence of South American monsoon outflow, easterly trade-wind moisture advection, and vigorous orographic uplift along an elevation gradient spanning 800–3800 m a.s.l. [10]. The La Argelia meteorological station (3.9999° S, 79.2000° W, 2160 m a.s.l.) provides one of the few continuous daily precipitation records available for the southern Ecuadorian Andes [11]. A single station cannot represent the full spatial variability of precipitation across this elevation gradient, and the projections presented here should accordingly be read as being representative of this specific site rather than of the surrounding region as a whole. Despite the hydrological significance of this location, no published study has simultaneously applied a multi-model CMIP6 ensemble with chronology-preserving statistical downscaling across multiple SSP scenarios at daily resolution using a homogenized observational baseline for this high-altitude station — a methodological gap that constrains science-based water management in southern highland Ecuador.

Statistical downscaling constitutes the principal methodological bridge between coarse-resolution GCM outputs (typically 100–250 km grid cells) and the point-scale observations required for local impact assessment [12,13]. Among available methods, quantile-mapping technique have become the standard approach for daily precipitation because they preserve the full distribution of observed precipitation intensity while correcting systematic GCM biases in both the mean and the tails [14,15]. A distinction of practical importance among quantile-mapping variants is whether the future daily sequence is drawn from the GCM's own projected trajectory or reconstructed from an analogue historical template: the former preserves genuine day-to-day persistence, interannual variability, and any trend embedded in the model's simulated future, while the latter — common in earlier delta-change implementations — can inadvertently substitute the historical record's own temporal structure for the model's projected one. Quantile delta mapping (QDM [16]) addresses this by computing each future day's bias-correction factor from its own rank within the model's future distribution, so that the corrected series inherits its chronology directly from the model's simulated trajectory rather than from a resampled historical analogue.

The Coupled Model Intercomparison Project Phase 6 (CMIP6) ensemble provides unprecedented scenario breadth through the Shared Socioeconomic Pathway (SSP) framework, enabling probabilistic assessments of future precipitation under optimistic (SSP1-2.6), intermediate (SSP2-4.5), and pessimistic (SSP5-8.5) emission trajectories [17]. Multi-model ensemble approaches applied to CMIP6 projections have demonstrated utility in quantifying future precipitation extremes and hydrological flood risk across a wide range of hydroclimatic settings [18,19].

Rigorous GCM performance evaluation prior to downscaling is a methodological prerequisite in applied climate science, as structural biases in raw model outputs directly propagate into projected extreme indices and return-period estimates [20]. The Kling–Gupta Efficiency (KGE [21]) and Nash–Sutcliffe Efficiency (NSE [22]) are widely adopted composite metrics that simultaneously account for correlation, bias, and variability ratio. Taylor diagrams [23] provide complementary visual diagnostics of model skill, enabling identification of structurally diverse model subsets. Independent out-of-sample validation — calibrating a bias-correction transfer function on part of the observational record and evaluating it against a held-out period that never entered calibration — is standard practice for establishing that a statistical downscaling procedure captures genuine predictive skill rather than simply reproducing the data used to fit it [13]. The analysis of ETCCDI extreme precipitation indices—particularly Rx1day, SDII, R95p, CWD, and CDD—is standard in regional climate impact studies because these indices map directly onto engineering design thresholds for drainage infrastructure, agricultural drought, and flood risk [24,25]. Return-period analysis using Gumbel extreme-value distributions, combined with bootstrap confidence intervals, provides decision-relevant probabilistic estimates of peak precipitation for infrastructure design [26]. Partitioning multi-model projection uncertainty into scenario- and model-driven components, following the variance-decomposition logic of Hawkins and Sutton [27], is routine in European and North American regional climate studies but has seen limited application in the tropical Andes, where the relative roles of GCM structural diversity and emission-pathway choice remain poorly quantified for daily precipitation.

To close this gap, this study applies a reproducible multi-model multi-scenario quantile delta mapping framework to five CMIP6 GCMs—ACCESS-CM2, CMCC-ESM2, CanESM5, MIROC6, and MPI-ESM1-2-LR—under SSP1-2.6, SSP2-4.5, and SSP5-8.5 for 2051–2080, calibrated against the SC-PREC4SA homogenized observational record over 1981–2010. The primary contributions are: (i) a chronology-preserving multi-model daily precipitation projection across three SSP scenarios for a high-altitude station in this climatically critical region; (ii) independent split-sample validation demonstrating genuine out-of-sample bias-correction skill; (iii) a comprehensive ETCCDI extreme-index assessment with inter-model confidence envelopes; (iv) Gumbel-fitted return-period estimates with bootstrap 95% confidence intervals suitable for infrastructure design; (v) a two-way variance decomposition quantifying the relative contributions of GCM structure, scenario choice, and their interaction; and (vi) a fully documented computational pipeline, described at algorithmic level in Appendix A, transferable to other high-altitude tropical stations across the Andes.

2. Materials and Methods

2.1. Study Area and Observational Data

La Argelia meteorological station (INAMHI code M0033; 3.9999° S, 79.2000° W; 2160 m a.s.l.) is located on the eastern slope of the Andes in Loja Province, southern Ecuador, within the Northern Andes and Sierra (NAS) ecoregion [28]. The station’s precipitation regime is bimodal, with primary maxima in March–April and a secondary wet period in October–November, driven by the seasonal migration of the ITCZ and easterly trade-wind incursions [10,11]. Daily precipitation observations were sourced from the SC-PREC4SA dataset (Huerta et al. [29]), which provides homogenized and bias-corrected daily precipitation series for 1981–2010 across the NAS ecoregion. The SC-PREC4SA product applies quality control, gap-filling, and homogenization following World Meteorological Organization protocols.

For station M0033, the 1981–2010 reference series comprises 10,957 daily values with 0.00% missing data, a mean of 2.554 mm day⁻¹ (standard deviation 5.418 mm day⁻¹), a maximum of 59.7 mm day⁻¹, 36.8% wet-day frequency (threshold: >1 mm day⁻¹), and a mean annual precipitation of 933 mm year⁻¹, consistent with published climatologies for southern Ecuador [8,11].

2.2. CMIP6 Model Ensemble

Five CMIP6 GCMs were selected based on established performance over the tropical Andes, structural independence, and availability of daily precipitation (variable pr) in the Pangeo CMIP6 cloud catalog: MPI-ESM1-2-LR (Max Planck Institute, Germany); T63 spectral resolution, ~200 km), CMCC-ESM2 (Euro-Mediterranean Centre, Italy; ~100 km), MIROC6 (Japan; ~250 km), CanESM5 (Canada; ~500 km), and ACCESS-CM2 (CSIRO, Australia; ~250 km). This set spans distinct modeling centers, dynamical cores, and convective parameterization schemes, providing structural diversity for the ensemble beyond a single institutional lineage. All models use realization r1i1p1f1. Historical data cover 1981–2010; future data cover 2051–2080 under ssp126, ssp245, and ssp585. A sixth candidate model, NorESM2-LR, was evaluated but excluded because its daily precipitation output was unavailable in the Pangeo catalog for this station and period at the time of analysis. ACCESS-CM2 future data were additionally unavailable for SSP1-2.6 and SSP5-8.5; the effective ensemble size is therefore n = 4 for those two scenarios and n = 5 for SSP2-4.5.

Global GCM output was downscaled directly from the Pangeo CMIP6 catalog rather than from a dynamically downscaled regional product such as CORDEX-CORE, for two reasons specific to this application. First, at the time of analysis the CORDEX-CORE South America domain did not provide a continuous daily precipitation archive at 25 km covering both the historical and all three SSP branches evaluated here for this specific grid cell. Second, the analysis is designed to quantify inter-model structural uncertainty across independent modeling centers (Section 2.8), which requires the broadest available multi-model CMIP6 ensemble rather than the smaller set of driving GCMs typically downscaled by any single regional climate model initiative. Extending this comparison to include CORDEX-CORE simulations as they become available for this domain is identified as a priority for future work (Section 4.5).

2.3. GCM Performance Evaluation

Nine performance metrics were computed by comparing monthly aggregated historical GCM precipitation against the SC-PREC4SA observational record at station M0033 over the 1981–2010 reference period. Negative NSE and sub-optimal KGE values in raw CMIP6 outputs evaluated against point-scale observations over complex terrain are a well-documented consequence of systematic coarse-resolution wet biases, not a reflection of methodological deficiencies [20,30]; composite ranking was therefore based on relative temporal correlation and bias magnitude rather than absolute efficiency thresholds. The metrics are: Pearson r; relative bias; RMSE; MAE; NSE [22]; KGE [21]; PBIAS; wet-day frequency error (WF_err); and the 95th-percentile ratio (P95_ratio [31]). GCMs were ranked by the sum of KGE, NSE, and r ranks (ascending). See Equations (1)-(3).

$$\mathrm{r}={\displaystyle \frac{\sum \left({\mathrm{O}}_m-\overline{O}\right)\left({\mathrm{S}}_m-\overline{S}\right)}{\sqrt{\sum {\left({\mathrm{O}}_m-\overline{O}\right)}^2}\hspace{1em}\sqrt{\sum {\left({\mathrm{S}}_m-\overline{S}\right)}^2}}}$$

(1)

$$\mathrm{NSE}=1 – {\displaystyle \frac{\sum {\left({\mathrm{O}}_{\mathrm{m}} -{ \mathrm{S}}_{\mathrm{m}}\right)}^2}{\sum {\left({\mathrm{O}}_{\mathrm{m}} - \stackrel{\mathrm{-}}{\mathrm{O}}\right)}^2}}$$

(2)

$$\mathrm{KGE}=1 – \sqrt{{\left(r-1\right)}^2+{\left(\alpha -1\right)}^2+{\left(\beta -1\right)}^2}$$

(3)

where α = σₛ/σₒ , β = μₛ/μₒ; Oₘ and Sₘ denote observed and simulated monthly precipitation.

2.4. Quantile Perturbation Downscaling

Before applying the bias-correction procedure to project future precipitation, its out-of-sample skill was tested independently of the future projections themselves. For each model, the transfer function described in Section 2.5 was calibrated using only the 1981–2000 sub-period of the historical record and observations and then applied to the model's own 2001–2010 historical data — a period that never entered calibration — to produce a corrected estimate that could be compared against the observed record for the same withheld years. NSE, KGE, and percent bias were computed for the raw and bias-corrected series in the validation period alike.

2.5. Quantile Delta Mapping Downscaling

Daily precipitation was downscaled using quantile delta mapping (QDM [16]), which corrects each future daily value using a bias-correction factor derived from its own quantile position within the future model distribution, so that the resulting series preserves the day-to-day sequencing, autocorrelation, and any embedded trend of the model's own projected trajectory. For each calendar month m and each model–scenario combination, empirical quantile functions were estimated at 1001 points q ∈ [0, 1] from the historical model series (Qhist), the observed series (Qobs), and the future model series (Qfut). For each future day with value x, the quantile rank τ is its own empirical percentile within Qfut for that calendar month; the historical-model and observed magnitudes at that same rank, Qhist(τ) and Qobs(τ), give the relative perturbation factor RF(τ) = x / [Qhist(τ) + ε] (ε = 10⁻⁶, clipped to the range [0.1, 3.0] to bound the correction near the dry/wet transition, where near-zero historical quantiles would otherwise produce unbounded scaling factors) and the absolute perturbation AF(τ) = x − Qhist(τ). The downscaled value is given by Equation (4).

$$P_{ds}=\max \left[Q_{obs}\left(\tau \right)·F_R\left(\tau \right)+F_A\left(\tau \right),0\right]$$

(4)

where FR(τ) and FA(τ) are blended relative and absolute perturbation factors. Blending follows a linear transition between a threshold T_low = 1.0 mm day⁻¹ and T_high = 5.0 mm day⁻¹: below T_low, the absolute factor dominates (preserving near-zero values); above T_high, the relative factor dominates (preserving intensity scaling); in the transition band, a linear weight w = (Q_obs(τ) − T_low)/(T_high − T_low) is applied. The result is indexed by the future day's own real calendar date, so the temporal order of the downscaled series reflects the model's simulated future trajectory rather than a resampled historical sequence. A fixed random seed (numpy seed = 42) was used throughout. All downscaled series cover 2051-01-01 to 2080-12-31 (10,958 daily values per model–scenario combination).

2.6. Extreme Climate Indices (ETCCDI)

Six ETCCDI indices were computed: Rx1day, Rx5day, SDII, R95p, CWD, and CDD, following standard definitions [24,25] (Table 1). All indices were computed per year and averaged over the 30-year analysis window before cross-scenario comparison. In addition to the annual-maximum consecutive wet/dry day counts (CWD, CDD) defined by the ETCCDI framework, the mean and 90th-percentile duration of individual wet and dry spells (threshold 1 mm day⁻¹) were computed across the full record to characterize the typical, rather than the annual-extreme, persistence of precipitation and dry periods.

2.7. Return-Period Analysis

Annual daily maxima were fitted to the Gumbel (Type I GEV) distribution by maximum likelihood, with annual maxima capped at their 99th percentile to limit the leverage of a single extreme outlier year on the fitted parameters; a sensitivity check comparing capped and uncapped fits across all model–scenario combinations showed a mean difference of +0.6% in the 100-year return level, indicating the cap has negligible practical effect at this station. Return quantiles Q(T) were estimated for T ∈ {2, 5, 10, 20, 50, 100} years. See Equation 11.

$$Q\left(T\right)=\mathsf{\mu } - \mathsf{\beta } · \ln \left[-\ln \left(1-{\displaystyle \frac{1}{T}}\right)\right]$$

(11)

where μ and β are the Gumbel location and scale parameters. The Gumbel (Type I) distribution was adopted rather than a three-parameter alternative such as the (log-) Pearson Type III or the generalized extreme value (GEV) distribution because the historical annual-maxima series available for direct fitting spans only 30 years; with this sample size, the additional shape parameter of a three-parameter distribution is estimated with substantially greater uncertainty than the two Gumbel parameters, and the two-parameter Gumbel model remains standard practice for daily precipitation extremes of this record length [26]. This choice is revisited in Section 4.4 as a limitation, since tropical precipitation extremes do not always conform to a Type I (light-tailed) shape. Bootstrap 95% confidence intervals (500 resamples with replacement, fixed seed) were derived for each model. Differences in return levels between the SSP2-4.5 and SSP5-8.5 ensembles were assessed with the two-sided Mann–Whitney U test at α = 0.05, computed two ways: first directly on the four-to-five point estimates (one per model), and second on the pooled bootstrap replicates of every model within a scenario (~2000 samples per scenario), which is markedly more informative given the small number of independent model realizations available for the first approach.

2.8. Uncertainty Decomposition

Total projection variance was partitioned into scenario, model, and model × scenario interaction components using a two-way analysis-of-variance sum-of-squares decomposition. For each calendar month m, the total sum of squares of the full scenario × model table of mean monthly values was computed around the grand mean and split into the between-scenario sum of squares, the between-model sum of squares, and the remainder attributed to their interaction. With one realization per scenario–model cell, the interaction term is confounded with any single-realization noise, so no independent residual/error term is separately estimable; this decomposition therefore reports scenario, model, and interaction contributions rather than a fourth, purely random component, and the interaction term should not be interpreted as an estimate of internal climate variability in the sense of Hawkins and Sutton [27], which requires long continuous simulations to isolate.

2.9. Signal-to-Noise Ratio and Model Consensus

The signal-to-noise (S/N) ratio was computed for each calendar month and scenario as |ensemble mean change| divided by the inter-model standard deviation. S/N > 1 indicates that the forced signal exceeds inter-model spread. Model consensus on the sign of change was defined as the fraction of models agreeing on the direction of change (positive or negative), with a robustness threshold of ≥80% adopted in place of a fixed fraction, given that the actual ensemble size varies between four and five models across scenarios (Section 2.2) rather than the six-model design against which a fixed 5-of-6 (83%) criterion would be calibrated.

2.10. Statistical Software and Reproducibility

All analyses were implemented in Python (pandas, numpy, scipy, matplotlib, intake-esm). CMIP6 data were accessed via the Pangeo CMIP6 catalog through Zarr streaming. The random seed was fixed at 42 for all stochastic operations. An algorithmic summary of the computational pipeline is provided in Appendix A.

3. Results

3.1. Observational Baseline and GCM Performance Evaluation

The SC-PREC4SA series (M0033, 1981–2010; N = 10,957; 0.0% missing) exhibits a mean daily intensity of 2.5540 mm day⁻¹ (σ = 5.4180 mm day⁻¹), a wet-day frequency of 36.8%, and a bimodal annual cycle with primary maximum in March (~150 mm) and a secondary wet period spanning October–December (~70–90 mm). See Figure 1.

GCM performance metrics are summarized in Table 2. MPI-ESM1-2-LR ranked first by composite KGE–NSE–Pearson r score. All five models yield negative NSE and sub-optimal KGE values when evaluated against point-scale observations without prior bias correction, reflecting systematic wet biases intrinsic to coarse-resolution GCMs over complex orographic terrain rather than methodological deficiencies. MPI-ESM1-2-LR achieved the strongest temporal correlation (r = 0.830) and the smallest overestimation bias (PBIAS = 23.6%). CanESM5 and ACCESS-CM2 exhibited severe the largest wet biases, exceeding 250%, consistent with coarser effective resolution and known overestimation of convective precipitation in complex terrain by these models [31].

Figure 2. Multi-metric evaluation of raw CMIP6 GCMs against SC-PREC4SA observations (1981–2010). (a) Normalized performance heatmap across seven evaluation metrics. (b) Mean annual precipitation cycle comparing the raw GCM outputs against the homogenized observational baseline.

Figure 2. Multi-metric evaluation of raw CMIP6 GCMs against SC-PREC4SA observations (1981–2010). (a) Normalized performance heatmap across seven evaluation metrics. (b) Mean annual precipitation cycle comparing the raw GCM outputs against the homogenized observational baseline.

3.2. Independent Validation

Split-sample validation, in which the bias-correction transfer function was calibrated on 1981–2000 and evaluated against the withheld 2001–2010 period, confirmed genuine out-of-sample skill for every model (Table 3). Raw historical GCM output showed strongly negative NSE and KGE in the validation window for all five models (NSE from −1.2 to −71.7; KGE from −0.08 to −5.04), reflecting the same systematic wet bias documented in Section 3.1. After bias correction, KGE improved to 0.60–0.76 and percent bias fell to within ±10.5% of observed for every model, none of which had contributed to the calibration period. This out-of-sample improvement — obtained on data never used to fit the correction — indicates that the transfer function captures a genuine, transferable statistical relationship between model and station precipitation rather than simply reproducing its own calibration data.

3.3. Downscaled Precipitation: Annual and Monthly Changes

Quantile delta mapping produced modest and directionally mixed changes in ensemble-mean annual precipitation relative to the 932.6 mm year⁻¹ historical baseline (Table 4; Figure 3 and Figure 4). Under SSP1-2.6, the ensemble mean declined by −12.1 mm year⁻¹ (−1.3%; inter-model range 825–1016 mm year⁻¹). Under SSP2-4.5, the ensemble mean increased by +47.2 mm year⁻¹ (+5.1%; range 830–1080 mm year⁻¹). Under SSP5-8.5, the change was negligible at −3.3 mm year⁻¹ (−0.4%; range 783–1077 mm year⁻¹). At the level of individual models, CMCC-ESM2 and CanESM5 consistently projected decreases across all scenarios (−4.3% to −17.0%), whereas MIROC6 and MPI-ESM1-2-LR consistently projected increases (+4.0% to +16.6%); ACCESS-CM2, available only for SSP2-4.5, projected an increase of +16.3%. This split between consistently drying and consistently wetting models, more than the choice of emission scenario, governs the sign of the ensemble mean in any given pathway.

Figure 5. Annual precipitation time series for raw GCMs versus observations and future downscaled projections. (a) Historical annual precipitation (1981–2010), raw CMIP6 outputs versus the homogenized SC-PREC4SA baseline. (b) Future downscaled ensemble mean precipitation with inter-model range under three SSP scenarios (2051–2080).

Figure 5. Annual precipitation time series for raw GCMs versus observations and future downscaled projections. (a) Historical annual precipitation (1981–2010), raw CMIP6 outputs versus the homogenized SC-PREC4SA baseline. (b) Future downscaled ensemble mean precipitation with inter-model range under three SSP scenarios (2051–2080).

3.4. Changes in Precipitation Distribution and Wet-Day Frequency

Quantile distributions of daily precipitation for each scenario, computed from each model's own distribution before aggregation, track the observed distribution closely across nearly the full probability range under all three scenarios, with a modest upward shift confined to the uppermost quantiles (Figure 6a). Wet-day frequency showed modest, non-significant adjustments (individual models ranging from about 28% to 46% against an observed 36.8%), with no directional consensus across models under any scenario (Figure 6b). The monthly change heatmap (Figure 7) shows that the largest relative anomalies occur across individual models in austral autumn and winter (April–September), with several models projecting changes exceeding ±30% in a given month; inter-model agreement on the sign of change is generally lower in this season than during the primary wet season.

3.5. ETCCDI Extreme Precipitation Indices

ETCCDI indices (Table 5; Figure 8) indicate precipitation tail intensification across scenarios, most consistently for the fixed-day maxima and mean intensity indices. Rx1day ensemble means increase relative to the observed 40.5 mm, reaching 46.4 mm under SSP1-2.6, 46.6 mm under SSP2-4.5, and 52.2 mm under SSP5-8.5 (+14.6% to +28.9%). Rx5day increases from an observed 80.3 mm to 88.5–92.7 mm (+10.2% to +15.4%) across scenarios. SDII rises modestly from an observed 6.7 mm day⁻¹ to 6.7–7.3 mm day⁻¹. R95p increases from an observed 399 mm to 407–421 mm. Consecutive-day indices are more mixed: the annual-maximum CWD (8.7 days observed) increases to 11.2–13.7 days across scenarios without a clear ordering by scenario, and annual-maximum CDD (19.5 days observed) increases to 26.8–29.3 days across all three scenarios, indicating a lengthening of the longest annual dry spell alongside intensified peak-day rainfall.

3.6. Return-Period Analysis

Gumbel maximum-likelihood fits to the observed series yield a RP-10 of 54.2 mm day⁻¹ [95% CI: 47.7–59.7] and RP-50 of 67.6 mm day⁻¹ [95% CI: 58.1–75.4] (Table 5; Figure 9). Future ensemble-mean projections show systematic upward shifts: the RP-100 increases from the observed 73.2 mm day⁻¹ to 91.9 ± 10.3 mm day⁻¹ under SSP1-2.6, 88.3 ± 6.7 mm day⁻¹ under SSP2-4.5, and 107.5 ± 15.7 mm day⁻¹ under SSP5-8.5—increments of 26%, 21%, and 47%, respectively. A sensitivity check comparing Gumbel fits with and without the 99th-percentile cap on annual maxima showed a mean difference of only +0.6% in RP-100 across all 13 available model–scenario combinations, indicating the cap does not materially affect the reported return levels.

Comparing SSP2-4.5 against SSP5-8.5 with the Mann–Whitney U test applied directly to the four-to-five model point estimates yields no statistically significant difference at any return period (p = 0.111–0.190), a result of limited value given the very small number of independent realizations available to the test. Pooling the 500 bootstrap replicates generated for each model's confidence interval into a single per-scenario sample (~2000 values) and repeating the comparison yields a statistically significant difference at all three return periods (RP-10, RP-50, RP-100; p < 0.05), with higher medians under SSP5-8.5 than SSP2-4.5 at every return period considered (e.g., RP-50 medians of 79.5 mm versus 95.3 mm). The two tests answer related but distinct questions — whether the four-to-five model realizations differ, versus whether the pooled distributions of plausible return-level estimates differ — and their divergent outcomes illustrate how sensitive small multi-model ensemble comparisons can be to the way the comparison sample is constructed

3.7. Temporal Trend Analysis (Mann–Kendall)

Mann–Kendall trend analysis of the 30-year annual precipitation time series (2051–2080) revealed predominantly weak and non-significant trends across models and scenarios (Figure 10). Individual-model Kendall's τ ranged from −0.18 to +0.19, and the ensemble-mean τ from −0.01 to +0.06 across scenarios, none statistically significant at p < 0.05. Because the future series preserves each model's own simulated chronology (Section 2.5), this analysis evaluates genuine projected temporal structure rather than a resampled historical sequence; the absence of significant trends indicates that projected changes over 2051–2080 primarily manifest as shifts in the precipitation distribution rather than as a systematic monotonic trend within this 30-year window.

3.8. Uncertainty Decomposition and Model Consensus

The two-way variance decomposition (Figure 11a) shows that inter-model structural uncertainty is the largest single component of total projection variance in every calendar month, ranging from 51% to 93% (mean 75%). Scenario uncertainty is generally small, ranging from under 1% to 19% (mean 8%), and is largest in December–January. The model × scenario interaction term ranges from 0% to 33% (mean 17%) and exceeds 30% in January, March, and April; because only one realization exists per model–scenario combination, this term is confounded with single-realization noise and should not be read as an independent estimate of internal climate variability (Section 2.8).

Signal-to-noise ratios exceed 1.0 in only a small minority of month–scenario combinations, most notably January under SSP1-2.6 (Figure 12a), indicating that inter-model spread exceeds the ensemble-mean signal through most of the annual cycle. Applying an 80% model-agreement threshold — set to reflect the actual four-to-five model ensemble size rather than a fixed fraction calibrated to a six-model design — six of the 36 month–scenario combinations qualify as robust: January under SSP1-2.6, and February, March, June, September, and December under SSP2-4.5 (Figure 13). No month–scenario combination under SSP5-8.5 reaches this threshold. The concentration of robust signals under SSP2-4.5, and their absence under SSP5-8.5, is consistent with the higher inter-model dispersion evident in the return-period and annual-total results for the pessimistic scenario (Section 3.3 and Section 3.6).

3.9. Additional Hydroclimatic Diagnostics

Figure 14, Figure 15, Figure 16, Figure 17, Figure 18, Figure 19 and Figure 20 present further diagnostic of model performance, seasonal dynamics, and extreme precipitation behavior. The Taylor diagram (Figure 14) and multi-metric radar/ranking chart (Figure 15) summarize the historical evaluation of Section 3.1 in a complementary graphical form. Seasonal disaggregation (Figure 16) shows that all three scenarios track the observed seasonal means closely for DJF, MAM, and JJA, with the widest inter-model spread in SON. Dry- and wet-spell duration statistics (Figure 17) show mean consecutive dry-day spells lengthening modestly relative to the observed 3.8 days across all scenarios, with wider inter-model spread in SSP2-4.5. The percentile-bin analysis (Figure 18) shows a consistent pattern across models and scenarios: the lowest precipitation percentiles (P0–20) decrease substantially (ensemble mean −57% to −68%), while the upper percentiles (P95–100) are comparatively stable to slightly positive (+2% to +9% at the ensemble level, with individual models reaching +16%), consistent with a shift of probability mass away from light precipitation rather than a uniform intensification across the full distribution. Exceedance-probability curves (Figure 19) show that the frequency of days exceeding 10–50 mm day⁻¹ under future scenarios remains close to the observed frequency, with small absolute increases at higher thresholds under SSP2-4.5 and SSP5-8.5. The comprehensive statistics matrix (Figure 20) consolidates annual totals, spell statistics, and 50-year return levels across all model–scenario combinations into a single reference table.

The composite ranking in Figure 15b uses the mean of five normalized metrics (r, KGE, NSE, and inverted RMSE and Bias) as a broader visual summary of overall model fidelity, and can order models differently from Table 2, whose ranking is based on the three-metric criterion (KGE, NSE, and r ranks) defined in Section 2.3 and used for model selection throughout the study; both orderings place MPI-ESM1-2-LR first and CanESM5 among the lowest-performing models, and differ only in the relative position of CMCC-ESM2 and MIROC6.

4. Discussion

4.1. Key Findings in Context

Three principal findings emerge. First, the direction of ensemble-mean annual precipitation change is small and scenario-dependent (−1.3% to +5.1%), driven mainly by a consistent split between drying models (CMCC-ESM2, CanESM5) and wetting models (MIROC6, MPI-ESM1-2-LR) rather than by emission pathway. Second, extreme-event indicators intensify more consistently than the annual mean: Rx1day increases by 15–29% and 100-year return levels by 21–47% across scenarios, in line with regional syntheses documenting intensification of upper precipitation percentiles across the tropical Andes under warming [32]. Third, inter-model structural uncertainty dominates total projection variance in every calendar month (51–93%), substantially exceeding scenario uncertainty (under 1–19%).

By preserving each model's own projected daily chronology and the full distribution of observed precipitation intensity, the quantile delta mapping framework captures tail intensification that would not be visible from mean-change metrics alone. The projected +29% Rx1day increase under SSP5-8.5 falls within the 15–30% range documented in regional assessments for the tropical Andes at comparable levels of global warming [32]. Recent CMIP6-based assessments in other data-scarce and topographically complex settings have reported a similar pattern of pronounced extreme-index intensification alongside comparatively modest mean-annual change, in the Yangtze River basin [33] and across the Korean peninsula [34], reinforcing that return-period and extreme-index analyses, rather than mean-change metrics alone, should anchor climate-risk communication in this type of setting. The joint lengthening of the annual-maximum dry spell (CDD, +37–50%) alongside intensified 1-day and 5-day maxima is consistent with the compound-hazard framing of Zscheischler et al. [35], whereby reduced wet-day frequency combined with elevated per-event intensity can amplify flash-flood and soil-erosion risk in steep mountain watersheds, even where the annual water balance itself changes little.

The chronology-preserving quantile delta mapping procedure applied here is one of several established bias-correction strategies for daily precipitation [13]; it was selected over simpler mean-scaling and delta-change approaches because it corrects the full distribution rather than a single central-tendency statistic, and over more data-demanding multivariate or stochastic weather-generator approaches because it requires only single-site daily observations and single-realization GCM output, both of which are the practical constraint at this station. A direct quantitative comparison against alternative bias-correction families (e.g., multivariate approaches conditioned on atmospheric circulation analogs [36]) was outside the scope of the present single-station application and is identified as a priority direction in Section 4.5; the independent split-sample validation reported in Section 3.2 was carried out specifically to provide evidence of predictive skill in the absence of such a comparison.

4.2. Uncertainty Decomposition Implications

Inter-model structural uncertainty is the dominant component of total projection spread in all 12 calendar months, in contrast with mid-latitude findings where scenario uncertainty becomes comparable to or exceeds internal/model variability by mid-century [27]. This is consistent with the documented characteristics of tropical projections, where ENSO-driven decadal variability and inter-model differences in convective parameterization are large relative to the forced scenario signal. The practical implication is that, at this station, expanding the diversity of the GCM ensemble is likely to narrow decision-relevant uncertainty more effectively than further disaggregating emission scenarios.

4.3. Return Periods and Engineering Implications

The RP-100 increase from 73.2 mm day⁻¹ observed to 88.3–107.5 mm day⁻¹ across scenarios represents a 21–47% intensification of daily rainfall at this station under the return-period framework evaluated here. These are daily-resolution, single-station statistics, and their translation into specific drainage-infrastructure design values is not direct: municipal and hydraulic design standards are typically expressed in terms of sub-daily (e.g., 1-hour or shorter) intensity–duration–frequency curves, require the disaggregation of daily totals informed by local storm structure, and depend on spatially distributed precipitation inputs and hydrological or hydraulic modeling rather than a point estimate at a single station. The results presented here are accordingly best read as evidence of the direction and approximate magnitude of change in daily extreme precipitation at this station, providing a quantitative starting point for such downstream engineering analysis rather than a substitute for it. The bootstrap 95% confidence intervals reported in Table 6 give a probabilistic range for return-level estimates that explicitly incorporates GCM structural uncertainty, for use in that downstream analysis. The contrast between the two Mann–Whitney comparisons in Section 3.6 — non-significant on four-to-five point estimates, significant once bootstrap replicates are pooled across models — is itself an important practical finding: with an ensemble this small, whether a scenario comparison appears “significant” can depend as much on how the comparison sample is constructed as on the underlying climate signal, and any downstream design values should be drawn from the full multi-model distribution rather than from a single point-estimate comparison.

4.4. Limitations

Seven limitations qualify interpretation. First, the effective ensemble size is only four models for SSP1-2.6 and SSP5-8.5 and five for SSP2-4.5, following the unavailability of ACCESS-CM2 future data for two scenarios and of NorESM2-LR for the full ensemble; this constrains the statistical power of any cross-model or cross-scenario comparison. Second, the interaction term in the variance decomposition (Section 3.8) is confounded with single-realization noise and cannot be interpreted as an independent estimate of internal climate variability, which requires long, continuously evolving simulations to isolate. Third, SC-PREC4SA constitutes a single-point measurement, and the presented projections are representative of this station rather than of the surrounding region as a whole. Fourth, quantile delta mapping assumes that the statistical relationship between model and observed distributions used for bias correction remains valid under the level of warming considered; while chronology is preserved from the model's own future trajectory, the correction itself is a statistical transfer function whose stationarity cannot be independently verified. Fifth, the two-parameter Gumbel distribution, chosen over a three-parameter alternative for the reasons given in Section 2.7, may underestimate return levels if the true tail behavior of precipitation extremes at this station departs from the light-tailed Type I shape, particularly at return periods beyond those evaluated here (T > 100 years) [26]; a comparison against the (log-)Pearson Type III and GEV distributions, not undertaken here given the 30-year record length available, is identified as a priority check in Section 4.5. Sixth, the daily-resolution, single-station return levels reported here require sub-daily disaggregation, spatial extension, and hydrological or hydraulic modeling before they can inform specific drainage-infrastructure design values, as discussed in Section 4.3. Seventh, NorESM2-LR, IPSL-CM6A-LR, and EC-Earth3 were not available for inclusion in the present ensemble; their addition is identified as a priority for improving statistical power and inter-model diversity in future work.

4.5. Future Research Directions

Ensemble expansion to include NorESM2-LR, IPSL-CM6A-LR, and EC-Earth3 would improve statistical power for cross-model and cross-scenario comparisons. Comparison against multivariate bias-correction schemes conditioned on atmospheric circulation analogs [36] would help quantify the sensitivity of the results to the stationarity assumption discussed in Section 4.4. Spatial extension of the present pipeline to the broader SC-PREC4SA station network would enable distributed, rather than single-point, projections across the region, and would additionally allow direct comparison against CORDEX-CORE regional simulations as their daily precipitation archive for this domain becomes available (Section 2.2). Adoption of the generalized extreme value family or the (log-)Pearson Type III distribution with L-moment estimation would extend return-period reliability beyond the 100-year horizon evaluated here and test the sensitivity of the results to the Gumbel shape assumption noted in Section 4.4. Integration with a calibrated hydrological model, together with sub-daily disaggregation of the daily projections presented here, would be required before translating these station-level results into specific drainage-infrastructure design values.

5. Conclusions

A multi-model, multi-scenario quantile delta mapping framework was applied to five CMIP6 GCMs under SSP1-2.6, SSP2-4.5, and SSP5-8.5 at a high-altitude tropical station (2160 m a.s.l.), yielding five principal conclusions.

(1) GCM performance evaluation: MPI-ESM1-2-LR ranked first by composite KGE–NSE–Pearson r score (r = 0.830). All models exhibited positive wet biases of 24–264% in raw form. Negative NSE values are consistent with known limitations of coarse-resolution GCMs applied to point-scale observations over complex terrain; independent split-sample validation confirmed that the resulting bias correction achieves genuine out-of-sample skill (KGE improving from negative values to 0.60–0.76) rather than simply reproducing its own calibration data.

(2) Mean precipitation change: Ensemble-mean annual changes are modest and directionally mixed (−1.3% to +5.1%), governed principally by a consistent split between drying and wetting models rather than by emission scenario. No scenario supports a robust determination of net drying or wetting for the 2051–2080 window at this station.

(3) Extreme precipitation: ETCCDI indices indicate tail intensification that is more consistent across scenarios than the annual mean. Rx1day increases by 15–29%; the 100-year daily return level rises from 73.2 mm day⁻¹ to 88.3–107.5 mm day⁻¹. Translating these station-level, daily-resolution results into drainage-infrastructure design values would require sub-daily disaggregation, spatial extension beyond this single station, and hydrological or hydraulic modeling (Section 4.3). Concurrent lengthening of the annual-maximum dry spell characterizes a compound hazard of more intense but less evenly distributed precipitation.

(4) Uncertainty decomposition: Inter-model structural uncertainty is the dominant component of total projection variance in every calendar month (51–93%), substantially exceeding scenario uncertainty (under 1–19%). Expanding GCM ensemble diversity is likely to be a more effective lever for reducing decision-relevant uncertainty at this station than further disaggregating emission scenarios.

(5) Return-period sensitivity to comparison method: a Mann–Whitney comparison of SSP2-4.5 against SSP5-8.5 return levels is non-significant when applied to the small number of available model point estimates, but significant once bootstrap replicates are pooled across models — a methodological finding relevant to any small multi-model ensemble study of this kind, and a reminder that infrastructure design values should draw on the full multi-model distribution rather than a single pairwise comparison.

The documented pipeline, including its independent validation step, constitutes a transferable methodological framework for probabilistic precipitation projection at other data-scarce, high-altitude tropical stations across the Andes.

Key Findings:

• Non-monotonic annual response: ensemble-mean annual precipitation change does not scale with emission-pathway severity — SSP2-4.5 (intermediate) shows the largest projected increase (+5.1%), while SSP1-2.6 (optimistic) and SSP5-8.5 (pessimistic) show small, opposite-signed changes (−1.3% and −0.4%), reflecting a consistent split between drying and wetting models rather than a scenario effect.
• Extreme intensification outpaces the annual mean: despite a near-neutral annual-mean change under SSP5-8.5 (−0.4%), this scenario shows the strongest intensification of extremes among the three (Rx1day +28.9%; 100-year return level +47%, from 73.2 to 107.5 mm day⁻¹), confirming that extreme-index and return-period analysis, not mean-change metrics alone, should anchor risk assessment at this station .
• Return-level differences depend on how the comparison is constructed: a Mann–Whitney test on the four-to-five model point estimates alone finds no significant difference between SSP2-4.5 and SSP5-8.5 (p = 0.111–0.190), but pooling each model's bootstrap replicates (~2000 samples per scenario) reveals a statistically significant difference at every return period evaluated, with higher return levels under SSP5-8.5.
• Model uncertainty dominates: inter-model structural differences account for 51–93% (mean 75%) of total monthly projection variance, far exceeding scenario uncertainty (under 1–19%, mean 8%), indicating that expanding GCM ensemble diversity would narrow decision-relevant uncertainty more than further disaggregating emission scenarios.
• Validated downscaling: independent split-sample validation (calibration 1981–2000, evaluation on the withheld 2001–2010 period) confirms genuine out-of-sample skill of the quantile delta mapping procedure, with Kling–Gupta efficiency improving from negative values for every raw model to 0.60–0.76 after correction.