Evaluation of the University of Victoria Earth System Climate Model version 2 . 10 ( UVic ESCM 2 . 10 )

The University of Victoria Earth system climate model of intermediate complexity has been a useful tool in recent assessments of long-term climate changes, including both paleo-climate modelling and uncertainty assessments of future 2 warming. Since the last official release of the UVic ESCM 2.9, and the two official updates during the last decade, considerable model development has taken place among multiple research groups. The new version 2.10 of the University of Victoria Earth 4 System Climate Model (UVic ESCM) presented here, and to be used in the 6th phase of the coupled model intercomparison project (CMIP6), combines and brings together multiple model developments and new components that have taken place since 6 the last official release of the model. To set the foundation of its use, we describe here the UVic ESCM 2.10 and evaluate results from transient historical simulations against observational data. We find that the UVic ESCM 2.10 is capable of 8 reproducing well changes in historical temperature and carbon fluxes, as well as the spatial distribution of many ocean tracers, including temperature, salinity, phosphate and nitrate. This is connected to a good representation of ocean physical properties. 10 For the moment, there remain biases in ocean oxygen, which will be addressed in the next updates to the model.

For the moment, there remain biases in ocean oxygen, which will be addressed in the next updates to the model.
In the recent special report on limiting global warming to 1.5 °C from the Intergovernmental Panel on Climate Change one of the key uncertainties for the assessment of the remaining global carbon budget was the impact from unrepresented Earth 4 system feedbacks. On the decadal to centennial timescales this specifically refers to the permafrost carbon feedback (Lowe and Bernie, 2018). Quantifying the strength and timing of this permafrost carbon cycle feedback to climate change has been a 6 goal of Earth system modelling in recent years (e.g., Burke et al., 2012;Koven et al., 2011Koven et al., , 2013MacDougall et al., 2012;Schaefer et al., 2011;Schneider Von Deimling et al., 2012;Zhuang et al., 2006).

8
For version 2.10 of the University of Victoria Earth system climate model, we combined version 2.9 with the new marine ecosystem model component as published in Keller et al. (2012), and the soil dynamics and permafrost carbon component as published by Avis et al. (2011) and MacDougall & Knutti (2016). For the 6th Phase of the Coupled Model Intercomparison Project (CMIP6) simulations, the merge of these two components will allow a more comprehensive 12 representation of the carbon cycle in the UVic ESCM, while incorporating the model developments that have been taken place in the context of the UVic ESCM. In addition to the structural changes, we also changed the spin-up protocol to follow CMIP6 14 protocols and applied the newly available CMIP6 forcing.
The objective of the new model development is to have a more realistic representation of carbon and heat fluxes in 16 the UVic ESCM 2.10 in agreement with the available observational data and with current process understanding to be used within the context of the next round of model intercomparison projects for models of intermediate complexity. To set the 18 foundation of its use, we will in the following describe the UVic ESCM 2.10 (Sect. 2.1.), the newly formatted historical CMIP6 forcing that has been and will be used (Sect. 2.2.), explicitly describe changes that have been implemented in the UVic ESCM 20 with respect to the previously published versions (Sect. 2.3.), and then evaluate results from transient historical simulations against observational data (Sect. 3.).

24
The UVic ESCM is a model of intermediate complexity , all model components have a common horizontal resolution of 3.6° longitude and 1.8° latitude and the oceanic component has a vertical resolution of 19 levels, with 26 vertical thickness varying between 50 m near the surface to 500 m in the deep ocean. The Modular Ocean Model Version 2 (MOM2) (Pacanowski, 1995) describes the ocean physics, it is coupled to a thermodynamic-dynamic sea-ice model (Bitz et 28 al., 2001) with elastic visco-plastic rheology (Hunke and Dukowicz, 1997). The atmosphere is represented by a twodimensional atmospheric energy moisture balance model (Fanning and Weaver, 1996). Wind velocities used to calculate 30 advection of atmospheric heat and moisture as well as the air-sea-ice fluxes of surface momentum, heat and water fluxes, are prescribed as monthly climatological wind fields from NCAR/NCEP reanalysis data (Eby et al., 2013). In transient simulations wind anomalies, which are determined from surface pressure anomalies with respect to pre-industrial surface air temperature, are added to the prescribed wind fields . In addition, the terrestrial component represents vegetation 2 dynamics including five different plant functional types (Meissner et al., 2003). Sediment processes are represented using an oxic-only calcium-carbonate model (Archer, 1996). Terrestrial weathering is diagnosed from the net sediment flux during spin-4 up and held fixed at the equilibrium pre-industrial value for transient simulations .
The new version 2.10 of the University of Victoria Earth System Climate Model (UVic ESCM) presented here combines 6 and brings together multiple model developments and new components that have taken place since the last official release of the model in the CMIP5 context. In the following the novel model components are described in detail.

(i) Marine biogeochemical model
Novel compared to the 2009 version of the model, is the ocean biogeochemistry model as published by Keller et al. (2012). It now includes equations describing phytoplankton light limitation and zooplankton grazing, a more realistic zooplankton growth and grazing model, and formulations for an iron limitation scheme to constrain phytoplankton growth. In 12 this context the ocean's mixing scheme was changed from a Bryan-Lewis profile to a scheme for the computation of tidally induced diapycnal mixing over rough topography (Simmons et al., 2004) (see ocean diffusivity profiles in Figure S3). In 14 addition, the air to sea gas parameterization was updated following the ocean carbon-cycle model intercomparison project updates for these numbers (Wanninkhof, 2014), which impacts the carbon exchange between the atmospheric and marine 16 components. Furthermore, we now apply the stoichiometry from Paulmier et al. (2009) to consistently account for the effects of denitrification and nitrogen fixation on alkalinity and oxygen.
(ii) Soil model The terrestrial component has also been updated relative to the latest official release of the UVic ESCM. It now includes 20 a representation of soil freeze-thaw processes resolved in 14 subsurface layers of which the thicknesses exponentially increase with depth: the surface layer having a thickness of 0.1 m and the bottom layer a thickness of 104.4 m, the total thickness of 22 the subsurface layers is 250 m. The top eight layers (to a depth of 10 m) are soil layers; below this are bedrock layers having the thermal characteristics of granitic rock. Moisture undergoes free drainage from the base of the soil layers and the bedrock 24 layers are hydrologically inactive (Avis et al., 2011). In addition, the soil module includes a multi-layer representation of soil carbon (MacDougall et al., 2012). Organic carbon from the litter flux is allocated to soil layers as a decreasing function of 26 depth, and is only added to soil layers with a temperature above 1ºC. If all layers are below this temperature threshold the litter flux is added to the top layer of soil. Soil respiration remains a function of temperature and moisture (Meissner et al., 2003) 28 but is now implemented in each layer individually. Respiration ceases if the soil layer is below 0ºC. Soil carbon is present in the top six layers of the soil column down to a depth of 3.35 m.

30
(iii) Permafrost model A representation of permafrost carbon has also been added to the model. Permafrost carbon is prognostically generated within 32 the model using a diffusion-based scheme meant to approximate the process of cryoturbation (MacDougall & Knutti 2016). In model grid-cells with perennially frozen soil layers soil carbon is diffused proportional to the effective carbon concentration of each soil layer. Effective carbon concentration is carbon concentration divided by porosity and a saturation factor (MacDougall & Knutti 2016). Carbon that is diffused into perennially frozen soil is reclassified as permafrost carbon and is 2 given different properties from regular soil carbon. Permafrost carbon decays with its own constant decay rate and is subject to an `available fraction', which determines the fraction of permafrost carbon that is available to decay. The available fraction 4 slowly increased if permafrost carbon becomes thawed and decreased if permafrost carbon decays. Using this scheme, the model can represent the large fraction of permafrost carbon that is in the passive soil carbon pool, while still allowing the 6 passive pool to eventually decay (MacDougall & Knutti 2016).

8
Anthropogenic forcing from greenhouse gases (GHGs), stratospheric and tropospheric ozone, aerosols and stratospheric water-vapor from methane oxidation are considered. Natural forcing includes solar and volcanic. All data used in the creation of this dataset can be accessed from Input4Mip on the Earth System Grid Federation (ESGF) unless otherwise specified. In the following we will shortly describe how the input data for our simulations with the University of Victoria Earth 12 system climate model (UVic ESCM) was created.
In the standard CMIP6 configuration, the UVic ESCM is forced with CO2 concentration data (ppm) (Meinshausen et 14 al., 2017) which then calculated the radiative forcing internally. These equations were updated to represent the newest findings from (Etminan et al., 2016). In contrast to that, radiative forcing for non-CO2 greenhouse gases (GHGs) was calculated 16 externally and summed up to be used as an additional model input, using concentration data of 43 GHGs (Meinshausen et al., 2017). We use updated radiative forcing formulations for CO2, CH4 and N2O following the findings of Etminan et al. (2016).

18
Radiative forcing of other GHGs was calculated using the formulations in Table 8.A.1 from the IPCC AR5 (Shindell et al., 2013). Meinshausen et al. (2017) introduced three options for calculating radiative forcing from GHGs concentrations. For 20 this study we chose to use the option in which one uses specific calculations for all available 43 GHGs, rather than treating some groups of GHGs in a similar manner.

22
The radiative forcing of stratospheric water-vapor from methane oxidation was calculated following the suggestion from Smith et al. (2018) by multiplying CH4 effective radiative forcing by 12%.

24
To calculate radiative forcing of tropospheric ozone, !"#$ , the equations from Smith et al. (2018) were used: where are transfer coefficients, %&' are methane concentrations, / are emissions of the respective species (NOx -nitrate 28 aerosols, CO -carbon monoxide, NMVOC -non-methane volatile organic compounds), /,)* are the respective preindustrial constants for the specific species, and T is temperature in Kelvin. Note that ( ) was not included in our calculations because the forcing is not calculated dynamically. Concentrations and emissions data were obtained from Input4Mip on the Earth System Grid Federation. Pre-industrial values were taken from Table 4 from Smith et al. (2018).
Halon 1202 data are not provided by Input4Mips and therefore was not included.
Three-dimensional aerosol optical depth (AOD) input for the UVic ESCM was created using a UVic ESCM grid and the scripts and data provided by Stevens et al. (2017), describing nine plumes globally which are scaled with time to produce 12 monthly sulphate aerosol optical depth forcing for the years 1850-2018 (for comparison see Fig. S1). The resulting AOD caused a too strong forcing in the historical period, a scaling routine was therefore implemented into the UVic ESCM, which 14 allows to scale the derived albedo alteration and accordingly aerosol forcing from AOD data. For transient simulations, the scaling factor was set to 0.7, which gives a globally average forcing of -1.4296 Wm -2 in 2011, consistent with the IPCC AR5 range estimate of between -2.3 and 0.2 Wm -2 .
Anthropogenic land-use changes (LUC) in the UVic ESCM are prescribed from standardized CMIP6 land use forcing 18 (Ma et al., 2019) that has been re-gridded onto the UVic grid. These gridded land use data products (LUH2), which contain information on multiple types of crop and grazing lands, were adapted for use with UVic by aggregating the crop lands and 20 grazing lands into a single "crop" type, which can represent any of 5 crop functional types, and a single "grazing" variable, which represents both pasture and rangelands. This forcing is used by the model to determine the fraction of each grid cell 22 that is crop or grazing land with those fractions of each terrestrial grid cell then assigned to C3 and C4 grasses and excluded from the vegetation competition routine of the TRIFFID component. CO2 emissions from LUC affect the model runs so that 24 when forest or other vegetation is cleared for crop lands, range lands or pasture, 50% of the carbon stored in trees is released directly in to atmosphere and the remaining 50% is placed into the short-lived carbon pool.

26
Historical volcanic radiative forcing data is provided by Schmidt et al. (2018). Following CMIP6 spin-up forcing recommendations (Eyring et al., 2016) volcanic forcing is applied as an anomaly relative to the 1850 to 2014 period in the 28 UVic ESCM.
Solar constant data for 1850 to 2300 was accessed from Input4Mips (Matthes et al., 2017). The available monthly 30 data was annually averaged. Following CMIP6 spin-up forcing recommendations (Eyring et al., 2016) spin-up values were set to the mean of 1850-1873, equal to 1360.7471 Wm -2 .
is shown in Fig. S2.

4
The model was spun up with boundary conditions as described in the CMIP6 protocol by Eyring et al. (2016) for over 10,000 years, in which the weathering flux was dynamically simulated and diagnosed. To diagnose the transient climate 6 response (TCR), equilibrium climate sensitivity (ECS), the ocean heat uptake efficiency ( _4x), and the transient climate response to cumulative emissions (TCRE), as given in Table 1) we ran 1,000 year simulations starting with a 1% per year 8 increase in CO2 concentrations until a doubling (2xCO2) and quadrupling (4xCO2) was reached after which the concentration was kept constant. Before switching from CO2 concentration driven simulations to CO2 emissions driven simulations, a 1,500 years drift simulation was run. Finally, the historical simulation is forced with fossil CO2 emissions, dynamically diagnosed land-use change emission, non-CO2 GHG forcing, sulphate aerosol forcing, volcanic anomalies forcing and solar forcing.

Fine tuning of the UVic ESCM 2.10
We tested version 2.10 of the UVic ESCM, with the main incentive to improve its skill in simulating carbon fluxes,

14
historical temperature trajectories and ocean tracers. While evaluating the model with available observational data, specific additional changes and updates were applied with respect to the UVic ESCM versions 2.9-02 ) and 2.9-CE 16 (Keller et al., 2014).
After merging the two model versions, the UVic ESCM's simulated historical cumulative land-use change emissions were close to zero, since its pre-industrial vegetation resembled closely the pattern of plant functional types of today. In order to get a good representation of deforested biomass, we updated the vegetation parameterization to ensure that diagnosed The new model version equilibrated with a rather low oceanic overturning strength, we therefore increased the ocean 30 background vertical diffusivity from previous value of 0.15 cm 2 s -1 in Keller et al. (2014) to 0.25 cm 2 s -1 , to increase ocean overturning (see Fig. S3 and S4). This caused ocean diffusivity to slightly increase in depths between 0 and 3500 m relative to the previous model version (Fig. S3), but to follow the tidal mixing profile very closely for higher depth. Global diffusivity 2 increased by about 4 %. This change enabled us to reach a very similar ocean overturning as found for the UVic ESCM 2.9-02, which uses the Bryan-Lewis mixing scheme ( Fig. S3 and S4). This stronger overturning then in turn also improved ocean physical properties (see Sect. 3.3. and supplementary material) and the global mean temperature and warming trends. However, it also causes the ocean heat content anomaly for the upper 700m to amount to 23.9 10 22 J, which is an overestimation of the 6 observed 700m OHC anomaly of 16.7 ± 1.6 10 22 J (Levitus et al., 2012) (Table 1). This seems to be a general feature of EMICs (Eby et al., 2013), but might still be problematic. 8

Evaluation of model components
In this section we evaluate the performance of the different components of the UVic ESCM version 2.10 based on observations.

Global key metrics -Temperature, Carbon Cycle, Climate Sensitivity and Radiation Balance
The emission-driven, transient historical climate simulation of the UVic ESCM version 2.10 forced with CMIP6 data 12 reproduces well the historical temperature trend in the 20 th century of 0.75 ± 0.21 K as derived from the Global Warming Index (Haustein et al., 2017) (Table 1, Fig. 1). Starting from the year 2000 the simulated global mean temperature increases at a 14 higher rate than previously, but the total temperature change since preindustrial remains within the uncertainty range of the estimate in the latest IPCC special report on 1.5 °C (Fig. 1, light grey cross) (Rogelj et al., 2018). This steep temperature 16 increase over the last twenty years of simulation amounts to a rate of temperature change of 0.21°C/decade, which is slightly higher than the best estimate from the infilled HadCRUT4-CW dataset of 0.17°C/decade (uncertainty range of 0.13-0.33°C/decade) (Haustein et al., 2017).  Eby et al. (2013). This might be caused by the rapid drop in Southern Hemisphere sea ice cover after about 430 years / 500 years in the 2x and 4x CO2 simulations, which in turn decreases the surface albedo and therefore causes 28 more warming. Since the ocean heat uptake efficiency is assessed at year 140 of the TCR4x simulation, it is as the TCR4x with 0.98 W m -2 °C -1 on the higher end, but still within the EMIC range, and does not yet reflect this process. In the same way, the transient climate response to cumulative emissions (TCRE) is with 2.2 K (1000 PgC) -1 higher than in previous model version but still within the likely range reported by the IPCC AR5 (Table 1) (Table 2). While the land sink is slightly higher than the best estimates the ocean sink is at the lower end of the given range. The simulated top of the atmosphere (TOA) short wave and long wave radiation of the UVic ESCM for the year 2005 4 lies well within the range of the CMIP5 models as reported by Wild et al. (2013) and agrees reasonably well with the observed estimates for both the solar and the thermal radiation fluxes (Table 3 or Fig. S5). The same is true for the simulated net surface 6 thermal flux, which is with -51 W m -2 , at the lower end of the CMIP5 range (Table 3). Now following the solar radiation through the energy moisture balance model, however, we find that the simulated atmospheric albedo with 0.227 is 8 comparatively low to the observed value of 0.250. This causes a rather low simulated absorption of solar radiation by the atmosphere of 68 W m -2 , where the observed estimate from Wild et al. (2013) is 79 W m -2 . Thanks to a rather high simulated surface albedo (0.171 compared to 0.13 from observations) the resulting absorbed solar radiation at the surface is still high but in contrast to the atmospheric absorption within the CMIP5 range. This results in a global mean surface net radiation of 117 12 W m -2 which is rather high compared to the observed best estimate of 106 W m -2 . This is the radiative energy available at the surface to be redistributed amongst the non-radiative surface energy balance components. Accordingly, the simulated sensible 14 heat flux in the UVic ESCM is with 31 W m -2 also too high compared to the CMIP5 range of (14.5, 27.7) W m -2 . Finally, the latent heat flux calculated from simulated evaporation is with 76 W m -2 on the very low end of observational and CMIP5 estimates, which is likely linked to the high transpiration sensitivity of plants in the UVic ESCM.

16
The simulated northern hemisphere summer sea ice extent with 3.4 million km -2 at the lower end of the CMIP5 18 estimates, and considerably smaller than the observed sea ice concentration (Table 1, and Fig. 3). This lower extent seems to be mainly due to a lack of simulates summer sea ice concentrations between 15 and 60 %, whereas higher concentrations show 20 a good agreement with the observed pattern (Fig. 3). The southward extension of winter sea ice concentration in the UVic is considerably smaller than the observations from passive microwave satellite missions. Concerning NH summer sea ice trends,

22
the UVic ESCM shows lower trends of -0.24 million km -2 dec -1 during the last 30 years, compared to what is observed (-1.07 to -0.73 million km -2 dec -1 ) (IPCC AR5, Chp.4, . The summer sea ice extent in the Southern Hemisphere 24 is with 1.3 million km -2 also smaller than the observed 3.1 million km -2 , and in contrast to the observed increasing trends in sea ice shows a decline of -0.25 million km 2 dec -1 (Table 1). While this is consistent with the simulated warming trend in the

4
Observed global mean terrestrial precipitation between 1961-1990 amounts to 818mm (Hulme et al., 1998). The 6 adjusted CO2 fertilization strength in the UVic ESCM 2.10 results in global mean terrestrial precipitation of 814mm for the same period (Table 1), bringing it close to the observed amount. Concerning terrestrial precipitation trends the UVic ESCM 8 2.10 shows a negative trend in terrestrial precipitation of -0.43 mm decade -1 for the period between 1930 to 2004 (Table 1).
This is in agreement with the range of terrestrial precipitation trends of [-4.2 -1.2] mm decade -1 given by Kumar et al. (2013).
The simulated pattern of annual mean precipitation flux for the last 30 years generally agrees well with the observed pattern ( Fig. 4). Similar to the seasonal temperature maps, the UVic ESCM slightly underestimates the most extreme amplitudes of 12 annual mean precipitation located in the tropical areas. The latitudinal mean values agree well in magnitude, but the tropical rain bands are extending too far north and south. Total terrestrial precipitation with 814 mm agrees well with the observed 14 1961-1990 mean value of 818 mm (Hulme et al., 1998). The terrestrial precipitation trend of -1.2 mm per decade also agrees well with the observed terrestrial precipitation changes for the recent historical period (1951 to 2005) of −7 to +2 mm per 16 decade, with error bars ranging 3-5 mm per decade (IPCC AR4) (Table 1).  14 The UVic ESCM overestimates vegetation carbon density of in tropical rainforest regions, such as in South America,

16
and Central Africa when compared to the revised estimates of Olson (1983,1985,2001) (Fig. 6). More recent biomass studies have challenged Olson's estimates for some regions of the world, but Olson (1983Olson ( ,1985 still provides the only globally-  (Table 2). In contrast, the simulated vegetation coverage of carbon densities of 2-5 kgC m -2 is lower than observations especially in Central Asia and at higher northern latitudes. This, however, does not imply that the dominant plant functional types, namely C3/C4 grasses, are underrepresented in this area. In the UVic ESCM 2.10 the representation of C3/C4 grasses, as well as needleleaf trees, in high northern latitudes improved compared to earlier 4 versions (see Fig. S7) thanks to the more complex soil module and the corresponding vegetation tuning. In summary, the UVic ESCM overestimates broadleaf tree cover in the tropics, but improved the representation of the vegetation cover at latitudes 6 north of 20 °N compared to previous model versions.

Ocean metrics -physical and biogeochemical
In the following section we will compare simulated ocean metrics with observations from the World Ocean Atlas The Taylor diagram for eight different ocean metrics illustrates that the UVic ESCM 2.10 improves ocean C14, and 6 slightly improves in ocean temperature, salinity, nitrate and phosphate distributions (dots in Fig. 9), relative to the UVic ESCM 2.9 (crosses in Fig. 9), given the same forcing. In contrast, mainly ocean alkalinity, but also dissolved inorganic carbon and 8 oxygen, show either a larger deviation or lower correlation compared to observations than the previous model version.
Generally, the model demonstrates skill in simulating these ocean properties, with correlation coefficients higher than 0.9 for all but the salinity and alkalinity fields, and root mean square deviation (rmsd) of below 50% of the global standard deviation of the observations, again with the exception of salinity and alkalinity. In the following we will discuss these features in more 12 detail.
global ocean (Fig. 10). The general good agreement for these ocean metrics is also seen in comparisons with vertical ocean section and ocean surface maps of these simulated fields with observations (see supplementary material). The only noteworthy 4 biases are too low C14 in the central Indian Ocean basin, indicating a too low overturning rate in this ocean basin (see Fig.   10 and 11), and small biases in simulated nitrate showing too high values in the Arctic Ocean compared to observation and 6 too low values in the Indian Ocean (see Fig. 10, Fig. 12 and supplementary material).

16
The apparent oxygen utilization (AOU) shows lower values (~15% lower) in the deep ocean, but otherwise reproduces the shape including the local maximum around 1000 m depth well ( Fig. 10 and 14). The main bias in AOU is found in the Southern Ocean, where the UVic ESCM 2.9 simulated too high values the UVic ESCM 2.10 now simulated too low oxygen utilization. This is especially true for the Atlantic and Indian Oceans. These biases in AOU are probably linked to biases in the ocean 20 biogeochemistry in the Southern Ocean. Comparing the simulated oxygen minimum zones (OMZ) of the UVic ESCM 2.9 and UVic ESCM 2.10 to observed OMZs, there is an improved representation of the asymmetry of the Pacific OMZ in the newer 22 model version, as well as a reduced bias in the Indian Ocean (Fig. 15).

Summary, Conclusion and Outlook
In order to obtain a new version of the University of Victoria Earth System Climate Model (UVic ESCM) that is to 14 be used in the 6th phase of the coupled model intercomparison project (CMIP6), we have merged previous versions of the UVic ESCM to bring together the ongoing model development of the last decades. In this paper we evaluated the model's performance with regard to a realistic representation of carbon and heat fluxes as well as ocean tracers in the UVic ESCM 2.10 in agreement with the available observational data and with current process understanding.
We find that the UVic ESCM 2.10 is capable of reproducing changes in historical temperature and carbon fluxes.
There is a higher warming trend in the southern hemisphere south of 40S compare to observations, which causes a bias in SH 4 sea ice trends. The simulated seasonal cycle of global mean temperature agrees well with the observed pattern, but has a lower amplitude. The air to sea fluxes of the UVic ESCM agree well with the observed pattern. The newly applied CO2 forcing 6 formulation has increased the model's climate sensitivity. Land carbon stocks concerning permafrost and vegetation carbon are within observational estimates, even though the spatial distribution of permafrost affected soil carbon and vegetation carbon 8 densities show regional biases. The top of the atmosphere radiation balance of the UVic ESCM is well within the observed ranges, but the internal heat fluxes are slightly off. The simulated precipitation pattern shows good agreement with observations, but are regionally too spread out especially in the tropics and as expected do not reach the most extreme values.
Terrestrial total precipitation and precipitation trends agree well with observations. Many ocean properties and tracers show 12 good agreement with observations. This is mainly caused by a good representation of the general circulation. Although problems remain, mainly the too low Southern Ocean oxygen utilization and the salinity bias in the AAIW.
14 These model data deviations, especially for the ocean tracers have not yet been fully addressed (note the misfits have been reduced relative to the previous model version forced with new forcing), since we are already planning for the next update of the UVic ESCM, which will incorporate more comprehensive biogeochemical modules that will require re-tuning of the oceans biogeochemical parameters as well. Model developments that have not been incorporated in the model version 18 described here, like carbon-nitrogen feedbacks on land (Wania et al., 2012), explicit representation of calcifiers in the ocean (Kvale et al., 2015), dynamic phosphorus cycle in the ocean (Niemeyer et al., 2017), and others, will in the following be 20 implemented and tested within this new model version.

24
including all necessary documentation and data upon final publication of the manuscript.