The compact Earth system model OSCAR v2.2: description and first results

This paper provides a comprehensive description of OSCAR v2.2, a simple Earth system model. The general philosophy of development is first explained, followed by a complete description of the model's drivers and various modules. All components of the Earth system necessary to simulate future climate change are represented in the model: the oceanic and terrestrial carbon cycles – including a book-keeping module to endogenously estimate land-use change emissions – so as to simulate the change in atmospheric carbon dioxide; the tropospheric chemistry and the natural wetlands , to simulate that of methane; the stratospheric chemistry , for nitrous oxide; 37 halogenated compounds; changing tropospheric and stratospheric ozone; the direct and indirect effects of aerosols; changes in surface albedo caused by black carbon deposition on snow and land-cover change; and the global and regional response of climate – in terms of temperature and precipitation – to all these climate forcers. Following the probabilistic framework of the model, an ensemble of simulations is made over the historical period (1750–2010). We show that the model performs well in reproducing observed past changes in the Earth system such as increased atmospheric concentration of greenhouse gases or increased global mean surface temperature.


Introduction
Simple biogeochemistry-climate models, also qualified as compact or reduced-form models, are widely used in the climate change research community. These models share several features. First, they are not spatially resolved and as such they can be referred to as box models, although the number of boxes -and therefore of state variables -may vary greatly: from a couple to several hundred. This limited number of state variables make them relatively intelligible, when compared to complex three-dimensional models. Second, the time step of analysis and of numerical solving is about 1 year, which implies they usually cannot include representations of seasonal processes. One consequence of these two features is their very high computing efficiency: one simulation typically takes about 1 min on a laptop. Compact models can therefore be used in a variety of setups, such as the following for instance: to translate a large number of pathways of greenhouse gases emissions into projected climate change (e.g., Clarke et al., 2014), to complement a study by a process-based model (e.g., Schneider von Deimling et al., 2012) or an economic model (e.g., Rogelj et al., 2013), to extend a given experiment or assess its uncertainty with a Monte Carlo analysis (e.g., Gasser et al., 2015), to attribute changes in a variable of the climate system to physical processes (e.g., Raupach et al., 2014) or to emitting countries (e.g., Höhne et al., 2011), or to easily discuss theoretical frameworks (e.g., Raupach, 2013) or policy-relevant indica-Published by Copernicus Publications on behalf of the European Geosciences Union.
tors (e.g., Huntingford et al., 2012). Because they are simple models of complex phenomena, these models can hardly be process-based. Quite the opposite: they usually consist of ad hoc parametric laws which require calibration or optimization using either observations (e.g., Ricciuto et al., 2008) or more complex models (e.g., Meinshausen et al., 2011).
Here, we present an important update of a simple Earth system model that has already been used for some time. The model is named OSCAR, and this paper provides a comprehensive description of version 2.2. OSCAR can be described as a non-linear box model of which the number of boxes, however, is fairly large. It is not spatially resolved (i.e., it is not gridded) but key processes such as land-use change or aerosol physico-chemistry are regionalized to account for the disparity in such processes that is observed in the real world. OSCAR does not endogenously simulate intra-or inter-annual variability. Consequently, although the time step of its inputs and outputs is 1 year, the main purpose of the model is to simulate trends in the Earth system change, and not year-to-year variations. OSCAR is also a parametric model in which a relatively large number of parameters are almost all calibrated on complex models. We call this approach meta-modeling: each module of OSCAR is designed to emulate the behavior of other more specialized models (e.g., global climate models, dynamical vegetation models, or chemistry-transport models). For most modules, we have access to several sets of parameters (one per complex model used to calibrate) and, rather than taking the average or arbitrarily choosing one, we adopt a probabilistic approach in which a given simulation with OSCAR is repeated many times with different sets of parameters picked at random. This paper is firstly a thorough description of OSCAR. Readers who are more interested in the model's applications should know that only a simulation covering the historical period is made here, for diagnostic and (partial) validation purposes. More concrete applications of the model or of its older versions can be found elsewhere in the literature (e.g., Gasser, 2014;Gasser et al., 2015Gasser et al., , 2016Li et al., 2016) and will be presented in future studies.
Section 2 provides the details of the mathematical formulation of the model, and it describes how the parameters are calibrated or derived. The first subsection provides general information about the model. The second subsection describes the forcings of the model: anthropogenic emissions of greenhouse gases, ozone precursors, aerosols and aerosols precursors, land-use and land-cover change, and some additional anthropogenic and natural radiative forcings. The next subsections describe the model's various modules dedicated to the prediction of carbon dioxide, methane, nitrous oxide, halogenated compounds, ozone, aerosols, surface albedo, and the climate response. The last subsection describes the numerical solving method. Section 3 then provides the first results from OSCAR v2.2 in the case of a simulation over the historical period from 1750 to 2010, as well as a brief discussion of these results. Section 4 provides some preliminary conclusions regarding the model, its performance and interest, and future potential development.

General points
Since version 2.0 (see Appendix A for an overview of the model's history), the development of OSCAR has been driven by three principles. First, we aim to embed in OSCAR as many components and processes of the actual Earth system as possible, with the overall idea of favoring the amount of processes, interactions, and feedbacks implemented over the modeling complexity of each of these elements. A diagram of the model's simplified causal structure is shown in Fig. 1. Second, OSCAR is built as a meta-model capable of emulating the sensitivities of models of higher resolution or superior complexity. To do so, in order to calibrate those sensitivities, we use outputs from complex models that participated in intercomparison projects whenever possible. Consequently, OSCAR is designed to be used in a probabilistic framework. The last table of the paper summarizes the various options available for this probabilistic setup. Third, the model is developed as a dynamic model of the Earth system with an assumed equilibrium corresponding to the preindustrial era. The reason for this is the original purpose of the model, to perform studies attributing the anthropogenic component of climate change (e.g., Ciais et al., 2013a;Gasser, 2014;Li et al., 2016).
This last point is also the reason why all the following equations are expressed as a difference to our so-called preindustrial equilibrium. In the model, all variables Z are therefore formulated as a constant term Z 0 , being the preindustrial value of the variable, to which a time-varying perturbation term Z is added, so that we always have the following: Z = Z 0 + Z. Only the forcings, described in the following subsection, are expressed without this -term, since per construction their constant term is zero. This -term is the actual state variable of the model, and it is consequently the one used to discuss the performance of OSCAR in Sect. 3. Also, all these state variables are nil at the beginning of the simulation (i.e., at t = 0). Because of this modeling approach, one may better describe OSCAR as being a "climate change" or "Earth system change" model.
Let us now introduce some of the main notations that are used throughout the description sections. The variables of the model are written with Roman letters, whereas its parameters are with Greek letters. Among the variables, some letters are consistently dedicated to a specific kind of variable: F for fluxes of matter; E for emissions, i.e., positive fluxes towards the atmosphere; C for carbon pools; A for surface areas; T for temperature; P for precipitation. Similarly, the parameters include the following: α for proportionality factors, and  Figure 1. Simplified causal chain of OSCAR v2.2. Each node of the graph corresponds to a module described in the section whose number is shown below the node's name. Colored lines show the forcings of the model, black lines show the natural cause-effect chain, and dashed lines show the climate feedbacks. "Halo" groups all the halogenated compounds, "Ocean" is the ocean carbon cycle, "Land" is the land carbon cycle, "Albedo" groups the surface albedo effects, "hν" is the stratospheric chemistry, "OH" is the tropospheric chemistry, "Cloud" is the indirect aerosol effect, and "Climate" groups the surface temperatures and precipitation. therefore also conversion factors; γ for relative sensitivities to a climate variable; for absolute sensitivities to a climate variable; χ for chemical sensitivities; τ for time constants; ν for rate constants, i.e., the inverse of time constants; κ for dimensionless constants; π for fractions, i.e., dimensionless constants within [0, 1]; and ω for weights, i.e., dimensionless constants whose reference value is 1. A few variables or parameters, however, do not follow these general rules. Additionally, we use the notation F for mathematical functions whose arguments are given in square brackets. Superscripts are used either to refer to a given atmospheric species or to denote the subdivision of a given variable or parameter along a given axis (e.g., a regional axis). The time variable t is kept implicit for legibility, unless it is required.

Anthropogenic emissions
Anthropogenic emissions of various active gases are the main drivers of climate change, and thus of OSCAR. In this version of the model, the anthropogenically emitted greenhouse gases directly prescribed as emissions over the 1750-2010 period are carbon dioxide from fossil-fuel burning and cement production (E FF ), methane (E CH 4 ), nitrous oxide (E N 2 O ), and many halogenated compounds (E X , see list in Sect. 2.7). The ozone precursors, similarly prescribed, are nitrogen oxides (E NO x ), carbon monoxide (E CO ), and non-methane volatile organic compounds (E VOC ), and the aerosols and aerosol precursors are sulfur dioxide (E SO 2 ), ammonia (E NH 3 ), organic carbon (E OC ), and black carbon (E BC ).
Because most of these anthropogenic drivers are actually poorly known, especially when going backward in time, we use various established inventories in order to introduce variation in our model's results, and thus explore the uncertainty in climate change reconstruction and projection. The historical emissions of fossil-based CO 2 can be based on the CDIAC dataset (Boden et al., 2013) or on EDGAR v4.2 (EC-JRC/PBL, 2011). Those of CH 4 can be based on EDGAR, ACCMIP ), or EPA (2012. Those of N 2 O can be based on EDGAR or EPA. Those of all hydrofluorocarbons (HFCs), perfluorocarbons (PFCs), sulfur hexafluoride (SF 6 ) and nitrogen trifluoride (NF 3 ) are based on EDGAR, while those of ozone-depleting substances (ODSs) are taken from Meinshausen et al. (2011). The emissions of NO x , CO, VOCs, SO 2 and NH 3 can be based on EDGAR or ACCMIP -this being set independently for each species. Those of OC and BC are based on ACCMIP. Note that the emissions from biomass burning of natural vegetation are removed from these datasets, since those emissions are endogenous to the OSCAR model (see Sect. 2.4.1), but all other sectors provided by the inventories are included, notably (but not only) agricultural waste burning, shipping, and aviation.
Similarly to what is done by Li et al. (2016, Supplement Sect. B.1), the time series of anthropogenic emissions are constructed as follows. First, we choose one of the datasets above as the reference dataset whose absolute values are taken as they are over its period of definition. Second, if needed, we extend the reference dataset beyond its period of definition by following the relative year-to-year variation of other datasets. The extension forward is needed in two cases: with the EDGAR dataset as reference, in which case the extension covers 2008-2010 and is made with the EDGAR-FT v4.2 dataset (EC-JRC/PBL, 2013) for greenhouse gases and the EDGAR-HTAP v2 dataset (EC-JRC/PBL, 2014) for other species, as well as with the ACCMIP dataset, in which case the extension covers 2000-2010 and is made with EDGAR first and then EDGAR-FT or EDGAR-HTAP. The extension backward is also needed in almost all cases. It is done with CDIAC for fossil-fuel CO 2 based on EDGAR, with EDGAR-HYDE v1.4 (van Aardenne et al., 2001) for N 2 O based on EDGAR or EPA and with ACCMIP for any other species based on EDGAR or EPA. For N 2 O, however, given that the EDGAR-HYDE dataset has global values that differ significantly from more recent estimates, the dataset is rescaled to the global estimates by Davidson (2009) before being used for extension. For non-CO 2 species, the obtained time series are then linearly extrapolated from 1850, or 1860 in the case of N 2 O: to be zero in 1500 for biogenic emissions, and in 1750 for fossil-related emissions. The year 1500 is taken to be consistent with the dataset we use for landuse change (see Sect. 2.2.2). For the HFCs and PFCs whose emissions are not zero in 1970 (first year of the EDGAR inventory) the backward extrapolation is slightly different.
From 1970, their emissions are extrapolated backward following a quadratic function of the time, to reach zero in a year taken from Meinshausen et al. (2011): 1930for HFC-23, 1950for SF 6 , 1922 for CF 4 , and 1889 for C 2 F 6 .
Finally, because of our assumption of a preindustrial equilibrium occurring before 1750, we have to offset the full time series of anthropogenic emissions obtained thus far by the level of emissions of 1750, thus assuming that everything before that point in time is part of the natural cycle -or at least negligible as an anthropogenic perturbation of this natural cycle. The obtained time series of anthropogenic emissions for all species are shown in Fig. 2, except for the halogenated compounds which are shown in Fig. S1 in the Supplement.

Land-use and land-cover change
As OSCAR embeds a book-keeping module to endogenously calculate CO 2 emissions from land-use change, it needs specific types of drivers related to land-use and land-cover change (LULCC). These are threefold. The first driver is for land-cover change itself (δA b 1 →b 2 , an area per unit time): it is the human-induced transitions from one biome b 1 to another b 2 . The second driver is for wood harvest (δH b , a mass of carbon per unit time): it is the extraction of woody biomass from a given biome b but without changing the land cover, and it can be seen as a coarse forestry driver. The third driver is for shifting cultivation (δS b 1 ↔b 2 , an area per unit time): it is the transition from one natural biome b 1 to another anthropogenic one b 2 which occurs simultaneously with the reciprocal transition. The latter driver is therefore a triangular matrix on the (b 1 , b 2 ) axes, and it is typical of (but not exclusive to) the agricultural practice happening in the tropics and known as "slash-and-burn".
In this version of OSCAR, only one LULCC dataset is available: the LUH v1.1 dataset (Hurtt et al., 2011) updated for the last TRENDY exercise . Given that this dataset provides information only for an aggregated "natural vegetation" biome, whereas OSCAR considers different natural biomes, we need to combine the dataset with the natural vegetation maps used for the terrestrial carbon cycle in Sect. 2.3.2, so as to disaggregate further the natural vegetation provided by the LUH project. Other than that, the dataset is used as is over the 1750-2010 period; it is shown in Fig. S2.

Radiative forcings
Finally, some remaining known climate forcings are prescribed to OSCAR directly as radiative forcings. This is the case of one anthropogenic forcing (aviation contrails and induced cirrus (RF con )) and of two natural forcings (volcanic aerosols (RF volc ) and solar irradiance (RF solar )). Those drivers are directly taken from IPCC (2013), except that we offset the volcanic forcing by its value averaged over 1750- 2011, thus assuming this value to be representative of volcanic preindustrial conditions.

Carbon dioxide
The carbon cycle in OSCAR is divided in two components: ocean and land. The ocean carbon cycle includes the socalled physical pump, i.e., the dissolution of anthropogenic CO 2 in the surface ocean and transport of this dissolved carbon to the deep ocean. No anthropogenic perturbation of the biological pump is implemented, however. The land carbon cycle includes the response to environmental perturbations such as change in atmospheric CO 2 or climate, but not change in nutrient availability. It also includes the response to human-made land-use and land-cover change, but no natural land-cover change -i.e., dynamical vegetation -is implemented. Notably absent from the model are fluxes linking the land to the ocean carbon pools (such as river discharges), the oxidation of non-CO 2 species into CO 2 , and permafrost carbon.

Ocean carbon cycle
The ocean carbon-cycle module is based on the mixed-layer impulse response developed by Joos et al. (1996), widely used among compact models (e.g., Meinshausen et al., 2011;Raupach et al., 2011), albeit with three important modifications. First, the convolution with the impulse response function is written as its equivalent box model, similarly to what Harman et al. (2011) did. Second, the ad hoc function used to emulate the carbonate chemistry is updated (Harman et al., 2011) and it now includes a dependency on sea surface temperature. Third, we extend the initial formulation to include a varying mixed-layer depth, assumed to vary with global sea surface temperature in order to represent the stratification of the upper ocean induced by global warming (e.g., Capotondi et al., 2012). With the last two modifications, key mechanisms of the global ocean uptake -such as carbonate saturation, warming-driven changes in solubility, and impact of ocean stratification (Ciais et al., 2013b) -are better accounted for in OSCAR.
Following Joos et al. (1996), we explicitly represent only one oceanic carbon pool which corresponds to the surface ocean (C surf ); other oceanic carbon pools are only implicitly considered. The two carbon fluxes -in and out -between this surface ocean and the atmosphere are then calculated separately. They both are proportional to the gas exchange rate (ν fg ) and to the atmospheric conversion factor for CO 2 (α CO 2 atm ). The latter is only used to express a partial pressure of CO 2 -in parts per million -into a quantity of carbon -in gigatons of carbon. The in-going flux (F in ) is a linear function of the atmospheric partial pressure in CO 2 (CO 2 ): and the out-going flux (F out ) is also a linear function of the sea surface partial pressure in CO 2 . This partial pressure is calculated thanks to an ad hoc function (F pCO 2 ) designed to emulate the non-linear oceanic carbonate chemistry. It depends on dissolved inorganic carbon concentration in the surface ocean (dic) and sea surface temperature (T S ): www.geosci-model-dev.net/10/271/2017/ Geosci. Model Dev., 10, 271-319, 2017 The dissolved inorganic carbon is deduced from the total carbon stored in the surface layer via a conversion factor (α sol ), the global area of the ocean (A ocean ), and the mixed-layer depth (h mld ).
where the mixed-layer depth is assumed to vary according to the following law, parameterized by the maximum relative intensity of the stratification (π mld ∈ [0, 1]) and its sensitivity to sea surface temperature change (γ mld < 0): We then represent the net effect of the oceanic circulation and mixing fluxes as a unique flux of carbon that goes from the surface ocean to an implicit deep ocean (F circ ). To do so, the surface ocean is subdivided into several boxes (superscript o ). Each box contains a fraction of the total surface carbon and is assigned an areal fraction (π o circ , o π o circ = 1) and a turnover time (τ o circ ), so that it works as a first-order model. So we have the following: This subdivision of the surface ocean is not geographical: it only corresponds to the different turnover times of mixing between the surface and deep oceans, accounted for by the response function of Joos et al. (1996), and the boxes are not distinguished otherwise as they, e.g., share the same carbonate chemistry. Finally, the global perturbation of the ocean carbon cycle is obtained by summing over the o-boxes, and by solving the carbon budget in each of the boxes, Note that this model of the ocean carbon cycle implicitly assumes no change in the biological pump -change that could be induced, e.g., by changes in temperature, ocean circulation, nutrient availability, or surface acidity (e.g., Ciais et al., 2013b). This is one of the several processes not implemented in this version of OSCAR. The atmospheric conversion factor is calculated following Prather et al. (2012, Table S2): a value of 0.1765 Tmol ppb −1 of dry air is assumed, and it is multiplied by the molecular mass of any species X to obtain α X atm . The conversion factor α sol is set to 1.722 × 10 17 µmol m 3 ppm −1 kg −1 (Joos et al., 1996). The function F pCO 2 can be either one of the two formulations (a Padé approximant or Power law fit) given by Harman et al. (2011). The parameters ν fg , A ocean , h mld,0 , π circ and τ circ , as well as the preindustrial sea surface temperature T S,0 , are taken from Joos et al. (1996) who provide four sets of values derived from four ocean models of various complexity.
We use the latter study in a way very close to what is done by Harman et al. (2011). We take the long-term response functions of Joos et al. (1996) which, analytically, are a weighted sum of exponential terms: π circ are taken as the weights and τ circ as the time constants. The number of oboxes is thus equal to the number of exponential terms. To ensure that the response is consistent in the short term, however, we add to these another box whose π circ is taken as the complementary fraction so that the sum of all fractions gives 1, and whose τ circ is arbitrarily set to 1/3 year in the case of the "HILDA" and "box-diffusion" models, and 1/2 year in the case of the "2D-Princeton" model. In the case of the "3D-Princeton" model, however, because the sum of the fractions provided by Joos et al. (1996) is greater than 1, we do not add that other box, we simply reduce the fraction of the fastest box so that the sum of the fractions is one.
The two parameters related to the mixed-layer depth, i.e., π mld and γ mld , can be calibrated on three CMIP5 Earth system models (see Appendix B for a list of the models used to calibrate OSCAR). To do so, we use the CMIP5 output variable named "omlmax" which corresponds to the maximum depth of the mixed layer over a given period of time. We then fit the parameters, on the basis of Eq. (4), using the relative variation of the "omlmax" variable over the historical period and the RCP8.5 up to 2300 with respect to the control simulation, and the sea surface temperature change. This fit is made with yearly data. Since for the "CESM1-BGC" model no value is available over 2100-2300, we use outputs from Randerson et al. (2015, data from Fig. 6b) instead. The CMIP5 fits are shown in Fig. S3.

Land carbon cycle: intensive perturbation
Before considering the extensive perturbation of the terrestrial carbon cycle, driven by land-use and land-cover changes, we first focus on its intensive perturbation, i.e., the perturbation that changes the areal properties of the terrestrial ecosystems. This intensive perturbation is driven by changes in environmental conditions such as atmospheric CO 2 , climate, or nutrient deposition -albeit not the latter in this version of OSCAR. Since only the areal properties of the terrestrial biosphere are affected, this section only describes the evolution of OSCAR's state variables per unit area. The intensive variables, i.e., the variables per unit area, are thereafter written in lowercase, in opposition to the extensive variables, written in uppercase, that we use in the next section.
The terrestrial carbon-cycle module is an upgrade of the previous versions (e.g., Gitz and Ciais, 2003;Gasser, 2014). The terrestrial biosphere is aggregated into several regions (superscript i ) and further divided into various biomes (superscript b ) in each region. Here, we note that the exact regional aggregation (both the number of regions and their def-Geosci. Model Dev., 10, 271-319, 2017 www.geosci-model-dev.net/10/271/2017/ inition), and to a lesser extent the way biomes are aggregated, can be chosen before every simulation, thus altering the numerical values of the related parameters. For this description paper, we set the regional aggregation to the nine broad world regions used by Houghton and Hackler (2001, see also our Fig. S4), and the biome list to bare soil, forests, mix of grasslands and shrublands, croplands, and pastures. Each doublet (i, b) represents the "average" biome b of the i-th region with assumed homogeneous biogeochemical characteristics. Each doublet (i, b) is then represented as a three-box model where each box exchanges carbon with the others and/or the atmosphere.
Contrarily to complex models that simulate separately the gross primary productivity, our terrestrial carbon cycle starts directly with the net primary productivity (NPP). The areal net productivity (npp) depends on a preindustrial intensity (η), and it responds to changes in atmospheric CO 2 via a fertilization function (F fert ), and changes in local surface temperature (T L ) and local yearly precipitation (P L ) with assumed linear sensitivities (γ npp,T and γ npp,P , respectively). This gives the following: The functional form of F fert can be either logarithmic or hyperbolic (Friedlingstein et al., 1995). If logarithmic, it is a function with one parameter that describes the intensity of the fertilization effect (β npp > 0).
and if hyperbolic, it is a function with two parameters: one also describing the intensity of the fertilization effect ( β npp > 1) and the other being the compensation point (CO 2cp ), i.e., the value of atmospheric CO 2 below which there is no NPP at all.
NPP fills a first carbon pool that corresponds to the vegetation's living biomass (c veg ). This biomass is partly oxidized by wildfires, creating a flux to the atmosphere (e fire ). This flux is proportional to the biomass available to be burnt and also depends on the preindustrial fire intensity (ι) and a function representing the variation of this fire intensity (F igni ): The change in fire intensity is a function of changes in atmospheric CO 2 -used as a proxy variable to encompass various effects such as change in leaf area index, which would help wildfires to spread, or change in evapotranspiration and thus in soil moisture that would reduce their intensity -in local surface temperature, and in local yearly precipitation. We arbitrarily choose a linear sensitivity for each of the three environmental factors (γ igni,C , γ igni,T > 0 and γ igni,P < 0, respectively). Here we note that our formulation implicitly assumes that there is no direct human intervention to, e.g., limit and control natural wildfires. So the function F igni is formulated as follows: The living biomass also partly dies at a fixed rate (µ), which generates a flux we call "mortality" (f mort ). The mortality rate does not depend on environmental conditions such as climate, but the lack of detailed outputs from complex models motivates this modeling choice. This is written as follows: The mortality flux goes into the litter carbon pool (c litt ), where heterotrophic respiration (rh litt ) occurs at a rate that depends on its own preindustrial value (ρ litt ) and a specific function of local climate conditions (F resp ): The functional form of F resp can be either exponential or Gaussian (Tuomi et al., 2008) regarding its sensitivity to temperature. It is always linear regarding sensitivity to precipitation. If exponential, it is therefore a function with two parameters (γ resp,T > 0 and precipitation γ resp,P ): and if Gaussian, it is a function with three parameters, two of which are the sensitivity to temperature split between a first-order term (γ resp,T 1 > 0) and a second-order term (γ resp,T 2 < 0), with the third being the sensitivity to precipitation ( γ resp,P ).
Part of the litter carbon is metabolized into soil organic carbon. This flux (f met ) is taken proportional to the heterotrophic respiration of the litter carbon pool (by a factor κ met ): Heterotrophic respiration (rh soil ) also occurs in the soil carbon pool (c soil ). It is a function of its preindustrial rate (ρ soil ) and of the same function F resp as for the litter: Finally, the terrestrial carbon cycle of a given doublet (i, b) is as follows: The equation system described above by Eqs. (19), (20), and (21) implies that our preindustrial equilibrium is as follows. fire,0 which, in terms of flux parameters and preindustrial carbon stocks, is equivalent to the following: Note that to obtain the global preindustrial terrestrial carboncycle equilibrium, one needs to multiply the above equilibrium by the preindustrial biome area extents (A 0 ), for instance NPP global 0 Note also that the extensive perturbation, described in the next section, alters the biome area extents so that we actually have The parameters for the preindustrial fluxes (i.e., η, µ, ρ litt , and ρ soil ) can be calibrated on nine TRENDY v2 dynamic global vegetation models Sitch et al., 2015). To do so, we use the first 30 years of the so-called "S2" simulation, in which changing climate and CO 2 are prescribed to the models but no land-use change happens. We assume the average fluxes and pools over that period are at steady-state, so that we can deduce the parameters from Eq. (23), taking κ met = 0.3/0.7 (Foley, 1995). For the few models that do not report the litter pool separately, we assume the total reported soil carbon pool is made at 5 % of litter carbon and 95 % of soil carbon. Also, to account for the harvest of croplands in OSCAR, we alter the parameters of this biome following some arbitrary rules: NPP is reduced by 80 %, thus we assume this fraction of the crops' productivity is harvested and oxidized within a year, and the mortality rate is set to 1 yr −1 , which corresponds to a yearly harvest.
Also, because the assumed preindustrial equilibrium based on TRENDY is 1901-1930 and not 1750, we scale down the NPP parameter η by a factor equal to the ratio of our preindustrial atmospheric CO 2 over the one for the TRENDY preindustrial period, i.e., by a factor of about 0.92.
The parameters for the transient response of NPP and heterotrophic respiration (i.e., β npp , β npp , CO 2cp , γ npp,T , γ npp,P , γ resp,T ,γ resp,T 1 ,γ resp,T 2 , γ resp,P ) can be calibrated on seven CMIP5 Earth system models (see, e.g., Arora et al., 2013). To do so, we use the outputs from three CMIP5 simulations: "1pctCO2", "esmFixClim", and "esmFdbk1", which correspond to simulations with an increase of atmospheric CO 2 of +1 % yr −1 in the case of a fully coupled configuration, a fixed climate, or a fixed carbon cycle, respectively. Depending on the functional form chosen, the fit for NPP is done on the basis of Eqs. (8) + (9) or (8) + (10). That of the heterotrophic respiration rate is done on the basis of Eq. (15) or (16). The calibration is done in two steps. A first fit is made with decadal moving averages of the relevant variables and for which the parameter related to local precipitation is set to zero. A second fit is then made with annual values to find the remaining parameter. This approach is used to avoid overfitting. The fit is made over the three simulations at the same time, using the "piControl" values to define the preindustrial equilibrium. In the case of the respiration rate, we also add a new term to Eq. (15) or (16) to calibrate the parameters. We multiply F resp by the term where β prim is a new sensitivity and F input is the input carbon flux from the vegetation pool to the soil pool, so as to account for the "false priming" effect observed in CMIP5 models  -that is, the effect of an increased respiration not because of increased respiration rate per se, but because of new carbon falling into a pool with faster turnover time than the average turnover time of the soil. This additional term is only used for calibration purposes and thus it is not added to OSCAR's formulation because the two soil boxes of OSCAR are expected to provide this "false priming" effect. The CMIP5 fits are shown in Figs. S5 to S11 for NPP, and Figs. S12 to S18 for heterotrophic respiration.
The fire-related parameters are similarly calibrated on TRENDY (for ι) and CMIP5 (for γ igni,C , γ igni,T , γ igni,P ) but this is done independently from the other parameters. Six models with wildfire emissions are available to calibrate on TRENDY, and four models are to calibrate on CMIP5. As was done previously, we alter the parameters obtained for croplands: we assume there is no wildfire at all within that biome. The CMIP5 fits are shown in Figs. S19 to S22 for fire intensity. Given how experimental it is to include fire processes in a model as simple as OSCAR, we also keep an option to turn off the preindustrial wildfire flux and/or its transient response.
Regarding the TRENDY and CMIP5 data processing, it has to be noted that none of the models provide biomespecific outputs. So we choose to deduce biome-specific data by weighting the biome-aggregated outputs of a model by its biome area fraction map, taken to the power 3. This approach is used to give more importance -in a given region -to the grid cells in which biomes are purer, without taking the risk of having too few of those grid cells if we were to set a threshold of biome area fraction instead. Also, some of the complex models used to calibrate OSCAR are lacking some of the biomes implemented in our model. Thus, we need rules to establish parameters for the lacking biomes on the basis of the available ones. When croplands are not in a model, we assume they have the same biogeochemical properties as grasslands, before any harvest or other human intervention. When pastures are not in a model, we assume their biogeochemical parameters are a mix of those of grasslands and bare soil, at 60 and 40 % respectively. In a configuration of OSCAR in which shrublands are separated from grasslands -which is not the case in this paper -and shrublands are not in a model, we assume they are made at 85 % of grasslands and 15 % of forests.
The preindustrial area extents A 0 are obtained by combining the preindustrial land-use map consistent with the LULCC drivers (see Sect. 2.2.2) for the anthropogenic biomes and 1 of 13 vegetation maps to distinguish between the natural biomes. The first map is used to know the fractions of "water/ice", croplands, and pastures in a given grid cell. The remaining fraction corresponds to natural vegetation, and this fraction is then subdivided into our different natural biomes following their proportions in each grid cell of the second map. Of the 13 possible vegetation maps, 2 are recent observations of land cover, MODIS (Channan et al., 2014) and ESA-CCI (2015), two are potential natural vegetation maps (Ramankutty and Foley, 1999;Levavasseur et al., 2012), and the other ones are the land-cover map of the same TRENDY models used to calibrate the preindustrial carbon fluxes and pools. In the first four cases, given that the maps provide land cover as "land-cover classes" and not as "plant functional types" -as used by TRENDY models -we use the cross-walking table developed by Poulter et al. (2011) to convert the former into the latter.

Land carbon cycle: extensive perturbation
Now we consider the extensive perturbation of the terrestrial carbon cycle, i.e., the one driven by changes in land use and land cover. This perturbation has a first-order effect that originates from the human-induced disturbance of a given biome which then transitions from its disturbed state to a new steady-state. When both extensive -change in biome extent -and intensive -change in areal properties -perturbations occur at the same time, their interaction creates a second-order effect, which is also included in the following equations. Here, we also note that in theory another extensive perturbation affects the terrestrial ecosystems: the migration of natural biomes caused by changes in environmental con-ditions (e.g., Jones et al., 2009). This is however not included in this version of OSCAR.
The book-keeping module used to estimate the carbon fluxes induced by the land-use drivers is very close to that of the previous version of OSCAR (Gasser, 2014). It is built on the approach developed by Gitz and Ciais (2003), although it now includes algorithms to treat not only land-cover change but also wood harvest and shifting cultivation. See Sect. 2.2.2 for a description of those drivers. Following the discussion and recommendation by Gasser and Ciais (2013), the bookkeeping is written so as to follow exactly the carbon fluxes and pools of transitioning ecosystems with regard to their expected but yet-to-be-reached new steady-state, so that the effect of the LUC perturbation tends toward zero in the case of infinitely old land-use disturbances. This corresponds to "definition 3" of Gasser and Ciais (2013) and to "definition B" of Pongratz et al. (2014).
For the book-keeping itself, we need to define a new series of extensive state variables for the terrestrial biosphere affected by LULCC (subscript luc ). These variables are defined following three axes: the region i axis, the biome b axis, and a new age-class a axis, so that the triplet (i, b, a) represents the "average" biome b of the i-th region that was originally disturbed at t = a. This implies that at any given time t, all the variables with a > t are nil.
The initialization of the book-keeping sequence, i.e., the initial disturbed state of a given triplet (i, b, a), depends on the kind of land-use disturbance. When land-cover change occurs, i.e., when there is conversion of a given land area from one biome b 1 to another biome b 2 (δA b 1 →b 2 ), we assume that all the living biomass of b 1 is taken away, and the living biomass of b 2 has yet to grow. When harvest occurs, we assume that the total amount of harvested biomass (δH b ) is taken from the living biomass pool of b, and this biomass will regrow in time. When shifting cultivation occurs, we assume it can be approximated by the harvest of all the living biomass over the shifting area (δS b 1 ↔b 2 ) of both biomes, b 1 and b 2 , except that the biomes are considered not to be fully grown. Their age is assumed to be equal to the shifting cultivation turnover rate (τ shift ), and thus their living biomass pool is taken equal to that of their fully grown counterpart multiplied by a factor π i,b So, the initialization of the LUC-disturbed living biomass (C veg,luc ) is as follows: www.geosci-model-dev.net/10/271/2017/ Geosci. Model Dev., 10, 271-319, 2017 In the case of land-cover change and shifting cultivation, the above-ground fraction (π agb ) of the biomass of the original biome b 1 is partly harvested and allocated to three harvested wood product pools (C hwp,luc ), following allocation coefficients (π hwp ). Each wood product pool (superscript w ) has a specific decay time (τ hwp ) that corresponds to a specific use (w = 1 is fuel wood, w = 2 is pulp-based products, w = 3 is hardwood-based products). In the case of harvest, all the harvested biomass follows the same allocation coefficients. It gives the following initialization of the harvested wood products: For all three kinds of land-use disturbance, the remaining fraction of the living biomass of the original biome b 1 is added to the litter carbon pool (C litt,luc ) of the new biome b 2 . This fraction is usually called "slash". Also, in the case of land-cover change, the soil carbon pool of b 1 has yet to transition to that of b 2 . This transition will lead to additional carbon fluxes. The initialization of the LUC-disturbed litter carbon variable is thus: Only in the case of land-cover change is the LUC-disturbed soil carbon pool (C soil,luc ) initialized by the difference in soil carbon density between the original and the new biomes. Therefore we assume harvest and shifting cultivation do not directly -i.e., at the initialization step -disturbed the soil carbon pool. So, it is simply given as follows: Here, it should be outlined that the initialization round is carbon-neutral to the atmosphere: carbon is moved between the three biospheric pools and the wood product pools, but none of it is emitted yet. Finally, the change in biome area extents is also calculated. It is by definition: Once the initialization round is done, the LUC-disturbed biospheric pools follow the same carbon cycle as the one described in the previous section for undisturbed biomes: Note that these equations are affected by environmental conditions through the F igni and F resp functions, but the arguments of these functions are not shown for legibility. There is no term for NPP in Eq. (29) because the cycle described here is the LUC-disturbed cycle (see Gasser and Ciais, 2013). Therefore, because in this version of OSCAR there is no difference in NPP between a disturbed biome and its undisturbed counterpart, the LUC-disturbed NPP is zero. As for the harvested wood products, they are oxidized at a varying rate that depends on the characteristic time of the pool (i.e., on τ hwp ) and also on a function (F hwp ) of the time passed since they were harvested (i.e., a function of t − a): The function F hwp is introduced to allow choice of the temporal profile of the wood product oxidation. For instance, if F hwp ≡ 1, the products are oxidized following an exponential profile (e.g., Houghton and Hackler, 2001). Alternatively to the exponential option, the profile can be linear (McGuire et al., 2001) or it can follow a gamma function (Earles et al., 2012). The oxidation profiles and the corresponding functions F hwp are shown in Fig. 3. The τ shift parameter is set to 15 years (Hurtt et al., 2011). The above-ground biomass fractions π agb can be calibrated on three TRENDY models, exactly as other preindustrial carbon-cycle parameters are (see Sect. 2.3.2). The allocation coefficients of the harvested wood products π hwp come from the work by Earles et al. (2012 , Table S1). Those being national, however, they are aggregated to obtain regional values by weighting them with the national estimates of aboveground biomass in forests assessed by FAO (2010, Table 2). To introduce more variation in our modeling, we have two options for processing the data. In the "low" biomass burning option, we assume all the "non-merchandable" biomass of Earles et al. (2012) becomes slash, while in the "high" biomass burning option, we assume 50 % of it is added to the fuel-wood pool (w = 1). Finally, the time constants of oxidation of the wood products τ hwp can come either from Earles et al. (2012) -which is based on the IPCC guidelines -with values of 0.5, 2 and 30 years for w = 1, 2, and 3, respectively, or from Houghton and Hackler (2001) with values of 1, 10, and 100 years, respectively.

Atmospheric CO 2 and RF
The incremental change in atmospheric CO 2 can be written as the balance between two sources (fossil-fuel and cement emissions (E FF , see Sect. 2.2.1) and land-use change emissions (E LUC )) and two sinks (the ocean sink (F ↓ocean ) and the land sink (F ↓land )). Note that despite being usually called "source" and "sink", since it is their historical role, each term of the budget can theoretically be of the opposite sign, thus changing from a source to a sink or vice versa. Mathematically it is written as follows.
where, on the basis of the three previous sections, we have the following: In Eq. (33) this version of OSCAR notably ignores the permafrost carbon that may be emitted under a warming climate (e.g., Ciais et al., 2013b). The radiative forcing (RF) induced by the increase in atmospheric CO 2 follows the logarithmic formula by Myhre et al. (1998): where α CO 2 rf = 5.35 W m −2 is given by Myhre et al. (2013b, Table 8.SM.1). For the preindustrial atmospheric concentration, we take CO 20 = 278 ppm (IPCC, 2013, Table AII.1.1a).

Non-CO 2 species
This intermediary section is dedicated to two elements which will be needed hereafter for non-CO 2 species: first, the endogenous estimate of the emission of a given species from biomass burning, and second, the estimate of the lagged concentration of a given species, assumed to be a proxy of its mid-stratospheric concentration.

Biomass burning
The atmospheric CO 2 budget above does not isolate the fluxes caused by biomass burning from those caused by all other sources of oxidation. But the biomass burning emissions are needed for non-CO 2 species in the next sections. Biomass burning emissions are altered by two aspects of the carbon cycle: one that relates to the land sink F ↓land , and one that relates to the land-use change emissions E LUC . The former can be isolated in Eq. (35) as being induced by changes in areal fire intensities and in land cover. The latter can be isolated in Eq. (36) as being induced by a change in living biomass stocks -itself induced by LULCC -and by the ox-T. Gasser et al.: Description of OSCAR v2.2 idation of the harvested wood product pool corresponding to fuel wood (w = 1). From these two CO 2 fluxes, we deduce the non-CO 2 ones by assuming that the emission of a species X from biomass burning (E X bb ) is proportional (by a factor α X bb ) to that of CO 2 , which gives the following: The α bb parameters come from the GFED v3.1 database (van der Werf et al., 2010). The biomass burning emissions of all species are averaged over the whole available timeperiod, and to each vegetation type or sector of GFED is associated a biome of OSCAR: "def" and "for" are forests, "woo" is shrublands, "sav" is grasslands, and "agr" is croplands; "pea", i.e., peatlands, are left alone. As in Sect. 2.3.2, pastures are assumed to be 60 % grasslands and 40 % bare soil. The parameters are then obtained by simply taking the ratio of the emissions of a given species over those of CO 2 .

Lagged concentrations
In the next sections, we need an estimate of the stratospheric concentration change of some species. For relatively long-lived species, we assume the stratospheric concentration change of this species can be approximated by its change in atmospheric concentration (X), albeit with a time lag (τ lag ). This change in "lagged" concentration (X lag ) is formulated as follows: This formula is a linearized form of the usual equation written with a delay (e.g., Newman et al., 2007): We opt for the linearized form because it is easier to implement in a numerical model, and because it allows the time lag to vary with time -although it is not the case in this version of OSCAR.
We set τ lag to a value of 3 years. That value corresponds broadly to the mean age of air in the mid-latitudes of the stratosphere (e.g., Newman et al., 2007). We also note that this approach to model stratospheric concentration, without an explicit representation of the stratosphere-troposphere exchange, does not hold for too short-lived species, i.e., for species with a lifetime lower than the time-lag parameter. This is one of the reasons why another approach is used for ozone in Sect. 2.8.

Methane
In OSCAR, all known atmospheric sinks of methane are included, and particular attention is paid to how the main tro-pospheric sink varies with anthropogenic and natural external factors. Amongst the natural sources of methane, only the emissions from wetlands are endogenously calculated in the model, implying that all other natural sources -such as freshwaters, termites, or permafrost -are assumed to remain constant.

Atmospheric sinks
The oxidation of atmospheric methane follows the same modeling approach as that of the previous version (Gasser, 2014) or of other simple models (e.g., Meinshausen et al., 2011). It is represented by a one-box model with one specific lifetime associated with each oxidative process. Those lifetimes may vary with time so that the resulting model is not linear.
The flux of oxidized CH 4 (F CH 4 ↓ ) is caused by four processes (e.g., Prather et al., 2012): tropospheric oxidation by the hydroxyl radical (preindustrial lifetime τ CH 4 OH ), stratospheric oxidation (τ CH 4 hν ), oxidation in dry soils (τ CH 4 soil ), and oxidation in the oceanic boundary layer (τ CH 4 ocean ). Transient change in the availability of hydroxyl radicals in the troposphere is a function (F OH ) of external factors: the atmospheric CH 4 concentration itself (CH 4 ); the stratospheric ozone concentration (O 3 s) which drives the actinic flux partially generating OH; global surface temperature (T G ) which is used to estimate changes in global atmospheric temperature and relative humidity; and emission of the three ozone precursors, represented in the form of another function (F prec ) for now. For the stratospheric sink, the lagged CH 4 concentration is used instead of the atmospheric one, and its actual lifetime is also a function (F hν ) which rationale and formulation are detailed in Sect. 2.6.1. Using also the atmospheric conversion factor α CH 4 atm defined in Sect. 2.3.1, we can write the following: The function F OH mostly follows the formulation by Holmes et al. (2013). It is parameterized with chemical sensitivities of OH to atmospheric CH 4 (χ OH CH 4 ), stratospheric O 3 (χ OH O 3 s ), global atmospheric temperature (χ OH T A ), and global atmospheric relative humidity (χ OH Q A ). The absolute change in global atmospheric temperature (T A ) is assumed to be proportional (by a factor κ T A ) to that in global surface temperature. The relative change in global atmospheric relative humidity is assumed to be proportional (by a factor κ Q A ) to that Geosci. Model Dev., 10, 271-319, 2017 www.geosci-model-dev.net/10/271/2017/ in saturation vapor pressure. The latter follows an empirical function of global atmospheric temperature change with two parameters (κ svp and T svp ). So far this can be written as follows: The functional form of F prec can be either linear (Ehhalt et al., 2001) or logarithmic (Holmes et al., 2013). In the linear case, it is a function parameterized with three absolute chemical sensitivities of OH to the ozone precursors: nitrogen oxides (χ OH NO x ), carbon monoxide (χ OH CO ), and nonmethane volatile organic compounds (χ OH VOC ): In the logarithmic case, it is parameterized by three relative chemical sensitivities ( χ OH NO x , χ OH CO , χ OH VOC ), and the preindustrial natural emissions of the three ozone precursors (E NO x nat , E CO nat , E VOC nat ). This gives the following equation: None of these two formulations shows regionalized chemical sensitivities of the OH sink, however, whereas in reality the sink is sensitive to where the ozone precursors are emittedespecially the NO x (e.g., Wild et al., 2001). The four lifetimes of methane are taken as the present-day lifetimes given by Prather et al. (2012 , Tables A1 and A2): 11.2, 120, 150, and 200 years for τ CH 4 OH , τ CH 4 hν , τ CH 4 soil , and τ CH 4 ocean respectively. The lifetime with regard to OH is then scaled down by an arbitrary factor of 0.80. We note that this does not follow the rescaling made by Prather et al. (2012) which was based on preliminary results from the ACCMIP models . The ACCMIP study is inconclusive about the change in methane lifetime between the preindustrial period and present day: some models predict an increase while others predict a decrease. Because our function F OH is based on a subset of the ACCMIP models (see below) which all find the methane OH lifetime increased, we scale down the preindustrial value of this lifetime, so that it roughly meets its present-day value during the simulation. Also, to introduce variation in this important parameter, we propose alternative values based on the ACCMIP chemistrytransport models   Table 1): optionally, the default lifetime can be rescaled by a factor equal to any of the 16 models' estimates of the lifetime over the multimodel mean estimate. Finally, the stratospheric lifetime is also scaled up by a factor 1.06, following Prather et al. (2015) (see also Sect. 2.6.1).
All the chemical sensitivities of the OH sink (i.e., χ OH CH 4 , χ OH CO and χ OH VOC ) are taken as one of the four sets of values from the study by Holmes et al. (2013, Table 2). Alternatively, for backward compatibility, these parameters can also be taken as the multi-model mean estimates from the OxComp project (Ehhalt et al., 2001, Table 4.11), in which case the sensitivities to temperature, humidity, and ozone are nil. The preindustrial global atmospheric temperature T A,0 is set to 251 K, and the proportionality coefficients are κ T A = 0.94 and κ Q A = 1.5 (Holmes et al., 2013). The saturation vapor pressure parameters are obtained from Jacobson (2005, Eq. 2.62), for which a small temperature perturbation is assumed, giving κ svp = 17.67 and T svp = −29.65 K. The preindustrial stratospheric ozone burden O 3 s 0 is set to 280 DU, roughly following Cionni et al. (2011).

Wetland emissions
Natural wetlands are the largest natural source of methane (Ciais et al., 2013b), and the future variation of this source could be significant for future climate change (e.g., O'Connor et al., 2010). We thus decided, since version 2.1 of OSCAR, to include a simple module describing the variation of this source of CH 4 . The current version is very close to the previous one (Gasser, 2014), except that a larger variety of parameterizations is now available. First, we estimate the regional change in CH 4 emission per unit area of wetlands (e wet ) as being proportional to its preindustrial value and to the relative change in heterotrophic respiration of the litter carbon pool in the same region. To this end, wetlands are considered to be a mix of the other biomes, with partition coefficients (π b wet , b π b wet = 1) having a nonzero value only for natural biomes. We note that this is an ad hoc assumption that we make because we lack detailed outputs from complex wetland models. The litter pool is chosen as a proxy of the changes in wetlands induced by more general changes in the carbon cycle. Therefore, here we implicitly assume that the sensitivity of areal wetland emissions to environmental conditions -e.g., CO 2 or temperature -is the same as that of heterotrophic respiration. So we have the following: Second, we assume that the regional change in wetland area extent (A wet ) depends on linear sensitivities to: atmospheric CO 2 (γ wet,C ), local surface temperatures (γ wet,T ), local yearly precipitation (γ wet,P ). This formulation is similar to that used for wildfire intensity in Eq. (12), and CO 2 is used as a proxy of changes in, e.g., evapotranspiration or vegetation species distribution. Mathematically it is written as follows: Consequently, the change in regional emission of methane by wetlands (E wet ) is calculated as follows: We calibrate two sets of parameters for wetlands. First, the preindustrial equilibrium of the wetlands can be calibrated on seven WETCHIMP models (Melton et al., 2013). We deduce the π wet parameters by combining the wetland map from the "exp 1" simulation, that is the control experiment of the WETCHIMP exercise, and the land-cover map used in Sect. 2.3.2 for natural vegetation. The preindustrial areal emissions e wet,0 are also taken from this "exp 1" simulation, but they are scaled down by a factor equal to the ratio of our preindustrial atmospheric CO 2 over the one used in WETCHIMP, i.e., by a factor of about 0.92, as we did with NPP in Sect. 2.3.2. Second, the parameters for the transient response of wetland extent (i.e., γ wet,C , γ wet,T , γ wet,C ) can be calibrated on six WETCHIMP models (reminder: see Appendix B for a list of those models). To do so, we use "exp4", "exp5", and "exp6": factorial simulations that separate the effect of temperature, precipitation, and atmospheric CO 2 , respectively. For the same reasons as with wildfires, we also keep an option to turn off the preindustrial wetland flux and/or its transient response.

Atmospheric CH 4 and RF
On the basis of the previous sections, the incremental change in atmospheric CH 4 follows the mass-balance equation: This equation implicitly assumes that all the natural sources of methane except for natural wetlands remain unchanged since the preindustrial period. Here, we also note that the anthropogenic emissions E CH 4 do include emissions from rice paddies -i.e., from anthropogenic wetlands.
The radiative forcing induced by the increase in atmospheric CH 4 follows a square-root formula to which an ad hoc function (F over ) is added to account for the overlap between the absorption bands of methane and nitrous oxide (N 2 O), following Myhre et al. (1998). It gives the following: where α CH 4 rf = 0.036 W m −2 ppb −0.5 and the analytical expression of F over are given by Myhre et al. (2013a, Table 8.SM.1). In addition to the RF induced by methane itself, we have to account for the RF induced by the increase in stratospheric water vapor caused by the oxidation of methane. To do so, as others have (e.g., Meinshausen et al., 2011), we assume it is equal to 15 % of the direct methane RF, but calculated with its lagged concentration: where α H 2 Os For the preindustrial atmospheric concentration, we take CH 40 = 722 ppb (IPCC, 2013, Table AII.1.1a).

Nitrous oxide
In OSCAR, the stratospheric sink of nitrous oxide is included, and particular attention is paid to how it varies with anthropogenic and natural external factors. However, no other natural processes are endogenous to the model, meaning that no change in natural sources or sinks of nitrous oxide (e.g., ocean, natural soils, biological fixation) is assumed.

Stratospheric sink
The oxidation of nitrous oxide follows the same modeling approach as that of methane, with only one sink in the stratosphere that has a varying lifetime. The law used to make the stratospheric lifetime vary, however, is recent and different from the previous version of the model.
is driven by the preindustrial lifetime of nitrous oxide with regard to stratospheric oxidation (τ N 2 O hν ). The transient change in this stratospheric lifetime is a function (F hν ) of the lagged N 2 O concentration (N 2 O lag ), the equivalent effective stratospheric chlorine (EESC, see Sect. 2.8.2), and global surface temperature change (T G ). The dependency on N 2 O and the EESC is meant to model the impact of a change in stratospheric ozone that changes the actinic flux, which in turn changes the stratospheric sink (e.g., Prather, 1998 The formulation of F hν is inspired by that used for methane and the study by Prather et al. (2015). It has three chemical sensitivities (χ hν N 2 O , χ hν EESC , and χ hν age ). This last parameter represents the sensitivity of the sink to a change in stratospheric age of air. This age-of-air change is itself driven by a changing Brewer-Dobson circulation which is induced by a changing climate (e.g., Butchart, 2014). In the following, we consider that the inverse of the relative change in age of air is a linear function of the absolute change in global surface temperature (parameterized by γ age , see also Fig. S23). This leads to the following formula: The preindustrial stratospheric lifetime τ N 2 O hν is taken as 123 years . As we do with methane, we introduce variation in the N 2 O lifetime by having the option to rescale the default value by a factor equal to the lifetime simulated by any of the eight models of Prather et al. (2015, Table 2) over the multi-model mean estimate. The first two chemical sensitivities of the stratospheric sink (i.e., χ hν N 2 O and χ hν EESC ) are taken as one of the four sets of values from the study by Prather et al. (2015). Three sets of value are given in their Table 3, and the fourth is the recommendation in their text. Also, to translate their Table 3 into our parameters, we assume that the preindustrial EESC in the models was 420 ppt -from IPCC (2013,  Table 1). Alternatively, for backward compatibility, these parameters can also follow (Prather et al., 2012, Table A1), in which case the sensitivity to EESC is zero.
Regarding the chemical sensitivity to the age of air, we assume it is not zero only when the other sensitivities are deduced from the "G2d" model, therefore following the results of Prather et al. (2015, Table 3) and their discussion pointing out the experimental aspect of such a parameterization. Nevertheless, in this specific case we need further information about the "G2d" model which we take from Fleming et al. (2011, Fig. 12), where one can see that the age of air at an altitude of 25 km changed from about 4.5 to 4.0 between the preindustrial and present-day periods. This is enough to deduce the χ hν age parameter. Then, the γ age parameter can be calibrated on seven CCMVal2 chemistry-transport models (Morgenstern et al., 2010). To do so, we use outputs from the "REF-B2" experiment which is a fully transient simulation over 1961-2099: we use the "mean_age" output at a pressure level of 25 hPa (∼ 25 km) and the temperature at the surface level. We then fit the parameter following our inverse linear relationship, defining the preindustrial conditions as the averaged first 10 years of the simulations. The CCMVal2 fits are shown in Fig. S23.

Atmospheric N 2 O and RF
The incremental change in atmospheric N 2 O is as follows: We note again that this implicitly assumes natural emissions remain unchanged since the preindustrial period. Similarly to methane, the radiative forcing induced by the increase in atmospheric N 2 O follows a square-root formula to which the ad hoc overlap function is added: where α N 2 O rf = 0.12 W m −2 ppb −0.5 and F over is given by Myhre et al. (2013b,

Atmospheric sinks
Conceptually, the modeling approach of the halogenated compounds' sinks is similar to that used for methane. Each of these species X is affected by three sinks, each sink with its specific preindustrial lifetime: a tropospheric oxidation by the hydroxyl radical (τ X OH ), a stratospheric oxidation (τ X hν ), and another sink which encloses all other processes such as oxidation in dry soils or in the oceanic boundary layer (τ X othr ). Note that a given oxidation process may not actually affect a given species; in this case the associated lifetime is set to a value of infinity (∞) -or equivalently its loss frequency (ν X = 1/τ X ) is set to zero. Mathematically, similarly to Eq. (40), we have the following for any species X being a HFC, PFC, or ODS: where the functions F OH , F prec , and F hν are the same as in Sect. 2.5.1 and 2.6.1. The lifetimes τ X OH , τ X hν , and τ X othr are taken from the compilation by Montzka et al. (2011, Table 1.3). However, the lifetimes with respect to the OH sink are rescaled using the same scaling factor as for methane (see Sect. 2.5.1), for consistency. Similarly, the lifetimes with respect to the stratospheric sink are scaled up by a factor 1.06, as done by Prather et al. (2015).

Atmospheric concentrations and RFs
The incremental change in atmospheric concentration of any species X being a HFC, PFC or ODS is as follows: With the exception of CF 4 , CH 3 Br, and CH 3 Cl, all the halogenated compounds are anthropogenic in nature, and thus no other natural fluxes need to be considered. For the three former species, however, their natural emissions are assumed to remain constant through time. The radiative forcing induced by the increase in atmospheric concentration of any of those species X is assumed to be proportional: where the values of α X rf are taken from Myhre et al. (2013b, Table 8.A.1). In the following, all these RFs will be combined into one: Finally, only the three species cited above have nonzero preindustrial atmospheric concentration: CF4 0 = 35 ppt (IPCC, 2013, Table AII.1.1a), CH3Br 0 = 5.8 ppt, and CH3Cl 0 = 480 ppt .

Ozone
In OSCAR, as in other simple models (e.g., Meinshausen et al., 2011), tropospheric and stratospheric ozone are estimated separately. As it is common in simple models, shortlived species such as ozone are not predicted using a dynamic model like long-lived species are. Rather, at each time step, the short-lived species are supposed to be at chemical steadystate with their drivers of change. This implies notably that no stratosphere-troposphere exchange is implemented in the model.

Tropospheric O 3 and RF
In our model, change in tropospheric ozone is driven by atmospheric methane, emissions of ozone precursors, and global climate change. We use a formulation close to that of the previous version of OSCAR, which was the formulation by Ehhalt et al. (2001). In version 2.2, however, the chemical sensitivities are updated and regionalized, and a sensitivity to climate change is added. The change in global tropospheric ozone burden (O 3 t) is a function of the following: atmospheric methane, with a logarithmic sensitivity (χ O 3 t CH 4 ); global surface temperature, with a linear sensitivity ( O 3 t ); and the three ozone precursors, with linear global sensitivities (χ that are regionalized thanks to region-specific weights (ω NO x , ω CO , ω VOC ). Here, we introduce a new regional axis (superscript r ) that is de facto different from the biospheric one (superscript i ). The regional axes are linked through parameters describing what fraction of a region i is actually included in a region r (π r,i reg , r π r,i reg = 1). So we finally have the following: The global chemical sensitivities (i.e., χ can be calibrated on four ACCMIP chemistrytransport models . To do so, we use their reference simulations for the year 2000, as well as the factorial simulations which were made so as to isolate each of the four drivers of the change in tropospheric ozone (namely "1850CH4", "1850NOx", "1850CO", and "1850NMVOC"). However, since we also have access to a simulation in which all of the four drivers vary at the same time -i.e., the difference between the experiments for 1850s and the 2000s -we can estimate the non-linearity of this chemical system. We account for this non-linearity by rescaling the individual sensitivities by a factor equal to the ratio of the ozone change in Geosci. Model Dev., 10, 271-319, 2017 www.geosci-model-dev.net/10/271/2017/ the all-varying simulation over the sum of the ozone changes in each of the factorial simulations. The sensitivity to global climate change O 3 t can be calibrated on eight models which participated in the same ACCMIP exercise and made simulations in which only climate varies. A simple linear fit is made over these simulations; we also keep an option to set this sensitivity to zero. The latter ACCMIP fits are shown in Fig. S24. The regional weights ω X can be deduced from the results of 11 HTAP chemistry-transport models (Fiore et al., 2009). To do so, we calculate regional ozone changes in the four HTAP regions thanks to Table S5 from Fry et al. (2012) and regional precursors emissions thanks to Table S1 from Fiore et al. (2009). Our weighting parameters are then deduced as the ratio of the regional ozone changes normalized by the precursors changes over the globally averaged normalized ozone change. A fifth region is then added, to account for areas of the globe that are not within the four HTAP regions, and for which the weighting parameter is set to exactly 1. The π reg parameters are logically defined as the fraction of the area of a region i that is inside a region r. Also, we keep an option to turn off that regionalization, i.e., setting all regional weights to 1.
Finally, the radiative forcing induced by the change in tropospheric ozone burden is assumed to be linear: where the value of α O 3 t rf is not unique -contrarily to what we do with greenhouse gases. This radiative efficiency can be the following: 0.042 W m −2 DU −1 , as reported by Myhre et al. (2013b); 0.032 W m −2 DU −1 , as reported by Forster et al. (2007); or one of the 15 radiative efficiencies given by the ACCMIP chemistry-transport models   Table 3).

Stratospheric O 3 and RF
With the same formalism as for tropospheric ozone, change in stratospheric ozone is driven by available stratospheric chlorine and bromine, stratospheric nitrous oxide, and global climate change. Compared to the previous version that was using only a linear dependency on chlorine and bromine, nitrous oxide and climate change are two new drivers in this module.
The first step to model stratospheric ozone is to estimate its first driver of change: the stratospheric chlorine and bromine available from the presence of the ODSs in the stratosphere. Those compounds release their chlorine and/or bromine atoms at various rates and thus interact differently with ozone. A proxy variable is thus created to lump together these various effects, namely the equivalent effective stratospheric chlorine (EESC). The EESC is calculated following Newman et al. (2007), on the basis of the fractional release of each ODS (π rel ), its numbers of chlorine atoms (n Cl ) and bromine atoms (n Br ), a parameter measuring the efficiency in destroying ozone of bromine relative to that of chlorine (α Br Cl ), and the lagged concentration of the ODS: Then, a change in stratospheric ozone burden (O 3 s) is assumed to happen with a change in EESC, with a linear sensitivity (χ O 3 s EESC ). To the effect of ODSs, we add the effect of nitrous oxide following the simple formulation by Daniel et al. (2010), which needs two additional parameters: one to quantify the linear sensitivity of stratospheric ozone to nitrous oxide (χ O 3 s N 2 O ), and one to account for the non-linear interaction between chlorine and nitrogen chemistries (EESC × ). As for tropospheric ozone, a linear sensitivity to global surface temperature change ( O 3 s ) is also added, which sums up as follows: Regarding the EESC parameterization, Newman et al. (2006, Tables A1 and A2) provide values of fractional release π rel for all our ODSs, assuming a mean age of air of 3 years taken equal to the time lag of Sect. 2.4.2. To introduce other possibilities of parameterization in the model, we can alternatively take fractional release values from Laube et al. (2013), either the values for the mid-latitudes or those for the high latitudes. In this case, if a value is missing for a given ODS we take that from Newman et al. (2006). The chemical formula of each ODS gives n Cl and n Br . We take α Br Cl = 60 (Newman et al., 2007). The chemical sensitivity of stratospheric ozone to EESC and that to global climate change (i.e., χ O 3 s EESC and O 3 s ) can be calibrated on 11 CCMVal2 chemistry-transport models studied by Douglass et al. (2014), using the results from their multi-linear regression. The sensitivity to nitrous oxide is calculated using the formula by Daniel et al. (2010) 10.4 ppt ppb −1 and EESC × 2642 ppt. Also, we add two extra options to this module: one for which this response of stratospheric ozone to nitrous oxide is simply turned off, and one for which it is assumed to be linear -instead of saturating -by setting the EESC × parameter to infinity.
Finally, the radiative forcing induced by the change in stratospheric ozone burden is assumed to be linear: www.geosci-model-dev.net/10/271/2017/ Geosci. Model Dev., 10, 271-319, 2017 where the radiative efficiency α O 3 s rf can be 0.004 W m −2 DU −1 , as reported by (Forster et al., 2007), or one of the four radiative efficiencies given by the ACCENT models (Gauss et al., 2006, Tables 4 and 6).

Aerosols
As for ozone, because the aerosols are short-lived, it is assumed that their global atmospheric burden reaches a steadystate with their respective drivers of change at each time step of the model. The direct and indirect radiative effects of five anthropogenic aerosols are considered here, with the notable caveat that their atmospheric physico-chemistry is only loosely coupled to the rest of OSCAR: for instance, the oxidation of SO 2 into SO 4 is implicitly independent of the OH availability, there is no interaction between the sulfate and nitrate chemistries, and the cloud effects are independent of any change in cloud cover induced by other species. Also notable is that this version of OSCAR assumes no change in natural aerosols and aerosols precursors such as mineral dusts, sea salt, dimethyl sulfide, or biogenic VOCs.

Direct effect
The direct effect of aerosols refers to the direct radiative forcing caused by the aerosol-radiation interactions, i.e., without consideration of any short-term adjustment of the climate system . This section describes how atmospheric burdens and the resulting RF of five anthropogenic aerosols, namely sulfate aerosols, primary organic aerosols, black carbon, nitrate aerosols, and secondary organic aerosols, are calculated within our model.
It must be noted that here we purposefully limit the number of these drivers of change: only two precursors are considered for each aerosol, to avoid over-fitting on data, which does not allow us to clearly separate the effect of each precursor, and we add the global surface temperature, used as a proxy of a changing climate. For the same reason -because of the calibration data -we keep the modeling simple with linear sensitivities. Note also that in this section every lifetime is said to be "apparent", because it corresponds to a globally averaged chemical sensitivity that has dimensions of time, and which results from several physical and/or chemical processes not explicitly modeled in OSCAR.
In the case of sulfate aerosols, their change in burden (SO 4 ) is parameterized by the apparent lifetime of sulfur dioxide (τ SO 2 ) with a regionalized weighting (ω SO 2 , analogous to that used for tropospheric ozone in Sect. 2.8.1), the apparent lifetime of dimethyl sulfide (τ DMS ), and their sensitivity to global surface temperature ( SO 4 ). So we have the following: The change in burden of primary organic aerosols (POA) is parameterized by the apparent lifetime of fossil-based organic matter (τ OM,ff ) also regionally weighted (with weights ω OM ), the apparent lifetime of pyrogenic organic matter (τ OM,bb ), and their sensitivity to global surface temperature ( POA ), as well as a factor used to convert organic carbon to organic matter (α OM OC ): The change in burden of black carbon (BC) is parameterized by the apparent lifetime of fossil-based BC (τ BC,ff ) also regionally weighted (with ω BC ), the apparent lifetime of pyrogenic black carbon (τ BC,bb ), and their sensitivity to global surface temperature ( BC ): In the case of nitrate aerosols, inspired by Shindell et al. (2009), we assume their formation is driven by nitrogen oxides and ammonia emissions, and therefore we uncouple the nitrate and sulfate chemistries even though they are coupled in reality . Hence, the change in burden of nitrate aerosols (NO 3 ) is parameterized by the apparent lifetime of nitrogen oxides (τ NO x ), the apparent lifetime of ammonia (τ NH 3 ), and their sensitivity to global surface temperature ( NO 3 ). So we have the following: Geosci. Model Dev., 10, 271-319, 2017 www.geosci-model-dev.net/10/271/2017/ Finally, the change in burden of secondary organic aerosols (SOA) is parameterized by the apparent lifetime of anthropogenic NMVOCs (τ VOC ), the apparent lifetime of biogenic NMVOCs (τ BVOC ), and their sensitivity to global surface temperature ( SOA ). Here, the dependency of SOA on other factors such as atmospheric NO x or POA  is neglected. So we have the following: Finally, here it must be noted that, despite being used for the calibration (see below) and being shown in Eqs. (63) and (67), DMS and BVOC emissions are constant in this version of OSCAR. In other words, in any simulation with OSCAR v2.2 we have E DMS = 0 and E BVOC = 0. In this version, we also do not model any change in natural aerosols, i.e., in mineral dust and sea salt. For SO 4 , POA, and BC, the apparent global lifetimes τ X and the climate sensitivities X can be calibrated on four CMIP5 or ACCMIP chemistry-climate models. To do so, we use the yearly outputs from the historical and RCP simulations, assuming the average of the first 10 years is our preindustrial equilibrium. We then fit the parameters on the basis of Eqs. (63), (64), or (65), and over all the simulations at the same time. For SOA, it is done in the same way, except that only two models are available. Additionally, because of our very low confidence in the SOA modeling, we also keep an option to turn it off. For NO 3 we use other simulations and models: we do the exact same fit with the input and output data from either Bellouin et al. (2011) or Hauglustaine et al. (2014. In the latter case, NO 3 is set to zero because climate does not vary in the available simulations. The conversion factor α OM OC -which is the same here for fossil-based and biomass burning emissions -can take three values: a default and widely used value of 1.4; 1.3 ); or 1.6 (Rotstayn et al., 2012). The CMIP5/ACCMIP fits are shown in Figs. S25 to S28.
The regional weights ω X can be deduced from the results of seven HTAP chemistry-transport models . To do so, for the four HTAP regions, we take the normalized aerosol-induced RF data from the detail of their Table 6. Our weighting parameters are then deduced as the ratio of the regional normalized RF over the globally averaged normalized RF. A fifth region is then added, to account for areas of the globe that are not within the four HTAP regions, and for which the weighting parameter is set to exactly 1. The π reg parameters are the same as in Sect. 2.8.1. Also, we keep an option to turn off that regionalization, i.e., setting all regional weights to 1.
For any of the five aerosols X described in this section, the direct radiative forcing induced by a change in atmospheric burden is assumed to be linear: where the radiative efficiencies α X rf are taken from the Aero-Com II intercomparison (Myhre et al., 2013a). This leads to 15 possible parameters for SO 4 (their Table 4), 15 for POA (their Table 6), 15 for BC (their Table 5), 8 for NO 3 (their Table 8) and 5 for SOA (their Table 7).

Cloud effects
Under this term, we group the so-called semi-direct and indirect effects -that is, the rapid adjustments in the atmospheric system induced by aerosol-radiation interactions and the adjusted aerosol-cloud interactions, according to the terminology by Boucher et al. (2013) (see also Sherwood et al., 2015). The formulation we propose here is new to the model.
For the semi-direct effect, the modeling approach is straightforward. According to Boucher et al. (2013) this effect can largely be attributed to absorbing aerosols, i.e., to BC in our model. We thus account for this effect simply by adding a RF term that is proportional (by a factor κ BC adj ) to the direct RF of BC. For the aerosol-cloud interactions, the modeling is done in two steps. First, we estimate the change in atmospheric burden of soluble aerosols (AER sol ) thanks to soluble fractions specific to each type of anthropogenic and natural aerosol (π X sol ). It gives the following: Second, inspired by several studies (Boucher and Pham, 2002;Hansen et al., 2005;Carslaw et al., 2013;Stevens, 2015), we assume the aerosol-cloud interaction effective RF varies with the logarithm of the change in this soluble aerosols burden, parameterized by the intensity of the effect ( ) and the preindustrial soluble aerosols burden. This logarithmic functional form represents a saturating yet not bounded capacity of the emitted hydrophilic aerosols to alter the clouds' albedo (see, e.g., Carslaw et al., 2013, Fig . 3). Finally, in OSCAR the RF of the two combined cloud effects is therefore formulated as follows: Note that the cloud effects are estimated independently from any change in cloud cover that is happening implicitly in the climate system module. One possible value for the coefficient used to account for the semi-direct effect is based on the fifth IPCC report : κ BC adj = −0.1/0.6. However, so as to introduce variation around this effect, we also add parameterizations based on the study by Lohmann et al. (2010). Using data from their Fig. 2, we multiply the IPCC-based value by one of the five models' estimate of the effect and divide it by the multi-model mean estimate, thus obtaining five alternative parameterizations. The derivation of the parameters for the aerosol-cloud interaction is done in three steps. First, we need the soluble aerosol fractions π X sol : they are taken either from the study by Hansen et al. (2005) or from that by Lamarque et al. (2011). When taken from Hansen et al. (2005, Sect. 3.3.1), we assume that the soluble fraction of BC is a mix in equal shares of that of fossil BC and biomass burning BC, SOA has the same soluble fraction as POA, and mineral dust -not modeled by OSCAR but necessary here to deduce AER sol,0 -has a soluble fraction of zero. When taken from Lamarque et al. (2011), all soluble fractions are equal to one, except for POA and BC whose solubility is taken as the percentage of hydrophilic aerosol provided by the study, and for mineral dust and sea salt whose solubility is taken as the percentage of aerosol with a diameter < 1 µm.
Second, we calculate the intensity parameter and a preliminary value of AER sol,0 using results from ACCMIP and CMIP5 models presented by Shindell et al. (2013). Using their Table 7, we can base our parameters on one of seven ACCMIP/CMIP5 estimates of the indirect aerosol RF over 1850-2000, or on their multi-model mean. However, because these estimates are far from the IPCC best guess , the chosen ACCMIP/CMIP5 value is rescaled by a factor equal to the IPCC best guess divided by the multimodel mean. We then extract from the ACCMIP or CMIP5 outputs the atmospheric burden of each aerosol type simulated by the chosen model. These burdens are then combined using our own solubility fractions to calculate the soluble aerosols burden in the years 1850 and 2000. These two points in time, combined with the previously rescaled RF estimate, are enough to deduce through the logarithmic formula. The soluble aerosols burden in 1850 is our preliminary value of AER sol,0 .
Third, because this preliminary value of AER sol,0 is for the year 1850 and not the year 1750, we rescale it by a factor from the study by Carslaw et al. (2013) and adapted to the logarithmic formula; its value is exp[(1.42 − 1.30)/ ] and it is named the "median" option. Again, in order to introduce variation in our modeling of the indirect effect, we also propose two other arbitrary rescaling options: one with no actual rescale, named "high", and one with the rescale factor applied twice, named "low". With these three steps, we expect to introduce enough variation for the model to cover a wide range of possible future evolutions of the aerosolcloud interactions, i.e., to span a large domain of the Fig. 3 of Carslaw et al. (2013) as it is illustrated in our Fig. 4.

Surface albedo
Anthropogenic perturbations of the Earth's energy budget through surface albedo change are difficult to model in a simple way, because they are local phenomena with significant seasonal variability. Moreover, they can involve nonradiative processes that are almost impossible to capture with simple models. The two OSCAR modules presented hereafter are first-order models of two surface albedo perturbations: BC deposition on snow, and land-cover change. As such, they are not coupled with one another, nor are they with the climate module.

Black carbon on snow
The radiative forcing induced by BC deposition on snow is taken directly proportional to the regional BC emissions. It is parameterized by a global radiative efficiency with respect to emissions (α BCsnow rf ), and further regionalized by regionspecific weights (ω BCsnow ). Mathematically this is written as follows: where the regionalization (superscript r ) is specific to this module, and therefore different from the regionalization based on HTAP seen in previous atmospheric chemistry modules (Sect. 2.8.1 and 2.9.1). The global radiative efficiency with respect to emissions α BCsnow rf can be taken from eight ACCMIP models   Table 3 and Fig. 15). The regional weights ω BCsnow are obtained from the study by Reddy and Boucher (2007 ,  Table 1). As in Sect. 2.9.1, the weighting parameters are deduced as the ratio of the regional radiative efficiencies over the globally averaged radiative efficiency. A 10th region is added, to account for areas of the globe that are not within the nine regions of Reddy and Boucher (2007), and for which the weighting parameter is set to exactly 1. The π reg parameters are logically defined as the fraction of the area of a region i that is inside a region r .

Land-cover change
The radiative forcing induced by changes in land cover is modeled following the first-order equation of Bright and Kvalevåg (2013). It is parameterized by the yearly averaged albedo at the biome and regional scale (α alb ), the regional radiative shortwave and downward flux at the surface (ϕ rsds ), and the global shortwave and upward transmittance (π trans ). Here we note that both the drivers and the regional disaggregation are the same as those of the terrestrial carbon cycle, which implies the i-axis is the same as in Sect. Here we show the simulated radiative forcing as a function of the total burden of soluble aerosols (left-hand side) or of the change in that burden since preindustrial times (right-hand side). In the former case, the grey dotted lines show the preindustrial burden we calculate; in the latter, the red area shows the 90 % range of RF provided by Myhre et al. (2013b), and therefore the associated change in burden implied by our formula.
2.3.3. So we have the following: where A Earth designates the surface area of the Earth. The upward transmittance is set to π trans = 0.854 (Lenton and Vaughan, 2009) We calculate the yearly-averaged biomespecific albedos by weighting the albedo climatology by one of two land-cover climatologies -either MODIS (Channan et al., 2014) or ESA-CCI (2015) -and by the radiation climatology used for ϕ rsds , in a similar fashion to He et al. (2014). This approach ensures that the yearly-averaged albedo accounts for the local seasonality, and especially that of snow cover. Also, regarding the deduction of biome-specific albedos, three more assumptions are made: we apply the same weighting method of the land-cover fraction as in Sect. 2.3.2; we remove the grid cells that see less than 1 % of their area changing over the historical period and the RCPs according to our LULCC dataset (Hurtt et al., 2011); and pastures are assumed to be made at 60 % of grasslands and 40 % of bare soils.

Climate
The climate module of OSCAR is relatively simple compared to other models, as the energy budget is done only on a global scale and no water budget is explicitly done. The two-box model used to estimate global temperature change is parameterized by -among other things -an equilibrium climate sensitivity. Therefore, all the feedbacks occurring within the climate system, such as changes in tropospheric water vapor, in ice cover, or in cloud cover are implicitly accounted for in a linear manner and with the assumption that all climate forcers induce the same level of feedback. Similarly, in this version of OSCAR, the regional patterns of temperature or precipitation changes are the same whatever the climate forcer. Also, no non-radiative biophysical or physiological effect relative to the land vegetation -e.g., a change in evapotranspiration induced by a change in land cover or in atmospheric CO 2 -is included in the model.

Radiative forcings
The first step to calculate global warming is to calculate global radiative forcing. So as to ease the notations, following Myhre et al. (2013b), we introduce two groups of anthropogenic forcings: the well-mixed greenhouse gases (WMGHGs) for which radiative forcing is defined as follows: The near-term climate forcers (NTCFs) for which radiative forcing is defined as follows: where the last three terms are the three drivers directly prescribed to OSCAR as radiative forcing and detailed in Sect. 2.2.3.
To estimate global warming, however, we have to account for the so-called "efficacy" of these forcings, i.e., we have to introduce new parameters (κ X warm ) that measure the relative efficiency at warming the Earth of a given RF when compared to the RF of CO 2 (see, e.g., Hansen et al., 2005;Forster et al., 2007). In OSCAR, we assume all efficacies are equal to 1 -although accounting for the semi-direct effect of BC could be defined as using an efficacy -except for the two surface albedo forcings and for volcanic aerosols. Therefore, the RF used to calculate warming (RF warm ) is as follows:  Table 7); and κ volc warm is set to an arbitrary value of 0.6 based on Gregory et al. (2016). However, regarding volcanic aerosols, we note that since the forcing is normalized to zero over the historical period in Sect. 2.2.3, its efficacy only influences the variability of our results and not the trend. Now, to estimate global precipitation change, we also need to estimate how much of this top-of-the-atmosphere RF is actually occurring within the atmosphere -thus creating a local energy imbalance -in opposition to the RF occurring at the Earth's surface. To do so, we introduce new parameters that quantify this atmospheric fraction for several groups of forcers: carbon dioxide alone (π CO 2 atm ); all the other longlived greenhouse gases, i.e., methane, nitrous oxide, and the halogenated compounds (π noCO 2 atm ); tropospheric ozone alone (π O 3 t atm ); stratospheric greenhouse gases, i.e., stratospheric water vapor and ozone (π strat atm ); scattering aerosols, i.e., sulfate, primary organic, nitrate, secondary organic, and volcanic aerosols (π scatter atm ); absorbing aerosols, i.e., black carbon (π absorb atm ); cloud-related forcings (π cloud atm ); forcings from surface albedo change (π alb atm ); and the solar forcing (π solar atm ). The atmospheric radiative forcing (RF atm ) consequently is as follows: We base our grouping of the forcers on Allan et al. (2013). This grouping assumes that the atmospheric fraction π atm of CH 4 applies for all non-CO 2 long-lived greenhouse gases and that of SO 4 applies for all scattering aerosols. Additionally, we assume that cloud, albedo-based and stratospheric forcers have a nil atmospheric fraction. Other than that, the atmospheric fractions are taken from Andrews et al. (2010, Table 3) or Kvalevåg et al. (2013, Table 2, case of highest perturbation), although in the latter case tropospheric ozone is also given a nil fraction.

Surface temperatures
Similarly to what is done in other simple models (e.g., Raupach et al., 2011, although not in MAGICC;Meinshausen et al., 2011), in OSCAR the global surface temperature change is based on an impulse response function (IRF) calibrated on more complex global circulation models. The impulse response function, however, is hereby coded as a twobox model, but theoretically speaking it is strictly equivalent (see Geoffroy et al., 2013). For regional temperatures, we use a simple linear approach. This is equivalent to pattern scaling, although an important limitation is that the pattern used is the same whatever the climate forcer.
The two-box model used to model the global surface temperature change has two state variables: the global surface temperature itself (T G ), and the temperature of the deep ocean (T D ). It is parameterized by the climate sensitivity (λ), the time-inertia of the surface box (τ T G ), that of the deep box (τ T D ), and a coefficient describing the exchange of energy between the surface and deep boxes (θ ). Mathematically, it is formulated as follows: So as to deduce the change in sea surface temperature (T S ) and in local surface temperatures (T L ) for each of our land regions (the i axis from Sect. 2.3.2), we use regional weighting coefficients (ω T S and ω T L , respectively): Geosci. Model Dev., 10, 271-319, 2017 www.geosci-model-dev.net/10/271/2017/ The first set of parameters of this module, for global temperature, can be calibrated on 25 CMIP5 global circulation models. First, using outputs from the "abrupt4xCO2" and "piControl" experiments, we estimate the steady-state temperature change at quadrupled CO 2 (T 4× ) following the methodology by Gregory et al. (2004). Second, we fit the temporal response of global surface temperature to this quadrupled CO 2 experiment using the typical formula for a two-box model: where π , τ 1 , and τ 2 are temporary parameters used for the calibration only. Third, we deduce our three dynamical parameters (i.e., τ T G , τ T D , and θ ) by using the correspondence between the temporary parameters and ours, given by Geoffroy et al. (2013, Table 1). Fourth, we deduce the climate sensitivity λ of the model by normalizing T 4× by the RF caused by a quadrupled CO 2 as quantified by the IPCC logarithmic formula given in Eq. (37).
The second set of parameters, for the pattern scaling, is calibrated on the same CMIP5 model chosen for the global temperature response. This pattern scaling can be based on the quadrupled CO 2 experiments, in which case the pattern is solely due to CO 2 -induced warming -although, depending on the CMIP5 model, part of the regional response may come from the physiological effect of CO 2 (Sellers et al., 1996). Alternatively, it can be based on the transient historical and RCP experiments (when those RCPs are available), in which case the pattern is induced by all anthropogenic and natural perturbations, and it is thus expected to be more "realistic" but without a clear distinction of the role of each forcing. The parameter ω T S is calibrated thanks to a linear fit between yearly values of global and sea surface temperatures, whereas in the case of ω T L the linear fit is made with decadal moving averages of global and local surface temperatures. The CMIP5 fits are shown in Figs. S29 to S39.

Precipitation
Changes in global yearly precipitation (P G ) -actually used as another climate change indicator and to deduce changes in local yearly precipitation -are calculated following the simple model of Allan et al. (2013) (see also Shine et al., 2015). In this model, global precipitation vary with global temperature change and with the atmospheric fraction of RF. Two parameters are thus needed: one for the first term (α G P > 0), which describes the long-term response of the hydrological cycle to global warming, and one for the second term (β G P < 0), which describes its short-term response to the local energy imbalance induced by radiatively active species: As for surface temperature, we use a pattern-scaling approach to deduce the local yearly precipitation (P L ) for each of our land regions, parameterized with regional weights (ω P L ): As for surface temperature, the pattern-scaling approach used here ignores the difference in effects that the various climate forcers may have on regional precipitation. The first set of parameters of this module, for global precipitation, can be calibrated on 25 CMIP5 global circulation models, chosen independently from the one used for the calibration of the temperature module. Using outputs from the "abrupt4xCO2" and "piControl" experiments, we calibrate the two parameters of Eq. (82) thanks to a linear fit with a constant term made between the global surface temperature and global precipitation. The constant term is assumed to correspond to the RF term, since the radiative forcing is actually constant in the quadrupled CO 2 experiment. α G P is the slope of the fit, and β G P is the y intercept, although the latter needs to be divided by the RF of a quadrupled CO 2 as per the IPCC formula of Eq. (37), and by the value of π CO 2 atm from Sect. 2.11.1.
The second set of parameters, for the pattern scaling, are also calibrated on the same CMIP5 model as the global precipitation response. The ω P L are fitted in the exact same way as the ω T L are in the previous section, but logically using the precipitation CMIP5 variable this time. These CMIP5 fits are shown in Figs. S40 to S49.

Ocean heat content
The ocean heat content (OHC) -a third climate change indicator -is simply deduced from the two-box model used for the temperature. However, we need to introduce a coefficient (π ohc ) to account for the extra energy received by the planet but that is taken up to heat the continents and the atmosphere and to melt the ice. We have the following: where we set π ohc = 0.94 (Otto et al., 2013, Supplement Sect. S1). Note also that by using RF instead of RF warm , we implicitly assume that the warming efficacies from Sect. 2.11.1 originate from non-radiative processes only, which is not fully the case for the volcanic forcing (Gregory et al., 2016).

Numerical solving
When put together, all previous equations from Eqs.
(1) to (84) form a system of ordinary differential equations of first order, for a subset of the variables of the model. These variables are the state variables -or prognostic variables -of the dynamical system described by the differential equations. They are compiled in Table 1, along with the drivers of the model. By definition, knowledge of both the drivers and the state variables, at any time step, gives knowledge of all the other variables of the system, at that time step. These other secondary variables -or diagnostic variables -are compiled in Table 2. The differential system is solved with the forward Euler method (Euler, 1768) with a time step (δt) that can be chosen before any simulation with OSCAR, although time steps greater than a quarter of a year systematically make the model diverge. This time step is usually set to δt = 1/6 year. We note that despite having a time step for solving that is less than 1 year, the model's results cannot be interpreted at a timescale shorter than the year, primarily because no seasonal process is implemented in the model.

Experimental setup
We make two series of historical simulations, with the goal of evaluating the performance of each module of OSCAR v2.2 separately and of the fully coupled model itself. The simulations are realized within a probabilistic framework: a set a drivers and parameters is drawn randomly, with equiprobability, from the pool of potential driving datasets and parameterizations that is summarized in Table 3. With the given drivers and parameters, two simulations are made: one in which the atmospheric concentrations of well-mixed greenhouse gases, the total and per component radiative forcings, and the various climate variables are prescribed to the model; and another in which nothing more than the drivers is prescribed. The first simulation is called "offline", and the second "online". The offline simulation has the interest of uncoupling the different modules of OSCAR, thus separating them from each other and allowing an easier diagnosis of any potential issue or bias in each module. The online simulation is meant to diagnose the behavior of OSCAR when it is used as a proper Earth system model, i.e., when it is driven only by the anthropogenic perturbations of the system. The Monte Carlo ensemble size is 10 000 simulations which are drawn from a pool of more than 10 44 potential combinations of parameters (see Table 3). As described in Sect. 2.3.2, the disaggregation of the terrestrial biosphere follows the nine regions of Houghton and Hackler (2001) and six biomes. The time step of solving is one-sixth of a year. For the atmospheric concentrations of well-mixed greenhouse gases, the forcing data used for the offline simulation are from the IPCC (2013, Tables AII.1.1a and AII.1.1b). For the component-based radiative forcings, the data are also from the IPCC (2013, Table AII.1.2), although we need a way to subdivide the two RFs that are kept aggregated by the IPCC: the one from non-CO 2 WMGHGs, and the one from aerosols (all effects). Regarding the former, we use the IPCC atmospheric concentrations which we combine with the data from Myhre et al. (2013b, Tables 8.A.1 and 8.SM.1) to have component-based RFs. Regarding the latter, we take the time series from Meinshausen et al. (2011) for each individual aerosol direct effect and for the indirect effect. To ensure consistency, we rescale the component-based RFs so that, first, their value in 2010 meets the value provided by Myhre et al. (2013b), and thereafter their sum meets the IPCC aggregated value every year. Finally, for the climate data used to force the offline simulation, we use the Had-CRUT4 data for global surface temperature (Morice et al., 2012), the HadISST1 for sea surface temperature (Rayner et al., 2003), and the CRU TS3.23 dataset for local temperature and precipitation . For these three datasets, we assume the preindustrial equilibrium is their respective average over the 1901-1930 period. Note also that, because the climate data are based on observation, the offline simulation will show natural variability, albeit not as a feature of OSCAR but as one of the driving data.

Results
The following sections are dedicated to discussing the results of the historical simulations for the main variables of the model. Each section refers to one of Figs. 5 to 12. In the case of the offline simulation, we show and discuss the "reconstructed" time series of those variables that are prescribed to the model. In other words, in the following, the offline atmospheric growth rate and concentration of a given WMGHG are reconstructed as the balance of the prescribed emissions and the simulated fluxes. The offline RFs are reconstructed on the basis of the reconstructed atmospheric concentrations. The climate variables, however, are reconstructed on the basis of the prescribed RFs, so that we can discuss the performance of the climate module alone, i.e., when it is not coupled to any other module.

Carbon dioxide (Fig. 5)
The median land-use change emissions simulated by the book-keeping module of OSCAR are of the same order of magnitude -though smaller than -the values reported by the global carbon project (Le  over the 1959-2010 period, be they for the online or offline simulations. The 90 % range of our simulated emissions, however, is much larger than the uncertainty range reported by Le , and its distribution is far from a regular distribution. It can be shown (see Gasser and Ciais, 2013, Appendix A) that these two results are a consequence of the biome-specific preindustrial carbon densities which are calibrated in Sect. 2.3.2 on the TRENDY models. The large differences in carbon densities are a feature of the dynamic vegetation models themselves, although it is possible that our way of processing their output data exacerbates this discrepancy. More investigation in the matter is required, for instance using observed biomass densities as constraints, especially as the non-constrained setup leads to negative emissions under some parameterizations. We also note that the  Carbon pool of the LUC-disturbed soil; in region i, biome b and age-class a.

2.3.3
Area of a biome b; in region i. offline and online land-use change emissions are almost the same, as a direct consequence of our choice of definition that makes the land-use flux only slightly sensitive to environmental changes such as atmospheric CO 2 or climate . The median land sink we simulate in the offline simulation is slightly smaller (in absolute value) than the estimate by Le , and more importantly smaller in the online simulation. The slightly smaller median value in the offline case can be explained by the weight of the four (out of 13) preindustrial land covers for which we use the crosswalking table of Poulter et al. (2011) to translate biomes into plant functional types (see Sect. 2.3.2). Using this table indeed gives a more important fraction of land covered by bare soil than is the case in most of the TRENDY models, explaining the bimodal distribution of the offline land sink. As for the online simulation, the reduced land sink is also a consequence of the warmer tropical climate simulated by OSCAR than the one prescribed with the CRU dataset in the offline simulation (see below). The inter-annual variability of the land sink simulated by OSCAR in the offline case, which is a result of the variability in the input climate data, does not match that from Le , but we do not expect our crude and aggregated approach to model the terrestrial biosphere's response to climate to be able to reproduce this variability, especially as some factors such as volcanoes do not directly influence the terrestrial carbon cycle of our model while they seem to do in reality (e.g., Raupach et al., 2014). The large spread in our estimated land sink has the same origin as that in our estimated land-use change emissions, i.e., the TRENDY models, although this time the distribution appears more regular.    γ resp,T ; γ resp,P ; γ resp,T 1 ; γ resp,T 2 ; γ resp,P    , and the dashed red lines show the 90 % uncertainty range (calculated as 1.645 times the 1σ -uncertainty). Reference 1 for the atmospheric growth rate and concentration is NOAA/ESRL . Reference 2 denotes Law Dome ice cores (Etheridge et al., 1996;MacFarling Meure et al., 2006).
The median ocean sink OSCAR simulates matches relatively well the estimate by Le , albeit it is slightly stronger (in absolute value) in the online case. This relative good performance of the ocean carbon-cycle module, given that no change in the biological pump is simulated by OSCAR, suggests that the physical pump is enough to satisfactorily simulate the (recent) past carbon uptake by the ocean, as noted by Prentice et al. (2001). Whether this would be enough to simulate future changes remains to be tested. Note that the discontinuous distribution of the ocean sink in the offline case reflects the fact that we only have 12 possible parameterizations for this module, when climate is fixed.
In both the online and offline simulations, the simulated atmospheric growth rate is very close, on average, to the one reported by NOAA/ESRL . In the online case, this happens in spite of the relatively small land sink discussed above, owing to the compensation of the reduced land sink by an enhanced ocean sink. This shows that there is a negative feedback loop occurring in the online setup. This loop occurs through the oceanic carbon cycle: when the land sink is too low, atmospheric carbon dioxide increases faster, which in turn increases the ocean sink. This kind of anti-correlation between two of the global carbon budget's fluxes is also found between the land sink and land-use change emissions: a high productivity configuration of the model simulates high emissions of land-use changebecause of high carbon density biomes -but also high terrestrial carbon sink. Consequently, the atmospheric CO 2 simulated by OSCAR under some parameterizations may appear to be correct because of biases compensation.
Finally, regarding excess atmospheric CO 2 , both median simulations follow the observations since 1959 fairly well, with a slight positive offset for the online case and a slight negative one for the offline one, of ∼5 ppm in both cases. In the online case, however, the simulated atmospheric CO 2 prior to the direct observations is very close to the estimates derived from ice cores (Etheridge et al., 1996;MacFarling Meure et al., 2006), at least until the simulation reaches the atmospheric plateau of the 1940s. Therefore, the offset we simulate over the recent period is a consequence of the model "missing" the plateau, as all complex models do (Bastos et al., 2016). The spread in the results from the two setups is high, but the spread in the offline simulation is much higher than in the online case, owing mainly to the spread in our simulated land sink. Some parameterizations in the offline setup even lead to negative atmospheric CO 2 , resulting from combined negative land-use emissions and a strong land sink. This unrealistic behavior of the model puts forward the need to use observational constraints to select only a subset of the parameterizations in future works. (Fig. 6) The emissions from biomass burning are shown and discussed here, despite being mainly a product of the carbon cycle in OSCAR, since they are part of the atmospheric balance of methane. One can see that our approach of calculating these emissions endogenously gives values of the same order of magnitude as those of Lamarque et al. (2010), albeit with a different temporal profile. This different profile of ours closely follows that of land-use change emissions in Fig. 5, which indicates that our emissions from biomass burning are mainly the product of land-use and land-cover change -or in other words that the second term of Eq. (38) dominates. In the offline simulation, however, there is a noticeable interannual variability, showing that the environmental conditions -and especially climate -also affect our biomass burning emissions; in other words the first term of Eq. (38) is not negligible. The peak at zero emission in the distribution is artificially created by the option we keep of turning off the wild-fire feature, which is drawn once out of seven in the Monte Carlo.

Methane
When compared to the multi-model mean of WETCHIMP (Melton et al., 2013), our offline predicted change in the emission of methane by natural wetlands is of the right order of magnitude, albeit without a good reproduction of the interannual variability simulated by complex models. We see this relatively good performance for the offline simulation, i.e., for an experimental protocol with OSCAR that is very close to the one used in WETCHIMP. For the online simulation, however, one can see that our simulated wetland emissions are much lower (by a factor 2) to that simulated in the offline case. This comes from the inability of OSCAR to simulate a regional climate change -and especially precipitation (see below) -close to the forcing data we use in the offline simulation, therefore affecting the wetland area extent predicted by the model.
The median lifetime of methane with regard to the OH sink which we simulate is very close to the best-guess value of Prather et al. (2012) for the present day. This is an expost justification to our arbitrary rescaling of the preindustrial lifetime τ CH 4 OH in Sect. 2.5.1. Here, we also note that our 90 % spread in methane's lifetime is greater than the corresponding uncertainty range provided by Prather et al. (2012), particularly in the online simulation. This stems from the large spread in our simulated emissions of biomass burning -which itself is a consequence of the spread in landuse change emissions -as the biomass burning emissions of NO x , CO, and VOCs impact the OH sink capacity.
In the online simulation the median atmospheric growth rate of methane we simulate is close to the observed one, over the short period of observation we have at our disposal. OSCAR manages to reproduce the slowdown of atmospheric increase around the year 2000; this slowdown is mainly driven by anthropogenic emissions in our model. After 2005, however, the atmospheric growth resumption is too fast when compared to observations. In the offline simulation the picture is completely different: the atmospheric growth rate -reconstructed as the balance between the concentrationdriven sinks and the anthropogenic emissions normally driving OSCAR in online mode -is systematically higher than in the online case, by 10 to 20 MtC yr −1 . If our wetland emissions can explain 5 MtC, the rest must come from the anthropogenic emissions of methane we use for reconstructing the growth rate. The remaining 5 to 15 MtC represent between 5 % (around 2000) and 30 % (in 1900) of the anthropogenic emissions. These relatively small percentages stress how sensitive to anthropogenic emissions predicted atmospheric methane is: the annual growth rate of ∼ 10 MtC yr −1 results from the balance between source or sink fluxes of ∼ 250 MtC yr −1 , and any small error in one of the two fluxes can have marked impact on the growth rate. Just as with CO 2 and the ocean sink (see Sect. 3.2.1), in the online configuration there is an obvious negative feedback loop that reduces  [w.r.t. 1901-1930] Change in emissions from natural wetlands: ∆E wet . Results of our simulations with OSCAR, for methane, with the same format as for carbon dioxide. References are ACCMIP  for biomass burning, WETCHIMP (Melton et al., 2013) for wetlands, Prather et al. (2012) for the lifetime and 90 % uncertainty range (calculated as 1.645 times the 1σ uncertainty), and AGAGE (Prinn et al., 2013) for the atmospheric growth rate and concentration. Reference 2 denotes Law Dome ice cores (Etheridge et al., 1998;MacFarling Meure et al., 2006, using the NOAA04 scale).
the importance of this: the sink is directly proportional to the atmospheric concentration; but this feedback loop is cut off in the offline configuration, leading to a much larger -and unrealistic -spread in the results, and again putting forward the need to use observational constraints.
Regarding atmospheric CH 4 , in the online configuration we simulate a concentration that is close to recent observations, albeit slightly lower. The distance between the median of our ensemble and AGAGE (Prinn et al., 2013) is ∼ 40 ppb over 1987-2005 and then decreases to be virtually zero in 2010. Before that, however, when compared to ice-core data (Etheridge et al., 1998;MacFarling Meure et al., 2006) the simulated atmospheric CH 4 is systematically higher by ∼ 100 ppb. With the offline configuration, as a direct consequence of the systematic overestimate of the reconstructed atmospheric growth rate, the reconstructed atmospheric concentration we simulate is completely off-track. This could be solved by using our own estimates of compatible methane emissions (see, e.g., Gasser et al., 2015) which would be 5 to 30 % lower than those used here (and described in Sect. 2.2.1), as explained above; but also by using con-straints to exclude unrealistic realizations of the Monte Carlo ensemble. (Fig. 7) The nitrous oxide emissions from biomass burning are shown here mainly to point out that they are strictly similar to that of methane in Fig. 6. This is true for all non-CO 2 species in OSCAR: given our modeling approach, their biomass burning emissions are roughly proportional by a factor equal to the ratio of two α X bb (see Sect. 2.4.1). Therefore, the comments on biomass burning made in Sect. 3.2.2 are valid for all non-CO 2 species.

Nitrous oxide
The median lifetime of nitrous oxide with regard to the stratospheric sink which we simulate is very close to the best-guess value of Prather et al. (2015) for the present day. Our simulated spread also reasonably covers the uncertainty range they give, although it does not reach the lowest values of the range. The distribution of the lifetime we simulate, however, is asymmetrical and somewhat discontinuous. Both features are direct consequences of the distribution in the model's estimates of the lifetime which we base our parameter on; however the latter one also indicates that we do not have enough available parameterizations to produce a proper uncertainty range.
On average, the median atmospheric growth rate we simulate is close to the observed one over 1979-2010, although slightly smaller for the offline simulation. The observed variability, however, is not reproduced by our model, be it in the online or offline setup. This suggests that a biological process related to nitrous oxide is missing in our model. Processes such as biological production in terrestrial or aquatic systems are viable candidates (Ciais et al., 2013b).
In the online simulation, the excess atmospheric concentration we simulate is lower than the one observed: the median is actually parallel to the observations with a distance of ∼ 4 ppb. This feature indicates that the growth rate simulated over the recent period is good -as we explained aboveand thus that the difference between simulation and observation originates from the earlier period. This is confirmed by the comparison with ice-core data (MacFarling Meure et al., 2006). Assuming that our estimate of the nitrous oxide sink is right, the difference could be explained by any phenomenon that would imply higher emissions in the past than we use as input here, be they anthropogenic or of natural origin. As for the offline configuration, the simulated atmospheric N 2 O is even lower, owing to the lower growth rate mentioned above, and its spread is larger because of the same reasons as for atmospheric CH 4 . (Fig. 8) While other species are shown in Fig. S50, here we show only the first compound of each group of halogenated compounds (i.e., HFC-23 for HFCs, CF 4 for PFCs, and CFC-11 for ODSs) to illustrate two points. First, OSCAR is able to reproduce relatively well the past evolution of the atmospheric concentration of these compounds, although not with very good performance in all cases. Second, the fact that we only have one set of preindustrial lifetimes and one dataset of anthropogenic emissions hampers our ability to produce a proper distribution of results with OSCAR. Hence, if any of the data is wrong, the simulation with OSCAR will also be wrong with respect to those species. Alternative parameters and/or input data should be used in future versions of the model, or -more importantly -in any future study that would focus on those compounds.

Halogenated compounds
If we look at the variables that summarize the two effects of the halogenated compounds within the climate system, i.e., effective equivalent stratospheric chlorine and radiative forcing, we can have an overview of the performance of this module. Regarding the EESC simulated by our model, it is lower than the one calculated on the basis of the IPCC (2013) atmospheric concentrations and the fractional release parameters from Newman et al. (2007) used by the WMO (Montzka et al., 2011). Note, however, that in OSCAR those fractional release factors can also take alternative values, as illustrated by the three lines in the distribution of the offline EESC. Regarding the combined radiative forcing of all halogenated compounds, the offline simulation gives a slightly higher value than the IPCC's (Myhre et al., 2013b), whereas the online simulation gives a slightly lower one. In both cases, the values remain within the 90 % uncertainty range assessed by the IPCC. Consequently, despite a relatively bad performance of the module for some individual halogenated species, we consider that the module performs well as a whole.
3.2.5 Ozone (Fig. 9) Regarding tropospheric ozone, the median change in burden simulated by OSCAR is very close to the only point in time we have from the IPCC (2013, Table AII.5.2) which is for the change in burden over 1850-2000. The corresponding RF, however, is higher in our simulation than the one provided by the IPCC (Myhre et al., 2013b) for the year 2010. Given that OSCAR seems to perform well over 1850-2000, the cause of the discrepancy between the IPCC RF estimate and ours can be a different estimate of change in burden before or after that period and/or a different radiative efficiency of tropospheric ozone. In any case, our estimate remains within the IPCC uncertainty range, but it must be noted that our 90 % range is almost systematically higher than the IPCC best guess, and that consequently there is probably a positive bias in our model for tropospheric ozone.
Regarding stratospheric ozone, our slightly underestimated EESC induces a slightly underestimated change in column burden (in absolute value), again over the reference period 1850-2000. Nonetheless, the estimate by the IPCC (2013) is well within our 90 % range -a range that is dis- 1 9 0 0 1 9 2 0 1 9 4 0 1 9 6 0 1 9 . Results of our simulations with OSCAR, for nitrous oxide, with the same format as for carbon dioxide. References are Prather et al. (2015) for the lifetime and 90 % uncertainty range (calculated as 1.645 times the 1σ uncertainty) and AGAGE (Prinn et al., 2013) for the atmospheric growth rate and concentration. Reference 2 denotes Law Dome ice cores (MacFarling Meure et al., 2006). continuous in the offline configuration, as could be expected from the discontinuity of the EESC seen in Fig. 8. The corresponding median RF we estimate is close to the IPCC best guess and its spread is also close to the uncertainty range provided by the IPCC, except that it does not go as far into the positive-value domain as the IPCC does. (Fig. 10) Regarding the direct effect of aerosols, OSCAR's ability to match the IPCC best guess (Myhre et al., 2013b) in 2010 varies with the aerosol considered. In the case of sulfates, the median RF we simulate is slightly smaller (in absolute value) than the IPCC reference, while the spread is larger than the reference and has a non-regular distribution. The non-regular distribution, as well as the model going too far into the negative values, can be explained by having only four models at our disposal to calibrate the atmospheric lifetime, among which one leads to a calibrated lifetime that is almost twice that of the others.

Aerosols
The cases of POA and BC are very comparable: our median RFs are significantly smaller (in absolute value) than the IPCC references, and the distributions are close to a lognormal one and with a relatively consistent spread. For both aerosols, however, if we remove the contribution of biomass burning aerosols to the IPCC best guesses, our median estimates are much closer. This odd feature does not greatly affect the overall performance of the model (see next section), as the IPCC best-guess estimate for combined biomass burning POA and BC is zero. It strongly suggests, however, that the way these biomass burning aerosols are treated in OSCAR can be improved.
In the case of nitrate, our median RF is relatively close to the IPCC best guess, whereas our distribution does not go as far in the negative values as the IPCC uncertainty range, as a result of having too few -only two -possible parameterizations for these aerosols. In the case of SOA, our median RF is very small, owing to the fact that one out of three simulations has the SOA turned off, and the distribution clearly shows that we only have three possible parameterizations for this aerosol. Also, because all the radiative efficiencies of SOA available to OSCAR are negative, the only way it could go into the positive-value domain would be to have varying biogenic emissions of NMVOCs, which is not the case in this version.
Regarding the cloud effect of aerosols, which includes both the so-called semi-direct and indirect effects, OSCAR performs well and its median estimate meets the IPCC best guess in 2010. This is mostly due to the way this effect is calculated in our model, as the main sensitivity parameter of the module (i.e., ) is rescaled using the IPCC estimate. Nonetheless, this and the shape of the distribution, which is close to a log-normal one, show that our simple formulation of the cloud effect is consistent. Note also that the online and offline simulations are very close, both for the direct and cloud effects, because of the limited role of climate in our aerosol module.  3.2.7 Radiative forcing ( Fig. 11) When we combine together the RF induced by all well-mixed greenhouse gases, we see that the median of both our online and offline simulations are slightly higher in 2010 than the estimate by Myhre et al. (2013b), albeit with a larger spread than the reference in the online case, and a much larger spread in the offline one. The latter feature is a direct consequence of the large spread in the offline simulations of atmospheric CO 2 and CH 4 discussed above. When the RF induced by all near-term climate forcers is combined, we see similarly that the median of both our online and offline simulations is close to the IPCC estimate for 2010. This time, however, our simulated spread is relatively consistent with the IPCC uncertainty range.
Regarding the two RFs induced by surface albedo, our two simple modules simulate values that meet the IPCC estimate for the year 2010. For black carbon deposition on snow, this could be expected from our rescaling of the global sensitivity parameter α BCsnow rf , although the spread in our results is smaller than the IPCC uncertainty. For land-cover change, however, no parameter was rescaled to meet the IPCC best guess, and the distribution of our simulated RF shows that this median result is actually the product of several parameterizations with very contrasted results. The contrast can be explained by how we process the LULCC dataset that is used as input of this module: we use the TRENDY models, or land-cover climatologies, to know the proportion of each natural biome within a given grid cell of the dataset. Therefore, the huge difference in preindustrial land covers among the TRENDY models or the climatologies leads to this large spread; and the climatologies in which there seem to be too much bare soil (see Sect. 3.2.1) lead to the very negative values of this RF. We also note that the offline and online simulations of this RF from land-cover change are strictly equal because the module is driven only by LULCC drivers, and it is therefore not coupled to any other module.
All in all, the total RF simulated by OSCAR -which is the sum of the above four RFs and the three drivers prescribed directly as radiative forcing -has a median value in the year . Results of our simulations with OSCAR, for aerosols, with the same format as for carbon dioxide. Reference 1 for the radiative forcing and its 90 % uncertainty range is IPCC (Myhre et al., 2013b). Reference 2 is the same except that contribution from biomass burning aerosols is removed.
2010 close to the IPCC best guess, but slightly higher. In the online case it has a relatively consistent spread, whereas in the offline one the spread is much larger. This large spread is dominated by the large spread in the RF of WMGHGs which itself is dominated by the large spread in offline atmospheric CO 2 (Sect. 3.2.1) and CH 4 (Sect. 3.2.2). . Results of our simulations with OSCAR, for radiative forcing, with the same format as for carbon dioxide. Reference for the radiative forcing and its 90 % uncertainty range is IPCC (Myhre et al., 2013b).
3.2.8 Climate (Fig. 12) Global mean surface temperature, which is our prime proxy of climate change, is relatively well simulated by OSCAR over the 1900-2010 period. We note, however, that the 1940s warmer period is not reproduced, and during the last 10 years of simulation the simulated temperature tends to be higher than the observations. Interestingly, OSCAR simulates a slowdown of the warming during these last 10 years -the so-called hiatus period. The fact that the slowdown is simulated in both the offline and online setups suggests it is a feature of our climate module alone, and thus can be explained by the RFs we use as inputs. However, the lack of inter-annual variability in OSCAR makes any further investigation on the topic virtually impossible. Note also that the offline simulation gives a narrower range than the online one because only one set of radiative forcings is prescribed in the former case. As for the global sea surface, one can see here the limits of our pattern-scaling approach: the single proportionality parameter makes the time series of sea surface temperature homothetic to that of global surface temperature. If the simulated temperature follows relatively well the observations over 1900-2010, the simulated temporal variability does not match the observed one. Similarly, the simulated local surface temperatures, shown in Fig. S51, are proportional to the   (Rayner et al., 2003) for sea surface temperature. Reference 2 is NOAA/NCDC (Smith et al., 2008) for global surface temperature, ERSST4 (Huang et al., 2015) for sea surface temperature, and NOAA/NODC (Levitus et al., 2012) for ocean heat content. Reference 3 is GISTEMP (Hansen et al., 2010) for global surface temperature.
global one, which gives temperature changes consistent with the CRU dataset  in most regions, with the notable exception of tropical regions. This suggests regional processes should be accounted for, especially as some anthropogenic activities, such as emission of short-lived species and land-use change, can have important regional impacts. This is discussed in the conclusion. Although we cannot compare our global yearly precipitation with a long enough time series of observation, we can note that OSCAR simulates a wide range of precipitation changes, with a non-negligible difference between the offline and online configurations. This is mostly caused by the difference between the simulated RF of aerosols in the online setup and the prescribed RF in the offline one. Regarding local yearly precipitation, shown in Fig. S52, OSCAR does not manage to capture the past variation of this variable, in any of our regions. This has limited impact on the model's results, since in Sect. 2.3.2 we calibrate the sensitivity parameters of NPP and heterotrophic respiration in two steps, the first of which being driven by temperature alone. It does, however, impact our simulated methane emissions from wetlands (see above). More work is needed to improve that aspect of the model.
Finally, the ocean heat content simulated with our model is of the right order of magnitude, logically owing to the good simulated RF and temperature. It follows relatively well the variations of the observations for both online and offline simulations, except over the last 10 years of simulation. This could be explained by our choice of a single value for π ohc , while this parameter should ideally be calibrated on each of the CMIP5 climate models we emulate. Alternatively, another explanation could be that our reference from NOAA/NODC (Levitus et al., 2012) actually estimates the ocean heat content down to a 2000 m depth, potentially creating a slight bias in our comparison. Note, however, that this variable is of limited importance in the model, as it is used for diagnostic purposes only.

Conclusions
In this paper, we have provided a complete description of the compact Earth system model OSCAR v2.2, and we have presented the model's results in the case of an historical simulation. Overall, despite some caveats discussed in the previous section, we conclude that the model performance is good, especially given its level of complexity. OSCAR manages to satisfactorily reproduce most of the past changes in the global Earth system, with an even better performance over the recent period for which better driving data are available. However, we note that a good performance of a simple model over the historical period does not warrant a good performance in any other simulation. In the case of OSCAR, since its parameters are generally calibrated on simulations that go relatively far from the historical conditions (e.g., under quadrupled CO 2 , or following the RCPs), we expect the model to provide reliable results over the plausible range of future climate change, in other words to cover all scenarios by Clarke et al. (2014). OSCAR's domain of validity is not as broad as that of complex models, however, and we would not recommend using the model to, e.g., perform paleoclimate studies. Ultimately, OSCAR's domain of validity should be investigated in future studies.
The fact that OSCAR has been developed to be used in a probabilistic setup is an additional strength of the model, although the spread in the model's results for some components may greatly differ from the uncertainty range assessed by studies based on more complex models and/or observations. In addition to the reasons discussed in the section above, there are two more general causes to that feature, owing to the principles underpinning OSCAR's development (expounded in Sect. 2.1). First, because all the modules of OSCAR interact with each other, the model's overall causal chain is fairly complex (as illustrated in Fig. 1) and it has many degrees of freedom -actually more than most CMIP5 complex Earth system models. These many degrees of freedom increase the odds of seeing a given simulation depart unreasonably from the plausible range of results. Second, OSCAR is not designed to emulate a given complex Earth system model as a whole: each of its modules is essentially an emulator, and OSCAR is the combination of these emulators. Consequently, in a given parameterization, two modules could emulate the sensitivities of two complex models that are physically inconsistent with one another (e.g., the implicit ocean transport of the climate module could be inconsistent with that of the carbon-cycle module), therefore potentially leading to unreasonable results. These two elements explain why OSCAR's average or median simulation can differ from the average of a model intercomparison exercise we used for calibration, and why the model's results can show very large spreads. A way to solve this and improve the probabilistic setup is to use observational constraints, either to rate a given parameterization and therefore give it a lower weight if its too far from the observations (e.g., Steinacher et al., 2013), or more abruptly to remove from the pool of the Monte Carlo experiment the parameterizations that lead to unrealistic results (e.g., Gasser, 2014). In any case, the observational constraints must be relevant to the study: for global climate change projections, atmospheric concentrations of greenhouse gases and global surface temperature could suffice, while for a study focusing, e.g., on land carbon cycle, additional constraints on NPP and carbon densities might be required.
To conclude, we want to suggest a few tracks for future development of the model. Despite its overall good performance, the model can indeed be improved, especially in terms of consistency of modeling. We see three broad aspects of the model for which such improvements would be advisable. First, the carbon cycle can be improved by inclusion of nutrient limitations for the land carbon cycle, and of the biological pump for the ocean carbon cycle. Inclusion of the nitrogen cycle would couple the carbon cycle and the atmospheric chemistry, as the carbon sinks would be affected by deposition of active nitrogen that would be induced by NO x and NH 3 emissions (e.g., Ciais et al., 2013b). This would also allow endogenous computation of some of the biogenic emissions of N 2 O, NO x or NH 3 , which would probably change our estimated past evolution of atmospheric nitrous oxide, e.g., by giving it more annual variability in the offline simulation. Second, the whole of OSCAR's atmospheric chemistry can be improved by making it consistent. In OSCAR v2.2, the atmospheric chemistry is a patchwork of many sensitivity studies. When we choose the parameters for, e.g., the stratospheric N 2 O sink, it should actually be coupled to the ozone stratospheric chemistry (e.g., Prather, 1998). Also, coupling of the tropospheric and stratospheric chemistries would be an improvement, especially for ozone, as would be a finer regionalization (or just a regionalization for the OH sink). We note however that a tremendous amount of factorial simulations by complex chemistrytransport models would be needed to make such an improvement. Third, the climate module can be improved, especially as it performs poorly at the regional scale. This would need to be done, however, by accounting for regional processes that affect temperature or precipitation, such as the physiological effect of CO 2 (e.g., Sellers et al., 1996) and thus the biophysical effect of land-use and land-cover change (e.g., Feddema et al., 2005), or the local effects of atmospheric pollution (e.g., Ramanathan et al., 2001). Implementing this in OSCAR would also require an important amount of factorial simulations, so as to be able to apply forcing-dependent patterns of climate change, or alternatively a complete rewriting of the climate module to explicitly model the local energy imbalance and water cycle. In addition to these three huge undertakings, we acknowledge that many smaller improve-ments could be made. But ultimately, the future development of OSCAR will depend on the data from complex models that will be made available.

Code availability
The source code of this version of OSCAR is available upon request to the corresponding author. Detailed information as to the complex model data processing can also be provided upon request. A brief user manual is provided with the code.

Data availability
All the data presented in this paper can be obtained upon request to the corresponding author.
Appendix A: Change log A1 OSCAR v2.1 Version 2.1 of OSCAR is completely described by Gasser (2014), although in French. Partial descriptions can be found in other studies (Cherubini et al., 2014;Li et al., 2016).
The main changes between v2.1 and v2.2 are the following: development of the ocean carbon-cycle module to include the stratification effect calibrated on CMIP5 models, extension of the terrestrial carbon-cycle module to be calibrated on many TRENDY and CMIP5 models, creation of a wildfire module, extension of the wetland module to be calibrated on many WETCHIMP models, development of the stratospheric sink module to include the effect of ozonedepleting substances and age-of-air change, development of the tropospheric ozone module to include a regionalization and the effect of climate change, development of the stratospheric ozone module to include the effect of nitrous oxide and climate change, development of the aerosols module to have explicit and regionalized parameterizations, creation of the surface albedo modules, and development of the climate module to include a global precipitation response.
Many other small and specific changes were also made during the development of the latest version.

A2 OSCAR v2.0
Version 2.0 of OSCAR is exactly the same as version 2.1, with two significant exceptions. First, non-CO 2 species were not modeled at all, which means that v2.0 was a carbonclimate model. Second, only one climate response was available, that developed by Hooss et al. (2001), instead of the CMIP5 responses available now. Version 2.0 was used by Gasser and Ciais (2013) and very briefly described therein.
It can also be noted that the main change between the previous versions of OSCAR and v2.0 is the computing language used to code the model. While previous versions were coded in Scilab, the following versions (i.e., from v2.0 onward) have been coded in Python.

A3 OSCAR v1.1
Version 1.1 of OSCAR is an update of version 1.0, described by Gitz (2004). The update is limited to the inclusion of a basic climate response and of a simple climate-carbon feedback, for the terrestrial carbon cycle only.

A4 OSCAR v1.0
Version 1.0 of OSCAR is described by Gitz and Ciais (2003). At that time, it was a simple carbon-cycle model designed to specifically focus on land-use change issues, as the bookkeeping module was already included in the model (albeit not exactly coded in the way it is now).

Appendix B: Complex models used for calibration
These models are those whose outputs we use to calibrate some of OSCAR's parameters. In other words, we do not list here the models for which we simply read OSCAR's parameter value in, e.g., a table of another study. Note that here we give the models' name as given by the study we base our calibration on. These names may vary across studies and from the official name itself.

B1 CMIP5
For the ocean carbon cycle, stratification effect (Sect.  For the atmospheric burden of sulfate, primary organic and black carbon aerosols (Sect. 2.9.1): GISS-E2-R.