The NASA-GISS ModelE2.1-CC2 ESM: development and evaluation

This paper describes the NASA Goddard Institute for Space Studies (GISS) ModelE2.1-Carbon Cycle version 2 Earth System Model (NASA GISS E2.1-CC2 ESM), assesses its skill against observations and the previous version of the same model (GISS‐E2.1‐G‐CC). NASA GISS E2.1-CC2 ESM includes the same physical climate model E2.1-G as its predecessor, updated ocean and land carbon cycle and longer equilibrium simulations. While the focus here is on the land and ocean carbon components and their interactions with the atmosphere and ice, we also describe in detail the physical coupled model for consistency and ease of reference. We detail parameterizations, tuning, and conservation diagnostics that are relevant to the global (land and ocean) carbon cycle. We also describe the pre-industrial control and the historical simulations with this model while future climate scenarios will be addressed in a companion paper. The model is an improvement to prior releases, has been better tuned and spun up (all reservoirs: the ocean, the ocean carbon but especially the soil carbon reservoir) and exhibits smaller drifts in all carbon components. However, some persistent biases remain, particularly with regards to low ocean productivity and land primary production. Low ocean productivity may be due to (i) low iron concentration in high-nutrient and low chlorophyll regions, and (ii) lack of photoadaptation, which would increase the $\mathrm{chl}/\mathrm{C}$ ratio in low-light environments, relieving light limitation beneath the surface of the ocean. Low land productivity may result from an insufficient gross primary productivity increase due to lack of prognostic LAI, or drought stress, or an underestimated CO₂ fertilization effect, or an overestimate of soil respiration due to biases in soil moisture, soil temperature, and/or how soil respiration responds to these environmental variables.

1 Introduction

Earth system models (ESMs) integrate the interactions of atmosphere, ocean, land and ice with the biosphere to estimate the state of regional and global climate under a wide variety of conditions (Heavens et al., 2013) and focus on comprehensiveness of Earth system interactions, feedbacks and interdependencies. ESMs include all the processes present in climate models, but they also simulate the interaction between the physical climate, the biosphere and the chemical constituents of the atmosphere and the ocean. They can also include the impact of human decision-making, e.g. changes in land use, in emissions mitigation etc. They are thus particularly useful for estimating remaining carbon budgets before anthropogenic warming exceeds certain thresholds, studying the feasibility, scalability and impacts of climate intervention practices and geoengineering activities, such as carbon dioxide removal (CDR) through soil carbon sequestration, biomass carbon removal and storage, enhanced mineralization, ocean alkalinity enhancement, afforestation/reforestation etc and studies of the Earth's climate evolution after curbing emissions.

The NASA Goddard Institute for Space Studies (NASA-GISS) ModelE version 2.1 with Carbon Cycle version 2 (ModelE2.1-CC2) is an Earth system model with a fully interactive land and ocean carbon cycle which is an update of the GISS submission to the 6th generation Climate Model Intercomparison Project (CMIP6) which was identified as GISS-E2-1-CC in the CMIP6 repository. We will use the shorthand GISS-ESM CMIP6+ alternatively with GISS-E2-1-CC2 to denote the difference from the previous model.

In the following sections we describe as comprehensively as possible the physical and biogeochemical components of the model, the relevant parts of the radiation, atmospheric chemistry and clouds. We then present results from the spinups and historical simulations and assess them against observations. Lastly, we lay out the model response to future emission and mitigation pathways. Lastly, we offer a summary of the model, comparison with other models across different experiments in the post-CMIP6 assessments and lay out our future modeling plans.

2 Model Description

The underlying physical model, ModelE2.1, has been described in Kelley et al. (2020) and includes updates to earlier versions of the NASA GISS model (Schmidt et al., 2006; Hansen et al., 2007; Schmidt et al., 2014). For ease of access and comprehensiveness of model presentation, we will describe key components of the model here as well, although complete information and evaluations of the physical model should be sought in these publications.

2.1 Model configuration

The atmospheric model has horizontal resolution of 2 by 2.5° in latitude and longitude respectively and 40 layers in the vertical from the surface to 0.1 hPa in the lower mesosphere; it is coupled to an ocean with 1 by 1.25 horizontal resolution in latitude and longitude respectively and has also 40 layers in the vertical.

Atmospheric dynamics uses Arakawa's B-grid scheme where eastward and northward velocity components are located at primary cell corners; vertical discretization follows terrain‐following (sigma) layers in the troposphere which transition to constant‐pressure layers in the stratosphere. The transition is smooth and centered at 100 hPa with a half‐width of approximately 30 hPa. The vertical discretization in the ocean is 40 vertical layers in mass with higher resolution of 5–20 m near the surface and within the thermocline and about 200 m at depth.

Prognostic variables in the atmosphere are dry air mass, potential temperature, water vapor mixing ratio, and horizontal velocity components. Virtual potential temperature is used for all density/buoyancy-related calculations. Ocean prognostic variables include seawater mass, potential enthalpy that accounts for variations in specific heat capacity, salt and horizontal velocity components on both the C-grid (perpendicular to primary cell edges) and D-grid (parallel to primary cell edges).

Subsurface reservoirs are split into four types: open water (including lakes and oceans), ice-covered water (including lake ice and sea ice), ground (including bare soil and vegetated regions), and land ice/glaciers. Within each type there are further subdivisions (fraction of plant functional types, fractional snow cover, melt pond fraction over sea ice, etc.), which are communicated to the atmospheric model as weighted mean quantities (for example the albedo in the radiation computations). The non-ocean reservoirs include prognostic variables for water mass and heat content, and the sea ice also includes salt content.

The model uses a 30 min time step for physics calculations (i.e. surface interactions, condensation etc). The radiation code is called every five physics time steps (every 2.5 h), however, the zenith angle and sun light intensity are updated every physics time step. The model uses a 225 s time step for the atmospheric dynamics (transports) and 112.5 s time step for ocean gravity wave computations.

Simplifications in atmospheric thermodynamics ignore the sensible heat of condensate. Although latent heat of water vapor is treated correctly, its mass is not added to atmospheric dry mass. Consequently, precipitation and evaporation occur at 0 °C, specific potential enthalpy equals potential temperature times dry specific heat capacity, and humidity affects atmospheric density only in the pressure gradient force computation.

2.2 The ocean component

The ocean model is a dynamic ocean model E2.1-G (Russell et al., 1995; Schmidt et al., 2006, 2014; Kelley et al., 2020), non-Boussinesq (i.e. includes density variations in the horizontal as well as vertical momentum equations and in the continuity equation, following Schmidt et al., 2006). It conserves potential enthalpy, freshwater mass, and salt mass, and uses natural boundary conditions (for heat and freshwater). It has a free surface and thus it is not volume conserving. The number of vertical layers varies from 2 to 40. After the dynamics, the mass ratio of layer n and its layer above is fixed in time and space, but this does vary with n. The thickness of layer 1 varies and is approximately 10 m. The vertical coordinate is mass per unit area, including the mass of air, ice and water above any level.

The ocean model solves the dynamical equations via an invariant form. The Coriolis term is computed via the C and D-grid velocities and has latitude dependence. The quadratic upstream advection scheme (Prather, 1986) uses mean, three-dimensional gradients, and second order moments to advect scalar tracers, thus effectively increasing the resolution, in both horizontal and vertical direction. Bottom friction is applied at the lowest layer of each ocean column. The ocean model employs the K-profile parameterization (KPP) for vertical mixing (Large et al., 1994) and the Gent-McWilliams parameterization for eddy induced tracer transports (Griffies, 1998a, b). The mesoscale diffusivity allows for surface‐intensified eddies, by employing an exponential decay of diffusivity with depth that depends locally on the horizontal gradient of potential density (Kelley et al., 2020) which permits eddies to restratify the Southern Ocean over a large depth range, consistent with the observed density structure there, while not overacting in other regions of the World Ocean (such as the North Atlantic). The Rossby radius is replaced by a geographically constant nominal length scale, and the diffusivity is a function of the Brunt‐Vaisala frequency and the slope of isopycnal surfaces.

Mesoscale transport is expressed in local quasi‐isopycnal layers, avoiding difficulties associated with skew-flux representation (1) in regions with weak vertical stratification, e.g. the surface mixed layer and at the ocean surface, where the GM parameterization assumes that eddies transport tracers along isopycnal surfaces. However, near the surface, these surfaces can be highly tilted or intersect the ocean surface, which violates the assumptions of the parameterization. Applying the GM scheme directly in these regions can lead to spurious mixing and incorrect representation of processes like deep convection and mixed-layer re-stratification, (2) near lateral boundaries and the ocean bottom, where the eddy fluxes must be parallel to these boundaries, but the GM scheme's isopycnal-following nature doesn't naturally enforce this, and (3) in areas of strong density fronts (e.g., the Gulf Stream or the Antarctic Circumpolar Current) where the “small-slope approximation” which assumes that the slope of isopycnal surfaces is small and is valid in the deep ocean, breaks down in areas where the isopycnals are steeply tilted. Ventilation of marginal seas through their connecting straits has been increased via two mechanisms in E2.1‐G, reducing salinity biases there. While there is no explicit overflow parameterization in the ocean model, for straits deep enough that density contrasts can drive strong opposing flows at the surface and at depth, the finer upper‐ocean layering in E2.1‐G resolves this structure, in conjunction with a slight tuning of strait depths. Second, horizontal diffusivity was increased in straits that are shallow or have weaker density contrasts. The first mechanism impacted the Red and Black seas, and the second the Baltic and Hudson (Kelley et al., 2020). The first is the sole ventilation mechanism for straits narrower than the nominal resolution, which are parameterized using the Russell et al. (1995) one‐dimensional channel scheme that lacks horizontal mixing.

In addition to KPP, the vertical diffusivity also includes a contribution from tidal dissipation. The ocean model Atlantic meridional overturning circulation (AMOC) is found to be sensitive to tidally driven mixing, which occurs in the shallow waters bordering the North Atlantic following dissipation distributions generated by Jayne (2009).

The Equation of State employed in the ocean model computes specific volume (m³ kg⁻¹) as a function of specific potential enthalpy (J kg⁻¹), salinity (psu), and pressure (Pa). Using constant uniform gravity, pressure at any level is gravity times mass per unit area above the level. The one-to-one relationship between temperature and specific potential enthalpy is computed using the adiabatic lapse rate of sea water as function of temperature, salinity and pressure (IES-1980 et al., 1981) and specific enthalpy of sea water as function of temperature and salinity at atmospheric pressure. Using these formulas, the model's Equation of State can be computed from the specific volume function of temperature, salinity and pressure (IES-1980 et al., 1981). Complex computations for Equation of State, freezing temperature, specific heat capacity, and others are evaluated offline and are stored in tables; the model interpolates the table values as functions of specific potential enthalpy, salinity and pressure (Chen and Millero, 1976). Table resolutions are 4000 J kg⁻¹ of specific potential enthalpy, 1 psu of salinity, and 2×10⁶ Pa of pressure. Under adiabatic flow, potential enthalpy is a conserved quantity, therefore the ocean model conserves energy (heat) as well as mass and salt.

2.2.1 Ocean Carbon Cycle

The ocean carbon cycle model is based on the NASA Ocean Biogeochemistry Model (NOBM), originally developed at the NASA-Global Modeling and Assimilation Office (Gregg and Casey 2007) and later coupled coupled to the GISS climate model under prescribed atmospheric concentrations (Romanou et al., 2013, 2014) and under prescribed emissions (Ito et al., 2020). The model includes four phytoplankton species (diatoms, chlorophytes, cyanobacteria, and coccolithophores), four nutrient species (nitrate, silicate, ammonia, and iron), three detrital pools (nitrate/carbon, silicate, and iron) and one heterotroph species (Fig. 1). Carbon cycling is represented through dissolved organic (DOC) and inorganic carbon (DIC), where carbonate chemistry is solved using MOCSY (Orr and Epitalon, 2015). Ocean carbon interacts with atmospheric CO₂ through gas exchange parameterization, following the CMIP6 protocol (Orr et al., 2017). The current version of NOBM (NOBMg) has been expanded to include oxygen and alkalinity (Lerner et al., 2021, 2024) and an updated treatment of iron scavenging (see below). The model equations and parameters used are laid out in Appendix A and in Tables A1–A11 therein. A short description of the model equations is provided below.

Figure 1 Schematic diagram of the ocean carbon cycle model (NOBMg). credit: Jay Friedlander, NASA-GSFC.

All tracers in NOBMg are changed via physical transport and mixing which are following parameterizations of these processes in the physical ocean model.

Phytoplankton species are also subject to sinking, growth via nutrient uptake, respiration and exudation of DOC, senescence and zooplankton grazing. Nutrient distributions are changed via uptake by phytoplankton species for growth, remineralization of particulate organic carbon or nitrogen, and detrital breakdown. Cyanobacteria can also grow via nitrogen fixation, where growth is limited only by light, iron, and temperature. Nitrogen assimilated into organic matter is obtained from N₂. Cyanobacterial growth via nitrogen fixation follows the same Michaelis-Menten formulation as cyanobacterial growth that is nitrate limited, with two exceptions. First, the iron half-saturation constant for nitrogen fixation is 33 % higher than that for nitrate-limited growth, as nitrogen fixation generally relies upon an increased availability of iron. Second, nitrogen fixation is inhibited by the presence of nitrate and ammonium, following laboratory culture experiments using Trichodesmium (Holl and Montoya, 2005). Ammonium is produced via the herbivore grazing and heterotrophic loss. Silicate is taken up only by diatoms. Iron is produced by zooplankton, atmospheric deposition of dust, and is lost to scavenging. NOBMg also includes an updated marine iron cycle, with regionally varying prescribed ligands, and a significant increase in the rate of ligand-bound iron scavenging, to be consistent with observed rates for trace metals with similar scavenging characteristics as iron (see Sect. 3.1.1 on Tunings). At each time step, the prognostic iron pool is split into iron that is complexed by organic ligands and a “free” iron pool that represents the sum of all inorganic iron species. Iron complexation chemistry is computed at each time step using vertical and regionally varying prescribed ligand concentrations (see Sect. 3.1.2 Forcings) and a conditional stability constant (K_FeL) = 10¹¹ L mol⁻¹. Then, free and ligand-complexed iron are computed by solving the system of equations for total iron, total ligand concentration, and the conditional stability constant (Aumont et al., 2015).

Herbivores grow via grazing of phytoplankton and are lost to heterotrophic consumption and excretion to DOC.

Carbon and nitrate particulates (detritus) are linked via the constant C:N ratio. Therefore, only one equation is used to describe both. Their change is due to advection and diffusion, sinking, remineralization and detrital breakdown. Sinking of particulates is expressed via an exponential function of depth (see Eqs. A15–A17). Silicate detritus includes an additional source term due to senescence of phytoplankton and grazing from herbivores. The iron detritus includes zooplankton production and scavenging production. DOC is produced during phytoplankton growth, herbivore excretion and detrital breakdown while the only DOC sink is bacterial degradation. DIC sources are phytoplankton and herbivore respiration, bacterial degradation of DOC, remineralization of carbon detritus and air sea exchange while the only sink is the growth of phytoplankton.

The model includes prognostic alkalinity which is advected and mixed like all other tracers and follows the Ocean Carbon Model Intercomparison Project parameterizations (OCMIP-2, Najjar et al., 2007), i.e. it is produced via conversion of nitrate to organic N and carbonate dissolution, while it is lost to remineralization of detrital organic carbon $/$ nitrogen and DOC and carbonate precipitation (see Eqs. A43–A47).

Oxygen interacts with the rest of the ocean biogeochemical tracers via remineralization of organic carbon and by determining when the model switches to slower rates of remineralization that are consistent with denitrification (Eqs. A48–A56; Lerner et al., 2024).

2.3 The Land component

2.3.1 Vegetation

The Ent Terrestrial Biosphere Model (Ent TBM) (https://www.giss.nasa.gov/projects/ent/, last access: 21 August 2026) is a demographic dynamic global vegetation model (DGVM) (Kiang, et al., 2006; Kim et al., 2015) with subgrid vegetation patches distinguished by plant functional type (PFT), currently with single-cohort structure for patch communities (Fig. 2). The 12 plant functional types are: evergreen broadleaf trees, evergreen needleleaf trees, cold deciduous broadleaf trees, drought deciduous broadleaf trees, deciduous needleleaf trees, cold-adapted shrubs, arid-adapted shrubs, C3 perennial grass, C3 annual grass, C3 arctic grass, C4 grass, and herb crops (with physics of C3 perennial grass). Soil biogeochemistry is a version of the Carnegie-Stanford-Ames' (CASA') model (Potter et al., 1996; Randerson et al., 2009; Doney et al., 2006), modified for Ent for an improved response of soil respiration to soil moisture, as described in Kim et al. (2015). Land ecosystem carbon pools that are calculated are: foliage, hardwood stem, sapwood stem, fine roots, coarse roots, non-structural carbohydrate (which we also call labile carbon), and 9 soil carbon pools.

Leaf biophysics is a coupled conductance/photosynthesis model that simulates vegetation stomatal control of both transpiration and uptake of CO₂. The coupled leaf photosynthesis/ stomatal conductance model consists of the photosynthesis model of Farquhar and von Caemmerer (1982) and the Ball-Berry model of stomatal conductance (Ball and Berry, 1987), with leaf boundary layer conductances following the approach by Collatz et al. (1991) but with the boundary layer conductance derived from the canopy surface layer conductance of the land surface model; for the cubic equations of coupled photosynthesis/conductance, a numerical solver for net assimilation of CO₂ and leaf internal CO₂ concentration was developed. C4 photosynthesis follows the scheme of Collatz and Berry (1992). Leaf-to-canopy scaling is achieved through the canopy radiative transfer scheme of Spitters et al. (1986) (introduced to ModelE in CMIP5 in Schmidt et al., 2006) combined with a numerical scheme for vertical layering of the canopy by converging the total canopy fluxes. Vegetation canopy absorbance and transmittance of photosynthetically active radiation thus depend on vertical layering of leaf area index, solar zenith angle, and the fractions of diffuse vs. direct radiation, providing vertical profiles of sunlit and shaded leaf area fractions. Physiological algorithms for water stress, cold hardening in evergreen needleleaf trees, and temperature sensitivity of photosynthetic capacity are detailed in Kim et al. (2015). To account for variation in parameters of the leaf biophysics by climate and seasonality of vegetation, temperature dependence of leaf photosynthetic capacity, V_cmax, is modeled with a Q₁₀ function (Hegarty, 1973), with a reference temperature of 25 °C and parameters from Collatz et al. (1991). In addition, photosynthetic capacity is subject to seasonal cold hardening in evergreen needleleaf trees, as a function of local air temperature trend, and is modeled based on the work of Repo et al. (1990); Hanninen and Kramer (2007); and Makela et al. (2004, 2006). Autotrophic maintenance respiration is subject to temperature acclimation, in which the reference optimal temperature in the equation of Levis et al. (2004) is allowed to vary between 10–30 °C, tracking a 10 d running average of surface air temperature; 10 d averages less than 10 °C or greater than 30 °C are considered outside the range of acclimation.

Water stress is modeled after the approach of Rodriguez-Iturbe et al. (2001), in which stress is a factor between 0 (full stress) and 1 (unstressed), as a function of relative extractable water (volume water − minus hygroscopic water)$/$(pore volume minus hygroscopic water). The factor modifies the leaf stomatal conductance, and different plant functional types have different critical soil moisture values for the onset of stress and wilting points. Because transpiration depends on the exposed leaf area index, and wet leaf fractions will not transpire, a scheme for canopy interception of precipitation was added following the algorithm presented in Koster and Suarez (1996) which is like the modified Shuttleworth scheme of Wang et al. (2007). The scheme primarily involves modification of the wet fraction of the canopy by computing the drip rate considering the fraction of precipitation that falls onto a previously wetted portion of the canopy with the grid cell. This scheme also distinguishes between stratiform and convective precipitation, because the stratiform type covers a greater portion of the grid cell than the convective type (e.g., Koster and Suarez, 1996; Wang et al., 2007).

For surface albedo, Ent's 12 different vegetation types have different spectral and snow masking depth properties. In E2.1 Ent PFTs are mapped to the 8 biome types of Matthews (1984) for seasonal albedo for radiative purposes.

In Schmidt et al. (2006) and henceforth, all surface conductance of water vapor for vegetated land was assigned to the vegetation, which meant that there was no soil evaporation under canopies. The rationale for this approach (Abramopoulos et al., 1988) was that vegetated land was simply a surface with a particular albedo and conductance. Numerous studies have demonstrated the importance of soil evaporation as a component of evapotranspiration, especially for semiarid ecosystems (e.g., Lawrence et al., 2007). While evaporation from bare soil was included in previous versions of the model, in E2.1 we have added a simple representation of soil evaporation based on the work of Zeng et al. (2005) and Lawrence et al. (2007), by including the evaporative flux from the vegetated soil into the total evapotranspiration scaled as a function of the current leaf area index. Total evapotranspiration never exceeds the potential evaporation.

2.3.2 Irrigation and Groundwater

In E2.1, irrigation is implemented as a prognostic field (Cook et al., 2011, 2014; Krakauer et al., 2016; Puma and Cook, 2010; Shukla et al., 2014) where water demand for irrigation (Wada et al., 2014) is estimated using prescribed irrigation areal extent (Siebert et al., 2015). The water is drawn first from rivers and lakes locally, and if that is insufficient, it is drawn from an external groundwater source (which is tracked diagnostically). Groundwater is assumed to have the same temperature as the soil and default tracer concentrations. Groundwater recharge is not accounted for, and so there is a small increase in total water mass (and eventually, sea level) associated with the net global groundwater draw in these simulations. These effects have a complex impact on freshwater delivery to the oceans (and hence sea level). Irrigation from local surface water sources leads to increased soil moisture and reduced river outflow, but this is dominated by net additions of groundwater which add freshwater to the climate system, about 0.2 mm yr⁻¹ of global sea-level-equivalent in 2010 (Miller et al., 2021). As discussed in Miller et al. (2021) this freshwater imbalance is driven by an underestimate of surface water in lakes and rivers (which means that more groundwater is needed) and also the lack of counteracting increase in reservoirs.

2.3.3 The Land Carbon Coupling

The land carbon dynamics in the fully coupled carbon cycle physics configuration is the same as described in Ito et al. (2020) for the ModelE submission to CMIP6 C4MIP. The leaf biophysics is sensitive to the surface atmospheric CO₂ concentration. The model tracks CO₂ exchanges between the land and atmosphere in diagnostics for GPP, autotrophic respiration, soil respiration, and carbon fluxes due to land cover change. Land carbon change associated with the annual land cover change is distributed as CO₂ fluxes to the atmosphere over the following year to avoid having to add a pulse of CO₂.

Figure 2 Schematic diagram of the Ent in GISS ModelE2.1-CC2.

2.4 The sea ice and land ice models

Sea ice and lake ice are computed on the atmospheric grid, and the information is mapped on to the ocean grid. Temperature and salinity profiles in the ice are computed prognostically, and sea ice distribution trends and seasonal cycles are realistic (Kelley et al., 2020).

Salinity in the ice affects basal ice formation/melt and heat diffusion (Schmidt et al., 2006) using brine-pocket (BP) thermodynamics (Bitz and Lipscomb, 1999; Schmidt et al., 2014) and updated in Kelley et al. (2020). Sea ice dynamics use the viscous-plastic formulation of Zhang and Rothrock (2000). Lead fraction is treated differently in the Arctic and Antarctic wintertime thick-ice growth to produce more realistic heat fluxes and sea ice growth rates. Processes relevant to the salt budget (e.g., gravity drainage and flushing of meltwater) are consistently treated with the BP physics.

The mass and energy associated with net snow accumulation over each hemisphere's ice sheets is used to adjust the implicit iceberg calving fluxes into the adjacent oceans over a 10-year relaxation time, and this diagnostic is preserved as glacial melt. This ensures that the model can reach a long-term equilibrium under changed climate forcings but does not imply realism in the response of ice sheets to climate, since the ice sheets are assumed to remain in mass balance.

2.5 Radiation, Clouds and Atmospheric Composition

The total solar irradiance at the top of the model's atmosphere is based on satellite calibrations (Kopp and Lean, 2011) with a base value of 1361 W m⁻². Spectral irradiance values are taken from Coddington et al. (2016). The radiation model includes explicit multiple scattering calculations for solar radiation (shortwave, SW) and explicit integrations over both the SW and thermal (longwave, LW) spectral regions (Schmidt et al., 2006; Schmidt et al., 2014; Kelley et al., 2020). SW radiation is absorbed by H₂O, CO₂, O₃, O₂, and NO₂ using the k-distribution approach (Lacis and Oinas, 1991) that utilizes 15 noncontiguous spectral intervals to model overlapping cloud–aerosol and gaseous absorption and is scattered by clouds and aerosols using size-dependent Mie scattering, ray tracing and T-matrix theory (Mishchenko et al., 1996).

Longwave radiation is absorbed by H₂O, CO₂, and O₃ using the correlated k distribution with 33 intervals (Lacis and Oinas 1991; Oinas et al., 2001), designed to match line-by-line computed fluxes and cooling rates throughout the atmosphere to within about 1 %. Weaker bands of H₂O, CO₂, and O₃, as well as absorption by CH₄, N₂O, CFC-11, and CFC-12 are included in an approximate fashion as overlapping absorbers, but with coefficients tuned to reproduce line-by-line radiative forcing over a broad range of absorber amounts.

Variations in incoming solar forcing in the 20th Century simulation follow Matthes et al. (2017) and include updated spectral solar irradiance. Orbital parameters that determine the seasonal and latitudinal variation in insolation are calculated each year (updated on 1 January), based on Berger (1978).

The radiation model also includes the effects of 3D cloud heterogeneity via the 3D cloud parameterization (Cairns et al., 2000) to get more realistic albedos from realistic water paths and particle sizes (Schmidt et al., 2006). The procedure rescales optical depth, asymmetry parameter, and single scattering albedo according to the relative variance of the cloud particle density distribution. The partitioning between convective precipitation that descends and has the potential to evaporate in the environment rather than in the downdraft is increased from 0 % to 50 %, thus increasing the sensitivity of humidity to convection. The downdraft buoyancy is based on virtual temperature including condensate loading. Glaciation in the mixed‐phase temperature range has been updated in the stratiform cloud parameterization.

Stratiform cloud formation is done more realistically because there is no restriction in coverage, no dependence on vertical motion, the altitude is taken to be relative to local planetary boundary layer (PBL) height rather than a fixed 850 hPa, better demarcating cloud‐topped boundary layers from the free troposphere where the threshold relative humidity is used for the TOA radiation balancing process (described in Sect. 3.1.1) to maintain top‐of‐the‐atmosphere radiative balance.

The advection of humidity (and other tracers) is done only once every physics time step (30 min) with iterative time stepping to avoid any Courant–Fredrichs–Levy violations. The smaller time steps are used primarily in the stratosphere and upper troposphere where zonal winds are strong or as a result of extremely high flow deformation.

The atmospheric chemistry configuration used in GISS ModelE2.1-CC2 is the basic non-interactive composition (NINT f2 used in CMIP6) in which the background fields of ozone and aerosol concentrations (dust, sea-salt, sulfate, nitrate, ammonium and carbonaceous aerosols e.g. black and organic carbon) are derived from simulations of the interactive one-moment aerosols (OMA) version of the model, run as an atmosphere only model with prescribed SST (Kelley et al., 2020). All anthropogenic and biomass burning emissions of short‐lived species were updated to CMIP6 specifications (Hoesly et al., 2018; van Marle et al., 2017) and were prescribed annually or monthly.

Specifically, with regards to dust, the six dust bins (0.1–0.2, 0.2–0.5, 0.5–1, 1–2, 2–4, and 4–8 µm particle diameter), are used for coatings on dust particles, the dust particle densities, and the weights that are used to partition the total clay which is advected as a bulk species in the model. The weights reflect the size distribution of dust. The effect of wind biases on dust transport over the dust source regions is alleviated in the model by employing an additional term (“gustiness”) to express the high frequency variability in the surface winds that lift dust into the atmosphere. This variability results from stirring of the atmospheric boundary layer by dry convective thermals or frictional eddies, along with downdrafts spreading outward from the base of convective events. The magnitude of the associated wind gusts is parameterized as described by Cakmur et al. (2004). The longwave optical depth for dust was increased by 30 % to account for the longwave scattering effect. The lensing effect of sulfate and nitrate coatings on BC was parameterized by increasing the shortwave optical depth for BC by 50 %. For more information see Pérez García-Pando et al. (2016).

Currently in NOBMg within GISS ModelE2.1-CC2, dust iron deposition is obtained from the historical annual cycle climatology extracted from the GISS E2.1-G online dust simulations which include eight externally mixed minerals (illite, kaolinite, smectite, carbonates, quartz, feldspar, iron oxides, and gypsum) plus internal mixtures between seven minerals and iron oxides (Perlwitz et al., 2015a; Perlwitz et al., 2015b). The masses of free and structural iron and their fractions of total iron have been evaluated using measurements from Izaña Observatory (Pérez García-Pando et al., 2016). For more on the dust deposition, see Sect. 3.1.2 on Forcings. In future versions of the model, dust deposition will be interactive with the global carbon cycle.

A second order scheme is used for the transport (Schmidt et al., 2006) of all the tracer fields in the atmosphere which are advected using the highly non-diffusive quadratic upstream scheme (QUS; Prather, 1986), which keeps track of nine subgrid-scale moments as well as the mean within each grid box. This increases the effective resolution of the tracer fields and allows the GISS model to produce reasonable climate fields with relatively coarse nominal resolution.

2.6 Surface albedo

The radiation model generates spectrally dependent direct/diffuse flux ratios for use in biosphere feedback interactions. The surface albedo utilizes six spectral intervals, in nanometers: nm: VIS (300–770), NIR1 (770–860), NIR2 (860–1250), NIR3 (1250–1500), NIR4 (1500–2200), NIR5 (2200–4000), and is solar zenith angle dependent for ocean, snow, and ice surfaces.

The spectral and solar zenith angle dependence of ocean albedo is based on calculations of Fresnel reflection from wave surface distributions as a function of wind velocity (Cox and Munk, 1956) taking also into account the effects of foam and hydrosols on ocean albedo (Gordon and Jacobs, 1977). The spectral and solar zenith angle dependence of snow and sea ice is modeled following Warren and Wiscombe (1980). Snow “ages” following Loth and Graf (1998) and has a different albedo for wet or dry snow. Ocean ice albedo is spectrally dependent and is a function of ice thickness and parameterized melt pond extent (Schmidt et al., 2006 and references therein).

Light profiles from the atmospheric radiation module are propagated underwater into the ocean and spectrally decomposed to 33 wavebands that are used to compute growth of the phytoplankton groups, sinking profiles, as well as changes in the vertical distribution of ocean temperatures due to biologically mediated absorption and scattering in the water column (Gregg and Conkright, 2002). NOBMg is coupled prognostically to the atmospheric radiation module which is part of ModelE2.1. Direct and diffuse sunlight in 6 wavebands (5 in the near IR and 1 in the visible part of the spectrum) that are computed by the radiation module are decomposed into the 33 wavelengths (250–3700 nm) needed for the absorption and scattering by the different phytoplankton species. The decomposition factors are based on the fractions of direct (diffuse) sunlight at the equator for each wavelength relative to the total radiation as derived from a monthly climatology provided by the Ocean–Atmosphere Spectral Irradiance Model (OASIM) and forced by MODIS data for aerosols and clouds (Gregg and Casey, 2007).

Vegetation spectral albedo is prescribed as seasonally varying by plant functional type (PFT) and is treated as Lambertian, assuming a flat, diffusely reflecting surface. Bare soil is represented as grey, with values from Matthews (1984). Seasonal spectral albedos for vegetation are specified by mapping Ent PFTs to the eight Matthews biome types: tundra (Ent cold-adapted shrubs), C3 grass (Ent C3 perennial, annual, and arctic grasses), shrub (Ent arid-adapted shrubs), deciduous forest (all Ent deciduous trees), evergreen needleleaf forest (Ent evergreen needleleaf trees), tropical rainforest (Ent evergreen broadleaf trees), and crops (Ent crops). Total land albedo is calculated as a weighted average across the corresponding land-cover fractions and linearly interpolated between prescribed seasonal values. When snow is present, snow albedo is included in the weighted average according to its masking fraction.

2.7 Interface exchanges

2.7.1 Surface fluxes

Separate atmospheric surface fluxes are computed for each subsurface reservoir: open water, sea or lake ice, land ice, snow free and snow-covered ground. ModelE calculates atmospheric turbulence over the entire vertical column. The surface fluxes are calculated using a submodule that is embedded between the surface and the midpoint of the first resolved model layer. In the atmospheric planetary boundary layer (PBL), a formulation for the temperature, moisture, and scalar fluxes is employed and consists of a local (diffusive) term and a countergradient term derived from large-eddy simulation (LES) data (Holtslag and Moeng, 1991) as well as a nonlocal vertical transport scheme for virtual potential temperature, specific humidity, and other scalars (Holtslag and Boville, 1993).

A level-2.5 turbulence model (Cheng et al., 2002; Yao and Cheng, 2012) is applied to this region (using eight sublayers) independently over each surface type where the mixing profile depends on the height of the PBL, which is closely related to the large eddy size and thus characterizes the nonlocality, and the buoyancy and shear effects at the surface. The turbulence length scale formulation is obtained from the large eddy simulation data by Nakanishi (2001).

Above the PBL a second-order closure (SOC) model developed by Cheng et al. (2002) is used, which is based on Mellor and Yamada (1982). The Reynolds stress and heat flux equations have been solved with a more advanced parameterization of the pressure–velocity and pressure–temperature correlations. In addition, a few turbulent time scales are determined by a new two-point turbulence closure model (Canuto and Dubovikov, 1996a, b). The more realistic “Richardson number criterion” rather than the “TKE criterion” to calculate the PBL height, following Troen and Mahrt (1986) and Holtslag and Boville (1993).

The lower boundary conditions for the wind and the potential temperature equations are determined assuming continuity of the respective turbulent fluxes at the surface. To calculate the fluxes through the surface layer, drag and transfer coefficients are set using similarity theory following Hartke and Rind (1997). The roughness length for momentum (z0m) over land is specified as in Hansen et al. (1983). The roughness lengths for temperature (z0h) and moisture (z0q) over land and land ice are taken from Brutsaert (1982) to be proportional to z0m. Ocean and ocean ice are treated as rough bluff surface, with z0m combining the smooth surface value (Brutsaert 1982) with the Charnock relation for the aerodynamic roughness length; as for z0h and z0q Eqs. (5.26)–(5.27) of Brutsaert (1982) with background values $\mathrm{1.4}\times {\mathrm{10}}^{-\mathrm{5}}$ and $\mathrm{1.3}\times {\mathrm{10}}^{-\mathrm{4}}$ m, respectively, are used.

We reduce the saturation specific humidity by 2 % over the oceans (to account for sea salt aerosol effects; Gill, 1982).

Over land, the surface evaporative flux is determined by the vegetation, but to improve estimates of surface humidity by the submodule described above, we calculate the maximum available evapotranspiration and do not allow the atmosphere to draw more water than is available. This requires a change in the humidity surface boundary condition from a variable to a fixed flux in such situations.

Continental runoff occurs at the mouths of rivers as well as small tributaries. Prognostic river runoff is computed within ModelE2.1-CC as part of the simulated global hydrological cycle, thus increases (decreases) due to changes in precipitation and evaporation over land, and delivers accordingly nutrients and carbon components to the coastal ocean (for more details see Sect. 3.1.2 on Forcings).

2.8 Coupling: Ocean, lake, ice, and land surface coupling

Coupling between the different components in ModelE is synchronous and at the frequency of the physics time step (30 min). It is always implemented via fluxes of the fundamentally conserved quantities (mass, energy, and tracer). Given the surface conditions, the atmospheric model calculates the precipitation, radiative, and other surface fluxes (surface wind stress, evaporation, and sensible and latent heat) over each surface type. Radiation and precipitation fluxes are first applied directly to the land surface (soils, vegetated ground, and glaciers) and the land surface model computes latent and sensible heat fluxes, as well as land carbon fluxes. Using the atmosphere–ice wind stress, the sea ice dynamics calculates the horizontal ice velocities and the resulting ice–ocean stress. The ocean has finer horizontal resolution than the atmosphere and thus area conserving schemes of horizontal interpolation of fluxes from the atmospheric and the ocean grid are designed to conserve mass, energy, and tracer fluxes.

River flow occurs on the atmospheric grid and is deposited into the ocean by the coupler (Schmidt et al., 2014). Ice-ocean fluxes involving melt, heat conduction, and momentum drag are computed on the ice grid, while ice formation in the ocean occurs on the ocean grid and is then interpolated back to the atmosphere-ice grid ensuring conservation of fluxes.

2.9 Conservation

More than 100 vertically integrated diagnostics for every process in the atmosphere and the ocean are accumulated on each horizontal surface cell. They include fresh and saltwater mass, atmospheric kinetic and static energy, ocean energy, and water and energy for each continental reservoir, as well as all ocean tracer fields. Offline programs utilize these monthly accumulated diagnostics, as well as water mass, tracer and energy amounts for each reservoir at beginning and end of each month and compute the change in reservoirs and compare to flux changes. For both global and regional domains, errors are zero within the restrictions of machine accuracy (real*4 for fluxes and real*8 for reservoirs).

2.9.1 Physical model conservation

In the atmosphere, the advection scheme is mass conserving for humidity and tracers, and potential enthalpy conserving for heat. All processes including dynamics, cloud schemes, gravity wave drag, and turbulence, conserve energy and air, water, and tracer mass to machine accuracy. All dissipation of atmospheric kinetic energy through various mixing processes is converted to heat globally. In the long-term mean, net top-of-atmosphere (TOA) radiative flux approaches zero in the pre-industrial control simulation.

2.9.2 Land carbon conservation

In atmosphere-coupled experiments, the Ent TBM is currently run in biophysics-only mode without prognostic leaf area index (LAI) or stem growth. Vegetation structure is prescribed from satellite-derived vegetation cover, type, height, and maximum LAI. Maximum LAI is used to compute vegetation demography and cohort densities based on Ent's allometric relations by PFT. For each PFT, we use observed monthly LAI from Moderate Resolution Imaging Spectroradiometer (MODIS) which is interpolated at a daily time step. For land cover/land use change (LCLUC), Ent prescribes only annual historical change in cropland cover, rescaling the natural vegetation cover fractions in a grid cell (Miller et al., 2021; Ito et al., 2020); land use change such as deforestation is included only implicitly as part of the total cultivated land cover change.

Currently, we run Ent in biophysics-only mode where vegetation structure is prescribed, and all carbon pools are prescribed except the soil pools and the labile (non-structural carbohydrate) vegetation carbon pool which are prognostic. The labile carbon pool is dynamic with uptake of carbon through photosynthesis, withdrawal for respiration and leaf growth, and retranslocation with senescence. Any carbon excess, that is since LAI is prescribed, and no stem growth is allowed, is redirected to the litter pool in the soil, and that ensures carbon conservation. As mentioned earlier, land carbon change associated with the annual land cover change is distributed as CO₂ fluxes to the atmosphere over the course of the following year. Because the prescribed vegetation structure may not be in equilibrium with the simulated climate due to climate model biases or scenarios of climate change, respiration is restricted when labile carbon balances approach zero or slightly negative, allowing recovery of the labile pool during the growing season. This is also done to ensure carbon conservation.

2.9.3 Ocean carbon conservation

NOBMg conserves carbon, iron, ammonium and silicate but does not conserve nitrate, due to the fact that the nitrogen cycle is not closed. For nitrogen fixation, we assume an unlimited reservoir N₂ is available for diazotrophs. Moreover, while there is an oxygen-dependent switch from faster (aerobic) to slower (anaerobic) remineralization rates, denitrification (the conversion of nitrate to N₂) is not represented in NOBMg. Alkalinity is also not conserved since we do not include an explicit calcium carbonate tracer. Hence, rather than exchanges between dissolved and solid-phase carbonate ions being present, dissolved carbonate ions are either removed in near-surface waters by implicit calcium carbonate production (assumed to be proportional to primary productivity), or released by implicit carbonate dissolution in deep waters (assumed to follow a fixed exponential profile).

3 Model setup

3.1 Model coupled-carbon components

The general procedure for setting up and running the GISS ModelE2.1-CC2 consists of the following steps:

1. running the physical climate model coupled to the one-moment aerosols (OMA) composition submodel in atmosphere with prescribed SST to obtain distributions of aerosols for the non-interactive aerosol (NINT) runs

2. running the fully coupled climate model with aerosols prescribed from (1) to equilibrium, i.e. till the top of the atmosphere radiative forcing is near 0 (approximately 3000 model years)

3. running the fully coupled climate model after the completion of step (2) + the land carbon to equilibrium i.e. the air-land flux of CO₂ is near 0 (within ±0.1 PgC yr⁻¹ from 0, approximately 50 years for radiative equilibrium only), this involves a preliminary soil carbon spinup as described in Sect. 3.1.3.

4. running the fully coupled model after the completion of step (2) + the ocean carbon cycle to equilibrium, i.e. the air-sea flux is near 0 (within ±0.1 PgC yr⁻¹ from 0, about 500 years).

5. running the fully coupled climate model (3)+(4) (including soil carbon) till carbon equilibrium (within ±0.1 PgC yr⁻¹ from 0, spinup run)

Tuning is done separately in each of steps 2–4 by controlling parameters within each submodel to ensure that CO₂ flux to the atmosphere is near zero. The tuning procedure for the ESM model is described below.

3.1.1 Model tuning

Tuning of the atmosphere-ocean coupled model

The value of the critical relative humidity for clouds is used for the TOA radiation balancing process described to maintain top‐of‐the‐atmosphere radiative balance (Schmidt et al., 2014).

Tuning of the land carbon

The Ent TBM photosynthetic capacity, slope of the Ball-Berry stomatal conductance equation, and respiration parameters were tuned offline on Fluxnet sites to predict observed net ecosystem exchange (NEE) and attain equilibrium behavior of seasonal labile carbon pools, as documented in Kim et al. (2015). The Fluxnet sites were: Tapajo National Forest, Brazil (evergreen broadleaf trees); Hyytiala, Finland (evergreen needleleaf trees); Morgan Monroe State Forest, Indiana, and Harvard Forest, Massachusetts, USA (cold deciduous broadleaf trees); Tonzi Ranch, California, USA (drought deciduous broadleaf trees); Vaira Ranch, Ione, California, USA (C3 annual grassland). Other vegetation parameters tuned to summaries from the literature were: relative soil moisture at water stress onset and wilting point; temperature at which frost hardening commences and begins; linear response of soil respiration to relative soil moisture and unstressed level. All allometric relations were also derived from literature sources. For applicability at the global scale, specific leaf area (SLA) by PFT was derived from the TRY database (https://www.try-db.org/TryWeb/Database.php, last access: 21 August 2026) of leaf traits (Kattge et al., 2011).

The SLA was further tuned for each PFT in the coupled model spinup using values within the tuning parameter range in the TRY database of leaf traits to ensure that the plant labile carbon pool (non-structural carbohydrate) remains stable at the given size of mature trees. This stable behavior of the labile carbon pool depends on prescribed plant size boundary conditions being at equilibrium with the simulated climatology. Carbon starvation, therefore, may occur locally in ESM simulations, in which case the plant respiration is artificially reduced or stopped to prevent negative carbon balances; such plants would otherwise experience mortality. Excess labile carbon can be allocated to prognostic stem growth, or when the model is run in biopysics-only mode with fixed canopy structure it is simply senesced as extra litterfall to the soil.

Tuning of the ocean carbon cycle model

Tuning is performed to ensure that the preindustrial air-sea CO₂ is at equilibrium and secondarily to produce realistic patterns and magnitudes of ocean primary production. Firstly, phytoplankton maximum growth rates for chlorophytes and coccolithophores were reduced to be within the range reported in previous studies. For chlorophytes, the growth rate (μ_chlo) was reduced from 2.52 to 1.89 d⁻¹ (Laws et al., 2011; Grant et al., 2013; Liefer et al., 2018). For coccolithophores, the growth rate (μ_cocc) was reduced from 2.27 to 1.70 d⁻¹ (Buitenhuis et al., 2008; Krumhardt et al., 2017) (Table A5). From this change, several tuning experiments were conducted where the maximum carbon remineralization rate (α_C) was perturbed between 0.06 and 0.6 d⁻¹. A maximum remineralization rate of 0.3 d⁻¹ resulted in the air-sea CO₂ flux being at equilibrium, as well as net primary productivity increasing from 15 to 23 Pg C yr⁻¹ (Lerner et al., 2024). To tune primary production, we first increased the stoichiometric ratio of iron to chlorophyll (b_Fe) from $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{5}}$ to $\mathrm{8}\times {\mathrm{10}}^{-\mathrm{5}}$ µg Fe µg chl⁻¹, the latter value falls within the upper fourth quartile of the compilation of stoichiometries for different phytoplankton groups (Finkel et al., 2010). We also increased the silica half-saturation constant (k_Si) for diatom growth from 0.2 to 2.85 µM, the revised value being a compromise between the much lower values found for low-latitude species (0.02–1.4 µM), and the much higher values found for species in the Southern Ocean (4.2–88.7 µM; Nelson and Tréguer, 1992). Finally, we tuned the remineralization rate of detrital iron and the scavenging rates of ligand-bound to increase NPP while achieving an equilibrated air-sea CO₂ flux. For detrital iron, we explored a range of maximum detrital iron remineralization rates from 0.01 to 0.4 d⁻¹, while for the ligand-bound iron scavenging rate (k_Fe), the range explored was 2.75 × 10⁻⁴ to 2.75 × 10⁻² d⁻¹. The best combination of parameters that both increased NPP while allowing the air-sea CO₂ flux to achieve equilibrium was k_Fe=2.75 × 10⁻³ d⁻¹ and α_Fe=0.25 d⁻¹, for which NPP increased from 23 to 30 Pg C yr⁻¹. Note that k_Fe found by this tuning is two orders of magnitude greater than that in the CMIP6 incarnation of the model (Romanou et al., 2013, 2014). However, such an increase is consistent with the observed rates of trace metals with similar scavenging behavior to that of iron (Quigley et al., 2002; Parekh et al., 2004; Lerner et al., 2017). In addition to tuning the iron cycle, we tuned the iron half-saturation constant of N₂-fixation so that the globally integrated rate would fall between 68 and 134 Tg N yr⁻¹. We explored a range between 0.12–12 µM and found that the lower limit provided the best globally integrated value (83 Tg N yr⁻¹), while still allowing the air-sea flux to be at equilibrium.

3.1.2 Forcing fields for carbon cycle

We use greenhouse gas concentrations from transient historical forcing prescribed for CMIP6 (Meinshausen et al., 2017). In the ESM simulations presented here, radiation experiences atmospheric emission‐based CO₂ abundance changes, as described below, while aerosol and ozone constituents are prescribed based on output of an ensemble of GISS‐E2.1 AMIP simulations, which used full atmospheric chemistry and aerosol schemes (OMA; Kelley et al., 2020). The GISS CCv1 model version is labeled GISS‐E2.1‐G‐CC in the CMIP6 repository.

Ocean carbon cycle forcing fields

Forcing (external) fields for the ocean carbon cycle include dust iron deposition, ligands, and riverine nutrients and carbon cycle components delivered to the coastal ocean. These are described below separately.

Dust Iron Deposition

Atmospheric dust falls to the surface through wet deposition, gravitational settling and turbulent removal. Wet deposition is the result of collisions with precipitation either within or below the cloud. As in many models (Adebiyi and Kok, 2020; Meng et al., 2022), dust falls out too fast in ModelE, so that its concentration, especially of larger particles, is unrealistically low downwind of sources (Castellanos et al., 2024), as revealed by ocean deposition measurements over the tropical Atlantic (van der Does et al., 2016, 2020).

The ocean iron cycle is forced with the historical annual cycle dust climatology extracted from the consistent (ModelE2.1) online dust simulations using the GISS dust model (Miller et al., 2006), which has been expanded to include eight externally mixed minerals (illite, kaolinite, smectite, carbonates, quartz, feldspar, iron oxides, and gypsum) plus internal mixtures between minerals and iron oxides (Perlwitz et al., 2015a, b). The masses of free and structural iron and their fractions of total iron have been evaluated using measurements for location at Izaña Observatory (Perez García‐Pando et al., 2016). From the PI simulation with the prognostic dust, we extract 50 years of daily dust fluxes and compute a monthly climatology that is used to force the model. In the ModelE dust module, there are 6 size categories, Clay, Silt1, Silt2, Silt3 and Silt4 and Silt5. The last two describe the largest particles which make small contributions for ocean deposition, since they fall out near their source. For each category, we compute dry and wet deposition fluxes, gravitational settling and ice-ocean flux (in units of 10⁻¹³ kg m⁻² s⁻¹). Silt1 particles contribute the largest proportion to the global average dust flux, the annual emission or the total deposition. This flux affects the timing of dust delivery to the ocean, since dust falling on ice may not be released to the ocean until the ice melts a few months later. Furthermore, we assume that the iron content varies within each dust size class, with clay containing 3.5 % iron and silt containing 1.2 % iron (Fung et al., 2000). Lastly, a constant iron solubility of 2 % is also assumed (Romanou et al., 2013). Figure S1 in the Supplement shows that the peak in dust deposition happens in the spring and summer and mainly in the tropical/subtropical regions of Africa, Saudi Arabia and south-eastern Asia.

Iron Ligands

The fraction of total iron that is complexed with organic ligands (see Appendix, Eqs. A27–A29) is computed using spatially varying but temporally constant prescribed ligand concentrations. Ligands are prescribed from a compilation of ligand measurements obtained at various depths from open ocean waters between 1992–2014 (Caprara et al., 2016). Due to the spatial sparseness of the dataset, we employ a simplistic method of translating the data onto the model grid. First, we split the model's ocean domain into thirteen distinct oceanographic regions, including the Antarctic, the Southern Indian ocean, the Southern Pacific Ocean, the Southern Atlantic Ocean, the Equatorial Indian ocean, the Equatorial Pacific Ocean, the Equatorial Atlantic Ocean, the Northern Indian Ocean, the Northern Central Pacific Ocean, the Northern Central Atlantic Ocean, the Northern Pacific Ocean, the Northern Atlantic Ocean, and the Mediterranean Sea. Then we combine all ligand samples present within each region into a single composite vertical profile. We then bin-average the composite profile using 100 m bins above 1000 m and using 500 m bins below 1000 m. The bin-averaged profile is then smoothed using LOWESS smoothing, where the weighted local linear regressions consider 40 % of the data, and three iterations of regressions are performed where successive iterations apply residual-based reweightings. Finally, the smoothed profile is linearly interpolated on the vertical coordinates of the model. Since a single composite profile is applied to an entire region, the resulting spatial distribution of ligands resembles “patches” (Fig. S2). In future model development we will remove the patchiness.

Riverine inputs of biogeochemical constituents

As described earlier, continental runoff is prognostic in ModelE2.1-CC and delivers chemical constituents to the coastal zone. The concentrations of particulate organic carbon (POC), dissolved organic carbon (DOC), dissolved inorganic carbon (DIC), nitrate, silicate, and iron at all the major and many minor river estuaries are obtained from an annual climatology (da Cunha et al., 2007), and they modulate the biogeochemical characteristics of the freshwater outflow into the ocean at these sites. Therefore, as continental runoff increases (decreases) due to changes in precipitation and evaporation over land, the delivery of these ocean biogeochemical tracers will increase (decrease). Particulate carbon burial into the sediment, although prescribed for all areas shallower than 150 m, takes place in the estuaries where detritus delivered with the riverine flow tends to accumulate. Riverine sources of silica, iron, DIC and DOC concentrations are larger at northern high latitudes (Fig. S3), whereas nitrate and POC show larger values in subtropical regions.

Land carbon cycle forcing fields

The prescription of crop cover change here was updated to the Land Use Harmonization Version 2 Global Carbon Budget 2020 dataset (GCB-LUH2 2020) (Friedlingstein et al., 2020), which provides land cover changes through 2020 (Fig. S4). All crop types in the data set were merged into one crop type and treated as a C3 crop. As in Ito et al. (2020), the natural vegetation boundary conditions and crop monthly LAI are provided from the Ent Global Vegetation Structure Dataset (Ent GVSD) v1.0 (Ent GVSD v.10 Technical Report). The Ent GVSD v1.0 is derived from satellite data sources and includes land cover types and monthly varying LAI from the Moderate Resolution Imaging Spectroradiometer (MODIS) (Gao et al., 2008; Myneni et al., 2002; Tian et al., 2002a, b; Yang et al., 2006) for the year 2004, and tree heights from Simard et al. (2011), who utilized 2005 data from the Geoscience Laser Altimeter System (GLAS) aboard the ICESat (Ice, Cloud, and land Elevation Satellite). Specific leaf area (carbon mass per leaf area) data from the TRY database of leaf traits (Kattge et al., 2011) was classified for the Ent TBM 12 plant functional types (PFTs). With historical crop cover change, the natural vegetation fractions are rescaled. The Ent GVSD also uses nearest-neighbor extensions of LAI for grid cells where the land cover change results in cover types that were not there in the 2004 MODIS land cover. These boundary conditions allow for geographic variation in subgrid fraction, height, and LAI for each PFT, such that biomes can be different mixes of PFTs with different demography and canopy structure.

The use of the natural vegetation boundary conditions and an older historical crop cover data set with the spinups in the ModelE C4MIP contribution (Ito et al., 2020) resulted in annual carbon fluxes and stocks well in the middle of CMIP5 predictions, as well as in the middle of emerging CMIP6 comparisons. However, the prescribed leaf area index and stem biomass do not allow for prediction of phenomena like boreal greening with climate warming, and the land model does not yet have deforestation or fire as major components of the land carbon budget. The CO₂ fertilization effect is manifested by enhanced photosynthetic uptake, boosts to litterfall to the soil, and increased water use efficiency.

3.1.3 Initial conditions

Land carbon initialization and spinup

At initialization, plant labile carbon has default values in relation to plant individual maximum foliage biomass by PFT, twice this amount for woody plants and four times this amount for grasses, sufficient to grow this foliage in the next growing season. Soil carbon may start from zero or be initialized from an input file for the upper 30 cm of soil derived from the International Soil Reference and Information Centre – World Soil Information (ISRIC-WISE) data set (Batjes, 1996). All coupled carbon cycle simulations require that these plant labile and soil carbon pools be spun up to equilibrate with climate.

The plant labile carbon pool typically spins up in 30 years. To spin up soil carbon, the same method is used as described in Ito et al. (2020): the equilibrated preindustrial control simulation is run for 9 years and the simulated meteorology (soil moisture and temperature) and litterfall is saved in memory at hourly time steps and then iterated over the soil 750 times for a total of 6750 simulated years. The whole cycle typically takes less than 24 world clock hours on a modern cluster using 88 CPU cores. This is typically repeated 2–3 times to achieve a balance an order of magnitude smaller than the 0.1 GtC yr⁻¹ in the protocol for C4MIP.

Ocean carbon initial conditions and spinup

For the ocean carbon cycle, macronutrients (nitrate and silicate), DIC, Alkalinity, and O₂ were initialized from the Global Ocean Data Analysis Project, version 2 climatology (GLODAPv2; Lauvset et al., 2016), which maps fully quality controlled and internally consistent ocean biogeochemical data collected between 1972–2013 onto a 1° × 1° lat/lon grid and 33 vertical layers. Phytoplankton were initialized using equations derived from the synoptic relationships found between satellite-derived chlorophyll a and diagnostic pigments (Hirata et al., 2011). These equations are:

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathrm{Diatoms}=\mathrm{1.3272}+\exp (-\mathrm{3.9828}\times +\mathrm{0.4003})\text{(2)}& {\displaystyle }\begin{array}{rl}\mathrm{Chlorophytes}& =(\mathrm{0.2490}/y)\\ & \exp (-\mathrm{1.2621}(x-\mathrm{0.5523}{)}^{\mathrm{2}})\end{array}\text{(3)}& {\displaystyle }\begin{array}{rl}\mathrm{Cyanobacteria}& =(\mathrm{0.0067}/(\mathrm{0.6154}\cdot y\left)\right)\\ & \exp (-\mathrm{19.5190}(x+\mathrm{0.9643}{)}^{\mathrm{2}}/\mathrm{0.3787})\\ & +\mathrm{0.1027}{x}^{\mathrm{2}}-\mathrm{0.1189}x+\mathrm{0.0626}\end{array}\text{(4)}& {\displaystyle }\begin{array}{rl}\mathrm{Coccolithophores}& =\mathrm{1}-(\mathrm{0.9117}+\exp (-\mathrm{2.7330}x\\ & +\mathrm{0.4003}){)}^{-\mathrm{1}}-(\mathrm{0.1529}\\ & +\exp (\mathrm{1.0306}x-\mathrm{1.5576}){)}^{-\mathrm{1}}\\ & +\mathrm{1.8597}x-\mathrm{2.9954}\\ & -(\mathrm{0.2490}/y)\exp (-\mathrm{1.2621}(x\\ & -\mathrm{0.5523}{)}^{\mathrm{2}})\end{array}\end{array}$$

where x is the base-10 logarithm of the chlorophyll a concentration and y is chlorophyll a concentration. Chlorophyll a used to compute initial phytoplankton functional groups is obtained from averaging the SeaWiFS seasonal climatologies between 1998–2010 (https://oceandata.sci.gsfc.nasa.gov/l3/, last access: 21 August 2026). Equations (1)–(4) above were used to initialize phytoplankton in the upper 75 m of the water column; below this depth, phytoplankton concentrations are initialized at 0 mg Chl m⁻³. Ammonium, zooplankton, and detrital material (carbon, silica, and iron) are initialized with spatially invariant values 0.5 µM, 0.05 mg chl m⁻³, and 0 µg L⁻¹, respectively. Finally, dissolved iron is initialized by splitting the domain into the same 13 geographic regions used to derive the prescribed ligand concentrations (see Forcings). Within each region, an iron:nitrate ratio is specified that is laterally and vertically invariant. The iron:nitrate ratios are obtained from a compilation of in-situ measurements (Fung et al., 2000), and are shown in Table S1 in the Supplement. After iron:nitrate ratios are specified, iron concentrations are initialised by multiplying, in each grid cell, this ratio by the initial nitrate concentration from GLODAPv2.

After the tuning, the ocean carbon model spin ups were carried out for approximately 300 years off an equilibrated simulation of the pre-industrial ocean-atmosphere-ice system which was itself spun up for 7000 years.

Table 1 ModelE2.1-CC experiments described in this paper.

3.2 Experimental design

For model simulations (Table 1), we conducted a preindustrial control with 1850 conditions, a transient historical run over 1850–2021, future scenarios for SSP2-4.5, SSP5-8.5, and SSP5-3.4 for 2022–2100 and a compound 1 % increase in CO₂. For the historical simulations and the SSPs we used CO₂ emissions (not atmospheric concentrations) from input4mips files (Hoesly et al., 2018; Gidden et al., 2019) for the period 1850–2014 and constant thereafter, while for non-CO₂ GHGs we follow CMIP6 specifications for prescribed concentrations from 1850–2014 (Meinshausen et al., 2017) and use concentrations from the NOAA Global Monitoring Laboratory from 2015 to 2021 (Lan et al., 2022b). For the eA0 (the compound 1 % run) we back-calculated implied emissions of CO₂ from a concentration driven compound 1 percent run and we used those to force the eA0 scenario. All other GHGs and aerosols are prescribed following the CMIP6 specifications. The datasets used for model evaluation are shown in Table 2.

3.3 Datasets used for model assessment

The datasets used for data evaluation are described in Table 2, where the fields used, the spatial and temporal coverage, the resolution and a key reference are also provided.

Table 2 Observational and reanalysis data used for model assessments.

4 Results

4.1 Preindustrial control simulation

Figure 3 shows that surface exchanges and properties are in equilibrium for at least the last 500 years of the pre-industrial control simulation (under 1850 conditions) and Table 3 lists their long-term mean, standard deviation and trend. Radiative balance at the top of the atmosphere (TOA) is −0.04 ± 0.34 W m⁻² (Fig. 3a and Table 3). Surface air temperature (Fig. 3b) has a small drift of about 0.02 °C per 100 years. Net land carbon flux is 0 ± 0.65 PgC (Fig. 3c), with a long-term trend of −0.05 PgC per 100 years (Table 3) while gross primary production (GPP), autotrophic and soil respiration (Fig. 3d, e, f respectively) have small trends, and so do net ecosystem production and net primary production over land (Fig. 3g and h). The largest variability arises from GPP while the largest trend from soil respiration (Table 3). Ocean carbon uptake (Fig. 3i) is at equilibrium with some slow variations of about 50-year scales and of the order of ±0.12 PgC yr⁻¹ (Table 3). Net primary productivity in the ocean (Fig. 3j) is about 30 PgC yr⁻¹ which is low compared to 50–60 PgC yr⁻¹ from observations (Westberry et al., 2023). Net heat flux (Fig. 3k) is out of the ocean by 0.07 W m⁻² with a large interannual variability of the order of 0.49 W m⁻² (Table 3) and a small positive trend of $\mathrm{2.6}\times {\mathrm{10}}^{-\mathrm{2}}$ W m⁻² per 100 years. Global mean precipitation is about 2.95 mm d⁻¹. Figure 3 also shows a comparison with the CMIP6 version of the model, denoted here as “CC” (pink lines). Due to the incomplete spinup of that version of the model, both NPP on land (including gross primary production, autotrophic and soil respiration) and in the ocean had significant trends, which imparted a 0.2 °C warmer temperatures.

Figure 3 Temporal evolution of preindustrial control fields: (a) net radiation at top of the atmosphere (TOA), (b) surface air temperature, (c) land carbon flux, (d) gross primary production (GPP), (e) autotrophic respiration, (f) soil respiration, (g) net primary production on land (NPP = GPP − autotrophic respiration) (h) net ecosystem production (NEP = NPP − soil respiration), (i) ocean air-sea carbon flux, (j) net ocean primary production (NPP ocean), (k) net heat uptake at the ocean surface, (k) precipitation. Annual mean fields for the CC2 model are shown in gray and for the previous model version (CC) are shown in pink. 20 year averages are shown in black. The horizontal axis is in years from the PI control run. Units are shown on the vertical axes. All fluxes are positive down into the surface (land or ocean).

Table 3 Mean, standard deviation and trend over the last 100 years of the PI control run.

Spatial distributions of the preindustrial fields are shown in the Supplemental Material (Fig. S5). Annual mean atmospheric CO₂ concentrations peak at the equator (Fig. S5a) with values around 284.5 ppm and lowest values at high latitudes about 283.7 ppm. The largest interannual variability (Fig. S5b) appears over tropical land regions and in the equatorial Atlantic and is associated with the large interannual variability due to ENSO (Kelley et al., 2020). Surface air temperature (Fig. S5c and d) shows a strong zonal variation with warmer temperatures over low latitudes and significantly colder at high latitudes, North America and Asia. The largest variability appears over the eastern tropical Pacific and in the Arctic Ocean, rim associated with natural modes of variability which are stronger than observed in our simulations (Kelley et al., 2020). The air-to-land flux is largest in the continental North America and the northern Eurasian continent (Fig. S5e). The Amazon is a weak source due to Model's warm and dry bias in this region, and the high latitudes are net sinks due slow soil carbon storage in these cold regions, while the largest variability is present in Europe and North America, tropical Africa and the Amazon regions (Fig. S5f). The ocean carbon flux has a strong sink in the North Atlantic and the Southern Ocean (Fig. S5g), particularly the Labrador Sea and the Greenland-Iceland-Norway (GIN) Seas. World ocean subtropical gyres are regions of carbon uptake whereas the tropical belt is outgassing CO₂ to the atmosphere with peak values in the eastern Equatorial Pacific and secondary peaks in the Benguela and Canaries current upwelling regions. The largest variability in the ocean sink is exhibited in the North Atlantic subpolar gyre (mostly in the GIN seas and the Gulf Stream separation region) as well as in the tropical upwelling regions (Fig. S5h). The net ocean heat flux peaks in the tropical belt (Fig. S5i) where the ocean takes up heat from the atmosphere and is lowest in the North Atlantic, specifically along the Gulf Stream, the Labrador and Irminger Seas and the GIN seas, where the ocean releases heat to the atmosphere. The largest ocean heat flux variability emerges in the eastern tropical Pacific Ocean, and the North Atlantic subpolar regions (Fig. S5j).

4.2 Historical simulations

Global mean atmospheric CO₂ concentrations agree with observations (Fig. 4a) with the largest discrepancies between 1900 and 1970. CC2 historical concentrations are improved compared to CC (Ito et al., 2020) due to fixing the bug with land-use change forcing files. The small underestimation of atmospheric CO₂ in CC2 stems from insufficient land uptake and is consistent with the higher energy imbalance at the top of the atmosphere (Fig. 4e). Both land and ocean carbon uptake underestimated in comparison to CarbonTracker (Figs. 4d, f and S4b) over 2000–2022 but may be overestimated over the longer historical period (more discussion on land and ocean carbon sinks below) while the ocean sink is overestimated compared to SOM-FNN. Spatially positive biases up to 1.5 ppm are found over the northern hemisphere land and ocean especially the high latitudes and similar negative biases are found in the Southern Hemisphere, especially over the tropical forest regions (Fig. 5a1, a2). Generally, CC2 carbon fields as well as surface air temperature are improved compared to the CC version of the model (pink lines in Fig. 4). This is due to better spin up of the soil carbon pool and the preindustrial ocean carbon flux. The carbon disequilibrium affects substantially the cumulative heat uptake in Surface air temperature agrees well with reanalysis products in terms of trend (Fig. 4b). The model overestimates the cooling associated with Pinatubo (1991–1992) and generally shows larger interannual variability than the observations after 1950 possibly related to the strong ENSO variability (Kelley et al., 2020). Surface temperature anomalies from preindustrial track the GISTEMP observations very well (Fig. 4c) particularly after the 1950s. Between 1875 and 1950 the model is warmer than the baseline than the observations, and this is even though there is somewhat less atmospheric CO₂ in the model than was observed. Spatially, the biases are more pronounced in the Arctic (by about 4 °C), the GIN Seas and the Weddell Sea (by about 2 °C), while the Eastern Pacific is warmer in the model by about 1–2 °C, the subpolar North Atlantic is colder by about 1–2°, the Southern Ocean in general is warmer by about 1 °C, and the Indian Ocean is colder by 1 °C (Fig. 5b1, b2).

Figure 4 Historical annual Time series of (a) atmospheric CO₂ concentration from the model (eHI) and observations from input4mips (Meinshausen et al., 2017) from 1850 to 2015, and from the NOAA Global Monitoring Laboratory from 2015 to 2021 (Lan et al., 2022a, b), (b) surface air temperature (annual mean timeseries) and superimposed estimates from ERA5 and MERRA2 (see Table 3), (c) change in surface air temperature from the average over the period 1951–1980 and GISTEMP, (d) air-to-land carbon flux from the model (eHI) compared to observations from the Carbon Tracker (shaded is the uncertainty of the annual mean based on the averaging of inversions that consider different prior estimates of fluxes and their uncertainties; Jacobson et al., 2023), (e) net radiation at the top of the atmosphere (TOA) compared to observations from CERES 2000–2026 (f) air-sea carbon flux from the model compared to observations from SOM-FFN and Carbon Tracker (2000–2022) (shaded is the uncertainty around the annual mean, stemming from a combination of uncertainties including (i) the gridding of SOCAT observations, (ii) the RMS error between the SOM-FNN estimates and the gridded observations, and (iii) the uncertainty stemming from the chosen air-sea flux parametrization; Landschützer et al., 2016), (g) cumulative heat uptake in the model and comparison to estimates from Argo, both timeseries are anomalies from the 1993–2022 mean, the model timeseries are dedrifted using the linear PI drift, (h) ocean heat flux and comparison to estimate from ECCO-V4r4, years 1992–2017, (i) Precipitation. Black lines refer to the current model version (CC2) and pink lines to the previous model version (CC). Observations are shown in color. It should be noted here that model forcings post-2014 remain constant at 2014 levels except for atmospheric non-CO₂ greenhouse gases which follow observations from NOAA-GML (Lan et al., 2022b), a fact to be considered in model-to-observation comparisons post 2014 shown in this figure.

Precipitation (Fig. 5c1, c2) in the historical eHI simulation compared to land surface precipitation observations by the Global Precipitation Climatology Centre (GPCC), National Oceanic and Atmospheric Administration (NOAA), Full Data Reanalysis Version 6.0, over 2001–2010, shows that the model has an extensive dry bias over the tropics, a persistent model bias also documented in Kelley et al. (2020) for the ModelE-2.1 contribution to CMIP6. This bias is particularly pronounced in the Amazon Basin, which is a persistent issue in GCMs in general (Respati et al., 2024). The ModelE dry bias is also strong in India and Indonesia. The model has a wet bias in Western China and the West Coast of South America along the Andes.

Average sea level pressure (Fig. S6a, b) is higher than observations by more than 10 mb over Greenland and Antarctica and the mountainous ranges in south-east Asia and South America. It is systematically underestimated over the ocean by about 5–10 mb. Surface winds are weaker than in observations (Fig. S6c, d) by more than 5 m s⁻¹ over Greenland and Antarctica, 3–5 m s⁻¹ over the mid- and southwest US and Patagonia, Europe and Northern Africa, and by 2–4 m s⁻¹ over ocean subpolar gyres, an issue with all previous modelE versions, and which led to the inclusion of the “gustiness” term (see Sect. 2.7).

Global mean air-to-land carbon (net biome productivity, NBP) flux (Fig. 5d) hovers at a net flux of zero in the 20-year moving average over 1850–1942, after which the trend becomes a net sink. Over 2000–2022, it is underestimated compared to Carbon Tracker as it lies below the mean and outside the envelope of interannual variability in this dataset. Interannual variability is also much higher amplitude in the model, due sensitivity of the land carbon to ModelE's higher amplitude interannual surface temperature variability from simulated El Niños. A more detailed breakdown of eHI time series of the land fluxes (GPP, autotrophic respiration, soil respiration, emissions due to land use change, NPP, NBP, NEP) and soil carbon stocks are provided in Fig. S7, which shows both the annual totals and their 20-year running means to show trends. These show all natural fluxes accelerating after ∼ 1960. Regionally, the underestimation (which exceeds 20 mol CO₂ m⁻² yr⁻¹) stems from weaker Amazon and African rainforest uptake as well as weaker boreal forest uptake (Fig. 5d). The cumulative land carbon sink 1850–2023 estimated by the Global Carbon Budget 2024 is 220 ± 60 PgC (Friedlingstein et al., 2025), whereas in the historical simulation the model reaches a cumulative sink of only 56.6 PgC by 2021 (Fig. S7b). CarbonTracker and other flux inversion models estimate generally lower NBP than DGVMs, which Bastos et al. (2020) attribute to disagreement on NBP interannual variability responses to El Niño and uncertainties in land use emissions, particularly in South America and Southeast Asia. Meanwhile, based on remote sensing observations of vegetation change, Randerson et al. (2025) find that forest growth is far less than predicted by the ensemble of DGVMs in the GCB2023; the remote sensing data yield an estimated net land carbon sink over 2000–2019 of 0.8 ± 0.7 PgC yr⁻¹, whereas the GCB2024 DGVMs estimate a high 1.32 ± 0.2 PgC yr⁻¹ due to unrealistic predictions of forest greening or densification. Our model estimates 1.12 ± 1.04 PgC yr⁻¹ over this period. The GCB2024 DGVM estimates of earlier preindustrial through 20th century deforestation emissions also rely on highly uncertain predictions of forest biomass and land use change. Observational constraints of vegetation biomass remain a challenge: inventory and satellite surface observations of global total vegetation biomass estimates still range over a wide 380–536 PgC (Erb et al., 2018), while more recent radar and lidar satellite data products of forest heights are still being processed to produce new global biomass estimates, which we detail more later.

In Fig. 5e, the emissions from land cover and land use change (LCLUC) are plotted together with the estimates of the Global Carbon Budget 2024 (GCB 2024) (Friedlingstein et al., 2025), which span 1960–2024 only. These show that our land cover change (LCC) emissions strongly underestimate the mean estimate of the GCB2024 by more than 1 PgC yr⁻¹, and is outside their low (less negative) estimate; the model underestimate is very likely because our model does not include deforestation/forest regrowth, wood harvest, or peat emissions, but only land cover conversion from natural vegetation to croplands. Therefore, our model's LCC emissions from the preindustrial period to 1960 likely also underestimate total LCLUC emissions by a similar magnitude, leading to our underestimate of atmospheric CO₂ (Fig. 4a).

Based on our model physics, we would expect the simulated underprediction of atmospheric CO₂ to be explained by the ModelE land model's lack of full representation of certain components of land use change (LUC). Figure 4e shows the eHI net carbon emissions due to LUC compared to the breakdown of components of LUC change emissions estimated by the Global Carbon Budget 2024 (GCB2024, Friedlingstein et al., 2025). Although the model represents deforestation that converts forest to croplands or other cover types, it does not represent wood harvest or deforestation that results in degraded forest; also, the model is not a fire-enabled model, so fire emissions associated with deforestation and peatlands are not represented. Therefore, the eHI emissions miss on average 1.0 ± 0.3 PgC yr⁻¹ over 1959–2021 compared to GCB2024 mean net LUC emissions. Since deforestation emissions were higher over the earlier historical period and the model misses these emissions, we would expect the underpredicted net emissions to the atmosphere to result in the model's slightly lower atmospheric CO₂ concentrations, which catch up to the observed as deforestation declined in recent decades. Given the much lower NBP estimated by Randerson et al. (2025), the high bias of GCB2024, and uncertainties in flux inversion models, despite the eHI simulation's comparatively low cumulative land carbon sink, it may be that ModelE's land carbon sink is overpredicted in the context of the model's simulated flux components (Fig. S7), explaining our low predicted atmospheric CO₂. These discrepancies highlight the continued large uncertainties in land carbon stocks and fluxes.

Model SST (averaged over 1991–2020) is colder than World Ocean observations by about 2–4 °C in the subArctic regions (Fig. 5f1, f2). The tropical eastern Pacific Ocean and Tropical Atlantic are warmer than observations by about 2–4 °C, whereas the Southern Ocean is warmer than observations by about 2–4 °C. The North Atlantic Current region and the Arctic are colder by about 2–5 °C. The model historical climatology 1991–2020 of the sea surface salinity (SSS, Fig. S6f) is fresher than observations (Fig. S6e) by about 4 psu in the Eurasian side of the Arctic Ocean, and about 2 psu in the tropical regions of the central Pacific Ocean, the Atlantic Ocean and the Indian Ocean. The model overestimates salinity by about 2 psu in the tropical eastern Pacific Ocean and in the Arabian Sea and the Bay of Bengal in the Indian Ocean. The SSS biases largely reflect biases in evaporation and precipitation patterns over the global ocean (Kelley et al., 2020; Miller et al., 2021) which are in turn related to biases in surface pressure and surface winds (Fig. S6a, b, c, d).

The ocean carbon flux (Fig. 5f) lies between the observation-based estimates SOM-FFN and Carbon Tracker and within both datasets' interannual variability. Spatially, ocean carbon uptake negative biases (up to 5 mol CO₂ m⁻² yr⁻¹) are associated with strong wind-driven circulation and fronts (e.g. over the Gulf Stream, the Equatorial trade winds and the Southern Ocean subtropical front (Fig. 5g1, g2) whereas positive biases (up to 6 mol CO₂ m⁻² yr⁻¹) are linked to convection sites in the North Atlantic (Labrador Sea, Irminger Sea and GIN Seas) and the Antarctic Circumpolar Current (ACC) region. Subtropical gyres and the tropical Atlantic and Indian Oceans also show positive biases. Tropical ocean outgassing is also stronger in the GISS model, despite the weaker than observed winds and warmer SSTs there, because the pco2 disequilibrium is larger than observed, due to a combination of low productivity (as we describe below) and also because the vertical sinking of organic carbon below the thermocline is too small, which is discussed in more detail in Lerner et al. (2024).

Net primary productivity is underestimated everywhere in the world's oceans compared to CAFE observations by about 100–200 mgC m⁻² d⁻¹ (Fig. 5h1, h2), except in the subtropical gyre in the North Atlantic and in the Benguela and Canaries current regions. The generally low bias in NPP results from at least two issues. First, there was an error in our ocean radiative transfer routines that compute the direct and diffuse radiation that enters and leaves each ocean layer. For each layer, we incorrectly used the layer depth (distance from the ocean surface to the layer), rather than the layer thickness (distance between upper- and lower-layer boundaries), when computing the amount by which direct radiation should be reduced from the top to the bottom of each layer. This resulted in too little shortwave radiation reaching subsurface waters, reducing NPP. Second, phytoplankton in our model currently do not acclimate in ambient light levels, so that chl:C ratios in the sun-lit surface are the same as those deeper in the water column. Acclimation is expected to increase the growth of phytoplankton and NPP under low-light conditions (Geider et al., 1998), so the absence of this process in our model likely also contributes to lower productivity levels. We also should note that our comparison only encompasses a single NPP product, that derived from the CAFE algorithm applied to the Ocean Colour Climate Change Initiative merged remote sensing data product (Ryan-Keogh et al., 2023). NPP products are known to vary widely depending on the algorithm applied, particularly in the Southern Ocean but also to a significant extent in the Northern high latitudes. For example, using the same dataset, the Eppley-VGPM algorithm yields roughly the productivity of the CAFE algorithm in the Southern Ocean (Ryan-Keogh et al., 2023). Therefore, despite improvements that need to be incorporated, how well our model performs in some ocean environments, particularly the Southern Ocean, remains unclear. Finally, we note that a few regions, in particular the Benguela and Mauritanian upwelling regions, display positive NPP biases. The cause of these biases remains a subject of future investigation, although a potential contributing factor is the fact that our model lacks water column (Lam and Kyupers, 2011) and benthic denitrification (Dale et al., 2014), which could result in too much upwelling of nitrate in these regions.

Figure 5 Observational and observation-based estimates and model spatial biases averaged over the period of the observations. (a1) Atmospheric CO₂ concentration anomalies from the global mean from Carbon Tracker averaged over 2000–2014. (a2) Atmospheric CO₂ concentration anomalies, model biases from Carbon Tracker, (b1) GISTEMP average 2000–2014, (b2) Surface air temperature model biases from GISTEMP v4 (anomalies from 1951–1980) averaged over 2000–2014, (c1) Precipitation average 2001–2010 from GPCPv4, (c2) historical precipitation biases over land fro GPCPv4, (d1) Land carbon flux taken from Carbon Tracker between 2000–2014, (d2) land carbon flux model biases in the model compared to Carbon Tracker, (e1) land GPP taken from FluxSat between 2000–2009, (e2) GPP model biases in the model compared to FluxSat, (f1) Sea surface temperature (SST) averaged from the period of 1991–2020 from World Ocean Atlas 2023 (WOA2023), (f2) historical SST biases in the model compared to WOA2023, (g1) ocean carbon flux observations from SOM-FFN averaged over the period 2000–2014, (g2) Ocean carbon flux biases in the model compared to SOM-FFN, (h1) Ocean net primary productivity (NPP) averaged from 2000–2014 from CAFE (Silsbe et al., 2016), (h2) ocean NPP biases in the model compared to CAFE, (i1) ocean heat content from ECCO-V4r4, averaged over the period of 1993–2017, (i2) Ocean heat uptake bias in the model compared to ECCO-V4r4.

Ocean heat uptake anomaly from a late-20th century baseline is well in agreement with estimates from Argo floats (cumulative; Fig. 4g) and with ECCO-V4r4 (non-cumulative anomaly from preindustrial; Fig. 4h) as well as variability in the latter dataset. Spatially, the patterns of average differences between the GISS model and ECCO-V4r4 for the period 1993–2017 are within ±50 W m⁻², with positive differences (GISS minus ECCO-V4r4) dominating the low latitudes and negative differences the high latitudes and the eastern Pacific Ocean (Fig. 5i1, i2). More specifically, ocean heat uptake averaged over 1993–2017 is underestimated in the model by more than 40 W m⁻² in the North Atlantic subpolar gyre and the ACC region (Fig. 5i1, i2) due to the lower wind speeds there, whereas it is overestimated over other regions with strong currents (e.g. the eastern Pacific and eastern Atlantic equatorial current regions, the Brazilian and Benguela current regions, the Agulhas retroflection region, the eastern part of the western Australian current region, and the Kuroshio region) because the sea surface temperatures are much higher (Fig. 5f1, f2).

The changes from CC to CC2 are synopsized in Fig. S8: atmospheric CO₂ increases (decreases) in CC2 by about 4 ppm in the Northern (Southern) Hemisphere. The largest changes in land carbon fluxes are an increase in Russia by ∼ 10 mol CO₂ m⁻² yr⁻¹, and declines in northwestern Europe and Southeast Africa by ∼ 5 mol CO₂ m⁻² yr⁻¹ (Fig. S8b). GPP is generally larger in CC2 than CC, showing increase of up to 40 mol CO₂ m⁻² yr⁻¹ in central Africa and South America, though there are some regions such as the Amazon rainforest that show small declines (< 5 mol CO₂ m⁻² yr⁻¹) instead (Fig. S8c). The ocean flux shows moderate increases in the subtropical regions and Southern Ocean (3–5 mol CO₂ m⁻² yr⁻¹), and a large decline in the equatorial Pacific ($<-\mathrm{5}$ mol CO₂ m⁻² yr⁻¹; Fig. S8d). Finally, NPP shows a large increase and decline of ∼ 200 mgC m⁻² d⁻¹ in the subtropical regions and Southern Ocean, respectively, as well as smaller decline in the subpolar North Pacific and Atlantic (∼ 50 mgC m⁻² d⁻¹; Fig. 9e)

Seasonality

The seasonal cycle of the air‐to‐sea flux of CO₂ and the partial pressure of CO₂ in the sea water have been extensively evaluated against several observational data sets (Lerner et al., 2021).

For the land carbon, seasonal GPP fluxes are shown in Fig. S9, compared to estimates by FluxSat v2 (Joiner and Yoshida 2021). ModelE-CC estimates annual GPP as 118.6 ± 1.5 PgC yr⁻¹ over 2000–2009, and 124.7 ± 2.0 PgC yr⁻¹ over 2014–2021, which is currently at the lower end of both multi-model and observation-based estimates ranging 110–170 PgC yr⁻¹ (Anav et al., 2015; Cheng et al., 2017; Xu and Chen, 2024). The latest Global Carbon Budget 2024 cites a mean annual GPP of 130 PgC over 2014–2023 (Friedlingstein et al., 2025).

The land seasonal regional GPP patterns qualitatively agree with FluxSat estimates (Fig S3). However, there are stark regional differences. ModelE-CC simulates much higher GPP in the Amazon and Congo Basins in all seasons compared to FluxSat. Since ModelE has a warm/dry bias in these regions, we would expect GPP to be low. The FluxSat data are expected to be a lower bound estimate when based on flux tower data and not utilizing solar-induced fluorescence (SIF) data. Therefore, the higher ModelE GPP may be appropriate. In arid and semi-arid regions outside the Amazon and Congo in South America and Africa, ModelE predicted GPP is much lower than FluxSat; these mostly mixed shrubland and grassland areas are strongly controlled by soil moisture availability, and ModelE has a dry bias in these areas, leading to lower GPP. The shrubland PFT is also one for which the Ent TBM lacked a Fluxnet site to evaluate simulated fluxes, so parameterizations of leaf photosynthetic capacity and respiration may be improved. In the northern high latitudes, the season MAM exhibits delayed spring greenup and the summer has low productivity due to ModelE's extended cold bias, while mid-latitude croplands in the Eastern U.S. and Europe are predicted to have much higher GPP than FluxSat estimates. The Northern Hemisphere differences between ModelE and FluxSat decline in the fall, SON, when deciduous trees senesce.

4.3 Heat and Carbon storage on land and in the ocean

Long term storage of heat and carbon in the ocean and on land is important for the control they exert on future climate change.

4.3.1 Ocean storage

A comparison of heat storage trends at different levels in the ocean in the GISS model and the ECCO-V4r4 ocean state estimate (Fig. 6) shows that, in general, the GISS ocean model stores more heat in the subtropical gyres in the North Atlantic as well as the Atlantic sector of the Antarctic convergence zone at all levels. ECCO-V4r4 stores more heat in the subpolar North Atlantic and in the deep Southern Ocean around Antarctica. The deep western boundary current which is clearly outlined in the GISS ocean model at depths below 600 m is not as clearly identified in ECCO-V4r4. The GISS model also stores more heat in the upper North Pacific Ocean compared to ECCO-V4r4, however the latter stores more heat in the tropical Pacific. Differences in the storage regions and depths will influence the response of the ocean model to future forcing changes, as heat stored closer to the surface and at higher latitudes will more easily be released to the atmosphere compared to heat stored at lower latitudes.

Figure 6 Ocean heat content trend (J m⁻² s⁻¹) for the period 1993–2017 in ECCO-V4r4 and GISSE2.1-CC2, vertically integrated: (a, b) from surface to 600 m, (c, d) 600–2000 m, (e, f) 2000 m-bottom.

Figure 7 Integrated DIC (kg m⁻²) from GLODAP and GISS eHI for the period 1993–2014 at different levels: (a) surface, and integrated over (b) 0–600 m, (c) 600–2000 m, (d) 2000 m – bottom.

Compared to GLODAP (Fig. 7), the model overestimates carbon stored in the upper ocean (0–600 m) of the tropical/subtropical regions, especially in the Atlantic (Fig. 7c, d). The model also underestimates the storage of carbon in intermediate waters (600–2000 m), particularly in the North Pacific (Fig. 7e, f). These findings are consistent with those of Lerner et al. (2024) that the flux of carbon beneath 100 m is too high and the transfer efficiency (ratio of carbon export at depth to carbon export at 100 m) is too small, although that study was limited to the equatorial Pacific. High export fluxes near the surface may be due to overgrazing by zooplankton, as zooplankton loss terms (via linear and quadratic mortality: Eq. A32) produce the bulk of the detrital organic carbon in the model. On the other hand, the lack of remineralized DIC in deeper waters below 600 m may be due to the absence of processes that could contribute to higher detrital sinking speeds, such as a lack of ballasting by inorganic material (carbonates and lithogenics; Armstrong et al., 2001), and our use of a single-particle class model rather than multiple size classes that would encompass a wide range of sinking speeds as found in the real ocean (Burd, 2015).

4.3.2 Land carbon storage Plant labile carbon and soil carbon distribution

With canopy structure and seasonal LAI otherwise prescribed in Ent, there are two prognostic carbon pools: labile carbon (also known as non-structural carbohydrate, NSC) and soil carbon. The labile carbon pool is the net of carbon uptake from photosynthesis plus retranslocation from senescing foliage minus autotrophic respiration. In the Ent TBM, autotrophic respiration for most PFTs is tuned to Fluxnet sites to achieve carbon balances with well-quantified site canopy structure (Kim et al., 2015). The global prescribed canopy structure in the simulations in this study may not be in climatological equilibrium with the simulated climate for several reasons. The LAI is from observations from the specific years 2004 for LAI and 2005 for canopy heights. These satellite observations have errors or simplifications, where the spikes in MODIS LAI are largely smoothed but the interpolated monthly LAI updates in ModelE will miss the shorter time scale of green-up and senescence. The forest canopy heights derived from ICESat/GLAS used in these runs are now known to be ∼ 6 m too tall, given the recent more accurate observations of the Global Ecosystem Dynamic Investigation (GEDI) mission (Dubayah et al., 2020); this results in excess stem biomass from tall trees in the prescribed canopy structure, resulting in 533 PgC in aboveground biomass (AGB) in the eHI experiment (Fig. 8a). This is nearly twice that estimated from harmonized inventory and satellite cover data by Spawn et al. (2020; Fig. 23b, e), which is 287 PgC AGB and total 409 PgC including belowground biomass for the year 2010, or from the GEDI lidar data (Fig. 8c, f; Dubayah et al., 2022), which is 283 PgC AGB and limited to south of 53° N. The GEDI data product predicts lower tropical biomass than Spawn et al. (2020) and higher in the low- to mid-latitudes. Liu et al. (2025) estimate global AGB to be 378 PgC from satellite radar. Figure 8e, f difference maps show that the tall bias of Ent TBM forests in the Amazon and Congo Basins and in the boreal zone results in high AGB compared to the Spawn et al. (2020) and GEDI estimates. Grassland areas in eHI (Great Plains, United States; Sahel, Africa) have higher AGB than the observational estimates due to a bug in labile carbon storage for grasses in Ent. Besides only recently available satellite data on vertical canopy structure at appropriate spatial scales for vegetation heterogeneity, plant allometry parameterization remains a large issue for all DGVMs for estimating biomass (Grant et al., 2025; Liu et al., 2025; Jucker et al., 2022; Massoud et al., 2019), resulting in large uncertainties in autotrophic respiration that is a function of the plant biomass. Finally, the GCM simulated climate has regional biases; in ModelE2, the tropics have strong warm/dry bias, and in the high latitudes, a strong cold bias causes prolonged frozen soil, hence prolonged simulated winter-spring water stress and low carbon uptake. In the interest of maintaining canopy structure as close to observations as possible rather than predicting wrong vegetation cover and mortality with climate biases, we opt for managing imbalances in the labile carbon pool, instead.

The map in Fig. 8a of plant labile carbon identifies regions where there is insufficient carbon uptake in both boreal evergreen needleleaf and tropical evergreen broadleaf forests. The carbon deficit in forests in these regions is due both to simulated climate biases and overestimated forest stem biomass. The boreal zone is 8–10 °C too cold in ModelE, with frozen soils persisting from the winter into what should be the greenup period, causing water stress and reduced GPP; the Amazon Basin has a warm/dry bias causing water stress where there should not be and also reducing GPP. Excess stem biomass from the overestimation of height from the ICESat/GLAS estimate results in high mass-based autotrophic respiration. Thus, the underpredicted GPP and overestimated autotrophic respiration in these forest types results in plant carbon deficits and labile carbon near zero. The label carbon map also highlights an error in parameterization of grasses that results in too high labile carbon by a few kg m⁻², particularly in the U.S. interior as was noted also for Fig. 8e, f.

In Fig. 8b, soil carbon distributions reflect high storage in the high latitudes where permafrost tends to be (ModelE does not explicitly simulate permafrost). The total soil carbon averaged over years 2002–2020 of eHIis 1693 PgC, which is virtually the same as the 1700 PgC from inventory estimates compiled by Canadell et al. (2021) and cited by the GCB2024.

Figure 8 Land carbon stocks. Aboveground plant biomass (AGB) from (a) the historical simulation (eHI), (b) harmonized inventory and satellite-based estimate by Spawn et al. (2020) (year 2010), (c) satellite observations from GEDI (25 March–22 May 2019), (d) difference between Spawn and GEDI, (e) the historical simulation eHI (year 2010) – Spawn dataset, (f) the historical simulation eHI – Spawn dataset (year 2010). Globally integrated AGB is 533 PgC for eHI, whereas it is 287 PgC for SPAWN and 283 PgC for GEDI.

4.4 Köppen-Geiger Climate Classifications

In Fig. 9, we show the Köppen-Geiger classifications of observed 2001–2009 (ERA5) climate and the simulated climates, following the classification algorithms of Kottek et al. (2006), which produce 31 climate classes, listed in Appendix Table B1. For the tropics, we implemented a slight modification. Where Kottek et al. (2006) defined a 6-month summer, we used 4-month seasons for tropical climates As and Aw (northern hemisphere, summer = June–September, winter = November–February). This modification identifies the rain shadow areas in East Africa, Sri Lanka, and maybe a little too much of northeastern Brazil.

The observed climatology was classified at 0.5° × 0.5°, because the coarser ModelE grid cannot capture a less common climate type, Dsd. We note also that because the Köppen-Geiger scheme uses only precipitation and temperature data, it cannot capture wetland areas like those in the southern end of the Amazon Basin, which are fed by stream- and riverflow; these areas classified as Aw (pink, equatorial savannah) are in reality Af or Am (dark red, equatorial rainforest, or bright red, equatorial monsoon). Compared to observations, the PI and historical runs show that ModelE2 has a warm/dry bias in the Amazon and Congo basins, producing the climates more suited to savannah than tropical forest. ModelE2 also has a cold bias in the northern high latitudes; where observed climate shows Dfb and Dfc (snow fully humid), ModelE2 simulates EF and ET (polar frost and polar tundra). ModelE simulates the Western U.S. as moister (Cs, warm temperature with dry summer) than observed (Bw and Bw, arid) and does not produce the northeast Asia Dwc climate (subarctic/boreal taiga/boreal forest). ModelE simulates Central and Western Australia as BSh (arid steppe hot), more conducive to shrubs and grasses than observed BWh (arid desert hot).

The simulated climate change between the PI and the end of the historical period shows a slight migration of high latitude climate type Dfc (bright purple, boreal forest) northward, with no changes in the mid-latitudes and tropics. This northward migration continues with the SSP future scenarios, along with encroachment of lower elevation classes to higher elevation around the Tibetan Plateau, and moistening of southern Chile.

Figure 9 Köppen-Geiger climate classifications of (a) observed 2001–2009 climatology (ERA5), (b) the pre-industrial and (c) historical runs. Colorscheme is from Kottek et al. (2006). Observed climate from ERA5 at 0.5° × 0.5° resolution.

5 Synopsis and Conclusions

In this paper we described the development, skill and characteristics of the NASA GISS Earth system model coupled to emissions of CO₂. With respect to the CMIP6 version of the emissions driven model (CC: Ito et al., 2020), GISSE2.1-CC2 has the same physical climate formulation with non-interactive aerosols which are prescribed from fully interactive runs. Additionally, GISSE2.1-CC2 includes interactive oxygen and improved iron distributions.

The main advantage of CC2 vs CC is that the CC2 model has been run to equilibrium more consistently than in CC (CMIP6). The land carbon physics are the same as in CMIP6, but the improved LUHv2 historical crop cover dataset has been used to replace older estimates used in CMIP6.

The model has improved pre-industrial carbon fields and lower carbon drift than the previous version of the model (Ito et al., 2020). It also includes a better crops and land use forcing product (GCB-LUH2 2020) and various improvements in the ocean carbon cycle (iron ligands, oxygen coupled to ocean biogeochemistry, improved nitrogen fixation). However, ocean primary production remains as low as in the CMIP6 version of the model. While some of the model biases (Latto and Romanou, 2018; Marshall et al., 2017; Romanou et al., 2017) are certainly attributed to biases in the underlying ocean model (e.g. low wind speed, Arctic circulation and hydrography) many of the biases in the carbon cycle stem from biases in different parameterizations (e.g. remineralization and nitrogen fixation), missing denitrification sinks for nitrate, sea ice effects on primary production, and a discrepancy between the iron distributions in the preindustrial control climate and the historical simulations which leads to a strong decline in productivity early in the 20th century simulation (Lerner et al., 2021).

The land carbon model overall predicts global land GPP, NPP, and soil carbon well within the mean of estimates from multi-model intercomparisons (Friedlingstein et al., 2025) of fluxes, and nearly exactly matches inventory compilations for soil carbon (Canadell et al., 2021). Predicted GPP (and NPP) is at the lower end compared to other current DGVMs. The model's global GPP is found to be 127 Pg yr⁻¹ in 2021 while uncertainties estimated by Cheng et al., 2017 (the primary source for IPCC Sixth Assessment Report Working Group, Canadell et al., 2021, and the GCB2024) and 146.1 ± 21.3 Pg yr⁻¹. The model grossly overpredicts vegetation biomass due to prescription of forest heights from early satellite lidar estimates from ICESat/GLAS (Simard et al., 2011), which we plan to replace with new more accurate estimates from GEDI and ICESat-2. There is a bug in the storage of plant labile carbon for grasses, causing excess accumulation in the plants instead of senescing as litterfall. Meanwhile, the model's configuration in biophysics-only mode with only prescribed vegetation structure does not allow for vegetation change through greening, mortality, or competition; however, the prognostic carbon pools of plant labile carbon and soil carbon provide useful diagnostics for vegetation responses to climate changes and CO₂ fertilization, and the assumption of mature forest canopy structure may be closer to reality than the overpredicted forest growth currently from DGVMs in the GCB2024. The model does not represent several LUC processes, notably deforestation, wood harvest, vegetation and peat fires, resulting in missing about 1 PgC yr⁻¹ in LUC-related emissions compared to GCB2024 models over 2009–2019, and that deficit is likely larger over the preindustrial through 20th century. The cumulative land sink since 1850 is still very poorly constrained, so it is better not to reference other estimates currently, which differ wildly, but rather to scrutinize self-consistency within the model. We conclude that the model's lack of these LUC processes causes a deficit in simulated cumulative fluxes to the atmosphere, relative to the model's other flux components, and may explain ModelE's slightly underpredicted atmospheric CO₂ over most of the historical period.

Overall, the ModelE2.1-CC2 predicts the atmospheric concentrations of CO₂ consistent with prescribed emissions given by Community Emissions Data System (CEDS) and much the mean inter-model difference of 7 ppm found in Hajima et al. (2025).

The current version (NASA ModelE2.1-CC2) has been included in several model intercomparison studies where its performance can be directly assessed against other state-of-the-art Earth system models. In coordinated experiments with a data-assimilated model of the historical ocean, the Ocean Circulation Inverse Model (OCIM; Holzer et al., 2021; DeVries, 2014) historical and future emissions of CO₂ were estimated from trawling disturbances to the ocean floor (Atwood et al., 2024). From the GISS model both concentration driven as well as emissions driven simulations were used and the results showed that the GISS model (ModelE2.1-CC2) estimates of the fraction of CO₂ that escapes to the atmosphere agree very well with OCIM and OCIM for both concentration and emissions driven runs, with the latter providing significant uncertainty associated with internal climate variability of CO₂ which is not included in concentration driven simulations.

In a 13-model intercomparison study and an ambitious effort to estimate properly the remaining carbon budgets (Silvy et al., 2024) till 1.5 and 2° warming levels, was shown to be consistently within the model spread, except perhaps the land carbon sink, which is however very poorly constrained by observations.

Similarly, in Sanderson et al. (2025) which investigated 10 comprehensive climate models' and 3 intermediate complexity models' response to positive, zero and negative emissions of CO₂, the GISS ModelE2.1-CC2 response was within the uncertainty range of the ensemble mean and only exhibited the strongest response in the later phase of the carbon dioxide removal experiments due to the weakened AMOC.

Future development in the ocean biogeochemistry will focus on adding diel migration as an improvement to the biological pump parameterization in the model, coupling ocean biogeochemistry with sea ice, the addition of denitrification and closing the nitrogen cycle, adding particulate inorganic carbon to improve alkalinity and the CaCO₃ pump, improvements in the underwater light distributions and photoadaptation of phytoplankton, and coupling to physical ocean via albedo and changes in stratification.

Future developments for the land model include releasing a geometric optical radiative transfer (GORT) canopy model for global production runs. This canopy radiative transfer model, the Analytical Clumped Two-Stream (ACTS) model (Ni-Meister et al., 2010; Yang et al., 2010) leverages the cohort structure of the Ent TBM analytically to solve for foliage clumping in heterogeneous canopies in two-stream scheme. This approach simulates vertical foliage profiles (VFP) as observed by satellite lidar, which will enable ingestion of these profiles from GEDI and ICESat-2. The combination of foliage clumping prediction and VFPs then allows prediction of the effects of zenith angle and gaps on canopy spectral albedo, vertical light profiles for light competition, and transmittance of radiation to the ground for predicting soil temperature and snowmelt. We have also completed an analysis (Grant et al., 2025) of the Tallo dataset of tree allometry (Jucker et al., 2022), in which we identified significant distinctions between plant functional types; these new allometric parameterizations will be introduced into the Ent TBM along with more recent satellite lidar observations of forest canopy heights. Combining these canopy physics and updated canopy structure and allometry will allow reduction in uncertainty in vegetation biomass; it will provide a fusion between vegetation demography for simulating plant competition for light with emerging satellite lidar constraints on vertical canopy structure. This framework will allow addressing persistent issues in vegetation-atmosphere coupling, including the snow-vegetation albedo feedback and heterogeneity of open tree-grass systems. In addition, fire-enabled dynamics coupling both the land carbon cycle and chemistry are under development, with fuel loads from Ent exchanged with an updated version of ModelE's fire chemical emissions model, pyrE (Mezuman et al., 2020), adding CO₂ as an emission. These new fire dynamics include natural fires and peat fires, and mortality of vegetation will be enabled in tiers of ecological dynamics, from simple stand thinning to full community dynamics. The patch dynamics implemented for fires will later enable implementation of deforestation.

With regards to the global carbon-climate interactivity future versions will include interactive dust utilizing modelE's state-of-the-art atmospheric composition model (MATRIX, Bauer et al., 2008) and emissions of CO₂ from deforestation from wildfires using GISS pyrE (Mezuman et al., 2020) as well as land use change, peats and wetlands. These future model development efforts will aim to develop more comprehensive simulations of the Earth System subject to assessments against pioneering NASA observational products such as from PACE and SWOT, as well as the hyperspectral Earth Surface Mineral Dust Source Investigation (EMIT) (Ochoa et al., 2025).

Appendix A

A1 Carbon Cycle Model Formulation

A1.1 NOBMg

Model Equations

Phytoplankton species (diatoms, chlorophytes, cyanobacteria, coccolithophores):

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}{\displaystyle \frac{\partial P_i}{\partial t}}& =\mathrm{\nabla }\left(K\mathrm{\nabla }{P}_i\right)-\mathbf{\nabla }\cdot \mathit{u}{P}_i-\mathrm{\nabla }{w}_{\mathrm{s},i}{P}_i\\ & +\left[{\mathit{\mu }}_i-\left(\mathit{\delta }+\mathrm{\Omega }\right)\right]P_i-\mathit{\gamma }H-\mathit{\kappa }{P}_i\end{array}\end{array}$$

The terms on the right-hand side represent mixing, advection, sinking, growth via nutrient uptake, exudation, respiration of DOC, zooplankton grazing, and senescence. The subscript i indicates one of the phytoplankton functional groups. Parameterization of the sinking rates (w_s,I), specific growth rates (μ_i), and the grazing rate (γ) are described below.  The constants δ, Ω, and κ are the exudation, respiration, grazing, and senescence rates, respectively, as defined in Gregg and Casey (2007).

The total specific growth rate (d⁻¹) is given by:

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}{\mathit{\mu }}_i& ={\mathit{\mu }}_{\mathrm{m},i}\left[\mathit{\omega }{\left(E_{\mathrm{T}}\right)}_i,\mathit{\omega }{\left(N_{\mathrm{T}}\right)}_i,\mathit{\omega }{\left(\mathrm{Si}\right)}_i,\mathit{\omega }{\left(\mathrm{Fe}\right)}_i\right]R{G}_i\\ & +{\mathit{\mu }}_{i,\mathrm{cyan}}min\left[\mathit{\omega }{\left(E_{\mathrm{T}}\right)}_i,\mathit{\omega }{\left(\mathrm{Fe}\right)}_{i,\mathrm{Nfx}}\right]R{G}_i{N}_{\mathrm{in},i}\end{array}\end{array}$$

where μ_m,i is the maximum growth rate at 20 °C (see Table A5). Maximum growth rates were adjusted as described in the section “Tuning”. ω(E_T) is the growth rate fraction that is function only of the total irradiance (µmol quanta m⁻² s⁻¹) and is given below

$$\begin{array}{}\text{(A3)}& \mathit{\omega }{\left(E_{\mathrm{T}}\right)}_i={\displaystyle \frac{E_{\mathrm{T}}}{E_{\mathrm{T}}+k_{E,i}}}\end{array}$$

where E_T is the total irradiance, and k_E is the irradiance at which ${\mathit{\mu }}_i=\mathrm{0.5}{\mathit{\mu }}_{\mathrm{m},\mathrm{i}}$ and $k_{E,i}=\mathrm{0.5}{I}_{k,i}$, where I_k is the light saturation parameter.

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}\mathit{\{}{I}_{k,i}& =I_{k,i,\mathrm{hi}},E_{\mathrm{T}}>\mathrm{250}{I}_{k,i}=I_{k,i,\mathrm{low}},E_{\mathrm{T}}\le \mathrm{25}{I}_{k,i},\\ & =\left(\mathrm{1}-{\displaystyle \frac{E_{\mathrm{T}}-\left(\mathrm{0.5}\left(\mathrm{250}-\mathrm{25}\right)+\mathrm{25}\right)}{\mathrm{250}-\left(\mathrm{0.5}\left(\mathrm{250}-\mathrm{25}\right)+\mathrm{25}\right)}}\right)I_{k,i,\mathrm{med}}\\ & +\left({\displaystyle \frac{E_{\mathrm{T}}-\left(\mathrm{0.5}\left(\mathrm{250}-\mathrm{25}\right)+\mathrm{25}\right)}{\mathrm{250}-\left(\mathrm{0.5}\left(\mathrm{250}-\mathrm{25}\right)+\mathrm{25}\right)}}\right)I_{k,i,\mathrm{hi}},\\ & E_{\mathrm{T}}\ge \mathrm{0.5}\left(\mathrm{250}-\mathrm{25}\right)+\mathrm{25}{I}_{k,i},\\ & =\left(\mathrm{1}-{\displaystyle \frac{E_{\mathrm{T}}-\mathrm{25}}{\left(\mathrm{0.5}\left(\mathrm{250}-\mathrm{25}\right)+\mathrm{25}\right)-\mathrm{25}}}\right)I_{k,i,\mathrm{low}}\\ & +\left({\displaystyle \frac{E_{\mathrm{T}}-\mathrm{25}}{\left(\mathrm{0.5}\left(\mathrm{250}-\mathrm{25}\right)+\mathrm{25}\right)-\mathrm{25}}}\right)I_{k,i,\mathrm{med}},\\ & E_{\mathrm{T}}<\mathrm{0.5}\left(\mathrm{250}-\mathrm{25}\right)+\mathrm{25}\end{array}\end{array}$$

Here, $I_{k,i,\mathrm{hi},\mathrm{med},\mathrm{low}}$ are light saturation constants for low, medium, and high irradiance categories and are set following Gregg and Casey (2007) (see Table A5).

Table A1 State variables.

Table A2 Atmospheric variables needed by NOBMg.

Table A3 Sea-ice variables needed by NOBMg.

Table A4 Physical parameters.

Table A5 Phytoplankton parameters.

Table A6 Zooplankton parameters.

The nutrient dependent growth rates are:

$$\begin{array}{}\text{(A5)}& {\displaystyle }\mathit{\omega }{\left({\mathrm{NH}}_{\mathrm{4}}\right)}_i={\displaystyle \frac{A}{A+k_{\mathrm{N},i}}}\text{(A6)}& {\displaystyle }\mathit{\omega }{\left({\mathrm{NO}}_{\mathrm{3}}\right)}_i={\displaystyle \frac{N}{N+k_{\mathrm{N},i}}}\text{(A7)}& {\displaystyle }f{\left({\mathrm{NO}}_{\mathrm{3}}\right)}_i=min[{\left({\mathrm{NO}}_{\mathrm{3}}\right)}_i,\mathrm{1}-\mathit{\omega }{\left({\mathrm{NH}}_{\mathrm{4}}\right)}_i]\text{(A8)}& {\displaystyle }\mathit{\omega }{\left(N_{\mathrm{T}}\right)}_i=\mathit{\omega }{\left({\mathrm{NH}}_{\mathrm{4}}\right)}_i+f{\left({\mathrm{NO}}_{\mathrm{3}}\right)}_i\text{(A9)}& {\displaystyle }\mathit{\omega }{\left(\mathrm{Fe}\right)}_i={\displaystyle \frac{I}{I+k_{\mathrm{Fe},i}}}\text{(A10)}& {\displaystyle }\mathit{\omega }{\left(\mathrm{Si}\right)}_i={\displaystyle \frac{S}{S+k_{\mathrm{Si},i}}}\end{array}$$

where $k_{\mathrm{N},\mathrm{Fe},\mathrm{Si},\mathrm{I}}$ are half-saturation constants for nitrogen (ammonium and nitrate), iron, and silica limitation, respectively. Note that while k_N remains set to the value used by Gregg and Casey (2007), k_Fe,Si have been adjusted as described in the section “Tuning”.

The temperature dependence of the growth rate is from Bissinger et al. (2008):

$$\begin{array}{}\text{(A11)}& R=e^{-b(T-\mathrm{20})}\end{array}$$

where the constant b, originally set to 0.0631 to fit phytoplankton data from the Liverpool Plankton Database, has been reduced to 0.031 to improve the regional distribution of net primary productivity (in particular, increase NPP in the subtropics). This also yields a Q₁₀ of 1.36, which is closer to the Q₁₀ obtained from a more recent compilation of phytoplankton community growth rates (1.47; Sherman et al., 2016) than the Q₁₀ value found by Bissinger et al. (2008)

Specifically for cyanobacteria there is an additional temperature dependence which reduces their growth rate in cold waters (colder than 15 °C, Agawin et al., 1998, 2000).

$$\begin{array}{}\text{(A12)}& \begin{array}{rl}{G}_i& =\mathit{\{}\mathrm{0.0294}T+\mathrm{0.558},i=\mathrm{cyanobacteria}\mathrm{1},\\ & i=\mathrm{other}\mathrm{PFTs}\end{array}\end{array}$$

Our treatment of nitrogen fixation, which amends the specific growth rate for cyanobacteria (second term in Eq. A2), has been revised from that of Gregg and Casey (2007), with two major differences. In Gregg and Casey (2007), the limitation of growth supported by nitrogen fixation included the same iron limitation (i.e., same half-saturation constant for iron-limitation of cyanobacteria) as non-nitrogen fixation supported growth. This has now been, i.e.:

$$\begin{array}{}\text{(A13)}& \mathit{\omega }{\left(\mathrm{Fe}\right)}_{i,\mathrm{Nfx}}={\displaystyle \frac{I}{I+k_{\mathrm{Fe},\mathrm{Nfix}}}}\end{array}$$

where the half-saturation for iron limitation for nitrogen fixation (k_Fe,Nfix) was tuned to obtain a global nitrogen fixation rate that is within observational uncertainity, as described in the section “Tuning”.

Additionally, we have replaced the cyanobacterial biomass inhibition of nitrogen fixation with inhibition by the presence of nitrate and ammonium, following findings from laboratory studies using Trichodesmium cultures Holl and Montoya (2005). The inhibition factor (N_in,i) is set to 0 for all non-cyanobacteria remaining phytoplankton groups as these groups do not fix nitrogen and hence the term representing growth due to nitrogen fixation (second term in Eq. A2) is set to zero.

$$\begin{array}{}\text{(A14)}& \begin{array}{rl}{N}_{\mathrm{in},i}& =\mathit{\{}{\displaystyle \frac{r_{k,\mathrm{ni}}}{N+A+r_{k,\mathrm{ni}}}},i=\mathrm{cyanobacteria}\mathrm{0},\\ & i=\text{other PFTs}\end{array}\end{array}$$

where r_k,ni is the inhibition constant for nitrogen and ammonium from Holl and Monotya (2005).

Phytoplankton sinking is treated as additional vertical advection term (following Gregg and Casey, 2007) but with additional exponential dependence on diatom concentration for the diatoms and limiting term for coccolithophores which depends on the coccolithophore growth rate. For the other species is a function of viscosity (++++)

$$\begin{array}{}\text{(A15)}& {\displaystyle }\begin{array}{rl}{w}_{\mathrm{s},i}& =\mathit{\{}\left(\mathrm{0.451}+\mathrm{0.0178}T\right)a_{\mathrm{diat}}{e}^{P_i{b}_{\mathrm{diat}}},\\ & i=\mathrm{diatoms}{w}_{\mathrm{sch}}\left(\mathrm{0.451}+\mathrm{0.0178}T\right),\\ & i=\mathrm{coccolithophores}\left(\mathrm{0.451}+\mathrm{0.0178}T\right)w_{\mathrm{s},i},\\ & i=\text{other PFTs},\end{array}\text{(A16)}& {\displaystyle }{w}_{\mathrm{sch}}=min(max\left(\mathrm{0.752}\mathrm{gcmax}+\mathrm{0.225},\mathrm{0.3}\right),\mathrm{1.4})\text{(A17)}& {\displaystyle }\mathrm{gcmax}=max({\mathit{\mu }}_{\mathrm{cocc}},\mathrm{gcmax})\end{array}$$

Note that gcmax, which is initialized to the value of the coccolithpohre growth rate, is updated each time step to either match the coccolithophore growth rate (if the growth rate is larger), or retain it's current value (if the growth rate larger).

Zooplankton grazing increases with temperature and the total biomass of phytoplankton.

$$\begin{array}{}\text{(A18)}& {\displaystyle }\mathit{\gamma }={\mathit{\gamma }}_{\mathrm{m}}{R}_{\mathrm{H}}(\mathrm{1}-e^{-\mathrm{\Lambda }\sum _i{P}_i})\text{(A19)}& {\displaystyle }{R}_{\mathrm{H}}=\mathrm{0.06}{e}^{\mathrm{0.1}T}+\mathrm{0.70}\end{array}$$

where Λ and γ_m are the Ivlev constant and maximum grazing rate at 20 °C, respectively (Gregg and Casey, 2007).

Nutrients (nitrate, ammonium, silicate, iron)

The governing equation for nitrate is:

$$\begin{array}{}\text{(A20)}& \begin{array}{rl}{\displaystyle \frac{\partial N}{\partial t}}& =\mathrm{\nabla }\left(K\mathrm{\nabla }N\right)-\mathbf{\nabla }\cdot \mathit{u}N-b_{\mathrm{N}}\left[\sum _i{f}_{{\mathrm{NO}}_{\mathrm{3}}}{\mathit{\mu }}_i{P}_i\right]\\ & +R{\mathit{\alpha }}_N{\displaystyle \frac{D_{\mathrm{c}}}{\mathrm{C}:\mathrm{N}}}+{\mathit{\lambda }}_{\mathrm{D}}{\displaystyle \frac{D_{\mathrm{c}}}{\mathrm{C}:\mathrm{N}}}\end{array}\end{array}$$

where the terms represent, from left to right, mixing, advection, uptake of nitrate via phytoplankton growth, detrital organic carbon $/$ nitrogen remineralization, and detrital organic carbon $/$ nitrogen breakdown to DOC. Here, b_N is the nitrogen to chlorophyll ratio, C:N is the stoichiometric ratio of carbon to nitrogen in phytoplankton, and λ_D is the rate of detrital breakdown to DOC. Each of these parameters are set to the values in Gregg and Casey (2007). Note that R follows Eq. (A11) except that we retain b=0.0631 (the original value from Bissinger et al., 2008). α_N, the remineralization rate of organic carbon $/$ nitrogen, follows:

$$\begin{array}{}\text{(A21)}& a_{\mathrm{N}}=r_{\mathrm{aer}}H(\mathrm{1}-e^{-{\mathrm{O}}_{\mathrm{2}}/k_{{\mathrm{O}}_{\mathrm{2}}}})\end{array}$$

where r_aer is the maximum rate of remineralization at 20 °C, H is the step function that equals 1 if O₂ > 2 µmol kg⁻¹, and 0 otherwise, and $k_{{\mathrm{O}}_{\mathrm{2}}}$ is the e-folding length scale for the oxygen-dependence of remineralization. The latter dependence was included following a study by Keil et al. (2016) that found the remineralization rate to decline exponentially with decreasing oxygen concentration (see also Lerner et al., 2024 for a discussion on how the value of $k_{{\mathrm{O}}_{\mathrm{2}}}$ was assigned).

$f_{{\mathrm{NO}}_{\mathrm{3}}}$, the fraction of nitrogen assimilated into phytoplankton biomass that is obtained from nitrate, follows:

$$\begin{array}{}\text{(A22)}& f_{{\mathrm{NO}}_{\mathrm{3}}}=N/(A+N)\end{array}$$

The governing equation for ammonium is:

$$\begin{array}{}\text{(A23)}& \begin{array}{rl}{\displaystyle \frac{\partial A}{\partial t}}& =\mathrm{\nabla }\left(K\mathrm{\nabla }A\right)-\mathbf{\nabla }\cdot \mathit{u}A-b_{\mathrm{N}}\left[\sum _i{f}_{{\mathrm{NH}}_{\mathrm{4}}}{\mathit{\mu }}_i{P}_i\right]\\ & +b_{\mathrm{N}}\mathit{ϵ}\left[\mathit{\gamma }H+{\mathit{\eta }}_{\mathrm{2}}{H}^{\mathrm{2}}\right]\end{array}\end{array}$$

where the terms represent, from left to right, mixing, advection, uptake of ammonium via phytoplankton growth, and zooplankton excretion. Here, ϵ is the zooplankton excretion parameter, and $f_{{\mathrm{NH}}_{\mathrm{4}}}$, the fraction of nitrogen assimilated into phytoplankton biomass that is obtained from ammonium, follows:

$$\begin{array}{}\text{(A24)}& f_{{\mathrm{NO}}_{\mathrm{3}}}=A/(A+N)\end{array}$$

The governing equation for dissolved silica is:

$$\begin{array}{}\text{(A25)}& {\displaystyle \frac{\partial S}{\partial t}}=\mathrm{\nabla }\left(K\mathrm{\nabla }S\right)-\mathbf{\nabla }\cdot \mathit{u}S-b_{\mathrm{S}}{\mathit{\mu }}_{\mathrm{1}}{P}_{\mathrm{1}}+R{\mathit{\alpha }}_{\mathrm{S}}{D}_{\mathrm{S}}\end{array}$$

where the terms represent, from left to right, mixing, advection, uptake of silica via phytoplankton growth, dissolution of detrital silica. Here, b_S is the silica to chlorophyll ratio, and α_s is the detrital silica dissolution rate b_S is set to the value in Gregg and Casey (2007), while α_s was adjusted as described in the “Tuning” section. For dissolved iron, the governing equation is:

$$\begin{array}{}\text{(A26)}& \begin{array}{rl}{\displaystyle \frac{\partial I}{\partial t}}& =\mathrm{\nabla }\left(K\mathrm{\nabla }I\right)-\mathbf{\nabla }\cdot \mathit{u}I-b_{\mathrm{I}}\left[\sum _i{\mathit{\mu }}_i{P}_i\right]\\ & +b_{\mathrm{I}}\mathit{ϵ}\left[\mathit{\gamma }H+{\mathit{\eta }}_{\mathrm{2}}{H}^{\mathrm{2}}\right]+R{\mathit{\alpha }}_{\mathrm{I}}{D}_{\mathrm{I}}+{\displaystyle \frac{A_{\mathrm{Fe}}}{L}}\\ & -{\mathit{\theta }}_{\mathrm{1}}I+{\mathit{\theta }}_{\mathrm{2}}{I}^{\prime }\end{array}\end{array}$$

where the terms represent, from left to right, mixing, advection, uptake of iron via phytoplankton growth, zooplankton excretion, remineralization of detrital iron, atmospheric depiction, and the scavenging of ligand-bound and free iron. Here, b_I is the iron to chlorophyll ratio, which is equal to the value used by Gregg and Casey (2007). α_I, θ₁, and θ₂ are the remineralization rate, ligand-bound and free iron scavenging rates, respectively, whose values are adjusted as described in the “Tuning” section. A_Fe is the atmospheric deposition of iron bound to dust, which is prescribed as detailed in the “Forcings” section, while L is the thickness of the surface layer. I^′ is solved for at each timestep from the system of equations that describe the complexation of iron with a single ligand class:

$$\begin{array}{}\text{(A27)}& {\displaystyle }\left[I\right]=\left[IL\right]+\left[I^{\prime }\right]\text{(A28)}& {\displaystyle }\left[L_{\mathrm{tot}}\right]=\left[IL\right]+\left[L^{\prime }\right]\text{(A29)}& {\displaystyle }{K}_{\mathrm{FeL}}=\left[IL\right]/\left(\left[I^{\prime }\right]\right[L^{\prime }\left]\right])\end{array}$$

L_tot is from the prescribed ligand pool described in “Forcings.” The stability constant K_FeL is set to a constant value of 10¹¹ nmol L⁻¹, taken from the average value from the compilation of Caprara et al. (2016).

Herbivores

The governing equations for herbivores is

$$\begin{array}{}\text{(A30)}& \begin{array}{rl}{\displaystyle \frac{\partial H}{\partial t}}& =\mathrm{\nabla }\left(K\mathrm{\nabla }H\right)-\mathbf{\nabla }\cdot \mathit{u}H+\left(\mathrm{1}-\mathit{ϵ}\right)\mathit{\gamma }H-{\mathit{\eta }}_{\mathrm{1}}\\ & H-{\mathit{\eta }}_{\mathrm{2}}{H}^{\mathrm{2}}-\mathit{\zeta }H-\mathrm{\Theta }H\end{array}\end{array}$$

where the terms represent, from left to right, mixing, advection, grazing, zooplankton linear and quadratic mortality, zooplankton excretion of DOC, and zooplankton respiration. Here, η₁ and η₂ are the linear and quadratic mortality rates, and Θ is the herbivore respiration rate. These constants are the same as those used by Gregg and Casey (2007). ζ is the excretion rate of DOC by herbivores, and follows:

$$\begin{array}{}\text{(A31)}& \mathit{\zeta }=r_{\mathrm{H}}{\displaystyle \frac{H}{H_{\mathrm{0}}+H}}R\end{array}$$

where r_H is the excretion rate of herbivores at 20 °C, and H₀ is the half-saturation constant of herbivore excretion, both of which are set equal to the values from Gregg and Casey (2007). R is the temperature-dependence of herbivore excretion, which is equal to R for detrital nitrogen/carbon remineralization (Eq. A20).

Particulate organic components (carbon/nitrate detritus, silica detritus, iron detritus)

The governing equation for detrital carbon $/$ nitrogen is:

$$\begin{array}{}\text{(A32)}& \begin{array}{rl}{\displaystyle \frac{\partial D_{\mathrm{C}}}{\partial t}}& =\mathrm{\nabla }\left(K\mathrm{\nabla }{D}_{\mathrm{C}}\right)-\mathbf{\nabla }\cdot \mathit{u}{D}_{\mathrm{C}}-\mathrm{\nabla }\cdot w_{\mathrm{d},\mathrm{C}}{D}_{\mathrm{C}}\\ & -R{\mathit{\alpha }}_{\mathrm{C}}{D}_{\mathrm{C}}+\mathrm{\Phi }\left[\mathit{\kappa }\sum _i{P}_i+{\mathit{\eta }}_{\mathrm{1}}H\right]\\ & +\mathrm{\Phi }\left(\mathrm{1}-\mathit{ϵ}\right){\mathit{\eta }}_{\mathrm{2}}{H}^{\mathrm{2}}-{\mathit{\lambda }}_{\mathrm{D}}{D}_{\mathrm{C}}\end{array}\end{array}$$

where the terms represent, from left to right, mixing, advection, sinking, remineralization of detrital iron, phytoplankton senescence, mortality of zooplankton, and breakdown of detrital carbon to DOC. Φ is the carbon to chlorophyll ratio, which is a constant value equal to the median (50 g C g chl⁻¹) of the three values used to represent three photoadaptation states in Gregg and Casey (2007). α_C is the remineralization rate of organic matter, defined as:

$$\begin{array}{}\text{(A33)}& \begin{array}{rl}{a}_{\mathrm{C}}& =\mathit{\{}{r}_{\mathrm{aer}}\left(\mathrm{1}-e^{-{\mathrm{O}}_{\mathrm{2}}/k_{{\mathrm{O}}_{\mathrm{2}}}}\right),{\mathrm{O}}_{\mathrm{2}}\ge \mathrm{2}\mathrm{\mu }\mathrm{mol}{\mathrm{kg}}^{-\mathrm{1}}{r}_{\mathrm{den}},\\ & {\mathrm{O}}_{\mathrm{2}}<\mathrm{2}\mathrm{\mu }\mathrm{mol}{\mathrm{kg}}^{-\mathrm{1}}\end{array}\end{array}$$

where r_aer is as defined in Eq. (A21), and r_den is rate of remineralization at 20 °C under anaerobic conditions, where oxygen is no longer used as the terminal electron acceptor.

For detrital silica, the governing equation is:

$$\begin{array}{}\text{(A34)}& \begin{array}{rl}{\displaystyle \frac{\partial D_{\mathrm{S}}}{\partial t}}& =\mathrm{\nabla }\left(K\mathrm{\nabla }{D}_{\mathrm{S}}\right)-\mathbf{\nabla }\cdot \mathit{u}{D}_{\mathrm{S}}-\mathrm{\nabla }\cdot w_{\mathrm{d},\mathrm{S}}{D}_{\mathrm{S}}\\ & -R{\mathit{\alpha }}_{\mathrm{S}}{D}_{\mathrm{S}}+b_{\mathrm{S}}\left[\mathit{\kappa }{P}_{\mathrm{1}}+\mathit{\gamma }H\right]\end{array}\end{array}$$

where the terms represent, from left to right, mixing, advection, sinking, dissolution of detrital silica, phytoplankton senescence and herbivore grazing.

$$\begin{array}{}\text{(A35)}& \begin{array}{rl}{\displaystyle \frac{\partial D_{\mathrm{I}}}{\partial t}}& =\mathrm{\nabla }\left(K\mathrm{\nabla }{D}_{\mathrm{I}}\right)-\mathbf{\nabla }\cdot \mathit{u}{D}_{\mathrm{I}}-\mathrm{\nabla }\cdot w_{\mathrm{d},\mathrm{I}}{D}_{\mathrm{I}}\\ & -R{\mathit{\alpha }}_{\mathrm{I}}{D}_{\mathrm{I}}+b_{\mathrm{I}}\left[\mathit{\kappa }\sum _i{P}_i+{\mathit{\eta }}_{\mathrm{1}}H\right]\\ & +b_{\mathrm{I}}\left(\mathrm{1}-\mathit{ϵ}\right){\mathit{\eta }}_{\mathrm{2}}{H}^{\mathrm{2}}+{\mathit{\theta }}_{\mathrm{1}}I+{\mathit{\theta }}_{\mathrm{2}}{I}^{\prime }\end{array}\end{array}$$

where the terms represent, from left to right, mixing, advection, sinking, remineralization of detrital iron, phytoplankton senescence, linear and quadratic herbivore mortality, and the scavenging of ligand-bound and free-iron.

For Eqs. (A32)–(A34), the detrital sinking rates follow and exponential dependence, similar to that for diatoms:

$$\begin{array}{}\text{(A36)}& w_{\mathrm{d},\mathrm{C},\mathrm{S},\mathrm{I}}=(\mathrm{0.451}+\mathrm{0.0178}T)a_{\mathrm{C},\mathrm{I},\mathrm{S}}{e}^{P_i{b}_{\mathrm{C},\mathrm{I},\mathrm{S}}}\end{array}$$

where the coefficients $a_{\mathrm{C},\mathrm{I},\mathrm{S}}$ and $b_{\mathrm{C},\mathrm{I},\mathrm{S}}$ have been previously described in Ito et al. (2020).

Table A7 Nutrient parameters.

Table A8 Detrital Parameters.

Table A9 Carbon Parameters.

Table A10 Alkalinity parameters.

Table A11 Oxygen parameters.

Table A12 Iron parameters.

Carbon components (dissolved organic carbon, dissolved inorganic carbon, alkalinity)

The governing equation for dissolved organic carbon (DOC) is

$$\begin{array}{}\text{(A37)}& \begin{array}{rl}{\displaystyle \frac{\partial \mathrm{DOC}}{\partial t}}& =\mathrm{\nabla }\left(K\mathrm{\nabla }\mathrm{DOC}\right)-\mathbf{\nabla }\cdot \mathit{u}\mathrm{DOC}\\ & +\mathrm{\Phi }\mathit{\delta }\sum _i{\mathit{\mu }}_i{P}_i+\mathrm{\Phi }\mathit{\zeta }H+{\mathit{\lambda }}_{\mathrm{D}}{D}_{\mathrm{C}}-\mathit{\phi }\mathrm{DOC}\end{array}\end{array}$$

where the terms represent, from left to right, mixing, advection, phytoplankton exudation, zooplankton excretion, breakdown of detrital carbon, and DOC remineralization. The DOC remineralization rate follows:

$$\begin{array}{}\text{(A38)}& \mathit{\phi }=H{\displaystyle \frac{N}{N+k_{\mathrm{1}}}}{\displaystyle \frac{\mathrm{DOC}}{\mathrm{DOC}+k_{\mathrm{2}}}}\mathrm{DOCR}\end{array}$$

where H is the same step function used in Eq. (A21). λ_DOC is the maximum remineralization rate at 20 °C, and k₁ and k₂ are half-saturation constants for limitation of DOC remineralization by nitrogen and DOC, respectively. The constants are the same as those used by Gregg and Casey (2007). The temperature-dependent factor R follows the same parameterization as that used for the remineralization of detrital carbon $/$ nitrogen (Eq. A20).

The governing equation for dissolved inorganic carbon (DIC) is:

$$\begin{array}{}\text{(A39)}& \begin{array}{rl}{\displaystyle \frac{\partial \mathrm{DIC}}{\partial t}}& =\mathrm{\nabla }\left(K\mathrm{\nabla }\mathrm{DIC}\right)-\mathbf{\nabla }\cdot \mathit{u}\mathrm{DIC}-\mathrm{\Phi }\sum _i{\mathit{\mu }}_i{P}_i\\ & +\mathrm{\Phi }\mathrm{\Omega }\sum _i{\mathit{\mu }}_i{P}_i+\mathrm{\Phi }\mathrm{\Theta }H+\mathit{\phi }\mathrm{DOC}\\ & +R{\mathit{\alpha }}_{\mathrm{C}}{\displaystyle \frac{D_{\mathrm{C}}}{\mathrm{C}:\mathrm{N}}}+F_{\mathrm{AO}{\mathrm{CO}}_{\mathrm{2}}}\end{array}\end{array}$$

where the terms represent, from left to right, mixing, advection, uptake of DIC by phytoplankton, phytoplankton respiration, zooplankton respiration, DOC remineralization, detrital carbon remineralization, and the air-sea exchange of CO₂. The last term follows the OMIP-BGC protocol for CMIP6 (Orr et al., 2017);

$$\begin{array}{}\text{(A40)}& F_{\mathrm{AO}{\mathrm{CO}}_{\mathrm{2}}}=k_{w,{\mathrm{CO}}_{\mathrm{2}}}\mathrm{ff}\left({\mathrm{xCO}}_{\mathrm{2}}-{\mathrm{CO}}_{\mathrm{2},\mathrm{s}}\right)\mathrm{fsice}/L\end{array}$$

Where $k_{w,{\mathrm{CO}}_{\mathrm{2}}}$ is the piston velocity of air-sea CO₂ gas exchange, ff is the solubility of CO₂, xCO₂ is the dry-air fraction of atmospheric CO₂ at the surface of the ocean, CO_2,s is the partial pressure of CO_2,s at the surface of the ocean, and fsice is the fraction of the surface ocean grid cell covered by sea ice. The piston velocity is computed following

$$\begin{array}{}\text{(A41)}& k_{w,{\mathrm{CO}}_{\mathrm{2}}}=a_{\mathrm{wan}}(s{c}_{{\mathrm{CO}}_{\mathrm{2}}}/\mathrm{660}{)}^{-\mathrm{0.5}}{w_{\mathrm{s}}}^{\mathrm{2}}\end{array}$$

Where a_wan is the Wanninkhof coefficient (Wanninkhof 2014, Orr et al., 2017), w_s is the surface wind speed computed by the model, and $s{c}_{{\mathrm{CO}}_{\mathrm{2}}}$ is the Schmidt number that follows:

$$\begin{array}{}\text{(A42)}& \begin{array}{rl}s{c}_{{\mathrm{CO}}_{\mathrm{2}}}& =\mathrm{2116.8}-\mathrm{135.6}\cdot T+\mathrm{5.2122}{T}^{\mathrm{2}}\\ & -\mathrm{0.10939}{T}^{\mathrm{3}}+\mathrm{0.00094743}{T}^{\mathrm{4}}\end{array}\end{array}$$

While the solubility of CO₂ is computed following

$$\begin{array}{}\text{(A43)}& \mathrm{ff}=e^{-A+\frac{B}{{\mathrm{tk}}_{\mathrm{100}}}+C\log \left({\mathrm{tk}}_{\mathrm{100}}\right)-D{{\mathrm{tk}}_{\mathrm{100}}}^{\mathrm{2}}+S\left(E-F{\mathrm{tk}}_{\mathrm{100}}+G{{\mathrm{tk}}_{\mathrm{100}}}^{\mathrm{2}}\right)}\end{array}$$

Where tk${}_{\mathrm{100}}=\mathrm{100}/T$, and the constant A–G are defined in Orr et al. (2017) (see also Table A9).

Alkalinity

The governing equation for alkalinity is:

$$\begin{array}{}\text{(A44)}& \begin{array}{rl}{\displaystyle \frac{\partial \mathrm{Alk}}{\partial t}}& =\mathrm{\nabla }\left(K\mathrm{\nabla }\mathrm{Alk}\right)-\mathbf{\nabla }\cdot \mathit{u}\mathrm{Alk}+b_{\mathrm{N}}\left[\sum _i{\mathit{\mu }}_i{P}_i\right]\\ & -R{\mathit{\alpha }}_{\mathrm{C}}{\displaystyle \frac{D_{\mathrm{c}}}{\mathrm{C}:\mathrm{N}}}-{\mathit{\lambda }}_{\mathrm{D}}{\displaystyle \frac{D_{\mathrm{c}}}{\mathrm{C}:\mathrm{N}}}+\mathrm{2}{J}_{\mathrm{Ca}}\end{array}\end{array}$$

where the terms represent, from left to right, mixing, advection, the increase in alkalinity due to the conversion of nitrate (a proton donor) to organic N, the decrease in alkalinity due to the conversion of organic N to nitrate during remineralization of detrital organic carbon $/$ nitrogen and DOC, and carbonate dissolution and precipitation. Carbonate dissolution and precipitation (J_Ca) follow the OCMIP-2 protocol:

$$\begin{array}{}\text{(A45)}& \begin{array}{rl}{J}_{\mathrm{Ca}}& =\mathit{\{}-\mathrm{rr}\mathrm{C}:\mathrm{P}\cdot \left(\mathrm{1}-{\mathit{\sigma }}_{\mathrm{Ca}}\right)PP,z\le z_{\mathrm{c}}{\displaystyle \frac{-\mathrm{d}{F}_{\mathrm{Ca}}}{\mathrm{d}z}},\\ & z>z_{\mathrm{c}}\end{array}\end{array}$$

Where rr is the rain ratio of CaCO₃ to organic phosphorous (Yamanaka and Tajika, 1996), C:P is the organic carbon to organic phosphorous ratio, z_c represents the compensation depth which we set to a constant value, and σ_Ca is the fraction of calcium carbonate that dissolved within the compensation depth (Table A10). The change in the downward flux of calcium carbonate with depth uses a fixed exponential profile:

$$\begin{array}{}\text{(A46)}& F_{\mathrm{Ca}}=\mathrm{rr}\mathrm{C}:\mathrm{P}{F}_{\mathrm{C}}{e}^{-\frac{(z-z_{\mathrm{C}})}{d}}\end{array}$$

where d is a fixed exponential length scale for carbonate dissolution, and F_c is the primary productivity integrated from the surface to the compensation depth that is associated with the net production of calcium carbonate within the same depth interval:

$$\begin{array}{}\text{(A47)}& F_{\mathrm{C}}=\left(\mathrm{1}-{\mathit{\sigma }}_{\mathrm{Ca}}\right)\underset{\mathrm{0}}{\overset{z_{\mathrm{C}}}{\int }}\mathrm{PPd}z\end{array}$$

PP is computed simply as the difference between the specific growth rate and the sum of respiration and exudation, integrated over all phytoplankton functional types:

$$\begin{array}{}\text{(A48)}& \mathrm{PP}=\sum _i{\mathit{\mu }}_i-\left(\mathit{\delta }+\mathrm{\Omega }\right)\end{array}$$

Oxygen

The governing equation for dissolved oxygen is:

$$\begin{array}{}\text{(A49)}& \begin{array}{rl}{\displaystyle \frac{\partial {\mathrm{O}}_{\mathrm{2}}}{\partial t}}& =\mathrm{\nabla }\left(K\mathrm{\nabla }{\mathrm{O}}_{\mathrm{2}}\right)-\mathbf{\nabla }\cdot \mathit{u}{\mathrm{O}}_{\mathrm{2}}+b_{\mathrm{O}}\sum _i{\mathit{\mu }}_i{P}_i-b_{\mathrm{O}}\mathrm{\Omega }\\ & \sum _i{\mathit{\mu }}_i{P}_i-b_{\mathrm{O}}\mathrm{\Theta }H-\mathit{\phi }{\displaystyle \frac{\mathrm{DOC}}{\mathrm{C}:\mathrm{O}}}\\ & -R{a}_{\mathrm{C}}{\displaystyle \frac{D_{\mathrm{C}}}{\mathrm{C}:\mathrm{O}}}+F_{{\mathrm{AOO}}_{\mathrm{2}}}\end{array}\end{array}$$

where the terms represent, from left to right, mixing, advection, production of O₂ by phytoplankton, phytoplankton respiration, zooplankton respiration, DOC remineralization, detrital carbon remineralization, and the air-sea exchange of O₂. b_O and C:O are, respectively, the Redfield ratio of O₂ produced or consumed by autotrophic and heterotrophic processes to chlorophyll and carbon in organic matter, respectively. As for CO₂, The air-sea exchange of O₂ follows the protocol described in Orr et al. (2017):

$$\begin{array}{}\text{(A50)}& F_{{\mathrm{AOO}}_{\mathrm{2}}}=k_{w,{\mathrm{O}}_{\mathrm{2}}}\left({\mathrm{O}}_{\mathrm{2},\mathrm{sat}}-{\mathrm{O}}_{\mathrm{2},\mathrm{s}}\right)\mathrm{fsice}/L\end{array}$$

Where $k_{w,{\mathrm{O}}_{\mathrm{2}}}$ is the piston velocity for O₂ exchange between the atmosphere and ocean, O_2,sat is the O₂ concentration at saturation, and O_2,s is the surface concentration of O₂. The piston velocity is computed following:

$$\begin{array}{}\text{(A51)}& k_{w,{\mathrm{O}}_{\mathrm{2}}}=a_{\mathrm{wan}}(s{c}_{{\mathrm{O}}_{\mathrm{2}}}/\mathrm{660}{)}^{-\mathrm{0.5}}{w_{\mathrm{s}}}^{\mathrm{2}}\end{array}$$

where $s{c}_{{\mathrm{O}}_{\mathrm{2}}}$ is the Schmidt number of O₂, which is computed following:

$$\begin{array}{}\text{(A52)}& \begin{array}{rl}s{c}_{{\mathrm{O}}_{\mathrm{2}}}& =\mathrm{1920.4}-\mathrm{136.25}T+\mathrm{4.7353}{T}^{\mathrm{2}}\\ & -\mathrm{0.092307}{T}^{\mathrm{3}}+\mathrm{0.000755}{T}^{\mathrm{4}}\end{array}\end{array}$$

The O₂ saturation concentration is also computed from the protocol of Orr et al. (2017):

$$\begin{array}{}\text{(A53)}& {\mathrm{O}}_{\mathrm{2},\mathrm{sat}}=K_{\mathrm{sol}}(\mathrm{SLP}-{\mathrm{pH}}_{\mathrm{2}}\mathrm{O}){\mathrm{xO}}_{\mathrm{2}}\end{array}$$

where SLP is sea level pressure from the model, xO₂ is the O₂ dry-air air fraction, and the solubility coefficient K_sol is computed following:

$$\begin{array}{}\text{(A54)}& K_{\mathrm{sol}}={\mathrm{O}}_{\mathrm{2},\mathrm{sat}\mathrm{0}}/\left({\mathrm{xO}}_{\mathrm{2}}\right(\mathrm{1}-{\mathrm{pH}}_{\mathrm{2}}\mathrm{0}\left)\right)\end{array}$$

where the partial pressure of H₂O at the surface of the ocean is computed from model temperature and salinity:

$$\begin{array}{}\text{(A55)}& {\displaystyle }\begin{array}{rl}& {\mathrm{pH}}_{\mathrm{2}}\mathrm{O}=\\ & e^{\mathrm{24.4543}-\mathrm{67.4509}{\mathrm{tk}}_{\mathrm{100}}-\mathrm{4.8489}\log \log \left(\frac{\mathrm{1}}{{\mathrm{tk}}_{\mathrm{100}}}\right)-\mathrm{0.000544}S}\end{array}\text{(A56)}& {\displaystyle }\mathrm{ta}=\log \left(\right(\mathrm{298.15}-T)/(\mathrm{Ts}+\mathrm{273.15}\left)\right)\text{(A57)}& {\displaystyle }{\mathrm{tk}}_{\mathrm{100}}=\mathrm{100}/(T+\mathrm{273.15})\end{array}$$

And the saturation concentration of O₂ at standard sea level pressure is computed following the equation of Garcia and Gordon (1992):

$$\begin{array}{}\text{(A58)}& \begin{array}{rl}& {\mathrm{O}}_{\mathrm{2},\mathrm{sat}\mathrm{0}}=\left(\mathrm{1}{e}^{-\mathrm{3}}\right)\cdot \\ & e^{A+B\mathrm{ta}+C{\mathrm{ta}}^{\mathrm{2}}+D{\mathrm{ta}}^{\mathrm{2}}+D{\mathrm{ta}}^{\mathrm{3}}+E{\mathrm{ta}}^{\mathrm{4}}+F{\mathrm{ta}}^{\mathrm{5}}-S\left(G-H\mathrm{ta}+I{\mathrm{ta}}^{\mathrm{2}}-J{\mathrm{ta}}^{\mathrm{3}}\right)-K{S}^{\mathrm{2}}}\end{array}\end{array}$$

A2 Ent Terrestrial Biosphere Model

The complete formulation of the land carbon dynamics is given in Kim et al. (2015). In this study, prognostic growth is not turned on, but the vegetation model is run in “biophysics-only” mode, as described in Sect. 3.

Appendix B

B1 Köppen-Geiger classes

Table B1 Köppen-Geiger biome/climate classes, with color scheme from http://koeppen-geiger.vu-wien.ac.at/ (last access: 21 August 2026).