Paper deep dive
Understanding Perturbed Parameter Ensemble Sensitivities Using A Contrastive Learning Approach
Da Fan, David John Gagne, Gregory S Elsaesser, Brian Medeiros, Addisu G Semie, Qingyuan Yang, Akila Sampath, Subashree Venkatasubramanian
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 94%
Last extracted: 9/28/2026, 3:01:31 AM
Summary
This study develops an explainable contrastive learning model to analyze Perturbed Parameter Ensembles (PPEs) from the Community Atmosphere Model version 6 (CAM6). The model maps monthly cloud and radiation fields into a shared representation space to distinguish between two microphysics schemes: the default KK2000 and the neural network emulator TAU-ML. Results show the model achieves over 94% classification accuracy, preserves seasonal variability, and reveals that TAU-ML simulations are closer to satellite observations than KK2000. Integrated Gradients attributions link regional differences to specific physical parameters like cloud microphysics and deep convection.
Entities (8)
Relation Signals (8)
CAM6 → uses → KK2000
confidence 98% · The first warm rain parameterization is the default two-moment bulk microphysics scheme (Khairoutdinov & Kogan, 2000, hereafter KK2000)
CAM6 → uses → TAU-ML
confidence 98% · The other warm rain parameterization is the machine learning emulator of the Tel Aviv University (TAU) bin microphysical model (hereafter TAU-ML)
Contrastive Learning → distinguishes → KK2000
confidence 95% · The learned representations separates two PPEs with over 94% linear classification accuracy
Contrastive Learning → distinguishes → TAU-ML
confidence 95% · The learned representations separates two PPEs with over 94% linear classification accuracy
TAU-ML → iscloserto → Observations
confidence 94% · TAU-ML PPE has a lower distance to observations compared to KK2000 in the representation space
Integrated Gradients → identifiesregions → Subtropical low-cloud regions
confidence 93% · Integrated Gradients attributions highlights the contributions in subtropical low-cloud regions
Integrated Gradients → identifiesregions → Storm track regions
confidence 93% · Integrated Gradients attributions highlights the contributions in ... Northern and Southern Hemisphere storm track regions
Regional Attributions → →
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Perturbed parameter ensembles (PPEs) reveal how physics parameters affect climate simulations, but interpreting parameter sensitivities across multivariate, spatially structured outputs remains challenging, particularly when calibrating models against observations. We develop an explainable contrastive learning model that maps 5 monthly cloud and radiation fields into a shared representation space. We train the model on the fields of two 100-member Community Atmosphere Model version 6 (CAM6) PPEs, spanning 34 parameters, that only differ in the warm rain microphysics scheme: KK2000, the default bulk microphysics scheme, and TAU-ML, a neural network emulator of a bin microphysics scheme. The learned representations separates two PPEs with over 94\% linear classification accuracy while preserving the seasonal variability and ensemble spread due to parameter perturbations. In the shared representation space, the representations of satellite observations occupy the same low-dimensional manifold as the PPEs but are displaced from them most strongly during boreal spring and autumn. TAU-ML PPE has a lower distance to observations compared to KK2000 in the representation space. Integrated Gradients attributions highlights the contributions in subtropical low-cloud regions, Northern and Southern Hemisphere storm track regions, and tropical convection regions to differences between PPEs and observations. Regional attributions correlate most strongly with parameters associated with cloud microphysics, boundary layer turbulence, and deep convection. These results demonstrate that explainable representations of climate fields can attribute model differences to specific variables, regions, seasons, and physical parameters.
Tags
Links
- Source: https://arxiv.org/abs/2609.30420v1
- Canonical: https://arxiv.org/abs/2609.30420v1
Trouble viewing inline? Open PDF directly →
Full Text
74,584 characters extracted from source content.
Expand or collapse full text
manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES) Understanding Perturbed Parameter Ensemble Sensitivities Using A Contrastive Learning Approach Da Fan 1,2 , David John Gagne I 3 , Gregory S Elsaesser 4,5 , Brian Medeiros 3 , Addisu G Semie 1,2 , Qingyuan Yang 1,2 , Akila Sampath 1 , Subashree Venkatasubramanian 1 , 1 Learning the Earth with Artificial Intelligence and Physics (LEAP) National Science Foundation (NSF) Science and Technology Center, Columbia University, New York, NY, USA 2 Department of Earth and Environmental Engineering, Columbia University, New York, NY, USA 3 NSF National Center for Atmospheric Research, Boulder, CO, USA 4 NASA Goddard Institute for Space Studies, NY, USA 5 Department of Applied Physics and Applied Mathematics, Columbia University, New York, NY, USA • Contrastive learning separates two perturbed parameter ensembles (PPEs) in em- bedding space while preserving seasonal variability and ensemble spread due to parameter perturbations. • Attribution of embedding distances between PPEs and observations highlights sub- tropical low-cloud regions, storm tracks regions, and tropical convection regions. • Correlations between regional attributions and perturbed parameters indicate that embedding distances encode regional sensitivities of cloud and radiative fields to physics parameters. Corresponding author: Da Fan, df2944@columbia.edu –1– arXiv:2609.30420v1 [physics.ao-ph] 24 Sep 2026 manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES) Abstract Perturbed parameter ensembles (PPEs) reveal how physics parameters affect climate sim- ulations, but interpreting parameter sensitivities across multivariate, spatially structured outputs remains challenging, particularly when calibrating models against observations. We develop an explainable contrastive learning model that maps 5 monthly cloud and radiation fields into a shared representation space. We train the model on the fields of two 100-member Community Atmosphere Model version 6 (CAM6) PPEs, spanning 34 parameters, that only differ in the warm rain microphysics scheme: K2000, the default bulk microphysics scheme, and TAU-ML, a neural network emulator of a bin microphysics scheme. The learned representations separates two PPEs with over 94% linear classifi- cation accuracy while preserving the seasonal variability and ensemble spread due to pa- rameter perturbations. In the shared representation space, the representations of satel- lite observations occupy the same low-dimensional manifold as the PPEs but are displaced from them most strongly during boreal spring and autumn. TAU-ML PPE has a lower distance to observations compared to K2000 in the representation space. Integrated Gradients attributions highlights the contributions in subtropical low-cloud regions, North- ern and Southern Hemisphere storm track regions, and tropical convection regions to dif- ferences between PPEs and observations. Regional attributions correlate most strongly with parameters associated with cloud microphysics, boundary layer turbulence, and deep convection. These results demonstrate that explainable representations of climate fields can attribute model differences to specific variables, regions, seasons, and physical pa- rameters. Plain Language Summary Climate models cannot directly resolve many cloud and rainfall processes, so these processes are represented using simplified equations containing uncertain choices. Test- ing many choices produces a large collection of simulations, but comparing all of their spatial patterns and variables with observations is difficult. We trained a machine-learning method to organize monthly cloud-water and radiation fields from two groups of 100 cli- mate simulations. The groups used the same uncertain parameter values but different representations of how cloud droplets combine to form rain. The learned organization preserved the seasonal cycle and the spread caused by parameter choices while distin- guishing the two rain-formation methods. Both simulation groups differed more from ob- servations than from each other, although the more detailed rain representation was gen- erally closer to observations. Cloud liquid water and reflected sunlight explained most of the differences, especially over marine low-cloud regions and midlatitude storm tracks. The method also connected these regional differences to parameters controlling cloud for- mation, atmospheric turbulence, and deep convection. This approach can help model de- velopers determine not only which simulations agree better with observations, but also which physical processes and regions produce that agreement. 1 Introduction Earth system models (ESMs) connect a growing number of Earth system compo- nent models to enable higher fidelity simulations of past climate variability and future climate change. At climate timescales, uncertainties and biases arise more from the com- plex web of modeling choices made when configuring each ESM run. Many microphysics, turbulence, convection, and aerosol processes occur at scales smaller than the grid spac- ing of current ESMs and must therefore be represented through sub-grid parameteriza- tions (Flato, 2011; Schmidt & Sherwood, 2015; Gettelman et al., 2022). These param- eterizations contain uncertain parameters and structural assumptions that influence the simulated cloud, precipitation, radiation, and ultimately climate sensitivities (Sherwood et al., 2020; Tapiador et al., 2019). As models continue to increase in complexity and –2– manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES) resolution, we need more scalable and systematic ways to evaluate and constrain uncer- tainty in physical parameters and structural bias (Morrison et al., 2020; Fletcher et al., 2022). Climate model development has traditionally relied on hand tuning to reduce bi- ases in energy and hydrological cycle metrics. Although this approach has produced sub- stantial improvements, it is time-consuming, subjective, and increasingly difficult as the number of uncertain parameters, model outputs, and observational constraints grows (Hourdin et al., 2017; Mauritsen et al., 2012; Schmidt et al., 2017; Mignot et al., 2021). Moreover, multiple parameter combinations can produce similarly acceptable mean climate states, a condition often described as equifinality (Beven & Freer, 2001; Her et al., 2019). Agree- ment with a small set of scalar metrics also does not guarantee that the model correctly represents the underlying physical processes. A model may compensate for errors across variables, regions, or seasons, producing apparently improved global metrics while re- taining important structural deficiencies (Hourdin et al., 2017; G. S. Elsaesser et al., 2025). Perturbed parameter ensembles (PPEs) provide a powerful framework for explor- ing these uncertainties and deficiencies. By systematically perturbing parameters in cloud microphysics, boundary layer turbulence, shallow and deep convection, and aerosol schemes, PPEs reveal how uncertain parameter choices affect simulated climate (Sexton et al., 2021; Eidhammer et al., 2024; Duffy et al., 2024). They also provide a basis for generating cal- ibrated physics ensembles (CPEs), in which parameter combinations are selected to im- prove agreement with observations while retaining physically plausible solutions (G. S. El- saesser et al., 2025). However, interpreting PPEs remains difficult because model out- put is high-dimensional, multivariate, spatially structured, and seasonally dependent. Dif- ferences among simulations, or between simulations and observations, may reflect pa- rameter uncertainty, structural error, internal variability, observational uncertainty, or compensating biases among fields (Braverman et al., 2021; Yang et al., 2024; G. S. El- saesser et al., 2025). A key challenge is therefore to develop diagnostics that compare models and ob- servations in a compact but physically meaningful space. Such diagnostics should pre- serve regional and seasonal structure, represent interactions among variables, distinguish model configurations, quantify discrepancies between the model and observations, and identify which variables, regions, and parameters control those differences. This need is especially important for CPEs, where the goal is not only to reduce a scalar error met- ric but also to determine whether calibration improves the cloud and radiation processes (Hourdin et al., 2023; Watson-Parris et al., 2021; Yarger et al., 2024; G. S. Elsaesser et al., 2025). Machine learning offers new tools for addressing this problem. Contrastive learn- ing is a self-supervised representation learning approach that trains a neural network to place similar samples closer and dissimilar samples farther apart in a compact embed- ding space (van den Oord et al., 2019; Hjelm et al., 2019; Chen et al., 2020; He et al., 2020). Related developments have shown that learned representations can capture sim- ilarities and differences among samples without requiring a manually specified distance metric (Khosla et al., 2021; Tian et al., 2020). For climate applications, this approach is attractive because it can learn compact representations of multivariate spatial or spa- tiotemporal fields while retaining information relevant to downstream prediction, clas- sification, bias correction, and physical interpretation (Ballard & Erinjippurath, 2022; L. Wang et al., 2022; Lv et al., 2022; Gong et al., 2024; Chatterjee et al., 2023; C. Wang et al., 2025; Bailey, 2026). If trained appropriately, the resulting embedding space might encode seasonal variability, ensemble spread, and differences among model configurations and observations. Distances in this space can then be interpreted as data-driven mea- sures of similarity among climate states. –3– manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES) In this study, we develop a contrastive learning framework to quantify sensitivi- ties to physical parameters and diagnose structural deficiencies in climate model sim- ulations and support the development of parameter calibration work (Hourdin et al., 2017; Mignot et al., 2021; Yang et al., 2024; G. S. Elsaesser et al., 2025). We apply the frame- work to two CAM6 PPEs that differ in their warm rain microphysics schemes: one uses the default K2000 scheme, and the other uses the TAU-ML bin scheme (described in further detail later). This paper addresses three questions specifically: 1. Can contrastive learning construct a representation space that captures seasonal variability, parametric uncertainty, and structural differences among TAU-ML PPE, K2000 PPE, and observations? 2. Can distances in this representation space provide useful measures of parametric uncertainty and structural error? 3. Which variables, regions, and physical parameters control the variability in the representation and physical space? The remainder of the paper is organized as follows. Section 2 describes the CAM6 PPEs, warm rain parameterizations, and observational datasets. Section 3 introduces the contrastive learning framework, embedding evaluation methods, and integrated gra- dient attribution approach. Section 4 presents the learned embedding structure, seasonal embedding distances, parameter sensitivities, attribution maps, and regional parame- ter sensitivities. Section 5 summarizes the main conclusions and discuss their implica- tions for parameter calibration. 2 Data Here we describe the CAM6 model used (Section 2.1), the PPE simulations (Sec- tion 2.2), and observational datasets (Section 2.3). 2.1 CAM6 model description The Community Atmosphere Model Version 6 (CAM6) is the atmospheric com- ponent of the Community Earth System Model version 2 (Danabasoglu et al., 2020). CAM6 features a two-moment stratiform cloud microphysics scheme (Gettelman & Morrison, 2015, MG2;) with prognostic cloud liquid, cloud ice, rain, and snow hydrometeors. For these hydrometeors, MG2 predicts both mass and number concentration. MG2 is cou- pled to a unified moist turbulence scheme, Cloud Layers Unified by Binormals (Golaz et al., 2002; Larson et al., 2002), implemented by Bogenschutz et al. (2013). CLUBB treats boundary layer moist turbulence, and grid-scale cloud processes. CAM6 features a 4-mode aerosol model (Liu et al., 2016) with an extra mode for primary carbon and a represen- tation of natural and anthropogenic aerosols. CAM6 also features an ensemble plume mass flux deep convection scheme (Neale et al., 2008; Zhang & McFarlane, 1995) with simple microphysics, a physically based mixed phase ice nucleation scheme (Hoose et al., 2010), and the Rapid Radiative Transfer Model for General Circulation Models (Iacono et al., 2000, RRTMG;) radiation scheme. All CAM6 (version 6 4046) experiments use a standard 0.9 ◦ ×1.25 ◦ (latitude and longitude) horizontal resolution with 32 vertical levels (minimum at 2 hPa) and are con- figured with present day climatological forcing, where Sea Surface Temperature (SSTs), sea ice, and boundary forcings are prescribed from climatological conditions centered on year 2000. –4– manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES) 2.2 PPE simulations We leverage two 100-member PPEs using CAM6 with the only difference being the warm rain microphysics parameterizations. The two microphysics parameterizations are more completely described in (Gettelman et al., 2021), but here we repeat key informa- tion. The first warm rain parameterization is the default two-moment bulk microphysics scheme (Khairoutdinov & Kogan, 2000, hereafter K2000), which represents warm rain formation process by autoconversion and accretion. K2000 defines autoconversion and accretion rate by fitting empirical power law to large eddy simulations (LES) with bin- resolved microphysics. The autoconversion rate components have been adjusted to bet- ter match observations (Gettelman et al., 2019). The other warm rain parameterization is the machine learning emulator of the Tel Aviv University (TAU) bin microphysical model (hereafter TAU-ML). The TAU model resolves the drop size distrbutions using 35 bins, which allows for mass transfer between bins and alleviates anomalous drop growth. The autoconversion and accretion rates in K2000 are replaced by explicitly solving the quasistochastic collection process. How- ever, the quasistochastic collection process is computationally expensive for practical global climate and weather model simulations. A neural network was used in Gettelman et al. (2021) to emulate the stochastic collection tendency from cloud and rain water variables, which shows consistent results with the TAU model and a variant of the neural network (Gagne et al., 2026) is used in our study as TAU-ML. Table 1.Perturbed parameter names, descriptions, default values, rangesm and units for the 100 member PPE. ParameterDescriptionDefaultMinMaxUnit micro_mg_max_nicons Maximum allowed ice number concentration1.00E+081.00E+051.00E+10#/m3 micro_mg_vtrmi_factorIce fall speed scaling1.000.205.00- micro_mg_iaccr_factor Scaling ice/snow accretion1.000.201.00- micro_mg_berg_eff_factorBergeron efficiency factor1.000.101.00- micro_mg_accre_enhan_factAccretion enhancing factor1.000.1010.00- micro_mg_homog_sizeHomogeneous freezing ice particle size2.50E-051.00E-052.00E-04meter micro_mg_dcsAutoconversion size threshold ice-snow5.00E-045.00E-051.00E-03meter clubb_c1Low Skewness in C1 Skw.1.000.403.00- clubb_C8Coef. #1 in C8 skewness Equation4.201.005.00- clubb_c11Low Skewness in C11 Skw0.700.200.80- clubb_c14Constant for u’2 and v’2 terms2.200.404.00- clubb_c_K10Momentum coefficient of Kh_zm0.500.201.20- clubb_gamma_coefLow Skw.: gamma coef. Skw0.3080.2500.350- clubb_wpxp_L_threshLscale threshold, damp C6 and C760.0020.00200.00- clubb_betaSet plume widths for theta_l and rt2.401.602.50- clubb_c2rtDamping on scalar variances1.000.202.00- zmconv_c0_lndConvective autoconversion over land0.00750.0020.1- zmconv_c0_ocnConvective autoconversion over ocean0.030.020.1- zmconv_keConvective evaporation efficiency5.00E-061.00E-061.00E-05- zmconv_dmpdzEntrainment parameter-1.00E-03-2.00E-03-2.00E-04- cldfrc_dp1Parameter for deep convection cloud fraction0.100.050.25- cldfrc_dp2Parameter for deep convection cloud fraction500.00100.001000.00- zmconv_num_cinAllowed number of negative buoyancy crossings115- zmconv_tiedke_addConvective parcel temperature perturbation0.500.002.00- zmconv_capelmtTriggering threshold for ZM convection70.0035.00350.00- seasalt_emis_scaleSeasalt emission scaling factor1.000.502.50- dust_emis_factDust emission scaling factor0.700.101.00- sol_factb_interstitialBelow cloud scavenging of interstitial modal aerosols0.100.101.00- sol_factic_interstitialIn-cloud scavenging of interstitial modal aerosols0.400.101.00- microp_aero_wsub_minMin subgrid velocity for liquid activation0.200.000.50- microp_aero_wsubi_min Min subgrid velocity for ice activation0.0010.0000.200- microp_aero_wsub_scaleSubgrid velocity for liquid activation scaling1.000.105.00- microp_aero_wsubi_scale Subgrid velocity for ice activation scaling1.000.105.00- microp_aero_npccn_scaleScale activated liquid number1.000.333.00- The two PPE simulations for K2000 and TAU-ML consist of 100 ensemble mem- bers with 34 parameters (Table 1) sampled by Latin hypercube sampling. These param- eters span the turbulence (CLUBB; Golaz et al., 2002), cloud microphysics (MG2; Get- –5– manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES) telman & Morrison, 2015), aerosol (MAM; Liu et al., 2012), and deep convection (ZM; Zhang & McFarlane, 1995). Their ranges are determined with expert elicitation (their Table 1, Eidhammer et al., 2024). The same parameter sets are used in both PPEs, but note that the accretion parameter (micro mgaccreenhanfact) is not activated when the TAU-ML scheme replaces K2000. We work with the CAM6 PPE present-day sim- ulations. Each PPE member simulation is run using near present day cyclic boundary conditions of the year 2000. The greenhouse gases and atmospheric oxidants take the average values of the 1995-2005 period, the average monthly sea surface temperatures (SSTs) during 1995-2010 are used, and the emission of aerosols and precursors is set to 1995-2005 in the simulations. Each simulation runs for a period of 3 years and the out- put is archived at monthly intervals for climatological fields. We focus on 5 variables, which characterize the cloud and radiative responses to changes in microphysics and pa- rameters, of the CAM6 PPE in this work: shortwave cloud radiative effect (SWCF), long wave cloud radiative effect (LWCF), upwelling longwave flux at top of model (FLUT), net shortwave flux at top of model (FSNT), and total (cloud plus rainwater) gridbox- average liquid water path (TLWP). These variables were selected because they exhib- ited substantial differences compared to observations (Gettelman et al., 2021, their Fig. 12) and strong sensitivities to parameter changes in the CAM PPE (Eidhammer et al., 2024, their Fig. 3). Four static variables, land-sea mask, surface height, latitude, and longi- tude, are also included in the inputs to provide geographic and surface context. 2.3 Observations Table 2. Monthly climatological observational products used in this work. CAM output name Observational product nameDescriptionReferences SWCFmonthly_swcrf_toaMonthly shortwave cloud radiative forcing at the top of atmosphere (TOA) from CERES-EBAF-Ed4.2, 2002-2023 climatology. Loeb et al. (2018); Kato et al. (2018) LWCFmonthly_lwcrf_toaMonthly longwave cloud radiative forcing at TOA from CERES-EBAF-Ed4.2, 2002-2023 climatology. Loeb et al. (2018); Kato et al. (2018) FLUTmonthly_olrMonthly outgoing longwave radiation at TOA from CERES-EBAF-Ed4.2, 2002-2023 climatology. Loeb et al. (2018); Kato et al. (2018) FSNTmonthly_swabs_toaMonthly absorbed shortwave flux at TOA from CERES-EBAF-Ed4.2, 2002- 2023 climatology. Loeb et al. (2018); Kato et al. (2018) TLWPmonthly_tlwpMonthly total (cloud + rainwater) liquid water path from the MAC-LWP product, based on passive microwave satellite observations during 1988-2016. Elsaesser et al. (2017) To evaluate PPE outputs, we use the observational monthly climatologies of LWCF, SWCF, FLUT, FSNT (Loeb et al., 2018; Kato et al., 2018), and TLWP (G. S. Elsaesser et al. (2017)), derived from products in Table 2, as the metrics of cloud and top-of-atmosphere radiation fields. Observational uncertainty has limited impacts on the results and is there- fore not considered in this study, but will be incorporated into future calibration work. 3 Methods 3.1 Data Processing The dataset comprises monthly samples from the K2000 and TAU-ML PPEs and observational climatologies between 65.5 ◦ S and 65.5 ◦ N. The polar regions are excluded due to the large observational uncertainties in cloud retrievals over snow and ice covered surfaces. Each sample includes anomalies of five monthly mean variables (SWCF, LWCF, FLUT, FSNT, and TLWP) and four static fields (latitude, longitude, land-sea mask, and surface height). At each grid point, we compute the anomalies by standardizing each vari- able separately, using three-year mean and standard deviation for each PPE member and the mean and standard deviation across 12 monthly climatologies for observations. We exclude the first month of each PPE simulation because ensemble spread is initially re- –6– manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES) duced as members diverge from the same initial conditions. The PPE samples are split by ensemble member into training (members 1–60; 4200 samples), validation (members 61–80; 1400 samples), and testing (members 81–100; 1400 samples) subsets. 3.2 Contrastive Learning framework Figure 1.Schematic of the contrastive learning framework and explainable AI workflow. (a) Contrastive learning architecture showing the training and inference pipelines and a visualization of the two-dimensional projections of the learned embeddings. The model uses five atmospheric variables and four static fields; only SWCF is shown for illustration. (b) Shared neural network encoder and projection head that map input from different data sources into a common embed- ding space. (c) Explainable AI workflow that uses the Integrated Gradient method to attributes cosine distances between embeddings to geographic regions of the input fields. We develop a self-supervised contrastive learning framework (Fig. 1) that maps cli- mate model outputs and observations to low-dimensional representations in a shared em- bedding space. Contrastive learning (Chen et al., 2020) is a discriminative approach that –7– manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES) rewards pulling similar samples closer in the embedding space while also pushing dis- similar samples apart. Fig. 1a shows the contrastive learning model architecture. The model is trained on mini-batches, with each containing two pairs of samples: (A t ,A t+1 ) from a TAU-ML ensemble member and (B t ,B t+1 ) from its corresponding K2000 member. The two mem- bers share the same perturbed physics parameters and both pairs span the same adja- cent months, t and t+1. The pairs from the same microphysics scheme, (A t ,A t+1 ) and (B t ,B t+1 ), are treated as positive pairs because the adjacent monthly outputs within each pair represent consecutive segments of the same dynamical trajectory and there- fore expected to share large scale physical characteristics. The pairs from different schemes, (A t ,B t ), (A t ,B t+1 ), (A t+1 ,B t ), and (A t+1 ,B t+1 ), are considered as negative pairs be- cause their samples represent segments from different dynamic trajectories generated with different schemes. In each positive pair, the two monthly samples represent distinct, par- tial views of the same dynamical trajectory from a single TAU-ML or K2000 member. By contrast, the two samples in each negative pair represent partial views of different trajectories. By contrasting the embeddings of positive pairs and negative pairs, the model learns a representation that captures similarities within a dynamical trajectory and dif- ferences between trajectories. We use a shared neural network (N) architecture (Fig. 1b), including a shared encoder f (·) and a projection head g(·), to project input data x into embeddings that represent radiation and liquid water characteristics. The shared encoder comprises a channel- aware masked convolutional layer followed by three same-padding convolutional layers, all with 3×3 kernels. The masked convolution is applied to the TLWP channel to ex- tract features only over water regions, with the output renormalized to remain compa- rable to the other channels. Hyperparameters were selected through manual tuning. The number of convolutional filters increases after every convolutional layer, from 64 to 128 and then to 256, with the last two convolutional layers both using 256 filters. Gaussian Error Linear Unit (GELU; Hendrycks & Gimpel, 2016) activation function is applied af- ter each convolutional layer. Two average pooling layers, each with a 3×3 kernel and a stride of 4, downsample the feature maps after every two convolutional layers. The fi- nal feature map is flattened into a one-dimensional feature vector, h = f (x) and passed through the projection head. The projection head applies three dense layers with 2048, 256, and 128 units, to reduce the dimension of the embedding space. GELU activations are applied after the first two dense layers and the final dense layer generates a 128-dimensional embedding z = g(h). We apply the shared encoder and projection head to map inputs from TAU-ML and K2000 datasets into a common embedding space. During training, the model is optimized using the normalized temperature-scaled cross-entropy (NT-Xent) loss, that groups embeddings of positive pairs together while distinguishing those of neg- ative pairs. Each minibatch contains N = 2 positive pairs and four negative pairs, for a total of 2N = 4 samples. Let sim(z i ,z j ) = z ⊤ i z j ∥z i ∥ 2 ∥z j ∥ 2 (1) denote the cosine similarity between embeddings z i and z j , equivalent to the dot prod- uct of their ℓ 2 -normalized representations. For each sample i, let p(i) denote the other sample in its positive pair. The NT-Xent loss for a minibatch is L batch =− 1 2N 2N X i=1 log exp sim(z i ,z p(i) )/τ 2N X k=1 I [k̸=i] exp[sim(z i ,z k )/τ ] ,(2) where I [k̸=i] equals 1 when k ̸= i and 0 otherwise, and τ is the temperature hyperpa- rameter. The loss is averaged over all four samples in the minibatch, with each sample serving once as the anchor. –8– manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES) We trained the contrastive learning model with minibatches. The loss was com- puted for each minibatch and gradients were accumulated over 100 minibatches before each optimizer update, corresponding to an effective batch size of 400. We used the Adam optimizer (Kingma & Ba, 2015) with a learning rate of 0.0001 and weight decay of 10 −5 for gradient descent. Early stopping was applied when the validation loss does not im- prove over 10 epochs. After training, we applied the model to extract embeddings from K2000, TAU- ML, and observation datasets. We then applied two approaches to evaluate the robust- ness of embeddings. First, we trained a logistic regression model using the embeddings to separate TAU-ML and K2000 classes. The skill of the logistic regression in separat- ing the classses provided initial evidence that the contrastive learning procedure was learn- ing a robust signal. Second, we applied principle component analysis (PCA) to reduce the 128-dimensional embeddings from all data sources to 3 dimensions for visualization. We then compared the projected embeddings across data sources to assess the robust- ness of the learned representation. 3.3 Explainable AI We employed the Integrated Gradients (Sundararajan et al., 2017) to explain the learned embedding space and identify differences among different datasets (Fig. 1b). The Integrated Gradients method generates heatmaps that represent the relative contribu- tions of geographic regions for climate applications (Mamalakis et al., 2022). Each heat map was normalized by its absolute maximum value and averaged across all ensemble members and months to highlight the structural differences. Integrated gradients pro- duces similar attribution patterns to the model-agnostic SHAP method but is much faster to compute. 4 Results 4.1 Embedding evaluation Figure 2.Linear (TAU-ML vs K2000) classification evaluation of the learned embeddings from different layers of the contrastive learning model. –9– manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES) We first evaluate whether the learned representations contain information that dis- tinguishes K2000 and TAU-ML samples, while also capturing variability arising from parameter perturbations and structural errors relative to observations. We trained lo- gistic regression classifiers using embeddings extracted by each network layer from TAU- ML and K2000 training samples and evaluated them on testing samples. Figure 2 shows that embedding-based classification skill remains high across all layers, with accuracies over 0.94. The highest skill occurs at the dense2 layer, whereas the lowest skill is at the last dense3 layer. Notably, embeddings from the first conv1 layer already achieve high classification accuracy, suggesting that the earliest learned representations contain in- formation for distinguishing between datasets. The limited skill improvement from the first conv1 layer to the deeper dense2 layer suggests that intermediate layers provide only incremental information for TAU-ML and K2000 separation. The accuracy decrease from dense2 to dense3 may reflect the role of the final dense layer in compressing the rep- resentation into a lower-dimensional latent space optimized for the contrastive loss rather than for classification. This compression may discard some discriminative information while retaining features that better preserve embedding similarity. In the following anal- ysis, we focus on embeddings from pool2 and dense3 as representative outputs of the en- coder and projection head, respectively. The learned embedding space provides a compact, low-dimensional representation of the input fields. Figures 3 and 4 show pairwise two-dimensional (2D) PCA projections of TAU-ML, K2000, and observation embeddings generated by the encoder and last projection layer respectively. For the encoder embeddings (Fig. 3), the PC1-PC2 pro- jection shows a circular structure with around 12 clusters, corresponding to the seasonal cycle. Within each cluster, TAU-ML and K2000 embeddings are generally well sepa- rated with some overlap, consistent with its high classification skill in Fig 2. The spread of TAU-ML and K2000 embeddings within each cluster reflects ensemble variability due to parameter perturbations. Across the PCA projections, and particularly in the PC1- PC2 panel, observation embeddings are distributed along the same low-dimensional man- ifold as the K2000 and TAU-ML embeddings, rather than forming an isolated cluster. This suggests that the learned representations provide a common embedding space in which model outputs and observations can be compared directly. In the PC1-PC2 and PC3-PC4 projections, some observation embeddings are displaced from the correspond- ing TAU-ML and K2000 clusters, indicating structural error in the model outputs rel- ative to observations. The last layer embeddings (Fig. 4) show a more compact structure than the en- coder embeddings. The first two PCs account for nearly all of the embedding variance, while PC3 and PC4 together explain only 0.23%. The PC1-PC2 projection forms a ra- dial pattern with around 12 arms extending from the center, corresponding to the sea- sonal cycle. The TAU-ML and K2000 embeddings are more densely mixed near the center, where embeddings from adjacent months are closer to each another. This sug- gests that the contrastive loss preserves similarity within positive pairs of samples from adjacent months. In this projection, TAU-ML and K2000 embeddings overlap substan- tially, and observation embeddings fall on the same manifold. In contrast, PC3 shows a clear separation between TAU-ML and K2000 despite accounting for a small frac- tion of the total variance. This separation is evident in the PC1-PC3 and PC2-PC3 pro- jections, as well as in the PC3 histogram, where TAU-ML projections are shifted towards higher values and K2000 projections toward lower values. These behaviors are consis- tent with the effect of the contrastive loss, which emphasizes features that distinguish negative pairs of K2000 and TAU-ML samples while preserving similarity within pos- itive pairs of samples at adjacent months. In the following analysis, we focus on the encoder embeddings because they pre- serve a more general representation of the input fields, whereas the last layer embeddings –10– manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES) Figure 3.Pairwise two-dimensional (2D) PCA projections of the learned encoder (pool2) embeddings for TAU-ML, K2000, and observations. Diagonal panels show distributions for each principal component (PC), while off-diagonal panels show pairwise projections among PC1, PC2, PC3, and PC4. Blue circles denote TAU-ML samples, red circles denote K2000 samples, and yellow triangles denote observations. emphasize discriminative and similarity features that are optimized by the contrastive loss. The monthly 2D PCA projections of the encoder embeddings in Figure 5 further show that the learned representation captures physically meaningful variability while re- taining temporal coherency. Embeddings from each month form a distinct cluster in the PC1-PC2 space, and these clusters rotate clockwise with month, reflecting the seasonal cycle. Within each monthly cluster, the TAU-ML and K2000 samples generally occupy similar regions, with their spread representing ensemble variability arising from param- eter perturbations. The observation embeddings generally follow the same seasonal pro- gression as the model embeddings, suggesting that the encoder places observations and model outputs on a shared low-dimensional manifold. However, observations are displaced from the TAU-ML and K2000 clusters in several months in spring and autumn, includ- ing March, April, September, and October. These offsets indicate months in which CAM exhibits structural error relative to observations. These results suggest that the encoder –11– manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES) Figure 4. Same as Fig 3, but for the last layer (dense3) embeddings. embeddings capture ensemble spread associated with parametric uncertainty, structural differences relative to observations, and seasonal variability in the input fields. 4.2 Embedding distance analysis To better understand differences in the embedding space, we analyzed the cosine distance among monthly climatology embedding for TAU-ML, K2000, and observations (Fig. 6). The monthly mean distances show a clear seasonal dependence. The distances between K2000 and TAU-ML are generally smaller than the distances between either K2000 or TAU-ML and observations. This indicates that differences between the two warm rain microphysics in CAM are smaller than the overall discrepancy between CAM and observations. The TAU-ML and observation distances are consistently lower than the K2000 and observation distances, suggesting that TAU-ML has a smaller bias relative to ob- servations than K2000. This improvement likely results from more complex represen- tations of cloud and rain water in TAU-ML. The K2000 and observation distances are largest in spring and autumn, especially in March, April, September, and October, in- dicating seasons when K2000 differs most from observations. TAU-ML also shows rel- atively large distances from observations during these months, suggesting persistent sea- –12– manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES) Figure 5.Monthly two-dimensional PCA projections of the learned encoder embeddings for TAU-ML, K2000, and observations. Blue points represent TAU-ML outputs, red points represent K2000 outputs, and yellow triangles represent observations. The colorbars indicate ensemble member index for TAU-ML and K2000. Figure 6.Monthly variation in mean cosine distance among encoder embeddings of monthly climatology for TAU-ML, K2000, and observations. Shading represents standard deviation of the ensemble. sonal model biases. In addition, the K2000 and observation distances exhibit larger en- semble spread than the TAU-ML and observation distances, indicating that K2000 bi- ases are more sensitive to physical parameters than TAU-ML. –13– manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES) Figure 7.Top-15 parameter importance for latent distances of K2000 a TAU-ML, K2000 versus observations, and TAU-ML versus observations. Red bars indicate positive correlation while blue bars indicate negative correlation. Figure 7 shows the top-15 parameter importance for embedding distances of K2000 versus TAU-ML, K2000 versus observations, and TAU-ML vs observations. For K2000 versus TAU-ML, the top-3 most important parameters are the accretion enhancement factor (micro mgaccreenhanfact), the convective available potential energy (CAPE) triggering threshold for deep convection (zmconvcapelmt), and convective parcel tem- perature perturbation zmconvtiedkeadd. This indicates the major difference between K2000 and TAU-ML are their different responses to the accretion and deep convection triggering processes, consistent with the major difference between the two warm rain mi- crophysics schemes. For K2000 versus observations, the top-3 most important param- eters are accretion (micro mgaccreenhanfact), the damping of scalar variances for liquid water (clubbc2rt), and the autoconversion size threshold of cloud ice to snow (micro mgdcs). This indicates the major biases of K2000 are in the representation of autoconversion and accretion process. For TAU-ML versus observations, the top-3 most important parameters are the ice fall speed scaling (micromgvtrmifactor), the au- toconversion size threshold of cloud ice to snow (micropmgdcs), and convective par- cel temperature perturbation (zmconvtiedkeadd). This suggests that the TAU-ML bi- ases results from the autoconversion and deep convection process. 4.3 Embedding distance explanation To explain the embedding distances, we applied integrated gradients to identify the variables and geographic regions that contribute to the separation among TAU-ML, K2000, and observations in the learned embedding space. Fig. 8 shows the mean attribution maps for pairwise embedding distances averaged across months and ensemble members. The attributions are dominated by TLWP and SWCF, indicating that cloud processes and their radiative effects dominate the embedding space differences among the datasets. For the TAU-ML and K2000 comparison (Fig. 8a), the TLWP attributions are concentrated over the northeast Pacific off western North America, the North Atlantic storm track region, and the Southern Hemisphere storm track regions, with weaker sig- nals over the subtropical South Pacific Ocean off the western South America. The SWCF attributions show strong localized contributions over the subtropical North Pacific, the southeastern Pacific off South America, the North Atlantic, and the southern Indian Ocean. FSNT only shows weak signals in the tropical ocean. LWCF also shows weak attribu- tions in Southern and Northern Hemisphere storm track regions. FLUT contributions are negligible over the globe. –14– manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES) Figure 8.Mean integrated gradient attribution maps for embedding distances among TAU- ML, K2000, and observations across months and members. Rows show (a) TAU-ML versus K2000, (b) TAU-ML versus observations, and (c) K2000 versus observations. Columns show TLWP, SWCF, FSNT, LWCF, and FLUT. Positive and negative values indicate regions where each input variable contributes to the latent distance between the datasets. The attributions for TAU-ML and observation comparison (Fig. 8b) and K2000 and observation comparison (Fig. 8c) show a generally similar structure but with stronger TLWP attributions in the Southern Pacific Ocean and Southern Hemisphere storm track regions and stronger SWCF attributions in the tropics. Overall, the attribution maps indicate that the embedding distances are dominated by cloud processes over oceans, es- pecially in the Southern and Northern Hemisphere storm track regions and subtropical ocean regions. We further analyzed the seasonal variation of the attributions and their connec- tions to cloud and radiative features. Figs. 9-11 show the seasonal integrated gradient maps for SWCF and TLWP, and the seasonal differences in SWCF, TLWP, and low cloud fraction for the TAU-ML and K2000, TAU-ML and observation, and K2000 and ob- servation comparisons, respectively. For the TAU-ML and K2000 comparison (Fig. 9), the SWCF attributions are con- centrated in subtropical low cloud regions with large low cloud fraction differences. In subtropical ocean regions, TAU-ML generates more low clouds than K2000, leading to stronger shortwave feedback, as shown in the SWCF difference and attribution. This suggests that SWCF attributions identify regions where changes in warm rain microphysics affect SWCF through their influences on low cloud processes. The locations of the strongest signals vary seasonally with the spatial pattern of low cloud differences. The southeast- ern Pacific stratocumulus region shows stronger positive attribution in DJF and JJA than in the other seasons, while the northeastern Pacific signal is stronger in DJF. The TLWP attributions show strongest signals in the midlatitude storm track re- gions, especially over the North Pacific, North Atlantic, and Southern Hemisphere storm track regions. The signals in the Northern Hemisphere storm track region are stronger in JJA and SON, while the signals in the Southern storm track region are strongest in MAM. This indicates that TAU-ML and K2000 embedding distance varies seasonally with TLWP differences in midlatitude storm track regions. The TLWP and SWCF dif- ferences over the tropical oceans only show moderate signals in TLWP and SWCF at- tributions. This suggests that the differences between TAU-ML and K2000, identified by the embedding distance, are dominated by low cloud systems, directly affected by warm rain microphysics, with a smaller contribution from tropical deep convection. –15– manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES) Figure 9.Seasonal integrated gradient attribution maps of SWCF and TLWP and differences of TLWP, SWCF, and low cloud fraction for the TAU-ML and K2000 comparison. Columns show December-January-February (DJF), March-April-May (MAM), June-July-August (JJA), and September-October-November (SON). Rows show, from top to bottom, SWCF integrated gradient (IG) attribution, SWCF difference, TLWP IG attribution, TLWP difference, and low cloud fraction difference. The differences are calculated as TAU-ML minus K2000. Positive and negative IG values indicate regions where each input variable contributes to the latent distance between the K2000 and TAU-ML. Figure 10. Same as Fig. 9, but for TAU-ML and observations comparison For the TAU-ML and observation comparison (Fig. 10), the SWCF attributions are concentrated in subtropical low cloud regions and tropical convective regions. TAU- ML shows low cloud biases relative to observations over subtropical and tropical oceans, leading to weaker shortwave feedback from low clouds than observations. However, SWCF –16– manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES) differences show TAU-ML generates more SWCF than observations in these regions, sug- gesting that there are other biases leading to negative SWCF biases while partially off- set by low cloud fraction biases. The SWCF attributions identify regions with signifi- cant SWCF biases, especially over tropical oceans. This suggests that SWCF attribu- tions identify regions affected by the biases in deep convection, a structural issue in the CAM model. The locations of the strongest signals show seasonal variations. The north- eastern Pacific signal is stronger in DJF and SON, while Indian Ocean signals are stronger in MAM and JJA. Tropical ocean signals are stronger in DJF and MAM. The TLWP attributions show strong signals in the midlatitude storm track regions and subtropical low cloud regions. The TLWP attributions show seasonal variations, con- sistent with the low cloud bias variation. It suggests that TLWP attributions reflect the low cloud bias. The Southern Hemisphere storm track regions show stronger signals in MAM and JJA, while the North Pacific and South Pacific signals are stronger in DJF. This suggests that the differences between TAU-ML and observations, indicated by their embedding distance, are dominated by both low cloud systems and tropical deep con- vection. Figure 11. Same as Fig. 9, but for TAU-ML and observations comparison For the K2000 and observation comparison (Fig. 11), the SWCF attributions are concentrated in subtropical low cloud regions and tropical convective regions. K2000 shows widespread low cloud biases relative to observations over tropical and subtropi- cal oceans, leading to much weaker shortwave feedback than observations, as shown in the SWCF differences and SWCF attributions. This suggests that SWCF attributions identify regions where cloud fraction biases affect shortwave radiation. The locations of the strongest signals vary seasonally with the spatial pattern of low cloud differences. The northeastern Pacific signal is stronger in DJF and SON, while the southeastern Pa- cific and southern Indian Ocean signals are stronger in MAM and JJA. The TLWP attributions show strong signals in the midlatitude storm track regions and subtropical low cloud regions. The seasonal variations of TLWP attributions in the subtropical and midlatitude regions are consistent with the low cloud fraction differences. It suggests that TLWP attributions mainly reflect the low cloud bias. The Southern Hemi- sphere storm track shows stronger signals in MAM and JJA, while the North Pacific and North Atlantic signals are stronger in SON. K2000 generally underestimates TLWP –17– manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES) relative to observations over broad tropical and subtropical ocean regions, likely due to the low cloud bias. This suggests that K2000 and observation differences, indicated by their embedding distance, are dominated by significant low cloud biases in the subtrop- ical ocean and storm track regions. –18– manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES) 4.4 Regional parameter importance We further examine the relationship between regional integrated gradient and the perturbed parameters. Fig. 12 shows the Pearson correlations between individual param- eters and the regional attributions for selected variables and regions. Positive correla- tions indicate that increasing parameter values would increase the attribution in the re- gion, while negative correlations indicate that increasing parameter values would decrease the attribution in the region. Figure 12 suggests that many parameters appear to have little impact for all se- lected regions and variables. Significant correlations are found in MG microphysics pa- rameters, CLUBB turbulence parameters, and several ZM convection parameters. For the TAU-ML and K2000 comparison (Fig. 12a), the accretion enhancement factor (micro mgaccreenhanfact) shows strong positive correlations with TLWP in the southern Ocean storm track regions, SWCF in the subtropical stratocumulus regions, and FSNT in tropical Indian Ocean. The autoconversion size threshold of cloud ice to snow (microp mgdcs) shows positive correlation with FSNT in the tropical oceans and FLUT in the North Pacific storm track regions. The CLUBB parameter clubbc14 shows positive cor- relations with TLWP and SWCF in the subtropical stratocumulus regions. The damp- ing of scalar variances for liquid water (clubb c2rt) is important in the southern storm track regions. The convective available potential energy (CAPE) triggering threshold for deep convection (zmconvcapelmt), and convective parcel temperature perturbation zmconv tiedkeadd have opposite impacts on TLWP in the subtropical stratocumulus regions and North Atlantic storm track region. This suggests that the TAU-ML and K2000 differences are dominated by parameters related to low cloud processes and deep con- vection. The TAU-ML and observation comparison (Fig. 12b) highlights the im- portance of microphysics and turbulence parameters. The scaling ice and snow accretion (micromgbergefffactor) and Bergeron efficiency factor (micro mgiaccrfactor) show negative correlations with TLWP in Southern and Northern storm track regions. The autoconversion size threshold of cloud ice to snow (microp mgdcs) shows positive correlations with FSNT in tropical oceans and LWCF and SWCF in the Northern and Southern storm track regions, and neg- ative correlations with SWCF in the Indian Ocean. The clubb C8 shows negative correlations with SWCF and TLWP in the subtropical stratocumulus regions. The damping of scalar variances for liquid water (clubb c2rt) shows negative correla- tions with TLWP and SWCF in the southeast Atlantic stratocumulus region, FSNT in the tropical oceans, and LWCF and FLUT in Northern and Southern storm track regions. The convective available potential energy (CAPE) triggering threshold for deep convection (zmconvcapelmt), and convective parcel temperature perturbation zmconvtiedkeadd show positive correlations with FSNT in the tropical Pacific ocean. This suggests that the bias of TAU-ML is dominated by parameter associated with autoconversion and accretion processes in storm track regions and subtropical stratocumulus regions and deep convective processes in the tropical oceans. The K2000 and observations comparison (Fig. 12c) shows a similar pattern of parameter sensitivity to the TAU-ML and observation comparison. The main difference is the accretion enhancement factor(micro mgaccreenhanfact), which is only used by K2000 microphysics. The accretion process shows positive cor- relations with TLWP in the Southern ocean storm track region and SWCF in the stratocumulus region, and negative correlations with TLWP in the subtropical stra- tocumulus regions and Northern Hemisphere storm track regions. This suggests that the bias of K2000 is also dominated by parameter associated with autoconver- sion and accretion processes in storm track regions and subtropical stratocumulus regions. –19– manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES) 5 Conclusions This study developed a contrastive learning framework for diagnosing paramet- ric uncertainties and structural deficiencies among climate model simulations and observations. We applied the framework to two CAM6 perturbed parameter ensem- bles in which corresponding members share the same perturbed parameter settings but differ in their warm rain microphysics schemes: the default K2000 parame- terization vs. TAU-ML, a neural-network emulator of the Tel Aviv University bin microphysics model. By training on monthly TLWP, SWCF, LWCF, FSNT, FLUT, and static geographic fields, the model learned a common embedding space in which K2000, TAU-ML, and observations can be compared using embedding distances. The learned embeddings captured differences and variability within the input data. The embeddings separated TAU-ML and K2000 samples with high linear classification accuracy, indicating that the network learned robust structural differ- ences between the two warm rain microphysics. At the same time, the encoder em- beddings retained physically meaningful organization: monthly samples formed co- herent clusters that followed the seasonal cycle, while the spread within each cluster reflected ensemble variability associated with perturbed parameters. Observations generally lie on the same low-dimensional manifold as the model simulations, show- ing that the learned representation provides a common basis for model-observation comparison. However, observations were displaced from the model clusters during the spring and autumn, indicating more systematic atmospheric model discrepancies in seasons with spatially-shifting storm tracks. Embedding distance analysis showed that the difference between TAU-ML and K2000 is smaller than the difference between either model configuration and obser- vations. TAU-ML is generally closer to observations than K2000, suggesting that the more detailed representation of warm-rain collection processes reduces the model bias. However, both configurations exhibit substantial discrepancies from observa- tions, particularly in spring and autumn. The larger ensemble spread in the distance between K2000 and observations further indicates that K2000 biases are more sensitive to perturbed physical parameters than TAU-ML biases. The latent distances are physically explainable. Integrated gradient attribution showed that TLWP and SWCF dominate the separation among TAU-ML, K2000, and observations, while FSNT, LWCF, and especially FLUT contribute more weakly. The strongest attribution signals occur over low cloud and deep convection regions, including subtropical low cloud regions, the northeast and southeast Pacific, the North Atlantic storm track, the southern Indian Ocean, tropical Ocean, and the broader Southern Hemisphere storm track region. Seasonal attribution maps further show that these signals vary coherently with differences in TLWP, SWCF, and low cloud fraction. This indicates that the learned embedding distances are not arbitrary statistical separations, but are tied to cloud liquid water, low-cloud structure, and their shortwave radiative effects. The comparison between TAU-ML and K2000 is dominated by low-cloud and storm-track regimes that are directly sensitive to warm rain microphysics. TAU-ML produces more low clouds in several subtropical ocean regions, leading to short- wave cloud-forcing differences that are reflected in the SWCF attributions. TLWP attributions are strongest in midlatitude storm track regions, especially the North Pacific, North Atlantic, and Southern Hemisphere storm tracks. In contrast, model and observation comparisons show broader and stronger attributions, including both subtropical low cloud regions and tropical convective regions. This suggests that the remaining model biases involve not only warm-rain microphysics but also coupled errors in low-cloud, storm-track, and deep-convective processes. –20– manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES) The sensitivity analysis of parameters and attributions provides an additional link between the learned representation and model physics. Significant correlations are concentrated in cloud microphysics, CLUBB boundary layer turbulence and shal- low convection, and ZM deep-convection parameters. For the TAU-ML and K2000 comparison, accretion and deep-convection triggering parameters strongly influence the latent separation, consistent with the structural difference between the warm- rain schemes. For the model and observation comparisons, parameters related to accretion, autoconversion, ice and snow processes, boundary layer turbulence, and convective triggering control regional attributions in storm-track, stratocumulus, and tropical ocean regions. Many aerosol-emission, solubility, and dust parameters show weaker and less systematic relationships, indicating that the dominant embedding distances in this analysis are primarily controlled by cloud and precipitation process sensitivities. Overall, these results show that contrastive learning can provide a compact and physically explainable diagnostic space for climate model evaluation. The framework captures seasonal variability, ensemble spread, structural differences between model configurations, and discrepancies between models and observations. When combined with integrated gradients, the embedding distance can be traced back to the vari- ables and regions responsible for the separation. This makes the approach useful not only for detecting that two simulations differ, but also for identifying where and why those differences arise. These findings have direct implications for future calibration. A scalar global error metric may identify parameter settings that appear closer to observations, but it does not necessarily reveal whether the improvement reflects physically meaningful changes or compensating errors. The embedding distance provides a multivariate metric that preserves spatial structure, seasonal dependence, and coupled cloud- radiation behavior. The attribution and parameter-correlation analyses then identify the cloud regimes and process parameters responsible for the discrepancy between model and observation. Future work should extend this framework to uncertainty- aware observational comparisons, and calibration strategies that use the learned embedding distance as an objective function. Acknowledgments This work was supported by NSF funding through the Learning the Earth with Artificial intelligence and Physics (LEAP) Science and Technology Center (STC; Award No. 2019625). This work was also supported by the National Science Foun- dation (NSF) National Center for Atmospheric Research, which is a major facility sponsored by NSF under Cooperative Agreement No. 1852977. The authors thank Dr. Michael Pritchard for contribution of ideas and William A. Manriquez for dis- cussions on the analysis. 6 Data Availability Statement The datasets used in this work are made available on Zenodo (https:// doi.org/10.5281/zenodo.22943628). Upon acceptance of the manuscript, the repository will be officially published and a permanent DOI will be provided. The CERES EBAF data were obtained from the NASA Langley Research Center CERES ordering tool at https://ceres.larc.nasa.gov/data/. MAC-LWP data are available at https://doi.org/10.5067/MEASURES/MACLWPM (G. Elsaesser et al., 2016). –21– manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES) References Bailey, N.(2026).Contrastive learning boosts deterministic and generative models for weather data. Retrieved from https://arxiv.org/abs/2603.24744 Ballard, T., & Erinjippurath, G. (2022). Contrastive learning for climate model bias correction and super-resolution. Retrieved from https://arxiv.org/abs/2211 .07555 Beven, K., & Freer, J.(2001).Equifinality, data assimilation, and uncertainty estimation in mechanistic modelling of complex environmental systems us- ing the GLUE methodology.Journal of Hydrology, 249 (1–4), 11–29.doi: 10.1016/S0022-1694(01)00421-8 Bogenschutz, P. A., Gettelman, A., Morrison, H., Larson, V. E., Craig, C., & Scha- nen, D. P. (2013, December). Higher-order turbulence closure and its impact on climate simulations in the community atmosphere model. J. Clim., 26 (23), 9655–9676. Retrieved from http://dx.doi.org/10.1175/JCLI-D-13-00075.1 doi: 10.1175/jcli-d-13-00075.1 Braverman, A., Hobbs, J., Teixeira, J., & Gunson, M. (2021). Post hoc uncertainty quantification for remote sensing observing systems.SIAM/ASA Journal on Uncertainty Quantification, 9 (3). doi: 10.1137/19M1304283 Chatterjee, D., Acquistapace, C., Deneke, H., & Crewell, S.(2023).Under- standing cloud systems’ structure and organization using a machine’s self- learning approach.Artificial Intel ligence for the Earth Systems.doi: 10.1175/AIES-D-22-0096.1 Chen, T., Kornblith, S., Norouzi, M., & Hinton, G.(2020).A simple framework for contrastive learning of visual representations. In Proceedings of the 37th in- ternational conference on machine learning (p. 1597–1607).Retrieved from https://arxiv.org/abs/2002.05709 Danabasoglu, G., Lamarque, J.-F., Bacmeister, J., Bailey, D., DuVivier, A., Ed- wards, J., . . . others(2020).The community earth system model ver- sion 2 (cesm2).Journal of Advances in Modeling Earth Systems, 12 (2), e2019MS001916. Duffy, M. L., Medeiros, B., Gettelman, A., & Eidhammer, T. (2024). Perturbing pa- rameters to understand cloud contributions to climate change.Journal of Cli- mate, 37 (1), 213–227. doi: 10.1175/JCLI-D-23-0250.1 Eidhammer, T., Gettelman, A., Thayer-Calder, K., Watson-Parris, D., Elsaesser, G., Morrison, H., . . . McCoy, D.(2024, 11).An extensible perturbed parameter ensemble for the community atmosphere model version 6.Geoscientific Model Development, 17 , 7835-7853.Retrieved from https://gmd.copernicus.org/ articles/17/7835/2024/ doi: 10.5194/gmd-17-7835-2024 Elsaesser, G., O’Dell, C., Lebsock, M., & Teixeira, J.(2016).Multisensor ad- vanced climatology mean liquid water path l3 monthly 1 degree x 1 degree v1. Greenbelt, MD, USA: Goddard Earth Sciences Data and Information Services Center (GES DISC).Retrieved from https://doi.org/10.5067/MEASURES/ MACLWPM (Dataset) doi: 10.5067/MEASURES/MACLWPM Elsaesser, G. S., O’Dell, C. W., Lebsock, M. D., Bennartz, R., Greenwald, T. J., & Wentz, F. J.(2017).The multisensor advanced climatology of liquid water path (MAC-LWP).Journal of Climate, 30 (24), 10193–10210.doi: 10.1175/JCLI-D-16-0902.1 Elsaesser, G. S., van Lier-Walqui, M., Yang, Q., Kelley, M., Ackerman, A. S., Fridlind, A. M., . . . Strong, J. D. O. (2025). Using machine learning to gener- ate a GISS ModelE calibrated physics ensemble (CPE). Journal of Advances in Modeling Earth Systems, 17 (4), e2024MS004713. doi: 10.1029/2024MS004713 Flato, G. M. (2011). Earth system models: An overview. WIREs Climate Change, 2 (6), 783–800. doi: 10.1002/wcc.148 Fletcher, C. G., McNally, W., Kuchera, E., Liu, Y., Xie, S., Klein, S. A., & Zhang, M. (2022). Toward efficient calibration of higher-resolution earth system mod- –22– manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES) els.Journal of Advances in Modeling Earth Systems, 14 (5), e2021MS002836. doi: 10.1029/2021MS002836 Gagne, I., David John, Gettelman, A., Chuang, W., Fan, D., & Medeiros, B. (2026). ”an updated machine learning warm rain process emulator for the community atmosphere model”. (”Manuscript in preparation”) Gettelman, A., Gagne, D. J., Chen, C. C., Christensen, M. W., Lebo, Z. J., Morrison, H., & Gantos, G.(2021, 2).Machine learning the warm rain process.Journal of Advances in Modeling Earth Systems, 13 .doi: 10.1029/2020MS002268 Gettelman, A., Geer, A. J., Forbes, R. M., Carmichael, G. R., Feingold, G., Pos- selt, D. J., . . . Zuidema, P.(2022).The future of earth system prediction: Advances in model-data fusion.Science Advances, 8 (14), eabn3488.doi: 10.1126/sciadv.abn3488 Gettelman, A., Hannay, C., Bacmeister, J. T., Neale, R. B., Pendergrass, A., Dan- abasoglu, G., . . . others(2019).High climate sensitivity in the community earth system model version 2 (cesm2).Geophysical Research Letters, 46 (14), 8329–8337. Gettelman, A., & Morrison, H. (2015). Advanced two-moment bulk microphysics for global models. part i: Off-line tests and comparison with other schemes. Jour- nal of Climate, 28 (3), 1268–1287. Golaz, J.-C., Larson, V. E., & Cotton, W. R. (2002). A pdf-based model for bound- ary layer clouds. part i: Method and model description.Journal of the atmo- spheric sciences, 59 (24), 3540–3551. Gong, Y., He, T., Chen, M., Wang, B., Nie, L., & Yin, Y. (2024). Spatio-temporal enhanced contrastive and contextual learning for weather forecasting.IEEE Transactions on Knowledge and Data Engineering, 36 (8), 4260–4274.doi: 10.1109/TKDE.2024.3362825 He, K., Fan, H., Wu, Y., Xie, S., & Girshick, R.(2020).Momentum contrast for unsupervised visual representation learning. In Proceedings of the ieee/cvf con- ference on computer vision and pattern recognition (p. 9729–9738). Retrieved from https://arxiv.org/abs/1911.05722 Hendrycks, D., & Gimpel, K. (2016, 6). Gaussian error linear units (gelus). arXiv Preprint. Retrieved from http://arxiv.org/abs/1606.08415 Her, Y., Yoo, S.-H., Cho, J., Hwang, S., Jeong, J., & Seong, C.(2019).Un- certainty in hydrological analysis of climate change: Multi-parameter vs. multi-GCM ensemble predictions.Scientific Reports, 9 , 3229.doi: 10.1038/s41598-019-41334-7 Hjelm, R. D., Fedorov, A., Lavoie-Marchildon, S., Grewal, K., Bachman, P., Trischler, A., & Bengio, Y.(2019).Learning deep representations by mu- tual information estimation and maximization.Retrieved from https:// arxiv.org/abs/1808.06670 Hoose, C., Kristj ́ansson, J. E., Chen, J.-P., & Hazra, A. (2010, August). A classical- theory-based parameterization of heterogeneous ice nucleation by mineral dust, soot, and biological particles in a global climate model.J. Atmos. Sci., 67 (8), 2483–2503.Retrieved from http://dx.doi.org/10.1175/2010JAS3425.1 doi: 10.1175/2010jas3425.1 Hourdin, F., Ferster, B., Deshayes, J., Mignot, J., Musat, I., & Williamson, D. (2023). Toward machine-assisted tuning avoiding the underestimation of uncer- tainty in climate change projections. Science Advances, 9 (29), eadf2758. doi: 10.1126/sciadv.adf2758 Hourdin, F., Mauritsen, T., Gettelman, A., Golaz, J.-C., Balaji, V., Duan, Q., . . . Williamson, D.(2017).The art and science of climate model tun- ing.Bul letin of the American Meteorological Society, 98 (3), 589–602.doi: 10.1175/BAMS-D-15-00135.1 –23– manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES) Iacono, M. J., Mlawer, E. J., Clough, S. A., & Morcrette, J.-J.(2000, June). Impact of an improved longwave radiation model, RRTM, on the energy budget and thermodynamic properties of the NCAR community climate model, CCM3.J. Geophys. Res., 105 (D11), 14873–14890.Retrieved from http://dx.doi.org/10.1029/2000JD900091 doi: 10.1029/2000jd900091 Kato, S., Rose, F. G., Rutan, D. A., Thorsen, T. E., Loeb, N. G., Doelling, D. R., . . . Ham, S.-H. (2018). Surface irradiances of edition 4.0 clouds and the earth’s radiant energy system (CERES) energy balanced and filled (EBAF) data product.Journal of Climate, 31 (11), 4501–4527.Retrieved from https:// doi.org/10.1175/JCLI-D-17-0523.1 doi: 10.1175/JCLI-D-17-0523.1 Khairoutdinov, M., & Kogan, Y.(2000).A new cloud physics parameterization in a large-eddy simulation model of marine stratocumulus.Monthly Weather Review, 128 (1), 229–243. doi: 10.1175/1520-0493(2000)128⟨0229:ANCPPI⟩2.0 .CO;2 Khosla, P., Teterwak, P., Wang, C., Sarna, A., Tian, Y., Isola, P., . . . Krishnan, D. (2021). Supervised contrastive learning. Retrieved from https://arxiv.org/ abs/2004.11362 Kingma, D. P., & Ba, J.(2015).Adam: A method for stochastic optimization. In Iclr 2015.San Diego, CA.Retrieved from http://arxiv.org/abs/1412 .6980 Larson, V. E., Golaz, J.-C., & Cotton, W. R.(2002, December).Small-scale and mesoscale variability in cloudy boundary layers: Joint probability density functions.J. Atmos. Sci., 59 (24), 3519–3539.Retrieved from http://dx.doi.org/10.1175/1520-0469(2002)059<3519:SSAMVI>2.0.CO;2 doi: 10.1175/1520-0469(2002)059⟨3519:ssamvi⟩2.0.co;2 Liu, X., Easter, R. C., Ghan, S. J., Zaveri, R., Rasch, P., Shi, X., . . . others (2012). Toward a minimal representation of aerosols in climate models: Description and evaluation in the community atmosphere model cam5. Geoscientific Model Development, 5 (3), 709–739. Liu, X., Ma, P.-L., Wang, H., Tilmes, S., Singh, B., Easter, R. C., . . . Rasch, P. J.(2016, February).Description and evaluation of a new four-mode version of the modal aerosol module (MAM4) within version 5.3 of the community atmosphere model.Geosci. Model Dev., 9 (2), 505–522.doi: 10.5194/gmd-9-505-2016 Loeb, N. G., Doelling, D. R., Wang, H., Su, W., Nguyen, C., Corbett, J. G., . . . Kato, S. (2018). Clouds and the earth’s radiant energy system (CERES) en- ergy balanced and filled (EBAF) top-of-atmosphere (TOA) edition-4.0 data product. Journal of Climate, 31 (2), 895–918. doi: 10.1175/JCLI-D-17-0208.1 Lv, Q., Li, Q., Chen, K., Lu, Y., & Wang, L. (2022). Classification of ground-based cloud images by contrastive self-supervised learning.Remote Sensing, 14 (22), 5821. doi: 10.3390/rs14225821 Mamalakis, A., Barnes, E. A., & Ebert-Uphoff, I.(2022, 10).Investigating the fidelity of explainable artificial intelligence methods for applications of con- volutional neural networks in geoscience.Artificial Intel ligence for the Earth Systems, 1 . doi: 10.1175/aies-d-22-0012.1 Mauritsen, T., Stevens, B., Roeckner, E., Crueger, T., Esch, M., Giorgetta, M., . . . Tomassini, L.(2012).Tuning the climate of a global model.Jour- nal of Advances in Modeling Earth Systems, 4 (3), M00A01.doi: 10.1029/ 2012MS000154 Mignot, J., Hourdin, F., Deshayes, J., Boucher, O., Gastineau, G., Musat, I., . . . Madec, G.(2021).The tuning strategy of IPSL-CM6A-LR.Jour- nal of Advances in Modeling Earth Systems, 13 (5), e2020MS002340.doi: 10.1029/2020MS002340 Morrison, H., van Lier-Walqui, M., Fridlind, A. M., Grabowski, W. W., Harrington, J. Y., Hoose, C., . . . Xue, L.(2020).Confronting the challenge of modeling –24– manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES) cloud and precipitation microphysics.Journal of Advances in Modeling Earth Systems, 12 (8), e2019MS001689. doi: 10.1029/2019MS001689 Neale, R. B., Richter, J. H., & Jochum, M.(2008, November).The impact of convection on ENSO: From a delayed oscillator to a series of events.J. Clim., 21 (22), 5904–5924.Retrieved from http://dx.doi.org/10.1175/ 2008JCLI2244.1 doi: 10.1175/2008jcli2244.1 Schmidt, G. A., Bader, D., Donner, L. J., Elsaesser, G. S., Golaz, J.-C., Hannay, C., . . . Saha, S.(2017).Practice and philosophy of climate model tuning across six u.s. modeling centers.Geoscientific Model Development, 10 , 3207–3223. doi: 10.5194/gmd-10-3207-2017 Schmidt, G. A., & Sherwood, S. (2015). A practical philosophy of complex climate modelling. European Journal for Philosophy of Science, 5 (2), 149–169. doi: 10 .1007/s13194-014-0102-9 Sexton, D. M. H., McSweeney, C. F., Rostron, J. W., Yamazaki, K., Booth, B. B. B., Murphy, J. M., . . . Karmalkar, A. V. (2021). A perturbed parameter ensemble of HadGEM3-GC3.05 coupled model projections: Part 1: Selecting the parameter combinations.Climate Dynamics, 56 (11–12), 3395–3436. doi: 10.1007/s00382-021-05709-9 Sherwood, S. C., Webb, M. J., Annan, J. D., Armour, K. C., Forster, P. M., Har- greaves, J. C., . . . Zelinka, M. D.(2020).An assessment of earth’s climate sensitivity using multiple lines of evidence.Reviews of Geophysics, 58 (4), e2019RG000678. doi: 10.1029/2019RG000678 Sundararajan, M., Taly, A., & Yan, Q. (2017). Axiomatic attribution for deep net- works. Proceedings of the 34th International Conference on Machine Learning, 3319–3328. Retrieved from https://arxiv.org/abs/1703.01365 Tapiador, F. J., Navarro, A., Moreno, R., S ́anchez, J. L., & Garc ́ıa-Ortega, E. (2019).Regional climate models: 30 years of dynamical downscaling.At- mospheric Research, 235 , 104785. doi: 10.1016/j.atmosres.2019.104785 Tian, Y., Sun, C., Poole, B., Krishnan, D., Schmid, C., & Isola, P.(2020).What makes for good views for contrastive learning? Advances in Neural Information Processing Systems, 33 , 6827–6839. Retrieved from https://arxiv.org/abs/ 2005.10243 van den Oord, A., Li, Y., & Vinyals, O.(2019).Representation learning with contrastive predictive coding.Retrieved from https://arxiv.org/abs/ 1807.03748 Wang, C., Yang, N., & Li, X.(2025).Advancing forecasting capabilities: A contrastive learning model for forecasting tropical cyclone rapid intensi- fication.Proceedings of the National Academy of Sciences, 122 (4).doi: 10.1073/pnas.2415501122 Wang, L., Li, Q., & Lv, Q.(2022).Self-supervised classification of weather sys- tems based on spatiotemporal contrastive learning.Geophysical Research Let- ters, 49 (15), e2022GL099131. doi: 10.1029/2022GL099131 Watson-Parris, D., Williams, A., Deaconu, L., & Stier, P. (2021). Model calibration using ESEm v1.0.0: An open, scalable earth system emulator.Retrieved from https://arxiv.org/abs/2108.10747 Yang, Q., Elsaesser, G. S., van Lier-Walqui, M., & Eidhammer, T. (2024). A simple emulator that enables interpretation of parameter-output relationships, ap- plied to two climate model PPEs.Retrieved from https://arxiv.org/abs/ 2410.00931 Yarger, D., Wagman, B., Shand, L., & Chowdhary, K.(2024).Autocalibration of the E3SM version 2 atmosphere model using a PCA-based surrogate for spatial fields. Retrieved from https://arxiv.org/abs/2308.06877 Zhang, G. J., & McFarlane, N. A. (1995). Sensitivity of climate simulations to the parameterization of cumulus convection in the canadian climate centre general circulation model. Atmosphere-ocean, 33 (3), 407–446. –25– manuscript submitted to Journal of Advances in Modeling Earth Systems (JAMES) Figure 12.Pearson correlations between perturbed physics parameters and regional inte- grated gradient attributions for selected high-attribution variable and region pairs. Rows show selected regions for TLWP, SWCF, FSNT, LWCF, and FLUT, including marine stratocumulus, subtropical cloud-transition, ITCZ, and storm-track regions. Columns show perturbed param- eters associated with cloud microphysics, CLUBB boundary layer turbulence and shallow con- vection, deep convection, cloud fraction, aerosol emissions, solubility, and aerosol microphysics. Panels show correlations for (a) TAU-ML versus K2000, (b) TAU-ML versus observations, and (c) K2000 versus observations. Colors indicate Pearson correlation coefficients between pa- rameter values and regionally averaged integrated gradients. Open circles indicate statistically significant correlations over 0.3. Positive correlations indicate that increasing parameter values would increase the attribution in the region, while negative correlations indicate that increasing parameter values would decrease the attribution in the region. –26–