Assessment of pre-industrial to present-day anthropogenic climate forcing in UKESM1

. Quantifying forcings from anthropogenic perturbations to the Earth System (ES) is important for understanding changes in climate since the pre-industrial period. In this paper, we quantify and analyse a wide range of present-day (PD) 15 anthropogenic climate forcings with the UK’s Earth System Model (ESM), UKESM1, following the protocols defined by the Radiative Forcing Model Intercomparison Project (RFMIP) and the Aerosol and Chemistry Model Intercomparison Project (AerChemMIP). In particular, by quantifying effective radiative forcings (ERFs) that include rapid adjustments within a full ESM, it enables the role of various climate-chemistry-aerosol-cloud feedbacks to be quantified. 2014) pre-industrial 1850) tropospheric ozone precursors (near-term climate forcers, NTCFs) exert a global mean ERF of -1.12 W m -2 , mainly due to changes in the cloud radiative effect. There is also a negative PD ERF from land use (-0.32 W m -2 ). It is outside the range of previous estimates, and is most likely due to too strong an albedo response. In combination, the net anthropogenic ERF is 35 potentially biased low (1.61 W m -2 ) relative to other estimates, due to the inclusion of non-linear feedbacks and ES interactions. By including feedbacks between greenhouse gases, stratospheric and tropospheric ozone, aerosols, and clouds, some of which act non-linearly, this work demonstrates the importance of ES interactions when quantifying climate forcing. It also suggests that rapid adjustments need to include chemical as well as physical adjustments to fully account for complex ES interactions.


Introduction
In order to have a quantitative understanding of past and future climate change, and attribute it and its impacts to different anthropogenic and natural drivers, it is necessary to have a detailed process-based understanding of all aspects of the pathway from anthropogenic (or natural) activity through to climate response and impacts. A key mechanism for addressing this overarching objective is the 6th Coupled Model Intercomparison Project (CMIP6; Eyring et al., 2016), which designs and 45 distributes data from multi-model simulations. These simulations, with state-of-the-art climate models or Earth System Models (ESMs), are aimed at addressing these important climate science questions directly or via dedicated CMIP6-endorsed model intercomparison projects (MIPs) such as the Aerosol and Chemistry Model Intercomparison Project (AerChemMIP; Collins et al., 2017). The first part of this cause-effect chain from activity to climate response, mediated through the atmosphere and the land surface, is quantifying changes to the Earth's radiation budget, often termed radiative forcing. 50 Successive assessment reports of the Intergovernmental Panel on Climate Change (IPCC) have used the concept of radiative forcing (RF) as a metric to quantify the effects of different anthropogenic and natural drivers on the Earth's radiation balance.
For this purpose, radiative forcing (RF), or more precisely, the stratospherically-adjusted radiative forcing (SARF) is defined at the tropopause (Myhre et al., 2013a) as: 55 where IRF is the instantaneous radiative forcing and Astrattemp is the additional change in the net downward radiative fluxes at the tropopause solely due to stratospheric temperature adjustment (Hansen et al., 1997), while holding all other variables fixed. 60 Including the stratospheric temperature adjustment can significantly affect the magnitude of a forcing (e.g. Smith et al., 2018) and even change the sign of the forcing in the case of ozone depletion (Shine et al., 1995). As a result, RF rather than IRF is a better predictor of the drivers of global mean temperature. RF is part of a framework, based around energy budget analyses, that has split forcing from climate response (Boucher et al., 2013;Myhre et al., 2013a;Sherwood et al., 2015). It's been used https://doi.org/10.5194/acp-2019-1152 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License. extensively to evaluate and compare the strength of the various mechanisms, both anthropogenic and natural, affecting the 65 Earth's radiation balance and hence, their contribution to climate change (e.g. Hansen et al., 1997;Shine and Forster, 1999).
However, despite the extensive use of RF as a metric for climate change, it is often calculated inconsistently between the different drivers of climate change (e.g. Myhre et al., 2013a) and participating models in the 5th Coupled Model Intercomparison Project (CMIP5; Taylor et al., 2012) can show large differences in forcing (Andrews et al., 2012a;Forster et al., 2013). For example, the RF attributed to long-lived greenhouse gases (LLGHGs) is typically based on changes in observed 70 concentrations between the pre-industrial (PI) and the present day (PD) and using line-by-line radiative transfer calculations and/or simple, yet justified, expressions for RF based on Myhre et al. (1998) and Ramaswamy et al. (2001). These expressions have been updated for some LLGHGs recently (Etminan et al., 2016). However, the observed concentrations themselves may be subject to biogeochemical feedbacks (e.g. Arneth et al., 2010;O'Connor et al., 2010) and observational-based estimates of forcing are also relatively uncertain (Skeie et al., 2011). 75 In contrast to the quantification of the LLGHG RF, the RF from tropospheric ozone (O3) changes at the PD relative to the PI (Stevenson et al., 2013) is based solely on models. It has been calculated using the ensemble of models (Young et al., 2013) participating in the Atmospheric Chemistry and Climate Model Intercomparison Project (ACCMIP; Lamarque et al., 2013) providing input to offline radiative transfer models (e.g. Edwards and Slingo, 1996). The simulations used the corresponding 80 sea surface temperatures (SSTs) and sea ice (SI) conditions for the time periods of interest (PI and PD), including some climate response and feedbacks at the PD, meaning that the resulting RF doesn't fit into the simple forcing-feedback framework, whereby feedbacks are related to global mean temperature change and forcings are not (Sherwood et al., 2015). There are also additional uncertainties associated with the estimate of tropospheric O3 RF due to the lack of a robust and reliable observational constraint for PI O3 concentrations (e.g. Stevenson et al., 2013), the diversity in modelled PD tropospheric O3 burden across 85 multi-model ensembles (e.g. Young et al., 2013;Young et al., 2018), uncertainties in historical emissions of tropospheric O3 precursors, and the apparent inability of current state-of-the-art chemistry models to replicate near-recent observed trends in tropospheric O3 (Parrish et al., 2014;Young et al., 2018) although recently, isotopic measurements seem to corroborate the modelled trends (Yeung et al., 2019). Additional uncertainties in Stevenson et al. (2013) arise from neglecting the change in O3 in the lower stratosphere attributable to changes in tropospheric O3 precursors and the contribution from stratospheric O3 90 depletion on the modelled changes in tropospheric O3. Despite these uncertainties in tropospheric O3 RF, the even larger uncertainty in aerosol forcing (Myhre et al., 2013a;Bellouin et al., 2019) accounts for the majority of the uncertainty in the total anthropogenic forcing.

100
Although RF has been calculated inconsistently between the different drivers of climate change to date, the fifth assessment report (AR5) of the IPCC also adopted an extension to the definition of radiative forcing to include fast feedbacks or rapid adjustments other than the stratospheric temperature adjustment. This updated definition, the effective radiative forcing (ERF), is now the preferred metric of choice for ranking the drivers of climate change. ERF is defined at the top of the atmosphere (TOA), following Chung and Soden (2015), as: 105 where IRF is the instantaneous radiative forcing at the TOA following a perturbation andAi is a rapid adjustment in the atmosphere or over the land that alters the net downward radiative flux at the TOA either positively or negatively. These 110 adjustments include changes in stratospheric temperatures (as included in the definition of RF or SARF above; Eqn. 1) as well as adjustments such as tropospheric temperatures, water vapour, clouds, and land surface albedo, as examples, but global mean surface temperatures or global ocean conditions remain unchanged. Although the error associated with an ERF tends to be larger than that of RF (or SARF), climate sensitivity parameters (i.e. the degree of warming per unit forcing) are less dependent on the forcing agent and it is more representative of the climate response than the traditional RF . For 115 example, rapid adjustments lead to the ERF for black carbon (BC) being half that of its IRF (Stjern et al., 2017;Smith et al., 2018). On the other hand, Xia et al. (2016) found that cloud and sea ice feedbacks driven by stratospheric O3 recovery included in the ERF definition lead to a significant adjustment that changes the sign and magnitude of the stratospheric O3 forcing from the SARF. Although the definitions of RF and ERF differ in where the change in radiative fluxes is diagnosed, ERFs at the tropopause are nearly identical to those at the TOA for tropospheric forcing agents; this is also the case for SARFs. An overview 120 https://doi.org/10.5194/acp-2019-1152 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License.
of the historical evolution of the RF concept, its quantification for different forcing agents, and applications of RF can be found in Ramaswamy et al. (2019).
The use of ERF as a metric also offers the advantage that it can be readily calculated using a pair of parallel simulations with standard model TOA radiative flux diagnostics  albeit with a requirement to run for relatively long periods 125 (~30 years) to reduce the uncertainty associated with meteorological variability (e.g. Shindell et al., 2013a). Despite this requirement, it has now become the metric of choice over RF (Boucher et al., 2013) to rank the drivers of climate change. Figure 1 illustrates the various definitions of radiative forcing, including those of IRF, SARF (or RF), and ERF.
The aim of the current study is to quantify the PI (Year 1850) to PD (Year 2014) effective radiative forcings (ERFs) from 130 anthropogenic drivers of climate change with a fully coupled Earth System Model (ESM). Using the experimental protocol recommended for the Radiative Forcing Model Intercomparison Project (RFMIP; Pincus et al., 2016), PD forcings will be quantified relative to the PI from changes in emissions, concentrations, and/or land use due to anthropogenic activities. This approach offers a consistent methodology for diagnosing the ERFs of all forcing agents or Earth System (ES) perturbations and, if applied consistently across all models as part of the 6th Coupled Model Intercomparison Project (CMIP6; Eyring et al., 135 2016), should help to address some of the deficiencies and uncertainties associated with previous estimates of forcing (e.g. Myhre et al., 2013a) and improve our understanding of how the ES responds to forcing . The paper is organised as follows. Section 2 provides a brief description of the UK's Earth System Model, UKESM1, used in this study.
Section 3 outlines the experimental setup and the simulations carried out. PD anthropogenic forcings relative to PI are presented in Sect. 4, while conclusions can be found in Sect. 5. 140

Model Description
The model used in this study is the atmospheric and land components of the UK's Earth System Model, UKESM1 . UKESM1 is based on the Global Atmosphere 7.1/Global Land 7.0 (GA7.1/GL7.0; Walters et al., 2019) configuration of the Hadley Centre Global Environment Model version 3 (HadGEM3; Hewitt et al., 2011), herein referred to as HadGEM3-GA7.1, to which terrestrial carbon/nitrogen cycles  and interactive stratosphere-troposphere 145 chemistry (Archibald et al., 2019) from the UK Chemistry and Aerosol (UKCA; Morgenstern et al., 2009;O'Connor et al., 2014) model have been coupled. The model resolution is N96L85; this is equivalent to a horizontal resolution of roughly 135 km and the 85 terrain-following model levels cover an altitude range up to 85 km above sea level. The physical atmosphere model, HadGEM3-GA7.1, already includes the UKCA prognostic aerosol scheme called GLOMAP-mode (Mann et al., 2010;Mulcahy et al., 2018), in which secondary aerosol formation is determined by prescribed oxidant fields. Here, the UKCA 150 chemistry and aerosol schemes are coupled, with oxidants from the stratosphere-troposphere chemistry scheme (Archibald et al., 2019) determining secondary aerosol formation rates. A full description and evaluation of the UKCA chemistry and aerosol https://doi.org/10.5194/acp-2019-1152 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License. schemes in UKESM1 can be found in Archibald et al. (2019) and Mulcahy et al. (2019), respectively. Mulcahy et al. (2018) also implemented a number of aerosol process improvements in HadGEM3-GA7.1 which helped to address the strong negative PD aerosol forcings from its predecessor model (i.e. HadGEM3-GA7.0; Walters et al., 2019). 155

Model Setup and Experiments
To calculate the pre-industrial (PI; Year 1850) to present-day (PD; Year 2014) effective radiative forcings (ERFs) due to a PIto-PD perturbation (e.g. change in emissions), timeslice experiments with fixed sea surface temperatures (SSTs) and sea ice (SI) were carried out, following the protocol defined by the Radiative Forcing Model Intercomparison Project (RFMIP; Pincus et al., 2016). The experimental setup was also consistent with recommendations from Forster et al. (2016)  simulation, only the land surface sees the PI-to-PD perturbation in CO2, the radiation scheme still sees the PI concentration. b The AerChemMIP experiment piClim-NTCF is also known as piClim-aerO3 in RFMIP.

170
Effectively, this involves running a PI timeslice experiment, called piClim-control here, in which SSTs and SI and all forcings were fixed at year-1850 levels. The SSTs and SI used in piClim-control were monthly time-varying climatologies derived from 30 years (i.e. Years 2156-2185 inclusive) of output from the UKESM1 pre-industrial coupled control experiment (piControl) characterised in  and one of an underpinning set of coupled experiments for CMIP6 (Eyring et al., 2016). It 175 also used 30-year monthly mean climatologies for the vegetation distribution, canopy height, leaf area index, surface seawater dimethyl sulphide (DMS) and chlorophyll concentrations derived from the same period of piControl. Fixing the vegetation distribution was not included in the RFMIP protocol and any potential vegetation rapid adjustments will be somewhat constrained. This is due to the simulations being based on the configuration of UKESM1 used for the Atmosphere Model Intercomparison Project (AMIP) simulation (which prescribed vegetation characteristics). Although the AerChemMIP 180 protocol (Collins et al., 2017) requested use of the maximum capability possible, interactive vegetation was not a model requirement. The extra RFMIP experiments carried out here were only done as a late addition and the same experimental setup was kept for internal consistency.
The model was initialised using output from the start of the 30-year period used to produce the PI climatologies (i.e. January 185 2156 of piControl). All the other experiments are perturbation experiments, parallel to piClim-control, in which selected https://doi.org/10.5194/acp-2019-1152 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License. emissions, concentrations and/or land use were changed from year-1850 to year-2014 values. Although AerChemMIP and RFMIP recommend 30 years for fixed SST timeslice ERF experiments, the perturbations to the LLGHGs took up to 15 years to propagate fully into the stratosphere due to the turnover timescale associated with the Brewer-Dobson circulation (e.g. Butchart, 2014). Therefore, all simulations were 45 years in length. Using the latter 30 years of the paired simulations, the 190 ERFs from the PI-to-PD perturbations were diagnosed as the time-mean global-mean difference in the TOA net radiative fluxes. Details on how the ERF was further decomposed can be found in Sect  (Arfeuille et al., 2014;Thomason et al., 2018;Matthes et al., 2017) using those specified for CMIP6 (Eyring et al., 2016). To test this, we created an ensemble of short runs on each machine by perturbing selected variables in their initial conditions, 215 using a perturbation with a numerical value comparable to the machines precision. The spread of results (at each point in time and space) on each platform was then used to determine whether they could each have been sampled from the same ensemble of results generated on either machine. A permutation method was used to ensure statistical independence between neighbouring points according to the work described by Wilks (1997). A paper, describing this protocol in more detail, is in preparation (Teixeira et al., 2019). Further to this, a number of perturbation experiments were carried out in duplicate to test 220 https://doi.org/10.5194/acp-2019-1152 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License. the sensitivity of the global mean ERFs to differences in HPC platform. These duplicate experiments are listed in Table 2

Present-Day Anthropogenic Effective Radiative Forcings
The ERF has been calculated from the difference ( ) in the net TOA radiative flux ( ) between the perturbed simulation (e.g. 230 piClim-CH4; Table 1) and the piClim-control simulation as follows: However, many of the experiments in this study either perturb aerosol emissions and/or alter aerosol concentrations via 240 chemical and dynamical feedbacks. Changes in aerosol can bias the diagnosed CRE as aerosol scattering and absorption typically reduce the contrast in shortwave (SW) reflection between cloudy and clear-sky scenes; a process termed "cloud masking" (e.g. Zelinka et al., 2014). In consideration of this, we have calculated the change in CRE from "clean" radiation calls that exclude aerosol-radiation interactions (ari), as recommended in Ghan (2013): The ERF is thus separated into a component due to cloud property changes ( ′) and the non-cloud forcing (ERFcs'). 250 Here, ERFcs' is the sum of the aerosol IRF and any non-aerosol changes in CS flux and differs slightly from ERFcs in equation (5), in that it can include the impact of aerosol scattering and absorption in the clear-air above or below clouds. One acknowledged limitation is that variations in gaseous absorption and emission between clear and cloudy scenes also lead to cloud masking effects (e.g. Soden et al., 2008). Although Ghan's method removes the very prominent influence of aerosols, cloud-masking from ozone (O3) and GHGs may still affect the separation of ERF into the CS and CRE components. 255

Overview of Present-Day Effective Radiative Forcings (ERFs)
The effective radiative forcing (ERF), clear-sky (CS), and cloud radiative effect (CRE) contributions, following equations (6) to (8) are listed in Table 3 for all perturbation experiments relative to piClim-control, and are further decomposed into the SW (solar), LW (terrestrial) and net (SW + LW) components. The global mean ERFs are also plotted in Fig. 2. Together, they show that the global mean PD ERF by greenhouse gases (GHGs) is +2.89 ± 0.04 W m -2 , which is offset by an aerosol ERF of 260 -1.13 ± 0.04 W m -2 . These estimates are consistent with ERFs of +3.09 and -1.10 W m -2 , respectively , from the HadGEM3 GC3.1 (Williams et al., 2017) physical model (herein referred to as HadGEM3-GC3.1) upon which UKESM1 is based and with the range of previous estimates (e.g. Myhre et al., 2013a). The net anthropogenic ERF is +1.61 W m -2 , again consistent with the range of estimates from AR5 (Myhre et al., 2013a) Figure 3 shows global distributions of the multi-annual mean PD ERFs from changes in GHGs, aerosols, tropospheric ozone 280 (O3) precursors, land use, and total anthropogenic sources since PI. It shows that the PD ERF from GHGs ( Fig. 3a) is strongly positive everywhere although there is some evidence of negative forcing (up to -2 W m -2 ) in the high latitudes, particularly in the southern hemisphere (SH). This is due to a negative ERF from the piClim-HC perturbation experiment (Table 3)  lifetime. However, the distribution of the aerosol forcing is still largely negative everywhere, except for over bright surfaces and regions where the positive forcing from black carbon (BC) emissions outweighs the negative forcing from scattering aerosols such as sulphate and organic carbon (OC). A breakdown of the aerosol forcing between constituents and between 290 aerosol-radiation interactions (ari) and aerosol-cloud interactions (aci) will be presented and discussed in Sect. 3.3. Figure 3c shows the global distribution of the PD ERF from emissions of tropospheric O3 precursors. It shows that the forcing is spatially heterogeneous, has regions of both weak positive and weak negative contributions, but on a global mean basis is weakly positive (+0.15 W m -2 ; Table 3) in comparison with other forcings (e.g. GHGs). Further analysis of the ERF from tropospheric O3 precursors and their contribution, along with methane (CH4), to forcing by tropospheric O3, can be found in 295 Sect. 3.4. The distribution of the land use ERF (Fig. 3d) is spatially heterogeneous although much of the negative forcing is concentrated over the northern hemisphere (NH) continental regions (e.g. North America, South East Asia). Together with aerosols, their combined ERF outweighs the positive ERF from GHGs, leading to a negative total anthropogenic ERF over much of the NH continents ( Fig. 3e). As was the case with the GHG forcing, the negative ERF over the SH high latitudes is still evident in the total anthropogenic forcing. 300 As discussed in Sect. 2, the UKESM1 fixed SST experiments were run on multiple High Performance Computing (HPC) platforms. Despite statistical methods ensuring that the model is not scientifically different on the different HPC platforms, such duplicate experiments still produce slightly different results. This raises the question of the impact of such differences on the quantification of PD ERFs with UKESM1. To address this question, the experiments described in Table 2 were used to 310 compare the difference in TOA radiative fluxes of equivalent realisation experiment pairs. The two-sample Kolmogorov-Smirnov test between the monthly TOA radiative fluxes from the realisation pairs (null hypothesis that the samples are drawn from the same distribution) shows that the TOA radiative fluxes between the two realisation pairs are statistically identical at a confidence level () of 5 %. Despite the fact that one cannot conclude that the distributions are identical, for each pair of experiments there is no evidence suggesting that the two distributions are different. 315

Carbon Dioxide (CO2)
Atmospheric carbon dioxide (CO2) concentrations have risen from 284.3 to 397.5 ppm between 1850 and 2014 (Meinshausen et al., 2017) resulting in a global mean ERF of +1.83 ± 0.04 W m -2 (Table 3) and is the largest individual contribution to the total historical forcing. Myhre et al. (2013a) report a stratospherically adjusted radiative forcing (SARF) or RF of +1.82 ± 0.19 325 W m -2 for the period from 1750 to 2011 which, coincidentally, has a near identical rise in CO2 of 113 ppm to that assessed here. An updated radiative forcing assessment based on line-by-line calculations (Etminan et al., 2016) increased SARF for 2015 relative to 1750 to 1.95 W m -2 but for a larger CO2 rise of 121 ppm. Our ERF estimate for 2014 is, therefore, consistent with the SARF from line-by-line calculations and confirms that for CO2, the SARF is a useful metric for evaluating the radiative https://doi.org/10.5194/acp-2019-1152 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License. effects of increasing CO2 . As expected for a GHG, the ERF is dominated by absorption in the clear-sky 330 (CS) longwave (LW) (+1.56 ± 0.02 W m -2 ) with additional CS shortwave (SW) forcing (+0.122 ± 0.02 W m -2 ) coming from a rapid adjustment in snow-cover across the northern latitudes and a net cloud radiative effect (0.144 ± 0.03 W m -2 ) from a reduction in cloud cover.
Rising atmospheric CO2 also exerts an indirect forcing through rapid changes in plant stomatal conductance (Doutriaux-335 Boucher et al., 2009;Richardson et al., 2018) enhancing plant water use efficiency and reducing evapotranspiration leading to an increase in sensible heating at the surface and corresponding drying of the boundary layer and reduction in low-clouds. This mechanism is known as a physiological forcing. In UKESM1 (piClim-CO2phys; Table 3), this forcing is small for the PD (0.008 ± 0.04 W m -2 ) but results from a balance between a negative CS LW component (-0.065 ± 0.03 W m -2 ) associated with surface warming and a positive SW cloud radiative effect (CRE) (0.112 ± 0.02 W m -2 ) associated with a reduction in low-level 340 cloud. The balance comes from smaller terms including a negative CS SW contribution from reduced water vapour. Previous Hadley Centre models at 4xCO2 found a much larger effect: 1.1 W m -2 in HadCM3LC (Doutriaux- Boucher et al., 2009) and 0.25 W m -2 in HadGEM2-ES (Andrews et al., 2012b). Assuming the physiological effect scales linearly with atmospheric CO2, we find that UKESM1 has a stronger CS LW component compared to HadGEM2-ES implying that UKESM1 has a more pronounced surface warming adjustment associated with the physiological effect. 345

Nitrous Oxide (N2O)
The global mean ERF due to changes in N2O from PI (Year1850 value of 273 ppbv) to PD (Year 2014 value of 327 ppbv) is calculated as +0.13±0.03 W m -2 (Table 3), following the equations (6)-(8) described above. The predominant contribution to the N2O ERF is the CS LW forcing (+0.25 ± 0.03 W m -2 ), with a small and non-significant cancellation by the CS SW forcing (-0.03 ± 0.02 Wm -2 ); the net CS forcing sums up to +0.22 ± 0.03 W m -2 . The net cloud-radiative effect (CRE) forcing is -0.10 350 ± 0.03 W m -2 , with the SW and LW contributions of +0.03 ± 0.03 W m -2 and -0.13 ± 0.02 W m -2 , respectively. In comparison, the net ERF calculated here (+0.13 ± 0.03 W m -2 ) is slightly lower than the SARF values of +0.17 ± 0.03 W m -2 from AR5 for 2011 (Myhre et al., 2013a) and of +0.18 W m -2 for 2014 based on the updated expression from Etminan et al. (2016). This is likely to be due to the effect of adjustments associated with changing N2O that weren't considered as part of the SARF in AR5 (Myhre et al., 2013a) or Etminan et al. (2016), including O3 depletion and fast cloud adjustments. Previously, Hansen et al. 355 (2005) calculated an IRF and a SARF of 0.15 W m -2 due to the change in N2O from 278 to 316 ppbv.
The global distribution of the individual components contributing to the N2O ERF (Fig. 4) shows that the SW forcing is negligible under CS conditions and is largely made up of noise in cloudy conditions which might be due to an aerosol effect.
CS LW forcing is predominantly positive, due mainly to N2O direct forcing. There are relatively large regional changes in the 360 LW due to clouds leading to a significant net negative LW CRE of -0.13 W m -2 . The positive LW CRE over the southern high latitudes appears to be correlated to the change in total column ozone (TCO) and the shift of the surface wind shown in Zonal mean ozone (O3) and temperature changes (Fig. 5) show the warming of the upper troposphere and lower stratosphere due to increased O3 in that region and the cooling of the stratosphere due to O3 depletion associated with the N2O increase. As a result, the atmospheric circulation has been modified, reflected in a shift of the upward velocity (w component of the wind) 365 ( Fig. 5d), which is correlated to the water content of the clouds. Such changes in clouds drive the LW forcing, which shows a pattern of significant regional changes with cancelling effects. The change in TCO is largely negative due to a reduction in stratospheric O3 (Fig. 5a). That would lead to a negative indirect 375 forcing by O3; indeed, the ERF calculated here is 0.05 W m -2 smaller than the direct N2O forcing in AR5 (Myhre et al., 2013a). increases there, resulting in a regional suppression of convection and a reduction of associated LW radiation coming from 380 cloud tops. Prather and Hsu (2010) also noted that coupling of N2O and methane (CH4), via stratospheric O3 depletion and photolysis, can reduce the global warming potential (GWP) of N2O by 5 %, based on PD concentrations. Here, the response of CH4 is constrained due to being concentration-driven. However, the whole-atmosphere CH4 lifetime changes from 8.1 years in piClim-390 control to 8.0 years in piClim-N2O, which would lead to a small offset to the N2O ERF of 0.01 W m -2 due to the direct radiative effect of CH4 if the model was emissions-driven; indirect CH4 forcings (e.g. O3) could enhance this offset even further .

Halocarbons
The global multi-annual mean ERF from the PI to PD change in ozone depleting substances (ODSs) is quantified, using piClim-395 HC relative to piClim-control, as -0.33 ± 0.04 W m -2 . A quantitative and process-based understanding of why the ERF is negative, despite our understanding to date of the relative roles of ODSs and O3 depletion in climate forcing, is presented in .

Methane (CH4)
The global mean methane (CH4) concentration changed from 808.3 ppbv in the PI (Year 1850) to 1831.5 ppbv in the PD (Year 400 2014), resulting in a PD ERF of +0.93±0.04 W m -2 (Table 3; Fig. 6), where the 0.04 W m -2 is the standard error following Forster et al. (2016). Most of the forcing (+0.73±0.02 W m -2 ) occurs in the CS LW, with an additional positive contribution (+0.12±0.01 W m -2 ) from the CS SW, consistent with the growing recognition of the importance of the SW absorption bands in CH4 forcing (Collins et al., 2006;Li et al., 2010;Etminan et al., 2016). There are additional SW and LW contributions from the cloud radiative effect (CRE) but these largely cancel out, leading to a small net positive forcing (+0.08 ± 0.03 W m -2 ) in 405 addition to the net CS forcing (+0.85±0.03 W m -2 ). Estimates of the direct CH4 ERF at the PD from HadGEM2 model simulations (Andrews, 2014) and the updated RF expression for CH4 based on line-by-line calculations (Etminan et al., 2016) are of the order of 0.50-0.56 W m -2 . But the forcing calculated here is higher by more than 0.4 W m -2 . However, it is consistent with other studies (e.g. Hansen et al., 2005;Shindell et al., 2009;Myhre et al., 2013a), who concluded that the total climate forcing by CH4 is almost double that of the direct forcing due to indirect effects. The UKESM1 estimate is also larger than the 410 +0.69 W m -2 radiative impact of an increase in CH4 concentration of 1800 ppbv above present-day levels quantified by

Total Greenhouse Gases (GHG)
The major driver of anthropogenic climate change is greenhouse gases (GHGs), which is offset by aerosols (Myhre et al., 2013a). Therefore, the total greenhouse ERF and the total aerosol ERF are key metrics in understanding observed and modelled changes in the climate system since the PI. As a result, a separate timeslice simulation with all GHGs (piClim-GHG; Table 1) was conducted following the RFMIP protocol .

425
The UKESM1 piClim-GHG experiment leads to an ERF of 2.89 ± 0.04 W m -2 (Table 3), which is dominated by a positive CS LW contribution (3.08 W m -2 ) that is partially offset by a negative CS SW component (-0.16 W m -2 ). There are significant positive and negative contributions from the CRE in the SW and LW but these largely cancel out and contribute little (-0.01 W m -2 ) to the total ERF (Table 3). This PD GHG ERF is lower than the ERF of 3.09 W m -2 estimated from the physical model The UKESM1 1850-2014 estimate, however, is consistent with the AR5 year-2011 SARF estimate of 2.82 W m -2 . The latter estimate has been adjusted to an 1850 baseline from 1750, and taking stratospheric O3 depletion, CH4-driven SWV, and half of the tropospheric O3 forcing (based on the attribution by Stevenson et al., 2013) into account, although some GHG 435 concentrations (e.g. CH4) have increased between 2011 and 2014 (Nisbet et al., 2016) while others (e.g. ODSs) have declined (Engel et al., 2018).
Although, in principle, one could use the combination of GHG simulations to test for linearity, this isn't feasible here; piClim-GHG perturbs non O3-depleting halocarbons whereas piClim-HC does not. However, there is a discrepancy in ERF of 0.35 W 440 m -2 between piClim-GHG and the sum of the individual GHG simulations which cannot be fully accounted for by the small positive RF from non O3-depleting halocarbons (HFCs, PFCs, SF6) of 0.02 W m -2 at 2011 (Myhre et al., 2013a). This suggests that there is a potential non-linearity in the GHG forcing, perhaps involving O3 chemistry and aerosol-mediated cloud feedbacks, which may warrant further investigation. Figure 7a summarizes the results from the anthropogenic aerosol experiments, including a breakdown of the ERF into instantaneous radiative forcing (IRF) and rapid adjustments (RA) for each of the anthropogenic aerosol experiments (piClim-SO2, piClim-OC, piClim-BC, and piClim-aer). The IRF was calculated using the double-call system where aerosol-radiation interactions (ari) (i.e. scattering and absorption) are withdrawn from the second call to the radiation scheme, as in Ghan et al. (2012). The rapid adjustments (RA) have then been derived as the residual between the ERF and IRF and includes all aerosol-450 induced changes in cloud radiative effects.

Aerosols and Aerosol Precursors 445
Uncertainties in aerosol ERF are likely to be driven by uncertainties in the PI aerosol state (Carslaw et al., 2013) as well as the oxidising capacity of the atmosphere (Karset et al., 2018). The additional interactive sources of natural aerosol in UKESM1 for marine dimethyl sulphide (DMS), terrestrial and marine biogenic sources, along with the inclusion of a fully interactive 455 chemistry scheme (Archibald et al., 2019) are generally found to improve the evaluation of PD aerosol in UKESM1 . This provides some confidence in the underlying physical processes driving the PI aerosol state in this model.
The total anthropogenic aerosol ERF evaluated from the "all" (piClim-aer) experiment is -1.13 W m -2 , which is almost identical to the ERF of -1.10 W m -2 from the physical model HadGEM3-GA7.1  and slightly lower in magnitude 460 than the ERF of -1.45 W m -2 derived from HadGEM3-GA7.1 with CMIP5 emissions (Mulcahy et al., 2018). The estimate fits well within the likely range of -1.60 to -0.65 W m -2 (16-84 % confidence level) provided by a recent major assessment of PI to PD aerosol ERF (Bellouin et al., 2019) and the -1.5 to -0.4 W m -2 likely range previously assessed by AR5 (Myhre et al., 2013a). The annual-mean global distribution of aerosol ERF, IRF and aci are shown in Fig. 8. The total aerosol ERF (Fig. 8a) is negative over most regions, and is strongest over the cloudy ocean regions of the Northern Hemisphere where the forcing is 465 dominated by the RA term (Fig. 7a) driven by aerosol-cloud interactions (aci) (Fig. 8c). Some areas of Asia and North Africa have positive total aerosol ERF due to BC-rich aerosol loadings that give locally positive aerosol IRF (Fig. 8b). However, globally, the IRF is rather small and negative (-0.16 W m -2 ) (Fig. 7a) due to scattering by sulphate and OC that is partially offset by absorption from BC.   (Table 3, Fig. 7a). The IRF component of the sulphate ERF is -0.48 W m -2 , which agrees well with the best estimate from AR5 (-0.40 ± 0.2 W m -2 ) and from AEROCOM II (-0.58 to -0.11 W m -2 ) 480 (Myhre et al., 2013b). The BC ERF was +0.32 W m -2 , coming mostly from the IRF (0.38 W m -2 ) and small negative offset of -0.06 W m -2 from rapid adjustments. As noted in Johnson et al. (2019), BC absorption leads to strong cloud adjustments but the SW and LW components of these almost cancel in HadGEM3-GA7.1 and UKESM1 (see Table 3). This contrasts with many other models where the combination of low cloud enhancements and reductions in upper-level cloud typically result in more substantial negative adjustments, making the BC ERF on average about half the magnitude of the IRF (Stjern et al., 485 2017). The BC ERF given by UKESM1 is however well within the range assessed by AR5 (0.05 -0.8 W m -2 ), which took into consideration the possibility that BC emissions and/or absorption efficiency were underestimated in CMIP5 models (Bond et al., 2013). The anthropogenic emissions of BC were specified as 5 Tg yr -1 in CMIP5 (Year 2000 as PD) but have increased to 8 Tg yr -1 in CMIP6 (Year 2014 as PD). With CMIP5 emissions, the BC ERF from HadGEM3-GA7.1 was found to be 0.17 W m -2  and is comparable to direct BC forcing from other CMIP5 model estimates (Myhre et al., 2013a). 490 It is worth noting that the aerosol absorption was in fairly good agreement with AERONET observations in HadGEM3-GA7.1 simulations that used the CMIP5 emission set (Mulcahy et al., 2018). The slightly higher CMIP6-based estimate of 0.32 W m -2 provided in the present study could therefore be an overestimate, although this is difficult to judge given the uncertainties in comparing models with absorption measurements. It is also worth noting that BC forcing can have proportionately larger impacts on surface radiation, precipitation and regional circulation due to the localised absorption of solar radiation in the 495 atmosphere.

500
The PD OC ERF relative to PI is -0.27 W m -2 , with -0.16 W m -2 from the IRF and -0.11 Wm -2 from RA (Fig. 7a, Table 3). The OC IRF estimate agrees well with AR5 that assessed the RF of primary and secondary OC to be -0.12 W m -2 (-0.4 to +0.1).
Note that OC is non-absorbing in UKESM1 and, as such, neglects the role of brown carbon. This potentially misses a small positive contribution to the aerosol forcing (e.g. Feng et al., 2013), although biomass burning emissions are the dominant global source of brown carbon and these do not change significantly from 1850 to 2014 (Saleh et al., 2014). The contribution 505 of brown carbon to aerosol ERF is, therefore, likely to be small compared to the overall uncertainty in modelling aerosol absorption by BC (Bond et al., 2013). Nitrate aerosols are also not represented in HadGEM3 emissions could not be evaluated here. The total aerosol ERF in UKESM1 would presumably be more strongly negative if the role of nitrate aerosol was included (e.g. Bellouin et al., 2011) and this should be borne in mind when making comparisons 510 with other models or estimates based on observational constraints.
The "all" aerosol forcing experiment combines the increases in BC, OC and SO2 emissions together in one simulation, and interestingly, the ERF is 0.27 Wm -2 weaker (less negative) than the sum of the ERFs from the experiments that perturb those emissions separately. The main reasons for this lack of linearity are: (i) Cloud droplet numbers do not increase linearly with 515 aerosol loading (e.g. Jones et al., 1994) and begin to saturate in the "all" experiment (piClim-aer), which means OC and BC emissions no longer contribute significantly to aci once co-emitted with year-2014 levels of SO2. (ii) The absorption of upwelling shortwave radiation by BC is enhanced by increases in aerosol scattering and cloud brightness due to aci, making the IRF less negative in the "all" experiment.

520
To further understand which processes contribute most to total aerosol ERF, a series of additional sensitivity experiments were conducted with aci processes selectively disabled. In these tests, the cloud droplet number concentrations (CDNCs) used for https://doi.org/10.5194/acp-2019-1152 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License. the calculation of cloud droplet effective radius (Reff) and/or autoconversion were prescribed via a 3D monthly-mean PI climatology constructed from the final 30 years of the piClim-control simulation. The resulting global-mean ERFs are summarized in Fig. 7b. In one pair of simulations, CDNCs were prescribed for the Reff calculation to disable the so-called 525 Twomey effect (Twomey, 1977). By comparison with the main piClim-aer experiment, this indicated a Twomey effect (ACI_Twomey) of -0.73 W m -2 . A similar pair with CDNCs prescribed only for the autoconversion process led to an estimate of -0.34 W m -2 for the cloud lifetime effect (ACI_lifetime) (Albrecht, 1989). By prescribing CDNCs for both Reff and autoconversion both microphysical aci processes are disabled and only aerosol-radiation interactions (ari) are included. A pair of simulations with this setup provided an estimate for ERFari of -0.14 W m -2 . The change in CRE in the ari-only experiment 530 was only +0.04 W m -2 , which is not statistically significant at the 95% confidence level and indicates that the semi-direct aerosol effect is small or approximately neutral in this model. To complete the breakdown, the method in Ghan (2013) was further used to derive the contribution from changes in "clean" (aerosol-free) CS radiation (RA_cs). This term was found to be +0.05 W m -2 , and arises due to changes in surface albedo and atmospheric temperature and humidity. This overall breakdown suggests an aerosol-cloud forcing of -1.03 W m -2 , with a roughly 70/30 split between the Twomey effect and 535 aerosol effects on cloud cover and water content (lifetime & semi-direct effects).
The estimated aerosol-cloud forcing sits well within the 90 % likelihood ranges assessed by AR5 (-1.2 to 0.0 W m -2 ) and of ERFari of -0.14 W m -2 is within the 5-95% confidence range (-0.45 +/-0.5 W m -2 ) from AR5 but slightly weaker than the

Tropospheric ozone (O3) precursor gases 550
The global mean ERF from emissions of tropospheric ozone (O3) precursors is weakly positive (+0.15 W m -2 ; Table 3), although spatially heterogeneous, with regions of positive and negative forcing. Tropospheric O3 precursor emissions affect the ERF both through changes to tropospheric O3, a GHG, and by changing tropospheric oxidants such as OH, which in turn affect aerosols (Karset et al., 2018) and CH4. Here, we explore the composition and forcings resulting from O3 precursor emissions changes between the PI and PD by comparing the piClim-O3 simulation with piClim-control and the separate effects 555 https://doi.org/10.5194/acp-2019-1152 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License. of nitrogen oxide (NOx) and volatile organic compound (VOC)/carbon monoxide (CO) emissions changes using the piClim-NOx and piClim-VOC simulations, respectively.

Tropospheric Ozone (O3) Changes
The multi-annual mean tropospheric O3 column for the piClim-control simulation, and tropospheric O3 column differences between piClim-control and the piClim-O3, piClim-NOx and piClim-VOC simulations are shown in Fig. 9. O3 levels in the troposphere are the result of competing production and loss processes. Production occurs in the presence of NOx and VOC/CO, with most regions of the troposphere being NOx-limited. PI tropospheric column O3 values (Fig. 9a) show maxima over central 565 Africa and the South Pacific, the result of relatively large emissions of NOx from soil, biomass burning and lightning in this region, and minima over regions remote from NOx sources, such as the Western Pacific (where O3 loss is efficient) (Young et al., 2018). O3 columns are also low over regions of high surface elevation where the atmospheric column is shallower.

Tropospheric O3 column values increase when the model is perturbed with the larger PD O3 precursor emissions in piClim-570
O3, with the largest increases in the northern hemisphere, particularly in southern and eastern Asia. The tropospheric O3 burden increases from 280.1 Tg in piClim-control to 356.1 Tg in piClim-O3. As can be seen in Fig. 9, the dominant driver of these changes is increased NOx emissions, and there is a similar pattern of changes in the piClim-O3 and piClim-NOx simulations, although the change in O3 is smaller in the latter; the burden in piClim-NOx is 335.7 Tg. Small increases in O3 are observed in piClim-VOC, with some hotspots in regions such as South East Asia and the O3 burden increases to 297.2 Tg. 575 The tropospheric O3 column difference between piClim-O3 and piClim-control of 76.0 Tg is slightly larger than the sum of the individual O3 burden changes in piClim-NOx (55.6 Tg) and piClim-VOC (17.1 Tg) which total 72.7 Tg. However, these https://doi.org/10.5194/acp-2019-1152 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License. differences due to the non-linear nature of tropospheric chemistry are small, on the order of < 10 %. Similar behaviour is seen in the patterns of the tropospheric O3 column difference between piClim-control and other experiments. 580

Tropospheric Ozone (O3) Radiative Forcing
As seen from earlier sections, it can be difficult to compare ERFs estimated from these experiments against estimates of RF (e.g. Stevenson et al., 2013;Myhre et al., 2013a), due to the inclusion of indirect forcings and/or rapid adjustments other than the stratospheric temperature adjustment although Shindell et al. (2013b) noted that for O3, SARF and ERF estimates are comparable. Nevertheless, in order to compare against Stevenson et al. (2013), we estimate SARF for the tropospheric O3 585 precursor experiments by adopting a radiative kernel approach (e.g. Soden et al., 2008). This involves applying the tropospheric O3 radiative kernel from Rap et al. (2015) to the diagnosed change in tropospheric O3 (using the 150 ppbv O3 isoline in piClimcontrol as a tropospheric mask, as used in Young et al., 2013, Stevenson et al., 2013, and Rap et al., 2015 to calculate a stratospherically-adjusted radiative forcing (SARF) or RF. In this way, we can directly compare against the best estimate of the 1850-2010 tropospheric O3 RF of 364 mW m -2 by Stevenson et al. (2013) and quantify the contribution of different 590 tropospheric O3 precursors (including methane) to the PI-to-PD change in tropospheric O3 and its RF. We can also compare with more recent estimates from Checa-Garcier et al. (2018) and Yeung et al. (2019). Figure 10 shows the global distribution of the PD tropospheric O3 RF from piClim-CH4, piClim-NOx, and piClim-VOC experiments and their sum using the kernel approach. It shows that the tropospheric O3 RF is strongest over the northern 595 hemisphere tropics and weakest over the southern hemisphere high latitudes. The strongest RF occurs in regions of warm surface temperatures and high albedo, coinciding with the largest tropospheric O3 change (Shindell et al., 2013a). As was the case in Stevenson et al. (2013), the forcing is weaker over regions of high altitude (e.g. Tibetan Plateau) due to there being less O3 column aloft to absorb in the longwave (LW). On a global annual mean basis, the tropospheric O3 RF from the kernel approach is 407 mW m -2 , higher than the best estimate of Stevenson et al. (2013) by approximately 11 %. However, the 600 Edwards and Slingo (1996) radiation scheme, which was used to generate the kernel, can lead to forcing estimates that are higher than other comparable radiation schemes (Stevenson et al., 2013). Other potential differences are due to experimental setup. The Atmospheric Chemistry and Climate Model Intercomparison Project (ACCMIP; Lamarque et al., 2013) simulations were timeslices and included historical climate change although Stevenson et al. (2013) found that to have a minimal impact on the RF. More importantly, the increase in tropospheric O3 at the PD (and its resulting RF) were offset by decreases due to 605 the downward transport of O3-depleted stratospheric air (e.g. Shindell et al., 2013a). Applying the kernel method to the diagnosed decrease in tropospheric O3 from the piClim-HC experiment (Sect. 3.2.3), we find a RF offset of -101 mW m -2 , larger in magnitude by nearly a factor of 2 than the estimates from Søvde et al. (2011) andShindell et al. (2013a), thus reducing our original estimate to 306 mW m -2 . This revised estimate is 15 % lower than the Stevenson et al. (2013)  In considering the attribution of the tropospheric O3 RF to its precursors, initial estimates ( Fig. 10; Table 5) indicate that it is 620 predominantly NOx driven (48 %), followed by CH4 (34 %), with the smallest contribution from VOCs and CO (18 %). This is qualitatively consistent with Stevenson et al. (2013). However, as outlined in detail in that paper, these estimates do not take account of potential CH4 (and O3) changes that would occur were these experiments driven by CH4 emissions rather than concentrations (e.g. Shindell et al., 2005). Taking the same approach as Stevenson et al. (2013), we calculate an equilibrium CH4 concentration for the perturbation experiments, using the total CH4 lifetime relative to piClim-control based on the 625 following (Fiore et al., 2009): where lifetimes in the piClim-control and piClim-X perturbation experiments, respectively. The CH4-OH feedback factor (Prather, 1996) is denoted by f and is defined as: .
(11) 640 Using the whole-atmosphere CH4 burden and its removal by OH, the whole-atmosphere CH4 lifetime is 8.1 and 9.8 years in piClim-control and piClim-CH4, respectively, when adjusted for stratospheric removal (160 yr lifetime) and soil uptake (120 yr lifetime). From Eqns. (10) and (11), the UKESM1 feedback factor f is 1.28, consistent with the range of other estimates (Prather et al., 2001;Shindell et al., 2005;Fiore et al., 2009;Stevenson et al., 2013;Voulgarakis et al., 2013;Turnock et al., 2018) and within 5 % of the observationally constrained best estimate of 1.34 (Holmes et al., 2013). Then, using the difference 645 between the global mean prescribed and equilibrium surface CH4 concentrations, we calculate the additional tropospheric O3 forcing in the piClim-CH4, piClim-NOx and piClim-VOC experiments by applying the tropospheric O3 radiative kernel (Rap et al., 2015) to the O3 response derived from scaling the (piClim-CH4 minus piClim-control) O3 response based on the relationship in Turnock et al. (2018); this additional contribution to the tropospheric O3 RF is shown in Table 5.

650
Changing from a concentration-based perspective to an emissions-based view increases the tropospheric O3 RF between PI and PD from 407 to 461 mW m -2 (an increase of 13 %), which agrees better with Stevenson et al. (2013) once the offset by O3 depletion is accounted for (-101 mW m -2 ). It also changes the relative contributions of the different tropospheric O3 precursors.
The contribution of CH4 now dominates (45 %), with NOx playing a smaller role (36 %), while the contribution from VOC/CO emission increases is relatively unchanged (19 %). These emissions-based contributions are well within the spread of estimates 655 from Stevenson et al. (2013) who quantified contributions from 6 of the ACCMIP models: CH4 (44 ± 12 %), NOx (31 ± 9 %), and VOC/CO (25 ± 3 %) although there are some differences with Shindell et al. (2005) and Shindell et al. (2009). In both Shindell et al. studies, the contribution from NOx is lower at 15 ± 8 and 11 %. The difference may be due to the strong sensitivity of the CH4 lifetime to NOx in the GISS model compared with other models  and/or could be due to differences in VOC chemistry; Archibald et al. (2010) showed that the response of OH to increasing NOx strongly depends 660 on the treatment of VOC chemistry. Nevertheless, this approach demonstrates the importance of an emissions-based view of climate forcing and is more directly relevant to policy makers than a concentration-based view Shindell et al., 2009)

ERF: Role of Other Oxidants and Aerosols
In addition to O3 being a GHG, it is also important for secondary aerosol formation along with other oxidants: hydroxyl (OH), hydrogen peroxide (H2O2) and nitrate (NO3). The OH radical in UKESM1  is involved in aerosol nucleation 670 of gas-phase sulphuric acid (H2SO4) via the reaction with sulphur dioxide (SO2), leading to new particle formation . O3 and H2O2 are important for SO2 oxidation in cloud and aerosol droplets, creating sulphate aerosol mass but not number. Likewise, oxidation of monoterpenes by O3, OH, and NO3 determines the rate of formation of secondary organic aerosol (SOA) in UKESM1 (Kelly et al., 2018;Mulcahy et al., 2019) although it doesn't lead to new particle formation. Thus, sulphate aerosol alone, through changes in O3, OH, and H2O2, has potentially important impacts on cloud and aerosol radiative 675 properties. Indeed, Karset et al. (2018) found that PD to PI oxidant changes in the CAM5.3-Oslo model alter the relative importance of different chemical reactions, leading to changes in aerosol size distribution, cloud condensation nuclei (CCN), and the PD aerosol ERF. Figure 11 (bottom row) shows a global distribution of OH in piClim-control and changes in the perturbation experiments 680 relative to piClim-control. The piClim-control experiment shows an OH maximum in the equatorial humid regions, where photolytic production of OH from excited oxygen atoms (O 1 D) and water vapour (H2O) is at a maximum. In piClim-O3, OH increases throughout the northern hemisphere due to increases in O3, which is the precursor of O 1 D. These increases in OH are driven largely by increases in NOx. However, the piClim-VOC experiment shows the opposite behaviour -while VOC and carbon monoxide (CO) emissions increases serve to increase O3, they also remove OH via direct reaction with OH. This latter 685 effect outweighs the small increase in O3 in piClim-VOC and there are large decreases in OH throughout the troposphere. When OH is lower, we anticipate a decrease in the number of CCN, a decrease in cloud droplet number concentration (CDNC), leading to larger cloud droplets (Twomey, 1977), and an increase in the effective radius of cloud droplets (Reff). The middle rows of Fig. 11 show these effects at work. In piClim-control, the distribution of CDNC shows large values in equatorial 690 regions, in regions of continental outflow and regions of deep convection. Large increases in CDNC are seen in piClim-NOx, with large decreases in piClim-VOC. These results reflect the changes in OH in these experiments -increases in OH lead to increases in CDNC, and vice versa, but it should be noted that the effect of OH on CDNC is seen over a larger region downwind, particularly in East Asia and over the North Atlantic. The impact of NOx emissions appears to dominate, given the similarity between piClim-O3 and piClim-NOx, but the eastern Pacific is a region where all experiments show a decrease 695 in CDNC. While tropospheric column O3 increases in both piClim-NOx and piClim-VOC, leading to positive contributions to the PD tropospheric O3 RF (Sect. 3.4.2), the global mean PD ERF is very different between the two simulations (+0.21 W m -2 for 705 piClim-VOC and -0.08 W m -2 for piClim-NOx relative to piClim-control). Further, the spatial pattern of the ERF does not match those regions of largest O3 changes. Instead, the dominant driver of the ERF differences is changes to OH and the subsequent impacts on aerosol particle formation and clouds. In particular, the positive tropospheric O3 forcing in piClim-NOx (+0.2 W m -2 ) appears to be outweighed by a negative aerosol forcing, largely from aerosol-cloud interactions (aci) driven by https://doi.org/10.5194/acp-2019-1152 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License. changes in oxidants and aerosol nucleation; the contribution to the global mean ERF from aerosol-radiation interactions (ari) 710 is only -0.02 W m -2 , despite the strong regional changes in aerosol optical depth (AOD; Fig. 11; top row). Similarly, the forcing in piClim-VOC due to ari is negligible (0.00 W m -2 ). However, the forcing in piClim-VOC due to aci enhances the positive tropospheric O3 forcing (+0.07 W m -2 ), leading to a global mean ERF of +0.21 ± 0.04 W m -2 .
A negative PD forcing attributable to NOx emissions has been found in other studies (Shindell et al., 2009;Collins et al., 2010) 715 from a balance between the direct O3 response (positive RF), the NOx-driven CH4 response (negative RF), and the subsequent O3 response to CH4 changes (negative RF). Inclusion of ari and aci from sulphate and nitrate aerosol further increases the magnitude of the net negative forcing or cooling (Shindell et al., 2009;Collins et al., 2010). However, a study by Fry et al. (2012) found that the chemistry response is sensitive to the location of the emissions, with so large an uncertainty that it is difficult to determine whether NOx emissions cause a warming or cooling. Indeed, other indirect effects such as NOx 720 deposition to the terrestrial biosphere leading to fertilisation  and/or NOx-driven O3 damage (Sitch et al., 2007;Collins et al., 2010) increase the uncertainty further through changes in carbon dioxide (CO2). Thornhill et al. (2019) also shows that the ERF from PI-to-PD changes in NOx emissions among the AerChemMIP models differ in both sign and magnitude. Here, the longer time-scale CH4 response to NOx emissions  is constrained, nitrate aerosol is neglected Mulcahy et al., 2019), and CO2 is concentration-driven . Nevertheless, in 725 UKESM1, the negative forcing due to aci from sulphate aerosol outweighs the positive NOx-driven O3 RF, leading to a net negative ERF.
Previous work has also found a positive PD SARF or RF from VOC and CO emissions, due to the combined indirect forcings by O3 and CH4 Forster et al., 2007); they estimate a global mean RF of +0.21 ± 0.10 W m −2 at the PD 730 (Year 1998) relative to PI (Year 1750). As was the case with NOx, the magnitude of this RF increases (to +0.25 ± 0.04 W m -2 for the year 2000) when additional indirect forcings from sulphate, nitrate, and CO2 are included (Shindell et al., 2009). More recently, Stevenson et al. (2013) found the RF from VOC/CO emissions to be marginally higher, at +0.29 W m -2 , excluding aerosols, with contributions of 0.09, 0.08, and 0.12 W m -2 from O3, CH4, and CO2, respectively. The RF contribution from O3 alone quantified here (88 mW m -2 ; Table 5) is consistent with the Stevenson et al. (2013) estimate. Despite excluding the 735 longer-term CH4 response and CO2 , the total global mean ERF from VOC/CO emissions of +0.21 ± 0.04 W m -2 is consistent with previous estimates of RF due to the additional positive contribution from aci driven by OH changes. However, Fry et al. (2014) found that the RF from VOCs is sensitive to the location of emissions and could influence the strength of the contribution from aerosols. Interestingly, other AerChemMIP models show a negative ERF from VOC and CO emissions (Thornhill et al., 2019); these differences in sign of the ERF warrant further investigation. 740 Despite the very different chemical response between NOx and VOC emissions, both in terms of the magnitude of the O3 changes and the different impacts on OH, aerosols and clouds, a comparison of the global ERF changes (Table 3)  that the PD ERF from piClim-O3 is a linear combination of that from piClim-NOx and piClim-VOC for the NET, clear-sky (CS) and the cloud radiative effect (CRE) components. Nevertheless, the results clearly suggest that Earth System (ES) 745 interactions, particularly chemistry-aerosol coupling, can strongly affect estimates of climate forcing. Here, these interactions alter the PD ERF from tropospheric O3 precursor emissions, while other studies (e.g. Shindell et al., 2009;Karset et al., 2018) show that they also affect estimates of PD anthropogenic aerosol forcing.

Near-Term Climate Forcers (NTCFs) 750
The PD anthropogenic ERF due to CH4, aerosols and O3 abundances (also known as near-term climate forcers, NTCFs) was identified as the main source of uncertainty in the total anthropogenic ERF since pre-industrial times (Myhre et al., 2013a). This is due to the interaction between individual forcings, as well as the non-linear response of climate feedbacks due to aerosol-cloud interactions (Feichter et al., 2004;Deng et al., 2016;Collins et al., 2017). In this section, three experiments related to NTCFs are discussed: the combined simulation (piClim-NTCF) is identical to piClim-control except that the aerosol 755 and tropospheric O3 precursor emissions (excluding CH4) are set to PD (Year 2014) levels. The single-forcing runs change aerosol and aerosol precursor emissions only (piClim-aer), and tropospheric O3 precursor emissions only (piClim-O3) from PI to PD levels. More details are described in Table 1.
The global mean ERF of NTCFs (excluding CH4) at the PD relative to the PI is -1.12 ± 0.03 W m -2 (Table 3), where the 0.03 760 W m -2 is the standard error . The negative ERF results from the combination of a weak negative forcing in CS conditions (-0.04 ± 0.03 W m -2 ) and a strong forcing due to the CRE (-1.08 ± 0.02 W m -2 ). The weak negative forcing in the CS (NETcs) is due to the negative CS SW forcing (-0.29 ± 0.02 W m -2 ) being largely offset by the CS LW forcing (+0.24 ± 0.02 W m -2 ). The negative CS SW forcing is closely correlated with changes in AOD (Figs. 12a,12b). The CS LW component of the ERF (LWcs) is positive (+0.24 ± 0.02 W m -2 ) due to both aerosols and O3. The spatial variations in the global distribution 765 of the LWcs component are closely related to land surface temperature (Ts) changes (Figs. 12c, 12d); the correlation coefficient is -0.81 with a statistical significance well over 99 %. Considering this good correlation, the increased LWcs may be in response to the Ts change due to NTCFs. The negative SWcre (-1.02 ± 0.03 W m -2 ) is largely correlated with changes in cloudiness, with the spatial pattern correlation of -0.59 between the SWcre and total cloud fraction (Figs. 12e, 12f). The LW CRE (-0.06 ± 0.02 W m -2 ) is the smallest in magnitude, with good correlation with the changes in high level cloud fraction (Figs. 12g,770 12h). This is because the impact of cloud in the SW primarily depends on the cloud optical thickness, but the impact in the LW depends on the optically thin cloud for which both the optical depth and the cloud top temperature matter.
https://doi.org/10.5194/acp-2019-1152 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License. This study also attempts to estimate the effects of the non-linear interactions between the aerosol and O3 forcings. When combined (aerosol and O3), their interaction may induce an effect that differs from the sum of the individual single forcings.
The ERFs do not add linearly, particularly in the net CS components (Table 3). Firstly, we calculate the aerosol IRFs. In piClim-NTCF, it is -0.25 W m -2 , which is more negative than the sum of the aerosol IRF in piClim-aer (-0.16 W m -2 ) and 785 piClim-O3 (-0.02 W m -2 ). Secondly, the non-linearity may also be due to O3. Photochemistry in the atmosphere is a wellknown source of tropospheric O3 and is determined by ambient level of tropospheric O3 precursors (i.e. NOx and VOC/CO), and photolysis rates, which are largely influenced by meteorological factors such as solar irradiance and temperature (Xing et al., 2017). Despite no direct coupling between aerosols and photolysis in UKESM1 (Archibald et al., 2019), aerosol-mediated cloud feedbacks result in SW reduction and surface cooling (Figs. 12a,12d,12e), which impact thermal and photochemical 790 reactions leading to reduced photolysis rates in the lower troposphere, while enhancing the photolysis rate of O3 in the upper troposphere ( Figure 13). This is consistent with the vertical distribution of O3 (not shown) and changes in surface O3 (Fig. 14).

Land Use
Land-use change causes radiative forcing primarily by changing surface albedo; croplands and pastures have higher albedos than forests and are less able to mask the high albedo of snow cover. In UKESM1, historical land-use change causes surface 805 cooling and this reduces outgoing LW radiation which acts to damp the negative SW ERF, reducing the magnitude of the ERF from -0.39 to -0.32 Wm -2 . The land-use ERF of UKESM1 is more negative than other estimates and outside the 'very likely' range (−0.25 to -0.05 Wm -2 ) of AR5 Myhre et al., 2013a). However, its magnitude is reduced relative to UKESM1's predecessor, HadGEM2-ES (Collins et al., 2011), which produced a land-use ERF of -0.40 W m -2 . Robertson (2019) showed that, even in the absence of snow cover, HadGEM2-ES had too strong an albedo response 810 to land-use change and it is likely that this bias still exists in UKESM1.

815
While the global mean land-use ERF is small, regionally it can be the dominant source of radiative forcing. Figure 15  contribution. The ERF is mostly constrained to regions of land-use change, with deforestation in North America and western Eurasia causing a negative ERF and increased tree cover in central Europe causing a positive ERF. Alongside the distribution of land-use change itself, the magnitude of the ERF is larger in the mid-latitudes than the tropics, because the masking of snow 820 cover by trees greatly increases the albedo response to changes in tree cover. The seasonality of leafiness, snow cover and insolation cause the land-use ERF to be largest in northern hemisphere summer (June to August mean of -0.53 W m -2 ) and smallest in autumn (September to November mean of -0.16 W m -2 ).
The land-use ERF is calculated by modifying three land-surface fields: land-cover, leaf area index (LAI) and canopy height. 825 The values of the modified fields are taken from a coupled UKESM1 simulation that is forced only by historical land-use change. In the coupled configuration of UKESM1 , the three land surface fields are prognostic fields calculated using a dynamic global vegetation model (DGVM), while in the atmosphere-only ERF configuration they are prescribed fields. This choice was made because the land-surface fields can be very slow to respond to changes in land-use forcing and radiative forcing. In the coupled simulation of historical climate change including all forcings (hist), the land 830 surface fields respond to changes in CO2, climate and land-use, but if we use the land surface fields from this simulation, we find no substantial change in the land-use ERF (-0.35 W m -2 ).
In addition to the radiative forcing caused by albedo changes, land-use change can alter climate via a number of other mechanisms. Land-use change alters the land carbon sink. In particular, deforestation emits CO2 to the atmosphere, and so 835 some fraction of the CO2 ERF is attributable to land-use change. Land-use change also affects surface climate via changes in roughness length and transpiration. These non-radiative mechanisms usually drive larger temperature changes than the albedo response and both the non-radiative mechanisms and the CO2 forcings tend to oppose the albedo response.

Total Anthropogenic ERF
As summarized above, historical climate change has been driven by a wide range of anthropogenic forcings that act together, 840 alongside natural forcings, to perturb the Earth's radiation balance. The total anthropogenic ERF is, therefore, a key metric in understanding observed and modelled changes in the climate system since the PI era. These various anthropogenic forcing mechanisms are not necessarily independent of each other and it is therefore worthwhile calculating the total anthropogenic ERF from a separate timeslice simulation including all forcings together (piClim-anthro; Pincus et al., 2016). We completed a dedicated timeslice simulation with all anthropogenic forcings set to 2014 levels, including all GHGs (CO2, N2O, CH4, HCs 845 both ODSs and non-ODSs), tropospheric O3 precursors (VOC/CO and NOx), land-use changes, and anthropogenic aerosol or aerosol precursor emissions (SO2, OC, BC). This experiment was not proposed in AerChemMIP (Collins et al., 2017) but is included here as part of RFMIP . The main difference is that atmospheric chemistry in UKESM1 is fully interactive whereas other models participating in RFMIP (e.g.  use the CMIP6 ozone dataset to represent changes in tropospheric and stratospheric O3. 850 https://doi.org/10.5194/acp-2019-1152 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License.
The UKESM1 piClim-anthro experiment leads to an ERF of 1.61 W m -2 (Table 3), which is dominated by a positive CS LW forcing due to GHGs and partially offset by a negative SW component due to aerosols and negative adjustments to LW and SW cloud-radiative effects (Table 3). This ERF is a little lower than the equivalent estimate from HadGEM3-GA7.1, which was 1.81 W m -2 . The UKESM1 estimate is also somewhat lower than the median estimate from CMIP5 855 models assessed in AR5, which equates to approximately 1.9 W m -2 after adjustment to the reference period of 1861-1880 to 2010-11 (Andrews and Forster, 2019; hereafter AF19). AR5 also provided an overall central estimate of 2.2 W m -2 and a 5-95 % confidence range of 1.0 -3.2 W m -2 (after adjustment to the same reference period as above in AF19) taking into account multiple streams of evidence. AF19 re-evaluated this as 2.3 W m -2 with a narrower range of 1.7 -3.0 [5 -95 % confidence] using a combination of atmospheric model outputs and observational constraints. The lower bound of this range was reduced 860 to 1.5 W m -2 if larger uncertainties were assumed for the climate feedback parameter or for the global-mean surface temperatures anomalies used to constrain the forcing. The UKESM1 estimate of 1.61 W m -2 is therefore within the original uncertainty range given by AR5 but on the edge of the range proposed in AF19.
There are several factors that contribute to the relatively low estimate of anthropogenic ERF in UKESM1. Firstly, the 865 anthropogenic aerosol forcing in UKESM1 (and GC3.1) is -1.13 W m -2 (well within the uncertainty range) and offsets a major portion of the positive GHG forcing (2.89 W m -2 ). Secondly, the ERF from piClim-HC is negative (-0.33 W m -2 ) due to a strong O3 response and connected aerosol-mediated cloud feedbacks . Thirdly, adjustments in vegetation lead to an appreciable negative ERF from land-use changes (-0.32 W m -2 ). The stronger negative forcing from piClim-HC and piClim-LU largely explain why the UKESM1 estimate is ~0.2 W m -2 lower than HadGEM3-GA7.1. As shown 870 in Fig. 3e, the negative contributions more than offset the positive GHG forcing in certain regions. For instance, the anthropogenic ERF is negative over large parts of North America and Asia (from a combination of land use change and aerosol forcing; see Figs. 3d and 3b), over the North Pacific (due to aerosol-cloud forcing) and at southern high latitudes (from O3 depletion due to HCs, see Fig. 3a). In contrast, the anthropogenic forcing is strongly positive over the tropics and southern hemisphere sub-tropics (Fig. 3e) where the direct radiative effect of GHGs dominates. 875 The couplings between chemistry, aerosol and land surface processes included in UKESM1 increase the possibilities for nonlinear feedbacks and interactions among the various anthropogenic forcing agents. However, these apparently have little net overall effect on the total forcing (1.61 W m -2 ), which is almost identical to the sum of the forcings (1.59 W m -2 ) from the four separate groups that it includes (GHG, aerosol, tropospheric O3 precursors, and land use). The two estimates are not statistically 880 different given the standard error on these is around 0.03 -0.04 W m -2 (Table 3). This does not imply that the forcings act independently as it is possible that competing non-linear interactions cancel in this particular case. As shown in Sect. 4.5.1, the aerosol and tropospheric O3 forcings did not add linearly, and neither did the sum of individual GHGs forcings (Sect.

Conclusions 885
Quantifying forcings from anthropogenic perturbations to the Earth System (ES) is important for understanding changes in climate since the pre-industrial (PI) period. In this paper, we have quantified and analysed a wide range of PD anthropogenic forcings with the UK's Earth System Model (ESM), UKESM1 , using the concept of effective radiative forcing (ERF). ERFs have been shown to be a more useful metric for evaluating and comparing the relative roles of diverse forcing agents due to the relationship to global-mean temperature and other impacts that scale with it. In particular, by 890 quantifying ERFs within a full ESM, this study addresses gaps in previous assessments in which rapid adjustments were neglected and enables the role of indirect forcings and various climate-chemistry-aerosol-cloud feedbacks to be quantified.
We find that carbon dioxide (CO2) exerts a global mean ERF of 1.82 W m -2 at the present-day (PD), consistent with previous estimates, making it the single largest contributor to PD climate forcing. However, UKESM1 appears to have a more 895 pronounced surface warming adjustment associated with the physiological forcing by CO2 than its successor, HadGEM2-ES.
The nitrous oxide (N2O) PD ERF quantified here (0.13 W m -2 ) is lower than previous estimates due to indirect effects from ozone (O3) depletion, and fast cloud adjustments; a shift of circulation results in a redistribution of clouds with clouds structure moving poleward. A warming associated with an O3 increase in the tropical upper troposphere/lower stratosphere (UTLS) region may result in a regional suppression of convection and a reduction of associated LW radiation coming from cloud top. 900 The PI-to-PD change in methane (CH4) concentration leads to a global mean ERF of +0.93 W m -2 , with the majority of the

905
The global annual mean ERF from the PI to PD change in O3 depleting substances (ODSs) is -0.33 W m -2 . A quantitative and process-based understanding of what's driving the negative ERF, despite our understanding to date of the relative roles of ODSs and O3 depletion in climate forcing, is presented in . Considering all greenhouse gases (GHGs) together, we quantify a global mean ERF of 2.89 W m -2 , less than the 3.09 W m -2 estimate from the physical model HadGEM3-GC31, due to indirect effects. There is also some evidence of non-linearity between the combined GHG ERF and the sum of 910 the individual GHG ERFs, potentially due to interactive O3 and aerosol-mediated cloud feedbacks, which may warrant further investigation.
interactions included in the new aerosol scheme mean that neither aerosol IRF nor aci are linear (sulphate, organic carbon (OC), BC interact with one another) making the total aerosol ERF less than the sum of the individual ERFs.
Examining tropospheric O3 radiative forcing alone, results from UKESM1 suggest that the contribution from CH4 dominates 920 (45 %), with nitrogen oxides (NOx) and volatile organic compound (VOC)/carbon monoxide (CO) contributing 36 and 19 %, respectively. These emissions-based contributions are well within the spread of estimates from ACCMIP (Stevenson et al., 2013): CH4 (44 ± 12 %), NOx (31 ± 9 %), and VOC/CO (25 ± 3 %) although there is disagreement with other studies. Changes in oxidants, driven by the PI-to-PD changes in tropospheric O3 precursor emissions, lead to an indirect aerosol ERF from aci, which either supplements or offsets the positive ERF from tropospheric O3, leading to global mean ERFs of 0.21 and -0.08 W 925 m -2 for VOC/CO and NOx emissions changes, respectively. However, there appears to be disagreement across the AerChemMIP models on the sign and/or magnitude of the tropospheric O3 precursor ERFs and further analysis to understand what's driving these differences is required.
The aerosol and tropospheric O3 precursors (near-term climate forcers, NTCFs) together exert a global mean ERF of -1.12 W 930 m -2 , which is mainly due to changes in the cloud radiative effect (CRE). There is also evidence of non-linearity in the global mean ERF between the combined piClim-NTCF experiment and the sum of the individual piClim-aer and piclim-O3 experiments; this is mainly evident in the CS and is driven by changes in aerosol optical depth (AOD) and O3. Land use (LU) change since the PI has also exerted a negative forcing at the PD, estimated to be -0.32 W m -2 on a global mean basis. However, this estimate is outside the likely range from previous estimates, and is most likely due to too strong an albedo response. 935 Historical climate change has been driven by a wide range of anthropogenic forcings that act together, alongside natural forcings, to perturb the Earth's radiation balance. As a result, the total anthropogenic ERF is a key metric in understanding observed and modelled changes in the climate system since the PI era. The estimate of the total anthropogenic ERF from UKESM1 is 1.61 W m -2 , which is relatively low compared to previous assessments; this is mainly due to an intermediate 940 negative aerosol ERF, a modest negative LU forcing and strong stratospheric O3 depletion. Although it may be biased low, that combined with high climate sensitivity means that UKESM1 reproduces historical global mean warming over the 1850 -2014 period .
In addition to quantifying anthropogenic PD ERFs with a fully coupled ESM, this study and other studies (e.g. Morgenstern 945 et al., 2019;O'Connor et al., 2019) show the importance of indirect forcings and climate-aerosol-chemical feedbacks, and quantified their role within the ERF framework. There are substantial feedbacks between GHGs, stratospheric and tropospheric O3, and aerosols, some of which act non-linearly. These effects demonstrate the importance of including Earth System (ES) interactions when quantifying PD climate forcing. In particular, we consider that rapid adjustments included in the definition of ERF should include chemical as well as physical adjustments, consistent with Ramaswamy et al. (2019). They concluded in their recent assessment that although the radiative forcing concept is simple, it needs to increasingly account for the complex relevant processes in the Earth System. The authors thank T. Andrews for constructive comments on the manuscript. Last but not least, the authors wish to acknowledge the huge effort from the UKESM1 core group in building and evaluating UKESM1 and making it available for use in RFMIP and AerChemMIP. 975

Data Availability
Data from all the simulations which underpin this paper are in the process of being post-processed in readiness for uploading to the Earth System Grid Federation (ESGF). The authors endeavour to use the precise data citation for each piece of data once all of it is citable.

Completing Interests
The authors declare that they have no conflict of interest.