UKESM1: Description and Evaluation of the U.K. Earth System Model

3.1 Snow-Vegetation Interactions

The most significant tuning of the coupled ESM concerned the interaction between shortwave (SW) radiation, snow, and vegetation. We found it necessary to tune the parameters governing this interaction in order to reduce a substantial SW radiation bias. The extent to which vegetation is covered by snow makes an enormous difference to surface albedo, and there is a large diversity in how models parametrize this process (Boisier et al., 2012; Bright et al., 2015; Ménard et al., 2014). JULES represents the bending and partial burying of vegetation by snow (Ménard et al., 2014) on a per-PFT basis, by reducing the leaf area index (LAI) used in the albedo calculation as snow depth increases (Walters et al., 2019). When the depth of snow (d) is less than the canopy height h_i of a given PFT, the partial burying is described by the equation

(1)

where i denotes the PFT, L_i is the true LAI for each PFT, and L_i′ is the “exposed” LAI used in the albedo calculation. The parameter n_i thus governs the rate at which vegetation is covered by snow as snow depth increases, but it is not well constrained by observations and must be tuned for a given LAI distribution. The LAI in GC3.1 is prescribed using a seasonal climatology derived from MODIS satellite observations, while the phenology in UKESM1 is simulated interactively using the scheme of Clark et al. (2011). Figure 1 shows LAI for C3 grass in the two models in JJA and DJF. GC3.1 shows an almost complete loss of LAI in winter in central Asia and North America. This seasonal reduction is more extreme than seen in other observational products (Garrigues et al., 2008) and may in part result from the obscuring of vegetation by snow as seen by the satellite data which is used in constructing the LAI climatology. Even the needleleaf tree PFT in this climatology sees LAI reductions of 80–95% in regions dominated by evergreen needleleaf trees (not shown). Therefore, in GC3.1 there is a little need for the parametrization in equation 1, since LAI is already very small and n_i is accordingly tuned to very low values.

In contrast, the LAI in UKESM1 is more similar between summer and winter for grasses and evergreen PFTs, with reductions of 20–35%. This wintertime LAI may be too high because UKESM1 does not fully represent the phenology of grasses, but the large difference from the LAI climatology used in GC3.1 requires a different tuning for the snow-vegetation parametrization. Accordingly n_i is given higher values in UKESM1, particularly for grasses, so that the radiative effect of snow partially covering vegetation is taken into account. The specific values of n_i were chosen to maximize the similarity in grid box mean wintertime surface albedo between UKESM1 and GC3.1, to avoid large differences in regional mean radiative fluxes. Using the GC3.1 values of n_i in UKESM1 produces a widespread negative bias in DJF clear-sky outgoing SW across northern mid-latitude land masses (Figure 2c). The tuning of n_i increases the wintertime albedo in seasonally snow-covered regions, increasing the TOA outgoing clear-sky SW (Figure 2b), which reduces the negative bias but does introduce positive bias in some regions (Figure 2d). The net impact of this tuning can be summarized by the second metric in Table 2, defined as TOA outgoing clear-sky SW mean over land between 30°N and 60°N for DJF, representing the clear-sky radiative impact of snow on TOA radiation. A crude observational range of 50.5–53.6 W m⁻² has been derived by calculating the same metric for two satellite data sets: CERES-EBAF Ed4.0 (Loeb et al., 2018) and ISCCP-FD (Zhang et al., 2004), using the 2000–2010 period common to both data sets. Table 2 shows that the impact of the tuning was to increase this metric from well below the observational range to just above it.

3.2 Bare Soil and Dust Emissions

In UKESM1, emissions of mineral dust are tightly coupled to the fraction of bare soil in a grid box and hence are very sensitive to biases in vegetation fraction in arid and semi-arid areas. In particular the model, prior to tuning, was prone to excessive emissions of dust in Australia, India, and central Asia. We tuned bare soil fraction in these regions by increasing the vegetation dynamics parameter lai_min (Clark et al., 2011) for grass PFTs. When LAI for a given PFT exceeds this lai_min, some excess productivity is allocated to horizontal spread of the PFT within the grid box. For the same productivity, smaller values of lai_min allow greater vegetation fractions and hence smaller bare soil area. Two satellite land cover data sets were used as a reference for this tuning: European Space Agency (ESA) Climate Change Initiative Land Cover data (CCI-LC; Poulter et al., 2015) and that from the International Geosphere-Biosphere Programme Land Use and Cover Change project (IGBP-LUCC; Loveland et al., 2000). The impact of increasing lai_min is illustrated by considering the total area of Australian bare soil in the model, which is reduced by nearly a factor of two to 2.8 ×10⁶km², covering around one third of the continent. The resulting bare soil area lies between the equivalent values for IGBP-LUCC and CCI-LC of 2.5 ×10⁶km² and 4.0×10⁶km², respectively. Bare soil areas in central Asia were similarly reduced, but the same is not true for India where the model suffers from a deficit in monsoon rainfall (Williams et al., 2018), which leads to such low plant productivity that the lai_min parameter has little or no impact.

Even with reduced biases in vegetation cover, there was also a need to tune the parameters of the dust emission scheme itself. These have been described by Woodward (2011), and in summary they control the emission response to soil moisture and friction velocity and include an overall scaling factor for the total emission. These were tuned to maximize agreement with in situ surface concentration measurements from the University of Miami network (18 sites), and Aerosol Robotic Network (AERONET; Holben et al., 2001) total aerosol optical depth (AOD) measurements at eight dust-dominated sites. While regional metrics are used in the tuning process (see section 4.6), overall performance is summarized by a single global dust error metric, defined as the RMS difference between and for all of these University of Miami and AERONET observations. The square root is used as a weighting function to increase the contribution of lower concentrations at remote sites to the global metric. The calculation uses seasonal means of the model and observations; each site receives equal weight; thus the concentration observations contribute more strongly than AOD: 18 versus 8 sites. The final UKESM1 configuration has a global dust error of 0.47 (unit-less), which is just over half the error of the untuned configuration; the reduction comes from a combination of the changes to dust parameters and the reduction in bare soil area described above. For context, the value of this metric is 0.51 for GC3.1 and was 0.98 for HadGEM2-ES.

3.3 Terrestrial Productivity

Carbon fluxes and stores are key performance indicators for an Earth system model, due to their importance in regulating carbon feedbacks. The global total GPP of plants in UKESM1 was tuned downward to better agree with observations by reducing the quantum efficiency of photosynthesis. Quantum efficiency is specified separately for each PFT, and the UKESM1 values were adjusted to the lower end of the ranges found by the literature review of Skillman (2008). To inform this tuning we combined GPP from our pre-industrial tests with an estimate of the likely pre-industrial-to-present-day change based on previous models and compared this against the present-day GPP estimated by the FLUXNET model tree ensembles (MTE) data product of Jung et al. (2011). The goal was to improve both the annual mean and seasonal cycle of the global total GPP. The global total GPP for the last 20 years of the resulting UKESM1 historical simulations is 130 GtC year⁻¹, a little higher than the Jung et al. (2011) 5–95% range of 119±6 GtC year⁻¹.

3.4 Surface Chlorophyll

The ocean-atmosphere coupling of PMOA emissions described in section A5 makes aerosol emission dependent upon the interactively simulated surface chlorophyll concentrations, which allows the model to capture potential feedbacks between ocean biogeochemistry and climate. However, ocean models typically have larger biases in ocean chlorophyll than other marine biogeochemistry properties (Kwiatkowski et al., 2014; Yool et al., 2013), posing the risk of introducing a large bias in aerosol load. We mitigate this risk by multiplying the chlorophyll concentration by a global scaling factor. A value of 0.5 is used for this factor, chosen to produce a similar global mean outgoing SW radiation flux at TOA to that in the atmosphere-only implementation of the PMOA emission scheme, which was driven by the observation-based GlobColour data set (Maritorena et al., 2010). This scaling factor reflects the fact that MEDUSA chlorophyll concentrations are higher than those of GlobColour. Note that this scaling of surface chlorophyll is performed only within the PMOA emission parametrization: There is no impact on the prognostic chlorophyll simulated by the ocean model.

3.5 Energy Balance

The physical core of UKESM1 was already well-tuned before adding ESM components (Kuhlbrodt et al., 2018; Williams et al., 2018), and with the exception of the snow-vegetation interaction above it was not deemed necessary to re-tune the physical model parameters. The TOA radiation balance was however altered by the addition of ESM components, through a combination of changes in cloud, surface albedo, and radiatively active gases. Without further tuning the net downward radiation at TOA in 1850 was −0.81 W m⁻², which would have resulted in a significant downward drift in model temperatures. We brought the TOA radiation into balance by tuning parameters in the Anderson et al. (2001) parametrization of DMS seawater concentration to permit lower minimum values of DMS. The standard configuration of this parametrization has a prescribed minimum value of 2.29 nM for seawater DMS concentration, while gridded observational DMS data sets, such as those of Kettle et al. (1999) or Lana et al. (2011), contain significantly smaller values than this over large regions of the ocean. Anderson et al. (2001) themselves point out that the data from which their parametrization is derived is likely to contain a sampling bias toward higher values, so a lower minimum is reasonable. We reduced the minimum from 2.29 to 1.00 nM, extending the Anderson et al. (2001) linear relationship (between DMS concentration and log₁₀ of chlorophyll concentration, surface SW, and nutrient availability) to lower values of DMS. The ensuing reduction in DMS gave widespread decreases in cloud droplet number (and thus cloud albedo) across the Southern Ocean and stratocumulus regions, resulting in a drop in reflected SW at TOA of 2 to 5 W m⁻² over large areas of the ocean.

Using the ocean DMS parametrization to balance the TOA, rather than for example cloud parameters, allows greater traceability between UKESM1 and its physical core GC3.1, by minimizing the number of physical model parameter changes between the two models. While the goal of this tuning was to balance the pre-industrial TOA, we note that the mean TOA over the last 20 years of the tuned historical run is within the Johnson et al. (2016) present-day observational estimate of 0.61–0.81 W m⁻² (5–95% confidence).

3.6 Effectiveness of Model Tuning

Despite the difficulties noted at the start of this section, the tuning exercise was successful in reducing the biases it sought to address. The metrics in table 2 which summarize these biases are all improved relative to the untuned equivalents; that is, the pre-industrial net TOA is closer to zero, and present-day values of the other metrics are closer to their observed values. Recall that the tuning decisions were based mainly on pre-industrial tests, so the present-day performance could not be assured until historical simulations were completed. In the case of the snow-vegetation interactions affecting the outgoing SW metric, the tuning possibly went too far since the present-day value of the “snow SW” metric changed from low-biased to high-biased, though still closer to the observations. This illustrates the imprecise nature of the tuning exercise.

Finally we note that, while the parallel runs presented in Table 2 give a good indication of the magnitude of the direct impact of tuning, they are likely to underestimate the full impact. This is because these tests are designed to be relatively short (25 years), in order to minimize confounding internal variability. While there is minimal drift after the first 5 years, the increased GPP and dust emissions in particular would produce slow responses in terrestrial carbon stores and ocean biogeochemistry respectively which could take more than a century to equilibrate. Assessing the full impact of these changes would require the untuned pre-industrial configuration to be run to equilibrium, before initializing a new ensemble of historical simulations; such an exercise is too costly to perform at this time.

We have not yet assessed the impact of this tuning on emergent properties such as radiative forcing or climate sensitivity (as done for example by Mauritsen et al., 2012); this will be considered in future work.

4.1 Pre-Industrial Trends and Variability

The pre-industrial control simulation is important for quantifying the model's natural variability and any long-term drifts in the model state. Significant drifts in climate or the carbon cycle would compromise the model's utility for many applications, particularly in assessing future Earth system responses to anthropogenic emissions. Figure 3 shows how five key properties of the model vary over the pre-industrial control. Generally drift is very small relative to internal variability.

There is a relatively constant outgassing of CO₂ from the ocean at a rate of ≈0.02GtCyear⁻¹, slowing slightly as the run progresses, and a release of soil carbon of a similar magnitude. Both of these, and the total carbon uptake, are well below the ±0.1GtCyear⁻¹ threshold set by C4MIP (Jones et al., 2016) as a guide for acceptable drift in carbon pools, shown as dotted lines in Figure 3.

Net TOA radiation is very close to zero, which ensures that there are no discernible drifts in surface temperature or sea ice extent. TOA radiation does drift very slowly over the course of the run with a mean of 0.009 W m⁻² over the first 500 years, and 0.038 W m⁻² over years 500 to 1,000.

Surface air temperature shows strong multi-century variability, most notably during years 400 to 900, with swings of up to 0.3 K over a 100-year period (Figure 3). This is linked to episodes of oceanic deep convection in the Southern Ocean which bring warm saline deep Atlantic water to the surface, leading to a significant reduction in sea ice and loss of accumulated ocean heat to the model atmosphere. The warm, convecting, phases of these cycles coincide with reductions in Antarctic sea ice cover.

Aggregated vegetation types show variability on a range of timescales. Grasses expand and retract into regions of bare soil on interannual and decadal timescales in response to local variations in rainfall and surface climate, resulting in a clear anti-correlation between grass and bare soil in Figure 3. Tree coverage evolves on a much slower timescale, with a small decrease of around 1% during the first 600 years of the run, followed by a slight recovery, mostly at the expense of shrubs.

The UKESM1 historical simulations are initialized using atmosphere-sea ice-land-ocean conditions taken at different points in the pre-industrial control simulation. Ideally, these initial conditions will sample the model's multi-decadal internal variability so that each historical run is initialized using well-separated states of the pre-industrial coupled model. To achieve this, we consider two large-scale modes of sea surface temperature (SST) multi-decadal variability, the Inter-decadal Pacific Oscillation (IPO; Power et al., 1999; Zhang et al., 1997) and the Atlantic Multi-decadal Oscillation (AMO; Kerr, 2000).

The IPO pattern is the largest multi-decadal quasi-periodic SST variability over the Pacific basin. It is represented by the first area-weighted Empirical Orthogonal Function (EOF) of deseasonalized, detrended, and low-pass filtered monthly SSTs. Figure 4 shows the IPO pattern for HadISST1.1 (1870–2017; Rayner et al., 2003) and for the UKESM1 pre-industrial control simulation. The IPO index is the principal component of this EOF. The AMO pattern is a quasi-periodic oscillation of North Atlantic SSTs. The AMO index is calculated as the area-weighted mean time series of the deseasonalized, detrended, and low-pass filtered North Atlantic SSTs. Low-pass filtering the data minimizes the influence of the El Nino/Southern Oscillation (ENSO). Observed and modeled AMO pattern are shown in Figure 5.

The observed and simulated patterns are similar, for both modes of variability. The most notable difference between the two IPO patterns is the westward extent of the equatorial warm region, which extends further toward the maritime continent in the simulation. The AMO patterns differ in terms of the northerly position of the warmest region. Nevertheless, the similarity of these climate mode patterns provides confidence that the simulation is, to first order, replicating the observed multi-decadal variability in the Pacific and North Atlantic basins.

The phase space spanned by the observed and modeled IPO and AMO indices is shown in Figure 6. The observations show no clear relationship between the IPO and AMO and populate all four sectors of phase space. The model has a much longer time series and hence explores the phase space more fully. Consistent with the observations, the IPO and AMO appear to be independent of each other in the model simulation. We use the simulation phase space to choose well-separated initial conditions for each historical run (marked in Figure 6), thus ensuring that the beginning of the historical ensemble samples the multi-decadal variability of the ocean. For the same reason, we also prescribed that initial conditions for historical members were separated by at least 30 years in the pre-industrial control.

The model also exhibits a realistic simulation of the El Nino/Southern Oscillation (ENSO). The addition of ESM components has made no significant difference to the ENSO behavior of the physical model GC3.1 documented in Kuhlbrodt et al. (2018). We have repeated the ENSO analysis of Kuhlbrodt et al. (2018), and the results are statistically indistinguishable (not shown).

4.2 Present-Day Surface Climate

In Figure 7 we evaluate climatological seasonal mean sea surface temperature against the HadISST1.1 observation data set (Rayner et al., 2003). As is common in ocean models of this resolution (Flato et al., 2013), UKESM1 exhibits a strong cold bias in the NW Atlantic, related to the too-zonal representation of the Gulf Stream and North Atlantic current. Again, similar to other low-resolution models, we see a year-round cold bias in the north Pacific, and a too-westward extension of the equatorial Pacific cold tongue during JJA. Overall the pattern of SST bias is similar to that of GC3.1 (Kuhlbrodt et al., 2018), including a warm bias of 0.5–2 K across much of the Southern Ocean.

Figures 8 and 9 make similar seaonal comparisons for climatological 1.5 m air temperature over land and precipitation, against the observation data sets of CRUTEM4.6 (Jones et al., 2012) and GPCP2.3 (Adler et al., 2003), respectively. UKESM1 shows a wintertime cold bias over North America, Europe, and NW Russia. One driver of this cold bias relates to the SST biases noted above, influencing land temperatures through atmospheric advection. These continental cold biases are larger than the mean upwind SST bias, suggesting a potential feedback involving snow cover, snow albedo, and surface temperatures whereby a cold surface temperature bias is maintained in the model through excess snow cover, both spatially and temporally through the winter season.

In common with the CMIP5 generation of models (Flato et al., 2013), UKESM1 exhibits a double intertropical convergence zone in the tropical Pacific, linked to the Equatorial cold tongue bias. Tropical Atlantic precipitation is too far south, and in DJF this pattern extends westward over South America. In JJA UKESM1 exhibits the same dry Indian monsoon bias described by Williams et al. (2018) for GC3.1. This dry bias leads to a vegetation deficit (described in section 4.3) and some of the Earth system feedbacks discussed by Martin and Levine (2012).

While our emphasis in this paper is on evaluating UKESM1 at larger spatial scales (i.e., continental or ocean basin scales) and longer time scales (multi-decadal to centennial), UKESM1 will be used to make a set of future projections following the CMIP6 scenarioMIP protocol (O’Neill et al., 2016). Output from these projections will be used to drive impacts models (e.g., to assess future risks to food security or water resources). Such models require high temporal and spatial resolution surface climate data as forcing. Furthermore, assessing a model's higher time frequency statistics increases confidence in the simulated long-term means, variability, and future changes, all of which are based on the underlying high time frequency behavior. As an initial evaluation of UKESM1, Figure 10 shows simulated and observed daily timescale precipitation for four distinct climatic regions, and Figure 11 shows a similar analysis for monthly mean near-surface air temperature.

For each of the regions, Figure 10 plots the fractional contribution from different daily precipitation rates to the total seasonal accumulation, calculated for the seasons specified in the figure labels. Observations include the GPCC/HOAPS v1 global daily precipitation (Andersson et al., 2010; Schamm et al., 2014), based on both gauge and satellite observations, the Global Precipitation Climatology Project (GPCPv1.1; Adler et al., 2003), and the Tropical Rainfall Measuring Mission 3B42v6 product (TRMM; Huffman et al., 2007). All observational data were first re-gridded onto the UKESM1 grid. TRMM data are only used over the SEA and AMZ regions, whereas GPCC is not used over the predominantly ocean area of SEA.

Over EUR we show winter (DJF, full lines) and summer (JJA, dashed lines) precipitation, with the nine UKESM1 members in black and the observations in different shades of blue. In both seasons the distributions are skewed toward a long tail of light precipitation (<1 mm day⁻¹) that contributes to the seasonal accumulation. UKESM1 captures this skewness quite well. In DJF the model slightly overestimates precipitation contributions from moderate to heavy daily rates (3 to 10 mm day⁻¹) but does capture the heavy rain tail of the distribution. Conversely, in JJA moderate to heavy rainfall is underestimated by the model while very heavy events (>20 mm day⁻¹) contribute too much to the seasonal accumulation. As the moderate to heavy daily rates contribute more to the seasonal accumulation (i.e., they occur more frequently) than the very heavy rates, the underestimate of the former results in a modest dry bias over EUR in JJA, primarily in central and south-east Europe (Figure 9).

Turning to the EAS region, where precipitation is dominated by the summer monsoon, we plot distributions for the summer (wet) season (JJAS) and the winter (dry) season (NDJFM). In JJAS UKESM1 significantly overestimates the contribution of rainfall events >30 mm day⁻¹, with some unrealistic daily rates in excess of 100 mm day⁻¹. Figure 9 indicates that this results in a significant wet bias, primarily over the ocean part of EAS (i.e., the west Pacific warm pool region). In NDJFM seasonal precipitation is significantly lower than JJAS in both observations and the model, although the model systematically overestimates precipitation from most of the bins. Again Figure 9 suggests this overestimate is mainly along the coastal ocean in EAS. To the south of the EAS region, over the maritime continent and west Pacific warm pool (SEA) the rainfall distribution in the main wet months (NDJFM) is very well captured, with only a modest overestimate in the model across most bins. Finally, the AMZ region shows precipitation over the Amazon rainforest. While the overall spread of the observed distribution is captured by the model, there is a significant underestimate of moderate to heavy rainfall days (approximately 8 to 30 mm day⁻¹) that leads to the seasonal mean dry bias seen over the central Amazon basin in Figure 9, which is surrounded to the south by a wet bias of equivalent magnitude.

In Figure 11 for the same four regions we compare spatially averaged 1.5 m monthly mean air temperatures between UKESM1 and CRUTS4.02 observations (Harris et al., 2014). From this we wish to see if UKESM1 captures the spread of cold/warm months around its climatological mean value, providing an evaluation of the inter-annual variability of near-surface climate in the model. For all regions we first interpolate the 0.5° CRUTS data onto the UKESM1 grid and then mask out ocean and coastal points, using the UKESM1 land-sea mask and a criterion that a grid point has to have fractional land amount >0.995 to be included in the spatial average.

Considering EUR, in winter (DJF) a systematic bias is evident in the model distribution, with mode values shifted relative to CRUTS by 3–4°C. Nevertheless, the spread of cold/warm months around the mode of the distribution, as well as the shape (skewness), are well captured in the model. Figure 8 indicates the bulk of the cold bias in EUR is over Scandinavia and European Russia, as noted above. In the summer (JJA) a smaller cold bias remains in the model distribution, which continues to be located in the northern part of EUR. A similar cold bias is seen over EAS during boreal winter, while the summer distribution is well simulated, albeit with a slight tendency for too frequent cold months (surface temperatures <18°C) not seen in CRUTS. Over AMZ, the mode of the UKESM1 distribution is quite accurate, although the model now has a bias toward extreme warm months (mean temperatures >27°C) not seen in CRUTS. These are likely associated with the dry precipitation bias over the region (Figure 9), shifting the surface energy balance from predominantly evaporative toward a larger sensible heat component. Finally, the SEA region shows a cold bias in the model distribution. It should be noted for SEA that due to the small number of land points used in the spatial means and the extreme topography of islands in this region, neither the model nor CRUTS data are likely a robust representation of true land temperatures in this region. We plan to extend our analysis of UKESM1 near-surface temperatures in a future study, considering daily temperature variability, as well as daily maximum and minimum values to better isolate model processes controlling biases in regional surface temperature.

Figure 12 compares the areal extent of snow with the GlobSnow observation product (Takala et al., 2011) for four different snow depths. In addition to the area for observed monthly mean snow depths, the equivalent for monthly maximum depth is shown as an indicator of inter-annual variability. The model has too much very thin snow compared with the observations, but a comparable amount of deeper snow. The deeper snow tends to thaw slightly later than the observations in the spring but increase at a similar rate to the observations in the autumn.

Figure 13 compares permafrost extent diagnosed from the model against the data set of Brown et al. (2002). The model has permafrost in the majority of the regions where the Brown et al. (2002) data set suggest there is continuous permafrost. There is permafrost for approximately 47% of the land surface where the air temperature is less than zero. This falls within the range of observations which is 40% to 72%.

Figure 14 compares September and March sea ice concentration/extent, in both hemispheres, against the HadISST.2.2.0.0 data set of Titchner and Rayner (2014) for the period 1990–2009. Also shown are the modeled and observed sea ice extent (the total area of ocean with at least 15% ice cover) and the corresponding Integrated Ice Edge Error (IIEE). The IIEE metric, introduced by Goessling et al. (2016), is a measure of the error in ice edge placement and is essentially the area integral of all grid cells where the model and observations disagree about whether sea ice concentration is above 15% (see Goessling et al., 2016, for further details).

In general there is very good agreement between the modeled (shading) and the observed (orange lines) location of the ice edge as depicted by the 15% concentration contour. In the Arctic the model sea ice is, overall, a little more extensive than the observations, particularly so in the north Pacific (Bering Sea and Sea of Okhotsk) in winter (March). The model also has too much winter sea ice north of Iceland and does not reproduce the Oden ice tongue, which is captured in the HadISST data set. In the summer the Arctic sea ice is slightly over extensive within the Arctic Ocean basin, except for in the Kara Sea where the sea ice extent is much smaller than observations. In the Antarctic the sea ice is slightly less extensive than observed, particularly in the Ross Sea and east of the Weddell Sea sector south of Africa.

In response to the cool northern hemisphere surface temperature reported above and in section 4.7, the present-day sea ice is much thicker than observational estimates. Mean September sea ice thickness for 1990–2009 is over 4.5 m throughout the central Arctic (not shown). However, despite the high modeled thicknesses, the spatial extent of sea ice in UKESM1 agrees well with the HadISST observational data set.

4.3 Terrestrial Biogeochemistry

To have confidence in the processes responsible for terrestrial carbon uptake, one needs to evaluate model simulated vegetation structure, surface carbon fluxes, and carbon stores. If a model performed well for one or two of these properties, but not all three, it would suggest that the model was getting the right answer for the wrong reasons. In this section we evaluate and characterize the behavior of each of these three aspects.

UKESM1 represents vegetation structure with nine natural PFTs and four crop/pasture PFTs, and the fractional coverage of each PFT within a grid box is determined by a competition hierarchy. We evaluate the vegetation distribution using two complementary approaches: first with spatial maps of aggregated PFTs (Figure 15) and second by aggregating fractions within biomes (Figure 16). The biomes are regions of the world defined by Olson et al. (2006) based on their common ecosystem characteristics. Within each of these regions Figure 16 shows the fractional coverage of aggregated PFTs for UKESM1 and the IGBP-LUCC and CCI-LC land cover data sets; the same region definition is used for all three data sets. Broadly the model does a very good job in representing vegetation structure of the biomes. Notable biases include an excess of C3 grass in tundra regions which the observations indicate should contain more bare ground; this can be seen in the northern Siberian C3 grass fraction in Figure 15. Some of this bias is a side effect of the lai_min tuning described in section 3, which improved the bare soil fraction in dust-producing regions but has allowed grass to spread too extensively in northern Russia and Canada. In contrast, Saharan bare soil extends too far south, with a deficit of grassland in the Sahel. Similarly there is excess bare soil in eastern India. Both of these biases result from precipitation deficits in these regions, associated with errors in the position and intensity of monsoon rains (see Williams et al., 2018). There is a small overestimation of tree fraction in savanna biomes, most notably on the southern edge of the Amazon forest. This is due to the lack of fire disturbance, the inclusion of which would be expected to improve vegetation structure in these regions (Burton et al., 2019).

Figure 17 shows the evolution of terrestrial carbon fluxes in the ensemble mean of the historical simulations. From 1850 to 1980, decadal-mean GPP slowly increases from 118 to 122 Pg year⁻¹, after which it increases rapidly to 134 Pg year⁻¹ in the last decade of the experiment. GPP in UKESM1 is slightly higher than observational estimates, which are around 120 Pg year⁻¹ in the last few decades (Jung et al., 2011). Net ecosystem productivity (NEP), the balance of GPP and respiration, also increases more rapidly after 1980. NEP is balanced by a small land use change emission flux and a much larger crop harvest flux; the sum of these three fluxes is the net biome productivity (NBP), the total carbon flux from the atmosphere to the terrestrial biosphere. The land is a net source of carbon prior to 1980: cumulative NBP is −40 Pg from 1850 to 1979. After 1980 the land becomes a sink of atmospheric CO₂: cumulative NBP between 1980 and 2015 is 22 Pg.

Generally the model has an excellent zonal distribution of vegetation carbon (Figure 18), lying between the two observational estimates. The exception is a high-bias between 25°S and 10°S, reflecting the excess of trees in savanna biomes noted above. Global total vegetation carbon is 551 Pg in year 2000 and 583 Pg in 1850 which is consistent with observation-based estimates (450-650 Pg in 1750; Ciais et al., 2013); 72% of vegetation carbon is stored in tropical forests and 26% in extra-tropical forests, with the remaining 2% being stored in shrubs and grasses.

At 1,780 GtC, the UKESM1 global total soil carbon is less than that estimated by the Carvalhais et al. (2014) data set for the whole soil profile (2,450 GtC) and slightly less than that estimated by Batjes (2016) for the top 2 m (1,970 GtC). The model represents the top 3 m of soil. There is not enough soil carbon simulated for the northern high latitudes (Figure 18); this is partly because turnover times are too fast (not shown) but also because there is not yet a representation of permafrost carbon included within UKESM1. The simulated peak of soil carbon in the tropics is slightly further south than in any of the observational data sets.

4.4 Ocean Biogeochemistry

Ocean biochemistry in UKESM1 plays a fundamental role in the carbon cycle, acting as a source or sink of atmospheric CO₂ and in the aerosol simulation via natural emissions of DMS and PMOA. In this section we focus on these two roles, evaluating the present-day surface flux of CO₂ and the productivity which underpins the aerosol emissions.

Figure 19 shows simulated net primary production, averaged over the period 2000–2009, compared with three satellite-derived observational estimates for the same period (Behrenfeld & Falkowski, 1997; Carr et al., 2006; Westberry et al., 2008, these three estimates are averaged here to a single field, as per Yool et al., 2013). In general, the model simulates the main geographical patterns of high productivity in upwelling regions, low production in subtropical oligotrophic areas, and seasonally high productivity in high latitude areas. However, modeled productivity is more tightly focused on the equator, especially in the Pacific Ocean. Oligotrophic areas exhibit lower production in the model, in part due to unrepresented mesoscale processes, and high latitude productivity is elevated above that observed, especially in the Southern Ocean.

Figure 20 shows simulated net air-sea CO₂ flux, averaged for the period 2000–2009, and compared with the observational estimate of Rödenbeck et al. (2013). The model captures the main characteristics of the broad pattern of CO₂, with outgassing (blue) in the equatorial and tropical band, especially in the Pacific Ocean, and ingassing (red) at higher latitudes, in particular the North Atlantic. The model does show some discrepancies, for instance excessive outgassing in the upwelling areas associated with the Humboldt and Benguela Currents, and excess ingassing close to Antarctica linked to a deficit in sea ice cover in spring and summer.

4.5 Atmospheric Chemistry

A thorough evaluation of the chemistry performance of UKESM1 can be found in Archibald et al. (2019); here we focus on those aspects which represent major enhancements relative to HadGEM2-ES, namely, stratospheric chemistry, interactive photolysis, and biogenic emissions of VOCs.

First we assess the quality of the simulation of total column ozone (TCO), which is dominated by the stratospheric ozone layer, versus the NIWA-Bodeker Scientific TCO database (version3.4; Bodeker et al., 2005) covering 1979–2016. We also use the TCO record measured at South Pole using a Dobson photospectrometer (Rosen et al., 1993; WOUDC, 2014).

Figure 21 shows the zonal-mean TCO averaged over 1995–2014, averaged over two historical UKESM1 simulations, and the difference versus the climatology. In the extra-polar regions, there is a general overestimation of ozone of around 20 to 40 Dobson Units (DU), whereas in the Antarctic, the model underestimates ozone relative to the climatology by up to −60 DU in November and December. Since satellite-based observation of TCO is particularly difficult at high latitudes, we compare this result directly to the Dobson long-term measurements at South Pole (Figure 22). The comparison shows a good agreement with the Dobson record in October, but in November and December, the model only tracks the most stable, deepest ozone holes and does not reproduce the observed interannual variability. Linear trends during the period of increasing ozone depleting substances (1979–1997) mirror this finding (Figure 23). In the tropics the model does not produce significant trends, in agreement with the climatology, but in both polar regions the model produces too large negative trends during this period during spring and summer in both hemispheres. A separate analysis of trends during the ozone recovery period (1998–2014) shows generally non-significant trends, in agreement with the observations.

An analysis of the vertical distribution of ozone and its trends (not shown) reveals a generally good correspondence with observations but too strong a downdraught of mesosphere air during polar winter which contributes to low TCO during spring and also the general overestimation of ozone in the extra-polar ozone layer discussed above.

The good correspondence in October suggests that the sensitivity of ozone to anthropogenic chemical ozone depletion is generally correct, broadly in agreement with scientific consensus on the underlying gas-phase chemistry (World Meteorological Organization, 2018), but the underestimation of variability in November and December suggest insufficient, and insufficiently variable, meridional heat transport toward both poles which would serve to inject variability into TCO and would shorten, on average, the lifetime of the polar vortex (Garfinkel et al., 2013). This is consistent with a stratospheric cold bias of 5–8 K over the Antarctic and a polar night jet that is too strong by more than 10 m s⁻¹ (not shown).

The general overestimation of ozone in the extra-polar regions would be partly addressed by a more comprehensive inclusion of heterogeneous reactions on sulphate aerosols involving bromine compounds (Dennison et al., 2019), which is not part of UKESM1. The nine CMIP5 models that did simulate stratospheric ozone interactively exhibited substantial biases (Eyring et al., 2013). Relative to this previous generation of models, UKESM1 compares well to the reference climatology.

Turning now to the troposphere, Figure 24 shows that for the historical simulation of UKESM1, the burden of ozone in the troposphere (defined as in Tilmes et al., 2016) has increased steadily over the period 1850–2014. The burden of ozone has recently been evaluated for the present-day (2014–2016) by Gaudel et al. (2018) as 300±12 Tg from 60°S–60°N. The UKESM1 tropospheric ozone burden, subset over the same latitude range for the year 2014 is 314 Tg, in excellent agreement with the observations. UKESM1 simulates an increase over the historical period (1850–2014) of more than 60 Tg of ozone in the troposphere, an increase of approximately 23% from a pre-industrial baseline of 260±4 Tg (1850–1875) to 320±7 Tg (1995–2010).

Methane lifetime is a core metric for atmospheric chemistry models due to the role methane plays in multiple chemical cycles, with impacts on ozone, aerosols (via oxidant concentrations), as well as exerting a strong direct radiative effect in its own right. The methane lifetime in UKESM1 is 8.4 years (calculated as total atmospheric burden divided by chemical loss and deposition for the period 1995–2010). This lies within the 9.1±0.9 year range estimated by Prather et al. (2012) based on constraints from methyl chloroform, and the 7.6 to 14 year range of SPARC (2013). Table 3 summarizes methane lifetime along with five other metrics of the global chemical state of the model.

Figure 25 shows the global distribution of isoprene and monoterpene emissions simulated by the interactive BVOC scheme in UKESM1. The emissions in this figure represent the mean over the last 10 years (2005–2014) of the historical run. Annual global total isoprene and monoterpene emissions amount to approximately 520 and 120 TgC year⁻¹, respectively. The spatial pattern of emissions is broadly consistent with state-of-the-art BVOC emission models (e.g., Arneth et al., 2008; Guenther et al., 2012; Messina et al., 2016; Müller et al., 2008; Sindelarova et al., 2014; Stavrakou et al., 2009; Wiedinmyer et al., 2006; Young et al., 2009).

4.6 Aerosols

Figure 26 shows the UKESM1 ensemble annual mean AOD at 440 nm for the period 1992 to 2008. Annual mean observations of AOD covering the same period from the ground-based, global AERONET (Holben et al., 2001) are overlaid on the simulated global AOD map. Overall, the model captures the spatial distribution of AOD well. This is reflected in the good correlation found in the point-wise comparison shown in Figure 26 (r²=0.63). AOD over continental polluted regions and tropical biomass regions is generally well-represented in the model, including the east-west contrast in continental North America. There are fewer observations in remote marine regions, but where observations exist, the model is able to capture the low AOD values in the Pacific Ocean and downwind of anthropogenic pollution source regions in the Atlantic Ocean.

The model AOD in Figure 26 is biased low at sites close to dust sources in West Africa and the tropical Atlantic. This is examined in more detail in Figure 27 which compares the model against seasonal mean AERONET AOD at dust-dominated sites and near-surface dust concentrations from the University of Miami observing network. The AODs in this figure are consistent with the annual mean result, with seasonal mean AOD lower than the observations by 39% on average. We note that these stations are all at similar latitudes in a region strongly affected by the West African Monsoon, which is challenging to simulate, and small biases in wind-speed or precipitation in this area would strongly affect local dust production. Thus the bias in AODs may not be representative of the wider dust simulation. Indeed, dust concentrations over the Atlantic are generally well-simulated (Figure 27). Further afield, there is some high bias in concentration over the Pacific and low bias over the Southern Ocean.

4.7 Emergent Behavior of the Coupled Model

We consider some emergent properties of UKESM1 using metrics of global climate sensitivity and carbon uptake. First, equilibrium climate sensitivity (ECS) estimates the equilibrium response of global mean surface air temperature to a doubling of CO₂. While this cannot be directly evaluated against observations, it is a useful metric for comparing the response of models to CO₂ perturbations across multiple generations of models. Figure 28 shows the temporal evolution of surface air temperature in response to an abrupt CO₂ quadrupling (the CMIP6 abrupt4xCO2 experiment; Eyring et al., 2016). We calculate effective climate sensitivity as in Gregory et al. (2004) by regressing the change in TOA net radiation against surface temperature over the first 150 years of the run. The projected equilibrium temperature change from the abrupt4xCO2 experiment is divided by 2 to correspond to a CO₂ doubling. The resulting ECS for UKESM1 is 5.4 K.

Transient climate response (TCR) is calculated from the CMIP6 1pctCO2 experiment in which CO₂ is increased at 1% per year. TCR is defined as the temperature change from the pre-industrial control after CO₂ has doubled, averaged over years 61–80. A four-member initial condition ensemble for UKESM1 (Figure 28) have TCR in the range 2.68–2.85 K. Transient climate response to cumulative emission (TCRE; Gillett et al., 2013) extends the TCR concept to use carbon feedbacks from an ESM in order to project the warming due to a cumulative CO₂ emission of 1 TtC. Following Gillett et al. (2013) we calculate TCRE as TCR (as defined above) divided by the cumulative compatible anthropogenic emission up to year 70 of the same 1pctCO2 experiment. CO₂ concentration is prescribed in these experiments, but the anthropogenic emission compatible with the model's carbon uptake can be calculated as the change in the total mass of carbon in the ocean, land, and atmosphere. To remove the (very small) background drift in ocean and land carbon stores (Figure 3) we calculate the compatible emission as the difference between the 1pctCO2 store at year 70 and that of pre-industrial control in the same year, yielding a range of compatible emissions from 1.07 to 1.08 TtC and TCRE from 2.49 to 2.66 K TtC⁻¹.

The UKESM1 values of ECS, TCR, and TCRE are all higher than those of CMIP5 models (respectively 2.1 K to 4.7 K, 1.0 K to 2.6 K, and 0.8 K TtC−1 to 2.4 K TtC−1; Andrews et al., 2012; Gillett et al., 2013). Elsewhere in this special issue, Bodas-Salcedo et al. (2019) analyze the increase in atmospheric climate feedbacks in HadGEM3-GC3.1 relative to the previous version of HadGEM3, whose ECS (3.2 K; Senior et al., 2016) is within the range of CMIP5 models. They find that the feedbacks have become more positive as a result of improvements to cloud microphysics and cloud-aerosol interactions. Both of these developments reduce errors in the present-day simulation of clouds relative to observations. These stronger feedbacks contribute to increases in all three of these metrics; TCRE is additionally increased relative to CMIP5 models by the inclusion of nitrogen limitation which reduces the terrestrial carbon sink for a given CO₂ atmospheric concentration and hence increases the CO₂ concentration and corresponding warming compatible with a 1TtC cumulative emission. More detailed analysis of the climate sensitivity and carbon feedbacks in UKESM1 will be reported in subsequent papers.

Carbon uptake can also be diagnosed from historical, rather than idealized, simulations and thus compared against estimates rooted in historical observations. In Table 4 we compare cumulative 1850–2014 carbon uptake and land use emissions against observation-based estimates from the Global Carbon Budget (GCB; Le Quere et al., 2018). The ocean and land sinks are both within the GCB range. UKESM1 land use emissions in Table 4 represent the carbon flux due to clearance of above-ground biomass, including crop harvest and deforestation. The model's total emission due to land use includes an additional contribution from decomposition of below-ground biomass which we have not diagnosed. The total land use emissions will hence will be closer to, or above, the GCB central estimate.

Surface temperature is one of the few variables for which reliable observations cover the full period of the 1850–2014 historical simulation, which allows us to evaluate the model's first-order climate response to the evolving forcing over this period. The UKESM1 global mean surface temperature anomaly in the historical ensemble is shown in Figure 29, alongside the HadCRUT4 observation data set (Morice et al., 2012). The observations represent only a single realization of the internal variability of the climate system, so one should not expect the model ensemble to be centered on the observations, but rather that the range of observational uncertainty overlaps the ensemble range (under the assumption that the model ensemble is large enough to sample the relevant internal variability). Most ensemble members begin to warm in the early 20th century, then cool strongly between 1950 and 1970 before warming rapidly through to the end of the simulation. The observations show a limited cooling of 0.1 to 0.2 K during 1940 to 1970, but the model ensemble mean cools by nearly 0.4 K over the same period. All ensemble members also show a stronger cooling response to the large volcanic eruptions of 1883, 1963, and 1991 than is seen in the observations.

Figure 29 also separates the mean surface temperature into northern and southern hemisphere time series. The stronger-than-observed cooling is restricted to the northern hemisphere which, together with its temporal evolution, points to either aerosol or land use forcing as the prime driver. Further investigation into this discrepancy will be the subject of future work. In the southern hemisphere the model ensemble overlaps with the observational uncertainty for the entire duration of the experiment, with the exception of the dip following the Mt. Pinatubo eruption in 1991 noted above.