A reconstruction of radiocarbon production and total solar irradiance from the Holocene 14 C and CO 2 records: implications of data and model uncertainties

Abstract. Radiocarbon production, solar activity, total solar irradiance (TSI) and solar-induced climate change are reconstructed for the Holocene (10 to 0 kyr BP), and TSI is predicted for the next centuries. The IntCal09/SHCal04 radiocarbon and ice core CO2 records, reconstructions of the geomagnetic dipole, and instrumental data of solar activity are applied in the Bern3D-LPJ, a fully featured Earth system model of intermediate complexity including a 3-D dynamic ocean, ocean sediments, and a dynamic vegetation model, and in formulations linking radiocarbon production, the solar modulation potential, and TSI. Uncertainties are assessed using Monte Carlo simulations and bounding scenarios. Transient climate simulations span the past 21 thousand years, thereby considering the time lags and uncertainties associated with the last glacial termination. Our carbon-cycle-based modern estimate of radiocarbon production of 1.7 atoms cm−2 s−1 is lower than previously reported for the cosmogenic nuclide production model by Masarik and Beer (2009) and is more in-line with Kovaltsov et al. (2012). In contrast to earlier studies, periods of high solar activity were quite common not only in recent millennia, but throughout the Holocene. Notable deviations compared to earlier reconstructions are also found on decadal to centennial timescales. We show that earlier Holocene reconstructions, not accounting for the interhemispheric gradients in radiocarbon, are biased low. Solar activity is during 28% of the time higher than the modern average (650 MeV), but the absolute values remain weakly constrained due to uncertainties in the normalisation of the solar modulation to instrumental data. A recently published solar activity–TSI relationship yields small changes in Holocene TSI of the order of 1 W m−2 with a Maunder Minimum irradiance reduction of 0.85 ± 0.16 W m−2. Related solar-induced variations in global mean surface air temperature are simulated to be within 0.1 K. Autoregressive modelling suggests a declining trend of solar activity in the 21st century towards average Holocene conditions.


Introduction
Solar insolation is the driver of the climate system of the earth (e.g. Gray et al., 2010;Lockwood, 2012). Variations in total solar irradiance (TSI) have the potential to significantly modify the energy balance of the earth (Crowley, 2000;Ammann et al., 2007;Jungclaus et al., 2010). However, the magnitude of variations in TSI (Schmidt et al., 2011(Schmidt et al., , 2012Shapiro et al., 2011;Lockwood, 2012) and its temporal evolution (Solanki et al., 2004;Muscheler et al., 2005b;Schmidt et al., 2011) remain uncertain and are debated. Solar activity and TSI were reconstructed from the Holocene radiocarbon record. These reconstructions relied on box models of the ocean and land carbon cycle or on a 2-D representation of 11 kyr BP), (ii) variations in Holocene and last millennium climate, (iii) changes in ocean sediments, and (iv) changes in vegetation dynamics and in anthropogenic land use.
The goal of this study is to reconstruct Holocene radiocarbon production, solar modulation potential to characterise the open solar magnetic field, TSI, and the influence of TSI changes on Holocene climate from the proxy records of atmospheric 14 C Reimer et al., 2009) and CO 2 and a recent reconstruction of the geomagnetic field (Korte et al., 2011). The TSI reconstruction is extended into the future to year 2300 based on its spectral properties. We apply the Bern3D-LPJ Earth system model of intermediate complexity that features a 3-D dynamic ocean, reactive ocean sediments, a dynamic global vegetation model, an energy-moisture balance atmosphere, and cycling of carbon and carbon isotopes. The model is forced by changes in orbital parameters, explosive volcanic eruptions, well-mixed greenhouse gases (CO 2 , CH 4 , N 2 O), aerosols, ice cover and land-use area changes. Bounding scenarios for deglacial radiocarbon changes and Monte Carlo techniques are applied to comprehensively quantify uncertainties.
TSI reconstructions that extend beyond the satellite record must rely on proxy information. 14 C and 10 Be are two proxies that are particularly well suited (Beer et al., 1983;Muscheler et al., 2008;Steinhilber et al., 2012); their production by cosmic particles is directly modulated by the strength of the solar magnetic field and they are conserved in ice cores ( 10 Be) and tree rings ( 14 C). The redistributions of 10 Be and 14 C within the climate system follow very different pathways, and thus the two isotopes can provide independent information. The 10 Be and 14 C proxy records yield in general consistent reconstructions of isotope production and solar activity with correlations exceeding 0.8 (Bard et al., 1997;Lockwood and Owens, 2011). Important caveats, however, and differences in detail remain. After production, 10 Be is attached to aerosols, and variations in atmospheric transport and dry and wet deposition of 10 Be and incorporation into the ice archive lead to noise and uncertainties. This is highlighted by the opposite, and not yet understood, 20th century trends in 10 Be in Greenland versus Antarctic ice cores (Muscheler et al., 2007;Steinhilber et al., 2012). Thus, taken at face value, Greenland and Antarctic 10 Be records suggest opposite trends in solar activity in the second half of the 20th century. 14 C is oxidised after production and becomes as 14 CO 2 part of the global carbon cycle. It enters the ocean, sediments, and the land biosphere and is removed from the climate system by radioactive decay, with an average lifetime of 5730 yr/ln(2) = 8267 yr, and to a minor part by seafloor sediment burial. Atmospheric 14 C is influenced not only by shortterm production variations, but also by variations of the coupled carbon-cycle-climate system, and their evolutions due to the long timescales governing radioactive decay and carbon overturning in the land and ocean.
The conversion of the radiocarbon proxy record Reimer et al., 2009) to TSI involves several steps. First a carbon cycle model is applied to infer radiocarbon production by deconvolving the atmospheric radiocarbon budget. Radiocarbon production is equal to the prescribed changes in the atmospheric radiocarbon inventory and decay in the atmosphere plus the modelled net air-to-sea and net air-to-land 14 C fluxes. The radiocarbon signature of a flux or a reservoir is commonly reported in the 14 C notation, i.e. as the fractionation-corrected per mil deviation of 14 R = 14 C/ 12 C from a given standard defined as 14 R std = 1.176 × 10 −12 (Stuiver and Polach, 1977). 14 C production depends on the magnitude of the shielding of the earth's atmosphere by the geomagnetic and the open solar magnetic fields. Reconstructions of the geomagnetic field can thus be combined with a physical cosmogenic isotope production model Beer, 1999, 2009;Kovaltsov et al., 2012) to infer the strength of the magnetic field enclosed in solar winds as expressed by the so-called solar modulation potential, . Finally, is translated into variations in TSI.
Each of these steps has distinct uncertainties and challenges. Deglacial carbon cycle changes are large, and atmospheric CO 2 increased from 180 ppm at the Last Glacial Maximum (LGM) to 265 ppm 11 kyr BP. These variations are thought to be driven by physical and biogeochemical reorganisations of the ocean (e.g. Brovkin et al., 2012;. They may influence the 14 C evolution in the Holocene and the deconvolution for the 14 C production, but were not considered in previous studies. Many earlier studies (Marchal, 2005;Muscheler et al., 2005a;Usoskin and Kromer, 2005;Vonmoos et al., 2006;Steinhilber et al., 2012) deconvolving the 14 C record relied on simplified box models using a perturbation approach (Oeschger et al., 1975;Siegenthaler, 1983) where the natural marine carbon cycle is not simulated and climate and ocean circulations as well as land carbon turnover are kept constant. In the perturbation approach, the ocean carbon and radiocarbon inventory is underestimated by design as the concentration of dissolved inorganic carbon is set to its surface concentration. 12 C and 14 C are not transported as separate tracers but combined into a single tracer, the 14 C/ 12 C ratio. These shortcomings, however, can be overcome by applying a spatially resolved coupled carbon-cycle-climate model that includes the natural carbon cycle and its anthropogenic perturbation, and where 12 C and 14 C are distinguished.
A key target dataset is the data of the Global Ocean Data Analysis Project (GLODAP) that includes station data and gridded data of dissolved inorganic carbon and its 14 C signature (Key et al., 2004). This permits the quantification of the oceanic radiocarbon inventory, by far the largest radiocarbon inventory on earth. Together with data-based estimates of the carbon inventory on land, in the atmosphere, and in reactive sediment and their signatures, the global radiocarbon inventory and thus the long-term average 14 C production can be reliably estimated. Further, the spatial carbon and 14 C distribution within the oceans is a yard stick to gauge the performance of any ocean circulation and carbon cycle model. 14 C permits the quantification of the overturning timescales within the ocean (e.g. Müller et al., 2006) and of the magnitude of the air-sea carbon exchange rate (e.g. Naegler and Levin, 2006;Sweeney et al., 2007;Müller et al., 2008).
The conversion of the radiocarbon production record into requires knowledge on the strength of the geomagnetic field, which together with the open solar magnetic field contributes to the shielding of the earth's atmosphere from the cosmic ray flux. Recently, updated reconstructions of the earth magnetic field have become available (Knudsen et al., 2008;Korte et al., 2011). The palaeo-proxy record of should be consistent with instrumental observations. The deconvolution of the atmospheric 14 C history for natural production variations is only possible up to about 1950 AD. Afterwards, the atmospheric 14 C content almost doubled due to atomic bomb tests in the fifties and early sixties of the 20th century and uncertainties in this artificial 14 C production are larger than natural production variations. This limits the overlap of the 14 C-derived palaeo-proxy record of with reconstructions of based on balloon-borne measurements and adds uncertainty to the normalisation of the palaeo-proxy record to recent data. How variations in cosmogenic isotope production and in are related to TSI is unclear and there is a lack of understanding of the mechanisms. This is also reflected in the large spread in past and recent reconstructions of TSI variability on multi-decadal to centennial timescales (e.g. Bard et al., 2000;Lean, 2000;Wang et al., 2005;Steinhilber et al., 2009;Steinhilber et al., 2012;Shapiro et al., 2011;Schrijver et al., 2011). Recent reconstructions based on the decadal-scale trend found in the TSI satellite record reveal small amplitude variations in TSI of the order of 1 W m −2 over past millennia (Steinhilber et al., 2009), whereas others based on observations of the most quiet area on the present sun suggest that TSI variations are of the order of 6 W m −2 or even more (Shapiro et al., 2011).
The outline of this study is as follows. Next, we will describe the carbon cycle model and methods applied. In the results Sect. 3, we apply sinusoidal variations in atmospheric 14 C with frequencies between 1/(5 yr) to 1/(1000 yr). The results characterise the model response to a given radiocarbon time series as any time series can be described by its power spectrum using Fourier transformation.
In Sect. 3.2, we discuss the radiocarbon and carbon inventory and distribution in the model in comparison with observations (Sect. 3.1), before turning to the time evolution of the carbon fluxes, the atmospheric carbon budget, and radiocarbon production in Sect. 3.2. In Sects. 3.3 to 3.5 results are presented for the solar activity, TSI, and simulated, solar-driven changes in global mean surface air temperature over the Holocene and extrapolated TSI variations up to year 2300. Discussion and conclusions follow in Sect. 4 and the Appendix presents error calculations for 14 C production, solar modulation, and TSI in greater detail.

Carbon cycle model description
The Bern3D-LPJ climate-carbon-cycle model is an Earth system model of intermediate complexity and includes an energy and moisture balance atmosphere and sea ice model (Ritz et al., 2011), a 3-D dynamic ocean (Müller et al., 2006), a marine biogeochemical cycle (Tschumi et al., 2008;, an ocean sediment model (Tschumi et al., 2011), and a dynamic global vegetation model   (Fig. 1). Total carbon and the stable isotope 13 C and the radioactive isotope 14 C are transported individually as tracers in the atmosphere-ocean-sediment-land biosphere system.
The geostrophic-frictional balance 3-D ocean component is based on Edwards and Marsh (2005) and as further improved by Müller et al. (2006). It includes an isopycnal diffusion scheme and Gent-McWilliams parametrisation for eddyinduced transport (Griffies, 1998). Here, a horizontal resolution of 36 × 36 grid boxes and 32 layers in the vertical is used. Wind stress is prescribed according to the monthly climatology from NCEP/NCAR (Kalnay et al., 1996). Thus, changes in ocean circulation in response to changes in wind stress under varying climate are not simulated.
The atmosphere is represented by a 2-D energy and moisture balance model with the same horizontal resolution as the ocean (Ritz et al., 2011). Following Weaver et al. (2001, outgoing long-wave radiative fluxes are parametrised after Thompson and Warren (1982) with additional radiative forcings due to CO 2 , other greenhouse gases, volcanic aerosols, and a feedback parameter, chosen to produce an equilibrium climate sensitivity of 3 • C for a nominal doubling of CO 2 . The past extent of Northern Hemisphere ice sheets is prescribed following the ICE4G model (Peltier, 1994) as described in Ritz et al. (2011).
The marine biogeochemical cycling of carbon, alkalinity, phosphate, oxygen, silica, and of the carbon isotopes is detailed by  and Tschumi et al. (2011). Remineralisation of organic matter in the water column as well as air-sea gas exchange is implemented according to the OCMIP-2 protocol Najjar et al., 1999). However, the piston velocity now scales linear (instead of a quadratic dependence) with wind speed following Krakauer et al. (2006). The global mean air-sea transfer rate is reduced by 17 % compared to OCMIP-2 to match observation-based estimates of natural and bomb-produced radiocarbon , and is in agreement with other studies (Sweeney et al., 2007;Krakauer et al., 2006;Naegler and Levin, 2006). As fractionation of 14 CO 2 is not taken into account during air-sea gas exchange, the air-sea fluxes are corrected in postprocessing using 13 C, for which 1882 R. Roth and F. Joos: A reconstruction of radiocarbon production and total solar irradiance the fractionation is explicitly calculated in the model. Prognostic formulations link marine productivity and export production of particulate and dissolved organic matter (POM, DOM) to available nutrients (P, Fe, Si), temperature, and light in the euphotic zone. Carbon is represented as tracers of dissolved inorganic carbon (DIC), DIC-13 and DIC-14, and labile dissolved organic carbon (DOC), DOC-13 and DOC-14. Particulate matter (POM and CaCO 3 ) is remineralised/dissolved in the water column applying a Martin-type power-law curve. A 10-layer sediment diagenesis model (Heinze et al., 1999;Gehlen et al., 2006) is coupled to the ocean floor, dynamically calculating the advection, remineralisation/redissolution and bioturbation of solid material in the top 10 cm (CaCO 3 , POM, opal and clay), as well as porewater chemistry and diffusion as described in detail in Tschumi et al. (2011). In contrast to the setup in Tschumi et al. (2011), the initial alkalinity inventory in the ocean is increased from 2350 to 2460 µmol kg −1 in order to get a realistic present-day DIC inventory (Key et al., 2004). The land biosphere model is based on the Lund-Potsdam-Jena (LPJ) dynamic global vegetation model with a resolution of 3.75 • × 2.5 • as used in Joos et al. (2001), Gerber et al. (2003, Joos et al. (2004) and described in detail in Sitch et al. (2003). The fertilisation of plants by CO 2 is calculated according to the modified Farquhar scheme (Farquhar et al., 1980). A land-use conversion module has been added to take into account anthropogenic land cover change (Strassmann et al., 2008;Stocker et al., 2011). 14 C is implemented following the implementation of 13 C  and taking into account radioactive decay using a mean lifetime of 8267 yr. Fractionation is twice as large for 14 C than for 13 C. Fractionation during assimilation depends on the photosynthesis pathways (C3 versus C4) and on stomatal conductance. No fractionation is associated with the transfer of carbon between the different pools in vegetation and soils.

Experimental protocol
The transient simulations are started at 21 kyr BP well before the analysis period for the production (10 to 0 kyr BP) to account for memory effects associated with the long lifetime of radiocarbon (8267 yr) and the millennial timescales of ocean sediment interactions.
The model is initialised as follows. (i) The oceanatmosphere system is brought into preindustrial (PI) steady state (Ritz et al., 2011). (ii) The ocean's biogeochemical and sediment component is spun up over 50 kyr. During this spinup phase, the loss of tracers due to solid material seafloor burial is compensated by spatially uniform weathering fluxes

Fig. 2. (a)
Reconstructed atmospheric 14 C of CO 2 (black, left axis) and the atmospheric radiocarbon inventory (blue, right axis) over the past 21 kyr. From 21 to 11 kyr BP, radiocarbon data are from IntCal09 (Reimer et al., 2009). Afterward, we use IntCal09 data for the Northern Hemisphere and SHCal04  for the Southern Hemisphere, thereby taking into account the interhemispheric gradient of approximately 8 ‰ and hemisphere-specific anthropogenic depletion of 14 C during the industrial period. The two records are not distinguishable in this figure because of their proximity -see (c) for a zoom-in. The atmospheric radiocarbon inventory is computed from the 14 C records and the CO 2 record from Monnin et al. (2004). (b) The 1σ error assigned to the 14 C record and used as input for the Monte Carlo uncertainty analysis.
to the surface ocean. These weathering input fluxes are diagnosed at the end of the sediment spin-up and kept constant thereafter. (iii) The coupled model is forced into a LGM state by applying corresponding orbital settings, GHG radiative forcing, freshwater relocation from the ocean to the ice sheets and a LGM dust influx field. Atmospheric trace gases CO 2 (185 ppm), δ 13 C (−6.4 ‰) and 14 C (432 ‰) are prescribed. The model is then allowed to re-equilibrate for 50 kyr. Next, the model is integrated forward in time from 21 kyr BP until 1950 AD using the following natural and anthropogenic external forcings (Figs. 2, 3 and 4): atmospheric CO 2 as compiled by Joos and Spahni (2008), 13 CO 2 (Francey et al., 1999;Elsig et al., 2009;Schmitt et al., 2012), 14 C of CO 2 Reimer et al., 2009), orbital parameters (Berger, 1978), radiative forcing due to GHGs CO 2 , CH 4 and N 2 O (Joos and Spahni, 2008). Iron fertilisation is taken into account by interpolating LGM  and modern dust forcing (Luo et al., 2003) following a spline-fit to the EPICA Dome C (EDC) dust record (Lambert et al., 2008). Shallow water carbonate deposition history is taken from Vecsei and Berger (2004). The ice sheet extent (including freshwater relocation and albedo changes) during the deglaciation is scaled between LGM and modern fields of Peltier (1994) using the benthic δ 18 O stack of Lisiecki and Raymo (2005) which was low-pass filtered with a cutoff period of 10 kyr. From 850 to 1950 AD, volcanic aerosols (based on Crowley, 2000, prepared by UVic), sulfate aerosol forcing applying the method by Reader and Boer (1998) detailed by Steinacher (2011), total solar irradiance forcing from PMIP3/CMIP5 (Wang et al., 2005;Delaygue and Bard, 2011) and carbon emissions from fossil fuel and cement production are taken into account (Andres et al., 1999).
The land module is forced by a 31 yr monthly Climatic Research Unit (CRU) climatology for temperature, precipitation and cloud cover. On this CRU-baseline climatology, we superpose interpolated anomalies from snapshot simulations performed with the HadCM3 model (Singarayer and Valdes, 2010). In addition, global mean temperature deviations with respect to 850 AD are used to scale climate anomaly fields obtained from global warming simulations with the NCAR AOGCM from 850 to 1950 AD applying a linear pattern scaling approach (Joos et al., 2001). Changes in sea level and ice sheet extent influence the number and locations of grid cells available for plant growth and carbon storage; we apply interpolated land masks from the ICE5G-VM2 model (Peltier, 2004). Anthropogenic land cover change during the Holocene is prescribed following the HYDE 3.1 dataset (Klein Goldewijk, 2001;Klein Goldewijk and van Drecht, 2006). In all transient simulations, the model's atmosphere is forced with the IntCal09 (21 kyr BP to 1950 AD, Reimer et al., 2009) and SHCal04 (11 kyr BP to 1950AD, McCormac et al., 2004 records for the Northern and Southern Hemisphere, respectively (Fig. 2a). For the equatorial region (20 • N to 20 • S), we use the arithmetic mean of these two records. Since SHCal04 does not reach as far back in time as the IntCal09 record, we use IntCal09 data for both hemispheres before 11 kyr BP. Between the 5 yr spaced datapoints given by these records, cubic interpolation is applied.
The Earth system underwent a major reorganisation during the last glacial termination as evidenced by warming, ice sheet retreat, sea level rise and an increase in CO 2 and other GHGs (Shackleton, 2001;Clark et al., 2012). Memory effects associated with the long lifetime of radiocarbon (8267 yr) and the long timescales involved in ocean sediment interactions imply that processes during the last glacial termination (ca. 18 to 11 kyr BP) influence the evolution of carbon and radiocarbon during the Holocene . Although many hypotheses are discussed in the literature on the mechanism governing the deglacial CO 2 rise (e.g. Köhler et al., 2005;Brovkin et al., 2007;Tagliabue et al., 2009;Bouttes et al., 2011;, it remains unclear how individual processes have quantitatively contributed to the reconstructed changes in CO 2 and 14 CO 2 over the termination. The identified processes may be distinguished into three classes: (i) relatively well-known mechanisms such as changes in temperature, salinity, an expansion of North Atlantic Deep Water, sea ice retreat, a reduction in iron input and carbon accumulation on land as also represented in our standard model setup; (ii) an increase in deep ocean ventilation over the termination as suggested by a range of proxy data (e.g. Franois et al., 1997;Adkins et al., 2002;Hodell et al., 2003;Galbraith et al., 2007;Anderson et al., 2009;Schmitt et al., 2012;Burke and Robinson, 2012) and modelling work (e.g. Tschumi et al., 2011); (iii) a range of mechanisms associated with changes in the marine biological cycling of organic carbon, calcium carbonate, and opal in addition to those included in (i).
The relatively well-known forcings implemented in our standard setup explain only about half of the reconstructed CO 2 increase over the termination . This indicates that the model misses important processes or feedbacks concerning the cycling of carbon in the ocean. To this end, we apply two idealised scenarios for this missing mechanism, regarded as bounding cases in terms of their impacts on atmospheric 14 C. In the first scenario, termed CIRC, the atmospheric carbon budget over the termination is approximately closed by forcing changes in deep ocean ventilation. In the second, termed BIO, the carbon budget is closed by imposing changes in the biological cycling of carbon.
14 C in the atmosphere and the deep ocean is sensitive to the surface-to-deep transport of 14 C. This 14 C transport is dominated by physical transport (advection, diffusion, convection), whereas biological fluxes play a small role. Consequently, processes reducing the thermohaline circulation (THC), the surface-to-deep transport rate, and deep ocean ventilation tend to increase 14 C of atmospheric CO 2 and to decrease 14 C of DIC in the deep. Recently, a range of observational studies addressed deglacial changes in radiocarbon and deep ocean ventilation. Some authors report extremely high ventilation ages up to 5000 yr (Marchitto et al., 2007;Bryan et al., 2010;Skinner et al., 2010;Thornalley et al., 2011), while others find no evidence for such an old abyssal water mass (De Pol-Holz et al., 2010). In contrast, changes in processes related to the biologic cycle of carbon such as changes in export production or the remineralisation of organic carbon hardly affect 14 C of DIC and CO 2 , despite their potentially strong impact on atmospheric CO 2 (e.g. Tschumi et al., 2011).
Technically, these two bounding cases are realised as follows. In the experiment BIO, we imply (in addition to all other forcings) a change in the depth where exported particulate organic matter is remineralised; the exponent (α) in the power law describing the vertical POM flux profile (Martin curve) is increased during the termination from a low glacial value (Fig. 3f). A decrease in the average remineralisation depth over the termination leads to an increase in atmospheric CO 2 (Matsumoto, 2007;Kwon et al., 2009;, but does not substantially affect 14 C. α is increased from 0.8 to 1.0 during the termination in BIO, while in the other experiments α is set to 0.9. In the experiment CIRC, ocean circulation is strongly reduced at the LGM by reducing the wind stress globally by 50 % relative to modern values. The wind stress is then linearly relaxed to modern values over the termination (18 to 11 kyr BP, Fig. 3f). This leads to a transfer of old carbon from the deep ocean to the atmosphere, raising atmospheric CO 2 and lowering 14 C of CO 2 (Tschumi et al., 2011). We stress that changes in wind stress and remineralisation depth are used here as tuning knobs and not considered as realistic.
To further assess the sensitivity of the diagnosed radiocarbon production on the cycling of carbon and climate, we perform a control simulation (CTL) where all forcings except atmospheric CO 2 (and isotopes) are kept constant at PI values. This setup corresponds to earlier box model studies where the Holocene climate was assumed to be constant.
As noted, the transient evolution of atmospheric CO 2 and 14 C over the glacial termination is prescribed in all three setups (CTL, CIRC, BIO). Thus, the influence of changing conditions over the last glacial termination on Holocene 14 C dynamics is taken into account, at least to a first order, in each of the three setups.

The production rate of radiocarbon, Q
The 14 C production rate Q is diagnosed by solving the atmospheric 14 C budget equation in the model. The model calculates the net fluxes from the atmosphere to the land biosphere ( 14 F ab ) and to the ocean ( 14 F as ) under prescribed 14 CO 2 for a given carbon-cycle/climate state. Equivalently, the changes in 14 C inventory and 14 C decay of individual land and ocean carbon reservoirs are computed. Data-based estimates for the ocean and land inventory are used to match preindustrial radiocarbon inventories as closely as possible. The production rate is then given at any time t by: Here, I and F denote 14 C inventories and fluxes, and τ (8267 yr) is the mean 14 C lifetime with respect to radioactive decay. Subscripts atm, ocn, sed, and lnd refer to the atmosphere, the ocean, reactive ocean sediments, and the land biosphere, respectively. Subscript data indicates that terms are prescribed from reconstructions and subscript model that values are calculated with the model. 14 F burial is the net loss of 14 C associated with the weathering/burial carbon fluxes. I ocn,data-model (t = t 0 ) represents a (constant) correction, defined as the difference between the modelled and observation-based inventory of DI 14 C in the ocean plus an estimate of the 14 C inventory associated with refractory DOM not represented in our model. In analogue, I lnd,data-model (t = t 0 ) denotes a constant 14 C decay rate associated with terrestrial carbon pools not simulated in our model. 14 F budget is a correction associated with the carbon flux diagnosed to close remaining imbalances in the atmospheric CO 2 budget.

Observation-based versus simulated ocean radiocarbon inventory
The global ocean inorganic radiocarbon inventory is estimated using the gridded data provided by GLODAP for the preindustrial state (Key et al., 2004) and in situ density calculated from World Ocean Atlas 2009 (WOA09) temperature and salinity fields Locarnini et al., 2010). Since the entire ocean is not covered by the GLODAP data, we fill these gaps by assuming the global mean 14 C concentration in these regions (e.g. in the Arctic Ocean). The result of this exercise is 3.27 × 10 6 mol of DI 14 C. Hansell et al. (2009) estimated a global refractory DOC inventory of 624 Gt C. 14 C of DOC measurements are rare, but data in the central North Pacific (Bauer et al., 1992) suggest high radiocarbon ages of ∼ 6000 yr, corresponding to a 14 C value of ∼ −530 ‰. Taking this value as representative yields an additional 2.9 × 10 4 mol 14 C. The preindustrial 14 C inventory associated with labile DOM is estimated from our model results to be 1.3 × 10 3 mol 14 C. This yields a data-based radiocarbon inventory associated with DIC and labile and refractory DOM in the ocean of 3.30 × 10 6 mol. The corresponding preindustrial modelled ocean inventory yields 3.05 × 10 6 mol for BIO, 3.32 × 10 6 mol for CIRC and 3.10 × 10 6 mol for CTL. Thus, the correction I ocn,data-model (t = t 0 ) is less than 8 % in the case of BIO and less than 1 % for CIRC and CTL.

Observation-based versus simulated terrestrial radiocarbon inventory
As our model for the terrestrial biosphere does not include carbon stored as peatlands and permafrost soils. We estimate this pool to contain approximately 1000 Gt C of old carbon with a isotopic signature of −400 ‰ (thus roughly one halflife old). Although small compared to the uncertainty in the oceanic inventory, we include these 5.9 × 10 4 mol of 14 C in our budget as a constant correction I lnd,data-model (t = t 0 ).

Closing the atmospheric CO 2 budget
The atmospheric carbon budget is closed in the transient simulations by diagnosing an additional carbon flux F budget : where the change in the atmospheric carbon inventory (dN atm /dt) is prescribed from ice core data, E are fossil fuel carbon emissions, and F as and F ab the net carbon fluxes into the ocean and the land biosphere. The magnitude of this inferred emission indicates the discrepancy between modelled and ice core CO 2 and provides a measure of how well the model is able to simulate the reconstructed CO 2 evolution. We assign to this flux (of unknown origin) a 14 C equal the contemporary atmosphere and an associated uncertainty in 14 C of ±200 ‰. This is not critical as F budget and associated uncertainties in the 14 C budget are generally small over the Holocene for simulations CIRC and BIO (see Appendix Fig. A1e).
The production rate is either reported as mol yr −1 or alternatively atoms cm −2 s −1 . The atmospheric area is set to 5.10 × 10 14 m 2 in our model, therefore the two quantities are related as 1 atom cm −2 s −1 = 267.0 mol yr −1 .

Solar activity
Radiocarbon, as other cosmogenic radionuclides, are primarily produced at high latitudes of earth's upper atmosphere due to nuclear reactions induced by high-energy galactic cosmic rays (GCR). Far away from the solar system, this flux is to a good approximation constant in time, but the intensity reaching the earth is modulated by two mechanism: (i) the shielding effect of the geomagnetic dipole field and (ii) the modulation due to the magnetic field enclosed in the solar wind. By knowing the past history of the geomagnetic dipole moment and the production rate of radionuclides, the "strength" of solar activity can therefore be calculated.
The sun's activity is parametrised by a scalar parameter in the force field approximation, the so-called solar modulation potential (Gleeson and Axford, 1968). This parameter describes the modulation of the local interstellar spectrum (LIS) at 1 AU. A high solar activity (i.e. a high value of ) leads to a stronger magnetic shielding of GCR and thus lowers the production rate of cosmogenic radionuclides. Similarly, the production rate decreases with a higher geomagnetic shielding, expressed as the virtual axis dipole moment (VADM).
The calculation of the normalised (relative to modern) Q for a given VADM and is based on particle simulations performed by Masarik and Beer (1999). This is the standard approach to convert cosmogenic radionuclide production rates into solar activity as applied by Muscheler et al. (2007), Steinhilber et al. (2008) and Steinhilber et al. (2012), but differs from the model to reconstruct TSI recently used by Vieira et al. (2011).
The GCR flux entering the solar system is assumed to remain constant within this approach, even though Miyake et al. (2012) found recently evidence for a short-term spike in annual 14 C. The cause of this spike (Hambaryan and Neuhäuser, 2013;Usoskin et al., 2013) is currently debated in the literature.
At the time of writing, three reconstructions of the past geomagnetic field are available to us spanning the past 10 kyr (Yang et al., 2000;Knudsen et al., 2008;Korte et al., 2011), shown in Fig. 11 together with the VADM value of 8.22 × 10 22 Am −2 for the period 1840-1990 estimated by Jackson et al. (2000). The reconstructions by Yang et al. (2000) and Knudsen et al. (2008), relying both on the same database (GEOMAGIA 50), were extensively used in the past for solar activity reconstructions (Muscheler et al., 2007;Steinhilber et al., 2008;Vieira et al., 2011;Steinhilber et al., 2012). We use the most recently published reconstruction by Korte et al. (2011) for our calculations.
For conversion from into TSI, we follow the procedure outlined in Steinhilber et al. (2009);Steinhilber et al. (2010) which consists of two individual steps. First, the radial component of the interplanetary magnetic field, B r , is expressed as a function of : where v SW is the (time-dependent) solar wind speed ; B IMF,0 , 0 , v SW,0 are normalisation factors; R SE is the mean Sun-Earth distance; ω the angular solar rotation rate; and the heliographic latitude. The factor 0.56 has been introduced to adjust the field obtained from the Parker theory with observations. The exponent is set to be in the range α = 1.7 ± 0.3 as in Steinhilber et al. (2009). Second, the B r -TSI relationship derived by Fröhlich (2009) is used to calculate the total solar irradiance: This model of converting B r (here in units of nanotesla) into TSI is not physically based, but results from a fit to observations for the relatively short epoch where high-quality observational data is available. Note that the B r -TSI relationship is only valid for solar cycle minima, therefore an artificial sinusoidal solar cycle (with zero mean) has to be added to the (solar-cycle averaged) before applying Eqs. (3) and (4). In the results section, we show for simplicity solar cycle averages (i.e. without the artificial 11 yr solar cycle).
A point to stress is that the amplitude of low-frequency TSI variations is limited by Eq. (4) and small. This is a consequence of the assumption underlying Eq. (4) that recent satellite data, which show a limited decadal-scale variability in TSI, can be extrapolated to past centuries and millennia. Small long-term variations in TSI are in agreement with a range of recent reconstructions (Schmidt et al., 2011), but in conflict with Shapiro et al. (2011), who report much larger TSI variations.

Sensitivity experiments
We start the discussion by analysing the response of the Bern3D-LPJ model to regular sinusoidal changes in the atmospheric radiocarbon ratio. The theoretical background is that any time series can be translated into its power spectrum using Fourier transformation. Thus, the response of the model to perturbations with different frequencies characterises the model for a given state (climate, CO 2 , land-use area, etc.). The experimental setup for this sensitivity simulation is as follows. 14 C is varied according to a sine wave with an amplitude of 10 ‰ and distinct period. The sine wave is repeated until the model response is at equilibrium. Periods between 5 and 1000 yr are selected. Atmospheric CO 2 (278 ppm), climate and all other boundary conditions are kept fixed at preindustrial values. The natural carbon cycle acts like a smoothing filter and changes in atmospheric 14 C arising from variations in Q are attenuated and delayed (Fig. 5a). That is the relative variations in 14 C are smaller than the relative variations in Q. Here, we are interested in inverting this natural process and diagnosing Q from reconstructed variations in radiocarbon. Consequently, we analyse not the attenuation of the radiocarbon signal, but the amplification of Q for a given variation in atmospheric radiocarbon. Figure 5a shows the amplification in Q, defined as the relative change in Q divided by the relative change in the radiocarbon to carbon ratio, 14 R. For example, an amplification of 10 means that if 14 R oscillates by 1 %, then Q oscillates by 10 % around its mean value. The amplification is largest for high-frequency variations and decreases from above 100 for a period of 5 yr to around 10 for a period of 1000 yr and to 2 for a period of 10 000 yr. High-frequency variations in the 14 C reconstruction arising from uncertainties in the radiocarbon measurements may thus translate into significant uncertainties in Q. We will address this problem in the following sections by applying a Monte Carlo procedure to vary 14 C measurements within their uncertainties (see Appendix A1).
The atmospheric radiocarbon anomaly induced by variation in production is partly mitigated through radiocarbon uptake by the ocean and the land. The relative importance of land versus ocean uptake of the perturbation depends strongly on the timescale of the perturbation (Fig. 5b). For annual-to decadal-scale perturbations, the ocean and the land uptake are roughly of equal importance. This can be understood by considering that the net primary productivity on land (60 Gt C per year) is of similar magnitude as the gross air-to-sea flux of CO 2 (57 Gt C per year) into the ocean. Thus, these fluxes carry approximately the same amount of radiocarbon away from the atmosphere. On the other hand, if the perturbation in production is varying slower than the typical overturning timescales of the ocean and the land biosphere, then the radiocarbon perturbation in the ratio is distributed roughly proportional to the carbon inventory of the different reservoirs (or to the steady-state radiocarbon flux to the ocean and the land, i.e. 430 vs. 29.1 mol yr −1 ). Consequently, the ratio between the 14 C flux anomalies into ocean and land, F as : F ab (with respect to an unperturbed state), is the higher the lower frequency of the applied perturbation (Fig. 5b). Note that these results are largely independent of the magnitude in the applied 14 C variations; we run these experiments with amplitudes of ±10 and ±100 ‰. Assuming a constant carbon cycle and climate, Q can be calculated by replacing the carbon-cycle-climate model by the modelderived Fourier filter (Fig. 5a) (Usoskin and Kromer, 2005).

Preindustrial carbon and radiocarbon inventories
The loss of 14 C is driven by the radioactive decay flux in the different reservoirs. The base level of this flux is proportional to the 14 C inventory and a reasonable representation of these inventories is thus a prerequisite to estimate 14 C production rates. In the following, modelled and observation-based carbon and radiocarbon inventories before the onset of industrialisation are compared (Table 1) and briefly discussed.
The atmospheric 14 C inventory is given by the CO 2 and 14 C input data and therefore fully determined by the forcing and their uncertainty. The ocean model represents the observation-based estimate of the global ocean 14 C inventory within 1 % for the setup CIRC and within 8 % for the setup BIO. These deviations are within the uncertainty of the observational data. Nevertheless, this offset is corrected for when Q is calculated (see Eq. 1).
The model is also able to represent the observation-based spatial distribution of 14 C in the ocean (Fig. 7). Both observation and model results show the highest 14 C concentrations in the thermocline of the Atlantic Ocean and the lowest concentrations in the deep North Pacific. Deviations between modelled and reconstructed concentrations are less than 5 % and typically less than 2 %. The model shows in general too-high  radiocarbon concentrations in the upper 1000 m, while the concentration is lower than indicated by the GLODAP data at depth. Modelled loss of radiocarbon by sedimentary processes, namely burial of POM and CaCO 3 into the diagenetically consolidated zone and particle and dissolution fluxes from/to the ocean, accounts for 52 mol yr −1 (or roughly 11 % of the total 14 C sink). This model estimate may be on the high side as the ocean-to-sediment net flux in the Bern3D-LPJ model of 0.5 Gt C yr −1 is slightly higher than independent estimates in the range of 0.2 to 0.4 Gt C yr −1 .
The total simulated carbon stored in living biomass and soils is 1930 Gt C. As discussed above, this is order 1000 Gt C lower than best estimates, mainly because peat and permafrost dynamics (Yu et al., 2010;Spahni et al., 2012) are not explicitly simulated in the LPJ version applied here. The model-data discrepancy in carbon is less than 3 % of the total carbon inventory in ocean, land, and atmosphere. It translates into 5.9 × 10 5 mol of 14 C when assuming a 14 C of −400 ‰ for this old biomass. This is well within the uncertainty range of the total 14 C inventory.

Transient results for the carbon budget and deep ocean ventilation
Next, we discuss how global mean temperature, deep ocean ventilation, and the carbon budget evolved over the past 21 kyr in our two bounding simulations (CIRC and BIO). Global average surface air temperature (SAT, Fig. 6a) and sea surface temperature (SST, not shown) are simulated to increase by ∼ 0.8 • C over the Holocene. In experiment CIRC, the global energy balance over the termination is strongly influenced by the enforced change in wind stress and the simulated deglacial increase in SAT is almost twice as large in experiment CIRC than in BIO. This is a consequence of a much larger sea ice cover and a higher planetary albedo at LGM in experiment CIRC than BIO. Circulation is slow under the prescribed low glacial wind stress and less heat is transported to high latitudes and less ice is exported from the Southern Ocean to lower latitudes in CIRC than BIO. In other words, the sea-ice-albedo feedback is much larger in CIRC than BIO. Deep ocean ventilation evolves very differently in CIRC than in BIO (Fig. 6b). Here, we analyse the global average 14 C age difference of the deep ocean (i.e. waters below 2000 m depth) relative to the overlying surface ocean. The surface-to-deep age difference is recorded in ocean sediments as 14 C age offset between shells of benthic and planktonic (B-P) species. Results from the CTL experiment with time-invariant ocean ventilation show that this "proxy" is not an ideal age tracer; the B-P age difference is additionally influenced by transient atmospheric 14 C changes and varies between 600 and 1100 yr in CTL.
The wind-stress forcing applied in experiment CIRC leads to an almost complete shutdown of the THC during the LGM and a recovery to Holocene values over the termination. Simulated B-P age increases from 2000 yr at LGM to peak at 2900 yr by 16.5 kyr BP. A slight and after 12 kyr BP a more pronounced decrease follows to the late Holocene B-P age of about 1000 yr. B-P variations are much smaller in simulation BIO, as no changes in wind stress are applied; B-P age varies between 1000 and 1700 yr. 14 C, i.e. the global mean difference in 14 C between the deep ocean and the atmosphere is −500 ‰ in CIRC until Heinrich Stadial 1 (HS1) (−350 ‰ in BIO), followed by a sharp increase of approximately 150-200 ‰ and a slow relaxation to late-Holocene values (∼ −170 ‰). This behaviour is also present in recently analysed sediment cores, see e.g. Burke and Robinson (2012) and references therein. The sharp increase in 14 C following HS1 is mainly driven by the prescribed fast atmospheric drop (Broecker and Barker, 2007).
The forcings prescribed in our bounding experiments CIRC and BIO are broadly sufficient to reproduce the reconstructed deglacial CO 2 increase. This is evidenced by an analysis of the atmospheric carbon budget (Fig. 6d and f). In the control simulation (CTL) an addition of 1700 Gt C is required to close the budget. In contrast, the mismatch in the budget is close to zero for BIO and about −100 Gt C for CIRC at the end of the simulation. In other words, only small emissions from unknown origin have to be applied on average to close the budget. Both experiments need a CO 2 sink in the early Holocene as indicated by the negative missing emissions. Such a sink could have been the carbon uptake of Northern Hemisphere peatlands (Yu et al., 2010), as peatland dynamics are not included in this model version. In summary, simulation CIRC corresponds to the picture of a slowly ventilated ocean during the LGM, whereas deep ocean ventilation changes remain small in BIO and are absent in the CTL simulation. The atmospheric carbon budget is approximately closed in simulations CIRC and BIO, whereas a substantial external carbon input is required in simulation CTL. These three simulations provide thus three radically different evolutions of the carbon cycle over the past 21 kyr and will serve us to assess uncertainties in inferred radiocarbon production rates due to our incomplete understanding of the past carbon cycle.

Time evolution of the radiocarbon production rate
Total inferred radiocarbon production, Q, varies between 350 and 650 mol yr −1 during the Holocene (Fig. 8d). Variations on multi-decadal to centennial timescales are typically within 100 mol yr −1 . The differences in Q in the early Holocene between the model setups CIRC, BIO, and CTL are mainly explained by offsets in the absolute value of Q, while the timing and magnitude of multi-decadal to centennial variations are very similar for the three setups. The absolute value of Q is about 40 mol yr −1 higher in CIRC than in BIO and about 60 mol yr −1 higher in CIRC than in CTL at 10 kyr BP. This difference becomes very small in the late Holocene and results are almost identical for CIRC and BIO after 4 kyr BP. Note that the absolute (preindustrial) value of the production rate is equal in all three setups per definition (see Eq. 1).
The inferred Q is assigned according to Eq. (1) to individual contributions from a net 14 C flux from the atmosphere to the ocean ( 14 F as ), a net flux to the land ( 14 F ab ), and atmospheric loss terms (Fig. 8). Variations in these three terms contribute about equally to variations in Q on decadal to centennial timescales, whereas millennial-scale variations in Q are almost entirely attributed to changes in 14 F as . This is in agreement with the results from the Fourier analysis presented in Sect. 3.1. Holocene and preindustrial mean fluxes and their temporal variance are listed in Table 2.
The oceanic component entering the calculation of Q is threefold (Fig. 9): (i) the compensation of the DI 14 C (and a small contribution of DO 14 C) decay proportional to its inventory, (ii) changes in the inventory of 14 C itself mainly driven by F as and (iii) the export, rain and subsequent burial  Table 2. The atmospheric radiocarbon budget (positive numbers: loss of 14 C) averaged over the last 10 kyr and for preindustrial times (1750-1900 AD). The statistical uncertainty becomes negligible due to the averaging over hundreds of datapoints and is not given here. Instead, we report √ VAR (where VAR is the variance in time) to quantify the temporal variability of the corresponding fluxes. The bulk uncertainty in the time-averaged fluxes can be estimated by the uncertainty in the corresponding radiocarbon reservoirs and is of the order of 15 % (1σ ), except for the atmospheric component which is tightly constrained by data. of Ca 14 CO 3 and PO 14 C. Thus, the mentioned offset in 14 F as between CIRC and BIO at the early Holocene are the result of differences in the dynamical evolution of the whole-ocean DIC inventories and its 14 C signature. During LGM conditions, the oceanic 14 C inventory is larger for BIO than CIRC as the deep ocean is more depleted in the slowly overturning ocean of setup CIRC. Accordingly the decay of 14 C and 14 F as is higher in simulation BIO than in CIRC (Fig. 9). The simulated oceanic 14 C inventory decreases both in CIRC and BIO as the (prescribed) atmospheric 14 C decreases. However, this decrease is smaller in CIRC than in BIO as the enforced increase in the THC and in ocean ventilation in CIRC leads to an additional 14 C flux into the ocean. In addition, the strengthened ventilation leads to a peak in organic matter export and burial while the reduced remineralisation depth in BIO leads to the opposite effect (as less POM is reaching the seafloor). In total, the higher oceanic radiocarbon decay in BIO is overcompensated by the (negative) change of the oceanic inventory. Enhanced sedimentary loss of 14 C in CIRC further increases the offset, finally leading to a higher Q at 10 kyr BP of ∼ 40 mol yr −1 in CIRC than in BIO. This offset has vanished almost completely at 7 kyr BP, apart from a small contribution from sedimentary processes. In general, the influence of climate induced carbon cycle changes is modest in the Holocene. This is indicated by the very similar Q in the CTL experiment. The biggest discrepancy between results from CTL versus those from CIRC and BIO emerge during the industrial period. Q drops rapidly in CTL as the combustion of the radiocarbon-depleted fossil fuel is not explicitly included.
In a further sensitivity run, the influence of the interhemispheric 14 C gradient on Q is explored ( Fig. 8e; dashed line). In simulation INT09 the Northern Hemisphere dataset IntCal09 is applied globally and all other forcings are as in BIO. Differences in Q between CIRC/BIO and INT09 are generally smaller than 20 mol yr −1 , but grow to 50 mol yr −1 from 1900 to 1950 AD. The reason is the different slopes in the last decades of the Northern and Southern Hemisphere records. This mid-20th century difference has important consequences for the reconstruction of solar modulation as described in Sect. 3.3. This sensitivity experiment demonstrates that spatial gradients in atmospheric 14 C and in resulting radiocarbon fluxes should be taken into account, at least for the industrial period.
In conclusion, inferred Holocene values of Q and in particular decadal to centennial variation in Q are only weakly sensitive to the details of the carbon cycle evolution over the glacial termination. On the other hand spatial gradients in atmospheric 14 C and carbon emissions from fossil fuel burning should be explicitly included to estimate radiocarbon production in the industrial period. In the following, we will use the arithmetic mean of BIO and CIRC as our best estimate for Q (Fig. 10). The final record is filtered using smoothing splines (Enting, 1987) with a cutoff period of 20 yr in order to remove high-frequency noise.
Total ±1 standard deviation (1σ ) uncertainties in Q (Fig. 10, grey band) are estimated to be around 6 % at 10 kyr BP and to slowly diminish to around ±1-2 % by 1800 AD (Appendix Fig. A1f). Overall uncertainty in Q increases over the industrial period and is estimated to be ± 6 % by 1950 AD. The difference in Q between BIO and CIRC is assumed to reflect the uncertainty range due to our incomplete understanding of the deglacial CO 2 evolution. Uncertainties in the 14 C input data, the air-sea gas exchange rate and the gross primary production (GPP) of the land biosphere are taken into account using a Monte Carlo approach and based on further sensitivity simulations (see Sect. A1 for details of the error estimation).
These include maxima in Q during the well-documented solar minima of the last millennium, generally low production, pointing to high solar activity during the Roman period, as well as a broad maximum around 7.5 kyr BP. However, the production estimates of Usoskin and Kromer (2005) and Marmod09 are about 10 % lower during the entire Holocene and are in general outside our uncertainty range. If we can rely on our data-based estimates of the total radiocarbon inventory in the Earth system, then the lower average production rate in the Usoskin and Kromer (2005) and Marmod09 records suggest that the radiocarbon inventory is underestimated in their setups.
In the industrial period, the Marmod09 production rate displays a drop in Q by almost a factor of two. This reconstruction is intended to provide the surface age of the ocean for the Holocene for calibration of marine samples from that time period. Data after 1850 AD are not to be considered as the authors purposely do not include a fossil fuel correction. Similarly, Usoskin and Kromer (2005) do not provide data after 1900 AD.

A reference radiocarbon production rate
Next, we discuss the absolute value of Q in more detail. Averaged over the Holocene, our simulations yield a 14 C production of Q = 472 mol yr −1 (1.77 atoms cm −2 s −1 ), as listed in Table 2 Fig. 10. (a, b) Total 14 C production rate from this study (black line with grey 1σ shading) and the reconstructions from Usoskin and Kromer (2005) (blue) and the output of the Marmod09 box model (red). The grey shading indicates ± 1σ uncertainty as computed from statistical uncertainties in the atmospheric 14 C records and in the processes governing the deglacial CO 2 increase, air-sea gas transfer rate and global primary production on land (see Appendix and Fig. A1). (c) Shows the last 1000 yr of the production record in more detail. The two dashed lines indicate the potential offset of our best guess production rate due to the ∼15 % (1σ ) uncertainty in the data-based radiocarbon inventory. and cosmogenic radionuclide production rates Beer, 1999, 2009) estimate Q = 2.05 atoms cm −2 s −1 (they state an uncertainty of 10 %) for a solar modulation potential of = 550 MeV. As already visible from Fig. 10, Usoskin and Kromer (2005) obtained a lower Holocene mean Q of 1.506 atoms cm −2 s −1 . Recently, Kovaltsov et al. (2012) presented an alternative production model and reported an average production rate of 1.88 atoms cm −2 s −1 for the period 1750-1900 AD. This is higher than our estimate for the same period of 1.75 atoms cm −2 s −1 . We compare absolute numbers of Q for different values of (see next section) and the geomagnetic dipole moment (Fig. 13). To determine Q for any given VADM and is not without problems because the probability that the modelled evolution hits any point in the (VADM, )-space is small (Fig. 13). In addition, the value depends on the calculation and normalisation of , which introduces another source of error.
For the present-day VADM and = 550 MeV (see Sect. 3.3), our carbon cycle based estimate of Q is ∼ 1.71 atoms cm −2 s −1 and thus lower than the value reported by Masarik and Beer (2009). Note that the statistical uncertainty in our reconstruction becomes negligible in the calculation of the time-averaged Q. Systematic and structural uncertainties in the preindustrial data-based ocean radiocarbon inventory of approximately 15 % (1σ ) (Key et al., 2004;R. Key, personal communication, May 2013) dominates the uncertainty in the mean production rate, while uncertainties in the terrestrial 14 C sink are of minor relevance. Therefore we estimate the total uncertainty of the base level of our production rate to be of the order of 15 %.

Results for the solar activity reconstruction
Next, we combine our production record Q with estimates of VADM to compute the solar modulation potential with the help of the model output from Masarik and Beer (1999) which gives the slope in the Q − space for a given value of VADM (we do not use their absolute values of Q). Uncertainties in are again assessed using a Monte Carlo approach (see Appendix A2 for details on the different sources of uncertainty in and TSI).
One key problem is the normalisation of , i.e. how is aligned to observational data (Muscheler et al., 2007). The period of overlap of the Q record with ground-based measurements is very limited as uncertainties in the atmospheric injection of radiocarbon from atomic bomb tests hinders the determination of Q from the 14 C record after 1950 AD. Forbush ground-based ionisation chamber (IC) data recently reanalysed by Usoskin et al. (2011) (US11) are characterised by large uncertainties and only cover roughly one solar cycle (mid 1936(mid -1950. The considerable uncertainties both in Q and in this period make it difficult to connect reconstructions of past solar activity to the recent epoch. Still, we choose to use these monthly data to normalise our reconstruction. Converting the LIS used for US11 to the one used by Castagnoli and Lal (1980) (which we use throughout this study) according to Herbst et al. (2010) yields an average solar modulation potential of = 403 MeV during our calibration epoch 1937-1950 AD. Accordingly, we normalise our Q record such that (the Monte Carlo mean) equals 403 MeV for the present-day VADM and the period 1937 to 1950 AD. In this time period, the uncertainty in the IC data is ∼ 130 MeV. In addition, the error due to uncertainties in Q is approximately 70 MeV. This makes it difficult to draw firm conclusion on the reliability of the normalisation. Thus, the absolute magnitude of our reconstructed remains uncertain.
We linearly blend the Q based with an 11 yr running mean of the monthly data from US11 in the overlap period 1937-1950 AD and extend the reconstruction up to 2005 AD.
Smoothing splines with a cutoff period of 10 yr are applied to remove high-frequency noise, e.g. as introduced by the Monte Carlo ensemble averaging and the blending. This blending with instrumental data slightly changes the average in 1937-1950 AD from 403 to ∼ 415 MeV. varies during the Holocene between 100 and 1200 MeV on decadal to centennial timescales (Fig. 12) with a median value of approximately 565 MeV (see histogram in Fig. 13). Millennial-scale variations of during the Holocene appear small (Fig. 12). This suggests that the millennial-scale variations in the radiocarbon production Q (Fig. 10) appear to be mainly driven by variations in the magnetic field of the earth (Fig. 11). The millennial-scale modulation of Q is not completely removed when applying VADM of Korte et al. (2011). It is difficult to state whether the remaining long-term modulation, recently interpreted as a solar cycle (Xapsos and Burke, 2009), is of solar origin or rather related to uncertainties in reconstructed VADM.
Values around 670 MeV in the last 50 yr of our reconstruction (1955( -2005 indicate a high solar activity compared to the average Holocene conditions. However, such high values are not exceptional. Multiple periods with peakto-peak variations of 400-600 MeV occur throughout the last 10 000 yr, induced by so-called grand solar minima and maxima. The sun's present state can be characterised by a grand solar maximum (modern maximum) with its peak in 1985 (Lockwood, 2010).
To check to which degree our results depend on the choice of the 14 C production model, we apply also the newer model of Kovaltsov et al. (2012) to our Q time series. The result is given for the last millennium in Fig. 12c (dotted line). Obviously, the two models yield very similar result as we normalise in both cases to the same observational dataset.
Our reconstruction of is compared with those of Usoskin et al. (2007) (US07, converted to GM75 LIS), of Muscheler et al. (2007) (MEA07) and with the reconstruction of Steinhilber et al. (2008) (SEA08), who used the ice core record of 10 Be instead of radiocarbon data (Fig. 12). Overall, the agreement between the the 10 Be and the 14 C-based reconstructions points toward the quality of these proxies for solar reconstructions. In detail, differences remain. For example, the 14 C-based reconstructions show a increase in in the second half of the 17th century, whereas the 10 Be-derived record suggests a decrease in during this period.
In difference to SEA08, which is based on 10 Be ice core records, our is always positive and non-zero and thus within the physically plausible range. Solanki et al. (2004) suggests that solar activity is unusually high during recent decades compared to the values reconstructed for the entire Holocene. Our reconstruction does not point to an exceptionally high solar activity in recent decades, in agreement with the conclusions of Muscheler et al. (2007). The relatively higher modern values inferred by Solanki et al. (2004) are eventually related to their application of a Northern Hemisphere 14 C dataset (IntCal98) only. Thus, these authors neglected the influence of interhemispheric differences in 14 C. We calculated from results of our sensitivity simulation INT09, where the IntCal09  (US11, violet). The latter was used for normalisation and extension of our reconstruction over recent decades. Note that all values of are normalised to the LIS of Castagnoli and Lal (1980) according to Herbst et al. (2010). The dotted line in the lowest panel (usually not distinguishable from the solid line) shows as calculated with the alternative production model of Kovaltsov et al. (2012).
Northern Hemisphere data are applied globally (Fig. 12, dashed line). Due to the lower Q in the normalisation period from 1937 to 1950 AD (see Fig. 8, dashed line), the record during the Holocene is shifted downward by approximately 140 MeV for INT09 compared to CIRC/BIO; the same normalisation of to the Forbush data is applied. Then, the solar activity for recent decades appears unusually high compared to Holocene values in the INT09 case.

Reconstructed total solar irradiance
We apply the reconstruction of in combination with Eqs. (3) and (4) to reconstruct total solar irradiance (Fig. 14). Uncertainties in TSI are again estimated using a Monte Carlo approach, and considering uncertainties in and in the parameters of the analytical relationship (but not the ±0.4 W m −2 in the TSI-B r relationship), TSI is expressed as deviation from the solar cycle minimum in 1986, here taken to be 1365.57 W m −2 . We note that recent measures suggest a slight downward revision of the absolute value of TSI by a few per mil (Kopp and Lean, 2011); this hardly effects reconstructed TSI anomalies. The irradiance reduction during the Maunder Minimum is 0.85 ± 0.16 W m −2 (1685 AD) compared to the solar cycle 22 average value of 1365.9 W m −2 . This is a reduction in TSI of 0.62 ± 0.11 ‰.
Changes in TSI can be expressed as radiative forcing (RF), which is given by TSI × 1 4 × (1 − A), where A is the earth's mean albedo (∼ 0.3). A reduction of 0.85 W m −2 corresponds to a change in RF of about 0.15 W m −2 only. This is more than an order of magnitude smaller than the current radiative forcing due to the anthropogenic CO 2 increase of 1.8 W m −2 (5.35 W m −2 ln (390 ppm/280 ppm)). This small reduction in TSI and RF is a direct consequence of the small slope in the relationship between TSI and the interplanetary magnetic field as suggested by Fröhlich (2009). Applying the relationships between TSI and suggested by Shapiro et al. (2011) yields changes almost an order of magnitude larger in TSI and RF (right scale in Fig. 14).
We compare our newly produced TSI record with two other recently published reconstructions based on radionuclide production, Vieira et al. (2011) (VEA11) and Steinhilber et al. (2012) (SEA12). For the last millennium, the data from Delaygue and Bard (2011) (DB11) is shown for completeness. During three grand minima in the last millennium, i.e. the Wolf, Spörer and Maunder minima, SEA12 shows plateau-like values, apparently caused by a truncation of the record to positive values. Also VEA11 suggests lower values during these minima compared to our reconstructions, but in general the differences in the TSI reconstructions are small, in particular when compared to the large TSI variations suggested by Shapiro et al. (2011).
Well-known solar periodicities in TSI contain the Hallstatt (2300 yr), Eddy (1000 yr), Suess (210 yr) and Gleissberg cycle (70-100 yr) as also discussed by Wanner et al. (2008) and Lundstedt et al. (2006). These periodicities are also present in our reconstruction (Fig. 15), as well as the long-term modulation of approximately 6000 yr (Xapsos and Burke, 2009) (not shown). A wavelet (Morlet) power spectrum indicates that the power associated with the different cycles fluctuated somewhat during the Holocene.
We test the predictability of TSI given by its periodic nature to extrapolate insolation changes up to 2300 AD. The extrapolation is calculated by fitting an autoregressive model (AR(p)) of order p, applying Burg's method (Kay, 1988) for parameter optimisation (Fig. 15c). The order of the model is the number of past time steps (i.e. years) used to fit the model to the time series. We note that the spectrum of TSI is not strictly stationary (Fig. 15b), limiting the confidence in forecasting future changes in TSI. We forecast TSI in the historical period (dashed lines in Fig. 15c) to test to which degree TSI changes can be considered as an AR(p) process. Obviously, not all changes are matched in the historic period, although the present high values of TSI seems to be predictable based on past TSI changes. Common to all extrapolations is the decreasing TSI in the next 10-20 yr to a magnitude comparable to 1900 AD. Towards 2200 AD, all three extrapolations show again increasing TSI.

Changes in global mean surface air temperature from total solar irradiance variability
We translate our new TSI record as well as two earlier reconstructions (VEA11, SEA12) into past changes in SAT using the Bern3D-LPJ model. We follow the same protocol as in Sect. 2.2. Two transient simulations were performed for each TSI reconstruction: in one simulation TSI is kept constant at its mean value, in the other simulation the actual TSI reconstruction is applied. SAT was then computed as the difference between the two simulations, normalised to SAT = 0 • C in 2005 AD. As shown in Fig. 16, the reconstructions yield rather small temperature changes attributed to solar forcing with | SAT| less than 0.15 • C at any time during the Holocene. Compared to VEA12 and SEA12, our TSI reconstruction results in less negative and more positive values in SAT before 1900 AD as expected from the TSI record. No considerable TSI-induced long-term temperature trend is simulated. Applying our TSI reconstruction based on the Shapiro et al. (2011) scaling (green line and right scale in Fig. 16) translates into SAT changes ≈ 6 times larger than the other reconstructions.
We note that our energy balance model does not take into account spectral changes and changes in ultraviolet (UV) radiation, and thus related heating or cooling of the stratosphere by the absorption of UV by ozone and other agents (Gray et al., 2010).

Summary and conclusions
In the present study the Holocene evolution of cosmogenic radiocarbon ( 14 C) production is reconstructed using a state of the art Earth system model of intermediate complexity and the IntCal09/SHCal04 radiocarbon records. Then this production record is used in combination with reconstructions of the geomagnetic field to reconstruct the solar modulation potential . In further steps, total solar irradiance is reconstructed using recently proposed relationships between and TSI, and variations in TSI were translated into temperature anomalies. The uncertainty arising from uncertainties in input data, model parametrisation for the glacial termination, and model parameter values are propagated through the entire chain from production to solar modulation to TSI using Monte Carlo and sensitivity simulations as well as Gaussian error propagation.
The Bern3D-LPJ model is used to estimate radiocarbon production from the IntCal09/SHCal04 radiocarbon records. The model features a 3-D dynamic ocean component with an ocean sediment model and an interactive marine biological cycle, a 2-D atmosphere and a dynamic global vegetation model component. The model is able to reproduce the observation-based distributions of carbon and radiocarbon within the ocean as well as approximate the evolution of atmospheric CO 2 during the past 21 kyr. The explicit representation of dissolved inorganic carbon and dissolved inorganic radiocarbon and the biological cycle in a 3-D dynamic setting is a step forward compared to the often applied perturbation approach with the box diffusion model of Oeschger et al. (1975) or the outcrop model of Siegenthaler (1983) Usoskin and Kromer, 2005;Solanki et al., 2004;Muscheler et al., 2005b). These models underestimate by design the inventory of dissolved inorganic carbon in the ocean due to the neglect of the marine biological cycle. Thus, the marine radiocarbon inventory and implied cosmogenic production is also underestimated. The Bern3D-LPJ carbonclimate model is forced with reconstructed changes in orbital parameters, ice sheet extent influencing surface albedo, and natural and anthropogenic variations in greenhouse gas concentrations as well as sulfate aerosol forcings and climate variations, thereby taking into account the influence of climatic variations on the cycling of carbon and radiocarbon.
Our production rate estimates can be compared to independent estimates obtained with models that compute production of cosmogenic nuclides in the atmosphere by simulating the interaction of incoming high energy particles with atmospheric molecules Beer, 1999, 2009). Masarik and Beer (2009) suggests a production of 2.05 ± 0.2 atoms cm −2 s −1 for a solar modulation potential of 550 MeV. The radiocarbon budget analysis yields a lower production of 1.71 atoms cm −2 s −1 and for the same . The Holocene mean 14 C production depends on variations in and the geomagnetic field and is estimated to 1.77 atoms cm −2 s −1 .
An open question was by how much uncertainties in the carbon cycle processes leading to the deglacial CO 2 rise affect estimates of the radiocarbon production. To this end, two alternative, bounding scenarios for the deglacial CO 2 rise were applied in addition to a control simulations with constant climate and ocean circulation. We show that uncertainties in the processes responsible for the reconstructed CO 2 and 14 C variations over the glacial termination translate into an uncertainty of the order of 5 % in the absolute magnitude of the production in the early Holocene, but only to small uncertainties in decadal to centennial production variations. This uncertainty in millennial average production due to the memory of the system to earlier changes vanishes over the Holocene and becomes very small (< 1 %) in recent millennia. Although a detailed process understanding of the termination, e.g. the cause of the "mystery interval" (Broecker and Barker, 2007), is not required to reconstruct radiocarbon production in the Holocene, the decreasing trend in atmospheric radiocarbon over the glacial termination must be taken into account.
Small differences between simulations that include Holocene climate change and the control run with no climate change are found. A caveat is that our 2-D atmospheric energy and moisture balance model does not represent important features of climate variability such as interannual and decadal atmospheric modes. However, the available evidence and the generally good agreement between 14 C-based and independent solar proxy reconstructions suggest that atmospheric radiocarbon and inferred production are relatively insensitive to typical Holocene climate variations.
A prominent excursion in the early Holocene Northern Hemisphere climate is the 8.2 kyr BP event with a decrease in Greenland air temperature by about 3 K within a few decades (Kobashi et al., 2007). A reduction in North Atlantic Deep Water formation related to a spike in meltwater input has been proposed to have occurred with this event (Ellison et al., 2006;Hoffman et al., 2012). This 8.2 kyr event is not represented in our model setup, and a possible atmospheric radiocarbon signal from changing ocean circulation (e.g. Matsumoto and Yokoyama, 2013) would be erroneously attributed to a change in production. Vonmoos et al. (2006) compared solar modulation estimated from 10 Be versus 14 C. Deviations between the two reconstructions are not larger during the 8.2 kyr event than during other early Holocene periods. Similarly, the different reconstructions shown in Fig. 12 do not point to an exceptionally large imprint of the 8.2 kyr event on inferred solar modulation and thus on radiocarbon production. Solar modulation potential is computed using the relationship between radiocarbon production, the geomagnetic dipole (Korte et al., 2011), and from the models of Masarik and Beer (1999) and Kovaltsov et al. (2012). A key step is the normalisation of the 14 C-derived record to the Forbush instrumental observations over the period 1937 to 1950 AD. Atmospheric 14 C is heavily contaminated after 1950 by atomic bomb tests, preventing a reliable extraction of the natural production signal. In the mid-20th century, atmospheric 14 C and its interhemispheric gradient is influenced by the early signals of fossil fuel combustion and land use. Our spatially resolved model allows us to treat interhemispheric concentration differences and land use explicitly. We show that reconstructions in that rely on the Holocene. To our knowledge, the record by Usoskin et al. (2007) is based on Q from Usoskin and Kromer (2005) and thus on the Northern Hemisphere radiocarbon record only. Their values are systematically lower than ours and their best estimates are in general outside our 1σ confidence band during the last 7 kyr. The record by Steinhilber et al. (2008) shows in general values below our confidence range during the period 7 to 3 kyr BP, while the two reconstructions agree on average during the more recent millennia. Agreement with respect to decadal to centennial variability is generally high among the different reconstructions. However, notable exceptions exist. For example, the reconstruction by Steinhilber et al. (2008) suggests a strong negative fluctuation between 8 and 7.8 kyr BP, not seen in the 14 C-derived records. We reevaluated the claim by Solanki et al. (2004) that the recent sun is exceptionally active. As earlier studies (Muscheler et al., 2007;Steinhilber et al., 2008), we conclude that recent solar activity is high but not unusual in the context of the Holocene. was estimated to be on average at 650 MeV during the last 25 yr by Steinhilber et al. (2008) for a reference Local Interstellar Spectrum (LIS) of Castagnoli and Lal (1980). We find an average for the Holocene (10 to 0 kyr BP) of 555 MeV. Solar activity in our decadally smoothed record is during 28 % of the time higher than the modern average of 650 MeV and during 39 % of the time higher than 600 MeV, used to define periods of high activity by Steinhilber et al. (2008) (Fig. 13). This may be compared to corresponding values of 2 and 15 % by Steinhilber et al. (2008). In contrast to earlier reconstructions, our record suggests that periods of high solar activity (> 600 MeV) were quite common not only in recent millennia, but throughout the Holocene.
Our reconstruction in can be used to derive a range of solar irradiance reconstructions using published methods. Shapiro et al. (2011) suggest that variations in spectral and thus total solar irradiance are proportional to variation in .  (10000) AR (2000) future past Wavelet power spectrum calculated with a Morlet wavelet with log 2 periodicity axis and colour shading. (c) Forecasting of the solar irradiance using an autoregressive model of order p (AR(p)) fitted using Burg's method. Before fitting the time series, a high-pass filter (1000 yr) has been applied in order to remove low-frequency trends. The dashed lines show forecasts with the AR(5000) model starting at year 1300, 1500, 1700 and 1900 AD (magenta dots). Coloured lines show forecasts starting in 2005 AD for different orders of the AR model. Their analysis suggests a high scaling of with TSI and a Maunder Minimum reduction in TSI compared to the present solar minimum by about 0.4 %. This is significantly larger than the reduction suggested in other recent reconstructions Vieira et al., 2011), but within the range of TSI variations explored in earlier climate model studies (Ammann et al., 2007;Jansen et al., 2007). Here, we further converted into TSI using the nonlinear equations used by Steinhilber et al. (2008). This yields small variations in TSI and a Maunder Minimum reduction in TSI compared to recent solar minima of only 0.62 ‰. The resulting overall range in TSI is about 1.2 W m −2 . Future extension of TSI using autoregressive modelling suggests a declining solar activity in the next decades towards average Holocene conditions.
We applied different TSI reconstructions to estimate variations in global mean surface air temperature (SAT) due to TSI changes only. We stress that these SAT anomalies are due to TSI variations only and do not include SAT variations due to other drivers such as orbital forcing, natural and anthropogenic greenhouse gas concentration variations, or changes in the extent and the albedo of ice sheets. The results suggest that there were several periods in the Holocene were TSI-related SAT anomalies were larger than for year 2000 AD. In particular, several warm anomalies occurred in the period from 5500 to 3600 BP.  16. (a, b) Simulated changes in SAT as a response to different TSI forcings: this study (black), VEA11 (red) and SEA12 (blue), left y axis. The green line (right y axis) shows the simulated SAT changes by a solar variation based on our history but converted to TSI using the method of Shapiro et al. (2011). The anomalies were calculated as the difference between two transient simulations, one with all forcings and the other with the same forcings but solar forcing taken to be constant. The last 1000 yr are shown again in (c). All lines show smoothing splines with a cutoff period of 30 yr.
are comparable to twentieth century anomalies in the period from 200 BC to 600 AD. In contrast, the TSI reconstruction of Steinhilber et al. (2012) implies that SAT anomalies due to TSI changes were below current values almost during the entire Holocene. In conclusion, our reconstruction of the radiocarbon production rate and the solar modulation potential, as well as the implied changes in solar irradiance, provide an alternative for climate modellers.

Estimation of uncertainties A1 Assessment of uncertainties in the 14 C production record
There are several potential sources of error entering the calculation of the production rate, summarised in Fig. A1. As already discussed in detail in the main text, the uncertain deglacial carbon cycle changes translate into an uncertainty in Q, especially in the early Holocene. To account for this, we use the arithmetic mean of BIO and CIRC as our best guess Q reconstruction and the difference between them as our uncertainty range, which is around 4.5 % in the early Holocene and drops below 1 % at 5000 kyr BP. In addition, also the uncertainty of the IntCal09/SHCal04 record translates into an error of the production record. 100 radiocarbon histories were randomly generated with Gaussian-distributed 14 C values at each point in time of the record, varied within 1σ uncertainty as (shown in Fig. 2b) on their original time grid (5 yr spacing). For each realisation of the 14 C history, the transient simulation with Bern3D-LPJ was repeated and the resulting Q smoothed. The 1σ range from the 100 smoothed records (Fig. A1b) is of the order of 4-5 %, decreasing towards preindustrial times to 1.5 %. During 1930During -1950, this error increases again rapidly (hardly visible in Fig. A1b) to 10 % because the atmospheric carbon inventory increases at a high rate. Therefore, rather  A1. Assessment of uncertainties in our production rate record in mol yr −1 . (a) Radiocarbon production rate (arithmetic mean of CIRC and BIO); (b) the uncertainty stemming from the pre-Holocene carbon cycle evolution calculated as the absolute difference between results from simulations BIO and CIRC; (c) the smoothed 1σ error from the 14 C record (as shown in Fig. 2b) is calculated by a Monte Carlo approach including 100 individual simulations; and (d) the uncertainty in air-sea and air-land fluxes deduced by sensitivity experiments. The total error in production (e) varies between 30 and 5 mol yr −1 , corresponding to a relative error of ∼ 6 to 1.5 %. The bulk uncertainty in the average level of Q resulting form an uncertain ocean (and land) inventory is not included in this calculation. We estimate this uncertainty to be ∼ 15 %. small errors in 14 C of ∼ 1 ‰ translate into rather big uncertainties in the change of the atmospheric radiocarbon inventory, and thus into the production rate.
The air-sea flux of radiocarbon is the dominant flux out of the atmosphere with approximately 430 mol yr −1 . We assess its variations by nominally changing the air-sea gasexchange rate within ±15 % . Two additional runs were performed in which the seasonally and spatially varying gas exchange rates of the standard setup were increased/decreased everywhere by 15 %. In analogy, uncertainties in the flux into the land biosphere were assessed by varying the gross primary productivity of the model by ±15 %. In these two experiments, differences in the total 14 C inventory were subtracted to get the uncertainties in the fluxes alone; the effect of the inventory is discussed separately. These errors in 14 F as + 14 F ab vary together with changes in the atmospheric 14 C signature. The error from the uncertainty in the combined air-to-sea and air-to-land ( 14 F as + 14 F ab ) fluxes is smaller than 2.5 % of the production at any time (Fig. A1d). It is important to note that the uncertainty in gas exchange and GPP do not alter the absolute value of Q, but change the amplitude in the variations of Q.
Uncertainties are also related to the 14 C signature associated with the carbon flux applied to close the atmospheric CO 2 budget ( 14 F budget in Eq. 1). For example, the missing carbon sink process could be due to a slow down of the oceanic thermohaline circulation (THC) or a growth of land vegetation: two processes with different isotopic signatures. We assume that this signature varies between ±200 ‰ relative to the contemporary atmosphere. The associated error in the production is in general small with peaks of 1-2 % (Fig. A1e). However, towards 1950 AD, where the anthropogenic influence leads to increasing CO 2 levels, the uncertainty reaches up to 4 % as the model simulates a too-weak oceanic/land sink and/or too-high emissions.
We estimate the total standard deviation by quadratic error addition, σ Q = σ 2 Q,i , to be ∼ 6 % during the early Holocene, decreasing to 1.5-2.5 % towards preindustrial times (Fig. A1f). The total uncertainty in the early Holocene is dominated by the uncertainty in the 14 C input data and the deglacial history. These two sources of error decrease towards preindustrial times. In the first half of the last century, uncertainties increase again mainly associated with the anthropogenic CO 2 increase. This increasing uncertainty towards 1950 AD (∼ 9 %) is important when normalising reconstructed solar modulation to the recent instrumental records.
As noted in the main text, the 1σ error is estimated to be 15 % for the ocean 14 C inventory (R. Key, personal communication, May 2013) and for the total 14 C inventory. This potential systematic bias is not included in the uncertainty band of Q as it would not change the temporal evolution.

A2 Propagation of uncertainties from Q to and TSI
Next, the uncertainty calculation of and of TSI is discussed (grey bands in Figs. 12 and 14). Three general sources of uncertainty in are considered: the uncertainty in Q, the uncertainty in VADM and finally the error of the production rate simulations (i.e. the slope of /Q) assumed to be 10 %. Using Monte Carlo calculations, we varied the different sources of uncertainty individually in every year to assess the propagation of errors. In Fig. A2b-d the individual 1σ errors of this calculation are shown. The total error in varies then between 300 MeV (early Holocene) and 70 MeV (preindustrial). All three components of error contribute to approximately equal parts to the total uncertainty. The possible bias introduced by the normalisation is not included in this figure. This would add another ∼ 130 MeV uncertainty to the average level of (but not on its temporal evolution). Note that any constant systematic offset in Q would not translate into an offset in the solar modulation potential as Q is normalised before converting to .
Similarly, the uncertainty in the change in TSI as reconstructed with the use of Eqs. (3) and (4) is calculated (Fig. A3b and c). The 1σ error of is taken into account as well as the uncertainties in Eq. (4) and of the exponent α in Eq. (3) (±0.3) (Steinhilber et al., 2009). The error in the reference value of TSI (1264.64 W m −2 in Eq. 4) does not affect relative changes in TSI, but its absolute level. The total 1σ uncertainty in TSI ranges between 0.2 and 0.8 W m −2 (Fig. A3d), to equal part stemming from the error in and the conversion model.
More than half of the uncertainty in TSI is of systematic nature and related to the -TSI relationship (Fig. A3c). More specifically, the slope of the linear relationship between the interplanetary magnetic field and TSI (Eq. 3) is highly uncertain with 0.38 ± 0.17 units. A lower slope would correspond to correspondingly lower amplitude variations in TSI and a larger slope to larger amplitudes in TSI, while the temporal structure of the reconstruction would remain unchanged. To distinguish between random and systematic errors, we provide the uncertainty related to the random uncertainty in the radiocarbon data and the total uncertainty arising from both systematic and random errors in panel d of Fig. A3.