The NEX-GDDP-CMIP6 developed by NASA as an improvement of CMIP6. The dataset offers 35 models with 0.25° spatial resolution (approximately 25 km) to provide finer details than the original CMIP6 models that have native resolution typically around 100 km or coarser (Thrasher et al., 2022). The data archive contains downscaled 1950–2,100 historical and future projections that can be accessed from the website of NASA at https://www.nccs.nasa.gov/services/data-collections/land-based-products/nex-gddp-cmip6. NEX-GDDP-CMIP6 utilizes a daily variant of the monthly Bias Correction/Spatial Disaggregation (BCSD) algorithm for its statistical downscaling. BCSD corrects and downscales CMIP6 data to bridge the gap in coarse-scale model outputs (Jain et al., 2019). BCSD method employs in two steps. First, quantile mapping bias correction used to align GCM’s distributions with historical observations, removing systematic errors. Second, spatial disaggregation interpolates into 0.25° grid for improved spatial resolution.

This study used the daily temporal resolution of NEX-GDDP-CMIP6 precipitation dataset with the historical period of 1985 to 2014, accessed in September 2023. The dataset includes all historical simulations of 35 NEX-GDDP-CMIP6 models as shown in Table 1. Variant of the models identified by code specifying key attributes given by letter and number of rNiNpNfN including realization index, initialization method, parameterization, and forcing variant (Moradian et al., 2023). The realization (r) stands for certain running simulation, the initialization (i) stands for initial conditions, and the physics (p) stands for the model physics version. The forcing (f) denotes the version of forcing index employed. The variant r1i1p1f1 represents the first version of realization, first version of initialization, first version of physics and first forcing index simulation (Gummadi et al., 2025). These simulations offer a reliable basis for studying precipitation dynamics, climate impacts and adaptation.

ENSMEAN, in this study, employed to reduce the uncertainty and variability of individual models to a certain degree (Guo et al., 2023; Nwokolo et al., 2023). ENSMEAN model was generated after land-sea masking step by calculating the unweighted ensemble mean of monthly climatological output from 35 NEX-GDDP-CMIP6 models. An unweighted (simple) average was chosen for its simplicity and as a robust baseline for comparison. This approach does not account for the varying skill/performance of individual models, which a performance-based or competence-weighted ensemble would address.

The limitations of long-term observational records in Indonesia led to the adoption of the MSWEP (v2), which covered 1985–2014 of daily temporal and 0.1° spatial resolution in this study. The MSWEP was used as the observational reference against the NEX-GDDP-CMIP6. The dataset covers 1979 up to the present day (Beck et al., 2017) and has 0.1° spatial resolution. The MSWEP exemplifies systematic data merging by assimilating the strengths of satellite, gauge, and reanalysis-based daily precipitation estimates. The adoption of the dataset in this study was based on the detailed information presented by Beck et al. (2019). The MSWEP provides high spatial resolution of 0.1° and up to 3-hourly temporal resolution considered well-suited for regional hydroclimate studies (Nazarian et al., 2024). The dataset is a reliable source of gridded precipitation data for climatological and hydrological studies related to Indonesia (Ferijal et al., 2025).

35 NEX-GDDP-CMIP6 and ENSMEAN model employed to conduct statistical evaluation in this study. The precipitation dataset was assessed against the MSWEP for the historical period of 1985–2014 (Lakew, 2020). Monthly mean climatology was calculated for both datasets to establish a baseline for comparison (Lee et al., 2013). Spatial performance was examined through seasonal climatological means of December–January–February (DJF), March–April–May (MAM), June–July–August (JJA), and September–October–November (SON). Temporal performance was analyzed using the annual climatological cycle as shown in Figure 2.

This dual scale evaluation framework of seasonal for spatial analysis and monthly for temporal analysis is based on the physical drivers of Indonesia. Spatial evaluation is conducted on a seasonal basis to capture monsoon-driven migration of Intertropical Convergence Zone (ITCZ) that affected precipitation distribution across Indonesia archipelago (Aldrian and Dwi Susanto, 2003). Temporal evaluation utilizes monthly data to assess the model’s ability in simulating annual cycle. Specifically, temporal evaluation use in capturing the timing and magnitude of wet-dry months transitions which are critical to agricultural adaptation (Jun-Ichi et al., 2002). By evaluating spatio-temporal dynamics, this study ensure that the model ranking reflect the ability to simulate Indonesia’s unique seasonal geography and monthly transitions of precipitation cycles.

The initial calculation was the adoption of the daily precipitation data to determine the monthly average with the aim of ensuring consistency. For each month m in year y, the monthly mean precipitation P m , y was computed for all daily precipitation values. P ( d , m , y ) represents the daily precipitation on the day d and n m , y is the number of days in the month. Therefore, the monthly mean precipitation was determined using the following equation: P m , y = 1 n m , y ∑ d = 1 n m , y P ( d , m , y ) (1)

Equation (1) was used to calculate 30-year climatological monthly mean by averaging all months in the period. For a given month m, the climatological monthly mean was determined in Equation (2): P ¯ m = 1 M ∑ y = 1 n P m , y (2)

M = 30 which is the total number of years and the Equation (2) produces long-term monthly climatology to represent the typical precipitation characteristics in a month. Furthermore, the annual climatology was computed by calculating the annual mean precipitation in each year. The annual mean for year y was determined based on the 12 monthly means as shown in Equation (3): P ¯ y annual = 1 12 ∑ m = 1 12 P m , y (3)

The 30-year climatology annual mean was subsequently obtained by averaging the annual means across the period. This was achieved using Equation (4): P ¯ annual = 1 M ∑ y = 1 N P y annual (4)

The performance of the models was quantified using statistical metrics for the purpose of detecting systematic deviations. The metrics used were Correlation Coefficient (CC), Normalized Standard Deviation (NSTD), and Mean Bias (MB) (Ngo-Duc et al., 2024). The Root Mean Standard Deviation (RMSD) evaluation metrics used in the Taylor plot were also adopted (Liu et al., 2024). Taylor diagrams were applied to show the combination of CC, NSTD, and RMSD that demonstrated the ability of each model to reproduce the variability observed. The adoption of Taylor diagram was associated with it’s advantages in evaluating model performance by simultaneously showing the correlation, variability, and errors in a single framework (Simão et al., 2020). The dataset were processed in Python programming language using a specific module of the Open Climate Workbench (OCW) (H. Lee et al., 2018). In addition, SkillMetrics (https://pypi.org/project/SkillMetrics/1.1.3/) and Matplotlib 3.2.2 libraries also employed in data processing. The steps for NEX-GDDP-CMIP6 and MSWEP included spatial cropping and land-sea masking over Indonesia. The observational MSWEP dataset was re-gridded to 0.25° NEX-GDDP grid to ensure a grid-to-grid comparison at daily time step. Following these steps, evaluation metrics were calculated to assess model performance as provided in Table 1.

Specifically, MSWEP was used as the reference precipitation dataset R_i and the NEX-GDDP-CMIP6 as simulated precipitation Si. Both were used to determine the average of the reference dataset and simulated precipitation, evaluated over 30 years. The average simulations and observations were provided as S ¯ and R ¯ . Furthermore, the variable N used in the metric computation formula is defined differently for each evaluation scale. For the spatial (seasonal mean) evaluation, N is the total number of grid cells across the domain (as the comparison is grid-to-grid on a single seasonal mean field). For the temporal (monthly climatology) evaluation, N is the total number of grid cells multiplied by the number of months in the annual cycle (12 months), as the long-term monthly climatology is compared across the whole domain.

The formula for each metric is presented in the following equation:

Correlation Coefficient (CC), in Equation (5), was used to evaluate the magnitude and orientation of the linear relationship between R_i and Si: CC = ∑ i = 1 N ( S i − S − ) ( R i − R − ) ∑ i = 1 N ( S i − S − ) 2 ∑ i = 1 N ( R i − R − ) 2 (5)

Normalized Standard Deviation (NSTD) was used to determine the variability differences between R_i and Si through the ratio of standard deviation as shown in Equation (6): NSTD = σ Si σ Ri (6)

NSTD was employed to assess variability differences between the simulations of Si and the reference data of R_i by calculating the ratio of their standard deviation. This metric provides a dimensionless measure of relative variability. NSTD value close to 1 indicates that the simulations variability is comparable to the reference. Values higher or lower than unity represent overestimation or underestimation of the precipitation spread, respectively (Taylor, 2001).

Root Mean Standard Deviation (RMSD) evaluates the average deviation between R_i and Si values, combining both systematic and random errors, as given in Equation (7): RMSD = 1 N ∑ i = 1 N ( S i − R i ) 2 (7)

Mean Bias (MB) was applied to measure the average deviations between R_i and Si values in Equation (8). This was used to determine when a model over- or underestimated. MB = 1 N ∑ i = 1 N ( S i − R i ) (8)

The monthly precipitation cycles of 35 NEX-GDDP-CMIP6 and ENSMEAN models were compared with the reference dataset for 1985–2014 in Figure 3. The pink shades indicate the inter-model distribution, the wider bands represent higher uncertainty, and narrower bands signal stronger agreement. MSWEP annual variability as observed in the maximum precipitation of ~280–300 mm in January–March, minimum in July–August at ~150–170 mm, and an increase in September–December with ~280 mm. The figure showed that NorESM2-LM and MPI-ESM1-2-LR consistently overestimated the rainy season with the peaks reaching 320 mm in January–March to suggest a monsoonal pattern of precipitation. Meanwhile, CMCC-CM2-SR5 and MIROC6 indicated underestimated values in the dry season (JJA) below 150 mm. Thus, reflected the existence of dryness over the months. CanESM5, CESM2, and HadGEM3-GC31-LL exhibited values closer to MSWEP by maintaining seasonal amplitudes that were consistent with the observation data. BCC-CSM2-MR and GFDL-CM4 reproduced the seasonal shape but showed stronger deviations in April–May and October when precipitation tended to shift between rainy and dry seasons.

A critical observation in the annual cycle in Figure 3 is the consistent relative minimum in precipitation during February across most models, referred to as the ‘February dip’. This represents a systematic discrepancy, as the MSWEP reference dataset depicts a continuous peak during DJF. This model behavior likely stems from common biases in simulating the peak intensity of monsoonal flow in DJF. This underestimation reflects the challenge of simulating convection over the Indonesia complex topography (Peatman et al., 2014). The discrepancies also attributed to the model’s challenges in simulating large-scale monsoon circulation, land-sea contrast interactions and complex tropical dynamics (Mulsandi et al., 2024). Notably, the reference data of MSWEP documented to cause and underestimation of peak precipitation values (Beck et al., 2017).

The bias spread in the Taylor diagram compares the annual precipitation from NEX-GDDP-CMIP6 against the MSWEP. To reflect inter-model differences in seasonal precipitation, CC was scaled to range of 0–1. Visualization utilizes distinct symbols to represent individual models alongside the ENSMEAN. Figure 4 shows that majority of the models clustered close to the reference point. The CC scores exceeded 0.9 with low RMSD and the NSTD was close to one. The clusters showed that the models could replicate annual precipitation variability in Indonesia with reasonable accuracy. The best performance was determined based on the closeness of the scores to the reference point of MSWEP which included high CC, low RMSD, and near ideal NSTD. The results showed that ACCESS-CM2, CMCC-ESM2, MRI-ESM2.0, TaiESM1, and ENSMEAN had good performance by having positions close to the MSWEP. Meanwhile, poor performance was associated with high disparities compared to the reference point, determined using NSTD far from 1, higher RMSD, or weaker correlation. While most models cluster closely to the MSWEP reference, certain models exhibit deviations in variability. Specifically, ACCESS-ESM1-5 and KACE-1-0-G emerge as outliers with NSTD values exceeding 1.25 and 1.20, respectively. This indicates an overestimation of the observed precipitation variability. In contrast, GFDL-CM4 demonstrate high fidelity to observation, maintaining CC above 0.90 and an NSTD near unity.

The Taylor diagram for DJF, as shown in Figure 5a, demonstrated that most clusters of the models are close to the reference dataset. The CC score exceeding to 0.9, reflects the good performance in simulating seasonal precipitation variability during the period. ACCESS-CM2, CMCC-ESM2, MRI-ESM2.0, and ENSMEAN were in proximity to the reference point by showing high correlation, low RMSD, and NSTD close to one, suggested minimal bias in the intensity. The clustering pattern for MAM which is presented in Figure 5b is quite similar but the variances among all models become clearer. For example, TaiESM1, ACCESS-CM2, CMCC-ESM2, and ENSMEAN achieved CC scores close to 0.95–0.99 and recorded low RMSD. Thus, reflected the reliability of the models in reproducing precipitation variability and spatial distribution during the season.

In Figure 5c, most models cluster near the observation point for JJA, showing high correlations (CC > 0.9) and low RMSD. This indicated that most models effectively capturing dry season precipitation variability. ACCESS-CM2, CMCC-ESM2, TaiESM1, and ENSMEAN showed high correlation and minimal bias. The Taylor diagram for SON presented in Figure 5d shows that the models have greater spread compared to JJA. For example, ACCESS-CM2, CMCC-ESM2, MRI-ESM2.0, and ENSMEAN had high CC in the range of 0.95–0.99 and low RMSD, confirmed the model’s robustness in capturing seasonal precipitation patterns and distribution. KACE-1-0-G showed deviations in both NSTD and RMSD, reflected persistent deficiencies as previously observed in Figures 5a–d.

The multi-criteria assessments are widely applied in climate science to identify the best GCMs to capture climate variabilities and patterns (Nahar et al., 2017; Shakeel et al., 2025). Section 3.3 provides information on the multi-criteria analysis to identify the models. Figure 6 presents the heatmap ranking of the experiments. The scoring method adopted and described in Section 3.3 offered a comparative evaluation of climate models based on spatial and temporal performance. Figure 6 shows the validity of the models in temporal and spatial dimensions with the warmer shades signaling a higher trend. The four columns on the left represent the normalized value of each evaluation metrics including r CC _ norm , r NSTD _ norm , r RMSD _ norm and r MB _ norm for the temporal dimension while the four on the right are for the spatial dimension. The row order is based on the ranking of the models. It was observed that ACCESS-CM2, CMCC-ESM2, TaiESM1, MRI-ESM2-0, CESM2-WACCM and KIOST-ESM has high performance accuracy against the reference dataset.

ENSMEAN also showed high CC as identified in the warmer colors. The best models improved the representation of variability and reduction of errors to show superior reliability compared to the others. The Figure 6 also shows the ability of the models in replicating spatial distribution of precipitation compared to temporal variability. The mean bias reflected the widest spread to emphasize the persistent errors due to overestimation and underestimation.

Figure 7 depicts ACCESS-CM2 as the top-performing model across all dimensions. Even though trade-offs remained evident between temporal and spatial performance with some excelling in one dimension but weak in another. Temporal simulations measured the capacity of models to reproduce climate variability and long-term trends. The results showed that ACCESS-CM2, CMCC-ESM2, and MRI-ESM2.0 had good accuracy in replicating observed time series compared to the others. TaiESM1, KIOST-ESM, and CMCC-ESM2 also had the best performance in terms of spatial accuracy. Additionally, ACCESS-CM2 and CMCC-ESM2 emerged as the top-performing models in the composite evaluation of spatial and temporal performance, followed closely by TaiESM1 and MRI-ESM2.0. The evaluation ranked ACCESS-CM2, CMCC-ESM2, TaiESM1, MRI-ESM2.0, and CESM2-WACCM as the top five models, followed by KIOST-ESM in sixth and ENSMEAN in seventh.

The spatial distributions of the best ranking models including MSWEP, ACCESS-CM2 and ENSMEAN are compared in Figure 8. It was observed that the highest value of precipitation intensity during DJF exceeded 600 mm/month over some areas of Papua. ACCESS-CM2 exhibit better spatial performance than ENSMEAN during the season by providing closer spatial agreement to MSWEP as presented in heatmap (Figure 6). As shown in Figure 8, ACCESS-CM2 and ENSMEAN demonstrated comparable spatial performance over western and northern Indonesia, that captured the eastward expansion of wet condition during SON. ACCESS-CM2 also captured a reduction in precipitation magnitude over Sumatra, Java and Nusa Tenggara during MAM, compared to MSWEP. The same trend identified during JJA with magnitude less than 150 mm over Java, Nusa Tenggara, and southern Sumatra. Whereas ACCESS-CM2 exhibited slight overestimation in western Sumatra during JJA.

The temporal variation of precipitation across Indonesia can be classified into three climatic patterns of monsoonal, equatorial, and local (Aldrian and Dwi Susanto, 2003). Figure 9 compares monthly mean precipitation for three representative regions which include Jakarta region for monsoonal, Pontianak region for equatorial, and Sorong region for local. While the models generally replicate the phase of Indonesia climatic types, systematic deviations from the MSWEP (red dashed line) persist throughout the annual cycle.

In Jakarta monsoonal regime, model trajectories deviate from MSWEP in almost every month. Thus, making a general characterization of “good match” difficult to maintain. During DJF wet season, ACCESS-CM2 and MRI-ESM2.0 exhibit magnitude deviations, underestimating peak precipitation. While during JJA dry months, MRI-ESM2-0 and CMCC-ESM2 appear visually closest to the MSWEP. Other models tend to overestimate dry season precipitation, suggesting a superior fit for specific models during this period compared to the higher uncertainty in transition months.

In the equatorial regime of Pontianak, all models captured the bimodal precipitation pattern, although quantitative deviations remain. ENSMEAN is closely in line with the MSWEP, reproducing both intensity peaks with minimal bias. Compared to other individual simulations, ACCESS-CM2 captures the dual-peak structure but exhibits lower precipitation intensity than MSWEP. CMCC-ESM2 underestimates peak rainfall, while CESM2-WACCM shows a slight overestimation.

In Sorong local regime, precipitation variability is nearly uniform. Model fidelity shows marked variation peaking in JJA. During this period, MRI-ESM2-0 and CMCC-ESM2 appear most closely to the MSWEP. While ACCESS-CM2 and CESM2-WACCM show consistent monthly variation, CMCC-ESM2 exhibits a negative bias during the local rainy season. Although ENSMEAN follows the MSWEP pattern, inter-model differences in this region remain relatively small.