A Pronounced Spike in Ocean Productivity Triggered by the Chicxulub Impact

There is increasing evidence linking the mass‐extinction event at the Cretaceous‐Paleogene boundary to an asteroid impact near Chicxulub, Mexico. Here we use model simulations to explore the combined effect of sulfate aerosols, carbon dioxide and dust from the impact on the oceans and the marine biosphere in the immediate aftermath of the impact. We find a strong temperature decrease, a brief algal bloom caused by nutrients from both the deep ocean and the projectile, and moderate surface ocean acidification. Comparing the modeled longer‐term post‐impact warming and changes in carbon isotopes with empirical evidence points to a substantial release of carbon from the terrestrial biosphere. Overall, our results shed light on the decades to centuries after the Chicxulub impact which are difficult to resolve with proxy data.

surface ocean cooling in the months to decades after the impact (Vellekoop et al., 2014(Vellekoop et al., , 2016, followed by longer-term warming (Vellekoop et al., 2014(Vellekoop et al., , 2018. Furthermore, proxy studies report rapid ocean acidification which likely started in the century after the impact and lasted for several thousand years (Henehan et al., 2019). The breakdown of the δ 13 C gradient from the surface to the deep ocean (Henehan et al., 2019;Sepúlveda et al., 2019;Zachos et al., 1989) probably indicates severe changes in the ocean's biological pump, but the detailed processes are still poorly understood. In particular, the immediate consequences of the impact on marine productivity during the first years to centuries after the impact cannot be reconstructed by proxy data due to a lack in temporal resolution.
Model experiments simulating the effects of sulfate aerosols from the impact on the climate indicate severe cooling resulting in strong ocean convection (Brugger et al., 2017). The resulting upwelling of nutrients from the deep ocean could increase marine primary productivity following the period of post-impact darkness, but it is unclear whether productivity remains limited by other critical nutrients (in particular iron) and how the marine biosphere is affected by the combined changes in light, temperature, ocean circulation, and nutrient distribution. Here, we address this question by modeling immediate impact-induced changes in climate and nutrient supply and their repercussions for the marine biosphere.

Model
We use the Earth-system model of intermediate complexity CLIMBER-3α+C (Hofmann et al., 2019;Montoya et al., 2005). The ocean model is a modified version of MOM3 (Hofmann & Morales Maqueda, 2006;Pacanowski & Griffies, 1999), run at a horizontal resolution of 3.75° x 3.75° with 24 vertical levels of variable thickness. Additionally, CLIMBER-3α+C includes a sub-model (Hofmann et al., 2019) which is similar to HAMOCC3.1 (Six & Maier-Reimer, 1996) and accounts for the dynamics of marine biogeochemistry and the carbon cycle. The stable carbon isotope fractionation of phytoplankton is parameterized as a function depending on the concentration of aquatic carbon dioxide (Hofmann et al., 1999), while for the inorganic carbonate system and the air-sea gas exchange the approach by Maier-Reimer (Maier-Reimer, 1993) is used. The sea-ice model (Fichefet & Maqueda, 1997) is a dynamic/thermodynamic model in two dimensions with the same horizontal resolution as the ocean model. The statistical-dynamical atmosphere model (Petoukhov et al., 2000) has a coarse resolution of 22.5° in longitude and 7.5° in latitude; it contains a module to simulate terrestrial sources, atmospheric transport and deposition of dust (Bauer & Ganopolski, 2010;Hofmann et al., 2019). The model version used for this study includes an improved lapse rate parametrization and more realistic ice and snow albedo values (Feulner & Kienert, 2014).

Boundary Conditions and Pre-Impact Climate State
An equilibrium simulation of the end-Cretaceous climate is the basis of all impact simulations. This equilibrium state is modeled using a Maastrichtian (70 Ma) continental configuration and vegetation distribution (Sewall et al., 2007), a solar constant of 1354 W/m 2 , based on the present-day solar constant of 1361 W/ m 2 (Kopp & Lean, 2011) and a standard solar model (Bahcall et al., 2001), and idealized orbital parameters (circular orbit, obliquity 23.5°). Proxy data for the Late Cretaceous give a wide range of CO 2 concentrations ranging from 240 ppm (Foster et al., 2017;Nordt et al., 2003) to 1,500 ppm (Royer, 2006). For the time just before the impact, there is evidence that atmospheric CO 2 concentrations were below 800 ppm (Foster et al., 2017;Hong & Lee, 2012;Nordt et al., 2002Nordt et al., , 2003Royer et al., 2012). We have adjusted the carbonate chemistry in the ocean to yield atmospheric CO 2 levels close to 500 ppm, a medium value within the range of proxy estimates: Starting from close to pre-industrial conditions (≈280 ppm), we have modified the ocean's alkalinity until the model approaches equilibrium at a value of 502 ppm after 11,400 model years. The equilibrium pre-impact state has a total alkalinity of 2,348 μeq L −1 and global annual mean amounts of dissolved inorganic carbon (DIC) of 2,328 μmol L −1 , of inorganic phosphate of 2.28 μ mol L −1 and silicate of 81.4 μ mol L −1 . For comparison, the measured present-day values are 2,353 μeq L −1 for the total alkalinity, 2,241 μ mol L −1 for DIC (Key et al., 2004), 2.1 μ mol L −1 for phosphate and 88.2 μ mol L −1 for silicate (objectively analyzed climatologies, Garcia et al., 2018). In addition to this baseline simulation with ≈500 ppm of CO 2 , we have performed an equilibrium simulation with ≈1,165 ppm for sensitivity tests (see supporting information).
The simulated global annual mean surface air temperature for the late-Cretaceous climate state with ≈500 ppm of atmospheric CO 2 is 19.0°C, about 4°C higher than the pre-industrial temperature. Further discussion of this pre-impact state including figures for the temperature distribution as well as a description of the sensitivity simulation with a higher atmospheric CO 2 concentration can be found in the supporting information. The modeled global mean marine net primary productivity (NPP) for the pre-impact climate state is 48 GtC/year. In a late-Holocene configuration of the model, the simulated global mean marine NPP is about 27 GtC/year, at the lower end of the range of CMIP5 values (Fu et al., 2016), but well within the broad interval (50 ± 28 Gt C/year) of observational estimates (Malone et al., 2017).

Modeling the Climatic Effects of the Impact
We explore the effects of sulfate aerosols, carbon dioxide (CO 2 ) and dust produced by the Chicxulub impact on the climate and the marine biosphere. To model the radiative effect of the sulfate aerosols we use a time series of visible transmission for 2.1 years stratospheric residence time and for a load of 100 Gt S with 80% of the sulfur ligated as SO 2 and 20% as SO 3 (Brugger et al., 2017;Pierazzo et al., 2003). Note that more recent investigations with geophysical impact models suggest higher sulfur masses based on new limits for the impact angle and the target composition (Artemieva et al., 2017). However, as the radiative effect does not increase for sulfur loads larger than 30 Gt (Pierazzo et al., 2003), we do not consider the radiative effect of larger sulfur masses. In our standard simulations it is assumed that sulfate aerosols are transferred directly to the ocean after leaving the stratosphere (Pierazzo et al., 2003), with 100% of the sulfur passed to the ocean. To derive a time sequence for the sulfur flux to the ocean we use data for the stratospheric sulfate aerosol concentration decreasing over 6 years (Pierazzo et al., 2003) to calculate the monthly flux to the ocean. We also test the influence of additional sulfur deposition to the ocean on shorter timescales of 10 days after the impact as it was suggested that sulfate aerosols are scavenged by larger and faster settling silicate particles (Ohno et al., 2014;Tyrrell et al., 2015).
Additionally, we consider the atmospheric and oceanic effects of an increased atmospheric CO 2 concentration due to the impact. For our standard impact simulation we use an additional carbon mass of 115 Gt. This represents the carbon release primarily from the vapourization of carbonate rocks (Artemieva et al., 2017). CO 2 is equally distributed over 10 days after the impact assuming instantaneous formation. As for the sulfate aerosols, a globally uniform distribution is used. We also test the sensitivity to C masses between 0 and 4,115 Gt (Artemieva et al., 2017;Pierazzo et al., 1998Pierazzo et al., , 2003Tyrrell et al., 2015), considering additional C from wildfires (additional 1,500 Gt) and decomposed soil organic carbon (additional 2,500 Gt assuming to be instantaneously released, Tyrrell et al., 2015). We use δ 13 C values of −3‰ for the C from the carbonate rocks (Kettrup et al., 2000) and −27‰ for C from organic material, assuming a plant mixture of C 3 and C 3 +CAM plants (Gilman & Edwards, 2020;Osborne & Sack, 2012).
As it was shown that the radiative effects of impact dust were dominated by the effects of sulfate aerosols (Pierazzo et al., 2003;Pope, 2002) and possibly soot (Bardeen et al., 2017), we only consider the effects of the impact dust from the projectile on the ocean, in particular from the important nutrients iron and phosphorus. We explicitly model the atmospheric transport and deposition of impact dust from the projectile to derive the dust distribution in the ocean. To estimate the amount of iron and phosphor distributed in the ocean, we assume that the impactor was a carbonaceous chondrite of type CV, CO or CR (Kyte, 1998). These chondrites have typical iron and phosphorus fractions of 22.5% and 0.104%, respectively (Braukmüller et al., 2018). The modeling of the dust distribution and transport as well as the increase of the bioavailability after the impact is described in detail in the supporting information. For sensitivity tests, we also run the standard impact simulation (115 Gt C and 100 Gt S) and an impact simulation with additional 1,500 Gt C for the pre-impact state with 1,165 ppm of CO 2 .

Temperature Changes in the Millennium After the Impact
The strong reduction in incoming solar radiation due to the stratospheric sulfate aerosols from the impact (Brugger et al., 2017;Pierazzo et al., 2003) has far-reaching implications for climate and life. This is most evident in the severe global cooling in the first decades after the impact. Global annual mean sea surface temperatures (SSTs) are reduced from 21.4°C by about ∼15°C in the fourth year after the impact ( Figure 1). Temperatures recover over the following centuries, approaching pre-impact values after 1,000 years. For a higher late-Cretaceous CO 2 concentration, the impact induces a cooling of very similar magnitude in the immediate aftermath of the impact. More detailed results and discussion for the impact simulations with higher CO 2 concentration can be found in the supporting information (Table S1, Figures S9-S11 and supporting information S5).
Tracing the decadal-scale cooling induced by the impact in geological records is challenging. However, high-resolution proxy data for the Brazos River region (Vellekoop et al., 2014) and for the New Jersey shelf region (Vellekoop et al., 2016) indicate cooling of 7°C (resolved on timescales of months to decades) and 3°C (resolved on timescales of decades to centuries), respectively, in good agreement with the cooling simulated in our model and averaged over appropriate timescales (see Table S1 in supporting information). Furthermore, the data for the Brazos region (Vellekoop et al., 2014(Vellekoop et al., , 2018 indicate a long-term increase of SST of 1−2°C in the early Paleogene starting in the centuries after the impact. In our simulations, the amplitude of SST changes at this location is very similar to the global value; therefore we can directly compare the global SST curves shown in Figure 1 with the proxy estimates. With carbon emissions from the impact only (and the response of the marine biosphere), our standard simulation shows no significant warming in the centuries after the impact. However, one would expect additional carbon emissions from the terrestrial biosphere (e.g., from dying vegetation, wildfires, or soil decomposition) which are not taken into account in this model experiment. Simulations with a series of larger carbon emissions (Figure 1) show that the best agreement with the Brazos River proxy data is achieved for additional 1,500 Gt C from terrestrial sources for the late-Cretaceous simulations with 500 ppm atmospheric CO 2 concentration. For a higher pre-impact CO 2 level, the long-term warming after the impact is smaller for the same amount of additional carbon (see Table S1, Figure S9 and supporting information S5), which is to be expected for a significantly higher baseline CO 2 concentration. Therefore, if the late-Cretaceous atmospheric CO 2 concentration was higher than 500 ppm, the proxy comparison for the pre-impact warming suggests the addition of even more carbon from terrestrial sources.   centuries and first millennia after the impact (Kaiho et al., 1999) also supports a warming after the impact, but the suggested strong warming of ∼5°C (Kaiho et al., 1999) is highly uncertain due to diagenesis, the superposition of different vital effects (e.g., Ziveri et al., 2003) linked to the complex composition of fine fraction carbonate, and the fact that δ 18 O of sea water is not constrained.

A Brief Peak of Net Primary Productivity After the Impact
The Chicxulub impact had profound short-term repercussions for the marine biosphere. In our simulations, the global annual mean NPP in the ocean drops to almost zero in the three years following the impact before briefly peaking at a maximum value about a factor of seven higher than the pre-impact NPP (Figure 2). Although the NPP peak lasts only for a few years, higher productivity persists for decades to centuries, returning to pre-impact levels only after ∼500 years.
The almost complete shutdown of marine NPP in the first years after the impact is caused by both the low light levels and the resulting temperature decrease. With the solar flux starting to recover already two years after the impact and reaching its pre-impact value after seven years (Brugger et al., 2017), global annual mean marine NPP exceeds pre-impact levels of 48 Gt/year in the fourth year after the impact and culminates in a sharp peak at a value of 368 Gt/year in year five. This increase in NPP is caused by two mechanisms: First, the severe surface cooling of the ocean following the impact induces strong ocean mixing and deep water formation, leading to a deep mixed layer, in particular in the mid-latitudes (see supporting information and Brugger et al. (2017) for figures and further discussion). This results in the transport of a large amount of nutrients from the deep ocean to the surface, thus increasing the productivity. Second, considerable amounts of iron and phosphorus are delivered to the ocean by the impact dust originating from the iron-rich projectile. These nutrients with increased bioavailability (see Section 2.3 and supporting information) are sufficient to sustain elevated levels of the NPP for several centuries. In contrast, the nutrients brought to the surface ocean by deep mixing are quickly consumed, as shown by a simulation without dust from the impact (see Figure 2).
The uptake of large amounts of carbon during the period of high ocean productivity results in lower atmospheric CO 2 concentrations and consequently decreases temperatures. Hence, earlier simulations without marine biogeochemistry show a very similar cooling, but a faster recovery where pre-impact temperatures are exceeded already after less than 100 years (Brugger et al., 2017).
Regionally, the increase in marine NPP starts in the tropics where light and ocean temperature exceed the critical thresholds necessary for photosynthesis earliest ( Figure S1   et al., 2017) bringing large amounts of nutrients to the surface ocean, the high primary productivity quickly expands to these regions. The nutrients from the impact dust are quickly spread in the ocean and hence do not cause significant local effects, but amplify the regional NPP patterns also observed in the simulation without impact dust. NPP stays below its pre-impact value for more than 10 years after the impact only in the very high Northern latitudes as a result of the very low ocean temperatures in this region. Note that productivity in coastal regions could have been even higher due to nutrient input from the terrestrial biosphere (Vellekoop et al., 2018) not taken into account in our model.
Blooms of dinoflagellates and cyanobacteria associated with regional maxima in the NPP could result in the production of toxins with detrimental effects for marine near-surface ecosystems. This mechanism described here for the first time could thus play a role in the context of the immediate effects of the impact on the marine biosphere. We note that the increased marine productivity might be dampened by self-limiting feedbacks of the algal bloom (Jackson & Lochmann, 1992;Shigesada & Okubo, 1981) which are not represented in our model. Another aspect not considered here is the formation of cloud condensation nuclei as a consequence of algal blooms. Planktonic algae produce dimethylsulfide (DMS) which form sulfate aerosols in the atmosphere acting as cloud condensation nuclei. Higher densities of cloud condensation nuclei influence cloud albedo and lifetime, and could therefore prolong the cooling after the impact (Charlson et al., 1987;Feulner et al., 2015).
Although there are no empirical data on the productivity immediately after the impact, it has been suggested that marine primary productivity could have been strongly limited only during the short period with direct environmental consequences of the impact, potentially only for a decade (D'Hondt, 2005), consistent with our results. A combined interpretation of proxy data, the fossil record and model results indicates that marine primary productivity recovered to normal or regionally high values (Hollis et al., 1995(Hollis et al., , 2003Lowery et al., 2020;Sepúlveda et al., 2009Sepúlveda et al., , 2019Vellekoop et al., 2017Vellekoop et al., , 2018 much earlier than the biological pump, but that the remineralization depth of organic matter was shallower, leading to a decrease in the export productivity ("Living Ocean" hypothesis, D'Hondt et al., 1998). Estimates for the recovery time for the marine primary productivity are constrained by the age models of the proxy studies and vary between a century (Sepúlveda et al., 2009), and tens of thousands of years (Henehan et al., 2019;Sepúlveda et al., 2019). Hence, our results suggest a convincing scenario for the immediate effect of the impact on primary productivity which cannot be resolved by proxy data and fits the observation of regionally heterogenous recovery of productivity (Sepúlveda et al., 2019). Nevertheless, there might have been some additional, prolonged effects limiting primary productivity in the centuries to millennia after the boundary, both as consequences of the impact as well as other changes across the boundary (Arthur et al., 1987), which are not represented in our model.

Carbon Isotopes Suggest Significant Release of Terrestrial Carbon
The strong perturbations at the K-Pg boundary induce severe changes in the marine carbon cycle leaving traces in the ocean's dissolved inorganic carbon δ 13 C values. This signal depends both on the amount of carbon emitted due to the impact and on its isotopic signature; it is therefore important to assign appropriate δ 13 C values to modeled carbon fluxes from different sources (see Section 2.3). We find an immediate reduction of δ 13 C in the surface ocean ( Figure 3 for the region around Shatsky Rise in the Northern Pacific). The decrease is caused by the enhanced ocean mixing transporting inorganic C with lower δ 13 C from the deeper ocean to the surface, combined with the reduced primary productivity in the first years after the impact. These processes dominate the simulation with only 115 Gt C: after a short decrease δ 13 C starts to rise during the productivity maximum before approaching a new equilibrium slightly above the pre-impact value. In contrast, δ 13 C values in the simulations with carbon released from the terrestrial biosphere are determined by these additional amounts of isotopically light carbon, decreasing the ocean's δ 13 C significantly. Surface δ 13 C starts to increase and approaches a new equilibrium value when the added carbon is exported to the deeper ocean. Hence, the reduction of δ 13 C in the deep ocean is smaller and delayed compared to the surface ocean, resulting in a reduced δ 13 C gradient from the surface to the deep ocean.
The immediate effect of the Chicxulub impact on δ 13 C in the first millennium after the event are resolved in detail in our simulations, but not covered by proxy data due to a lack of temporal resolution. However, if we interpret the earliest surface proxy data after the impact as a smeared signal for the time since the K-Pg boundary, again our simulation with an additional amount of 1,500 Gt C from the terrestrial biosphere shows the best agreement with the proxy data (see Figure 3), confirming the release of significant quantities of carbon from the terrestrial biosphere already suggested by the temperature evolution. This finding also holds for a significantly higher atmospheric CO 2 concentration before the impact (see Figure S10 in the supporting information). We note that a comparison of modeled δ 13 C data for the deep ocean with benthic proxy data is challenging because of the higher variability and lower temporal resolution of the benthic δ 13 C data for this period (Henehan et al., 2019;Hull et al., 2020).

Moderate Ocean Acidification Followed by a Fast Recovery
The paleontologic record of the end-Cretaceous shows a significantly higher extinction rate of marine calcifying organisms (Tyrrell et al., 2015). This could be a consequence of ocean acidification caused by the fast release of CO 2 and H 2 SO 4 during the impact. However, the occurrence of an acidified ocean is still debated, as the observed extinction patterns differ strongly among different calcifying species, making it difficult to get a clear understanding of the extent of acidification and its effects on the biosphere ( D'Hondt, 2005;Dishon et al., 2020;Kiessling & Baron-Szabo, 2004;Kiessling & Simpson, 2011;Tyrrell et al., 2015).
To explore the influence of impact-related carbon and sulfur emissions, we model their effects on the ocean's carbonate saturation state and pH. Consistent with our findings from the changes in temperature and δ 13 C, we assume an additional carbon amount of 1,500 Gt C; for sulfur, we use the most recent estimate of 325 Gt S (Artemieva et al., 2017). Model experiments with higher C and S masses are discussed in the SI (Section 3 and Figures S2, S4 and S5).
We find an immediate reduction of surface global mean pH, with a maximum surface pH decrease of 0.36 units four years after the impact (see Figure S2 in supporting information). The subsequent emerging peak in ocean productivity and the accompanying high carbon uptake leads to a short increase in pH before it decreases again. These fast changes immediately after the impact are followed by a transient recovery over the next 500 years, equilibrating at 0.23 pH units below the pre-impact value. Data from the Geulhemmerberg succession show a pH reduction of ∼0.25 for the millennium following the impact, but possibly even for a shorter time scale of ∼100 years (Henehan et al., 2019). This is in excellent agreement with our modeled pH reduction of 0.28 pH units in this region averaged over the century after the impact.
To assess the influence on calcifiers, Figure 4 shows the surface ocean aragonite saturation state Ω a . The critical value for the saturation state Ω varies for different calcifying species (Ries et al., 2009)   . Changes of surface carbon isotope ratios of dissolved inorganic carbon at Shatsky Rise, Pacific. Colored lines show modeled changes in δ 13 C for the 1,000 years after the impact relative to the 100-year pre-impact mean for simulations with different amounts of carbon. The model results are embedded in changes in surface δ 13 C from proxy data from bulk carbonate ( Table S11 in Hull et al., 2020), assuming an absolute age for the K-Pg boundary of 66.022 Ma (Table S3 in Hull et al., 2020) relative to the mean of the pre-impact proxy values. on marine carbonate chemistry (Bach et al., 2015). Here we assume that calcifying organisms begin to dissolve for Ω < 1. Globally, Ω a is strongly reduced in the 10 years after the impact, and it falls below the critical value of 1 at latitudes higher than ∼40°. Although Ω a quickly starts to recover and the line defining Ω a = 1 retrenches to higher latitudes, Ω a globally stays below the pre-impact values in the 1,000 years after the impact. The calcite saturation state Ω c , which is generally higher than Ω a due to the lower solubility of calcite in ocean water, shows similar patterns (see Figure S3 in supporting information). Ω c falls below the critical value of 1 for latitudes greater than ∼60° averaged for the 10 years after the impact before Ω c = 1 retrenches to higher latitudes.
In summary, we find rapid surface ocean acidification which, however, recovers quickly. The potential habitat for calcifiers after the impact is restricted to lower latitudes. For corals and rudists favoring warm waters, however, it is hard to distinguish whether the strong temperature decrease or the aragonitic undersaturation was more critical in limiting their habitat. As the surface ocean is undersaturated in calcite only at high latitudes, the global extinction of planktic foraminifera (D'Hondt, 2005) is difficult to explain primarily with ocean acidification. Our results thus indicate that surface ocean acidification after the Chicxulub impact caused regionally the extinction of some groups of calcifiers. However, the extinction of a large fraction of globally distributed calcifiers, as observed in the fossil record, cannot be caused by undersaturation alone, because it was not widespread enough.

Discussion
Our modeling results for the immediate consequences of the Chicxulub impact on the climate and the marine biosphere indicate rapid cooling, a sharp and short-lived increase of marine primary productivity following a brief productivity collapse, and moderate ocean acidification. The longer-term warming trend observed in proxy data and the comparison of surface δ 13 C data with our model results suggest an additional carbon release of ∼1,500 Gt from the terrestrial biosphere. If the carbon is the result of wildfires caused by the infrared radiation of the reentering ejecta (Robertson et al., 2013), the produced soot would have further reduced the incoming surface short wave radiation with significant effects on photosynthesis (Bardeen et al., 2017;Tabor et al., 2020). In addition, nutrients brought into the ocean by the wildfire ash could have intensified and prolonged the increase in marine primary productivity and the algal bloom (Abram et al., 2003). Proxy data indicate an additional nutrient flux resulting from land-derived material during the terrestrial mass extinction, leading in particular to increased productivity in the neritic zone during the centuries to millennia following the K-Pg boundary (Vellekoop et al., 2018). Future modeling studies will have to explore the full range of interactions between climate and the terrestrial and marine biosphere. BRUGGER ET AL.  Recent studies on the interaction of Deccan volcanism and the Chicxulub impact support the notion that the impact was the major driver of the end-Cretaceous mass extinction Sprain et al., 2019). From these findings and our results we conclude that the combination of several abrupt and vigorous changes caused by the Chicxulub impact triggered the end-Cretaceous mass extinction. In addition, life-history traits (D'Hondt, 2005), geographic distribution (Dishon et al., 2020), the respective ecosystem (D'Hondt, 2005) and chance survival (D'Hondt, 2005) could have determined the selectivity during the extinction. The interaction of these ecological and evolutionary factors with the modeled direct consequences of the impact would have to be further investigated, for example with ecological niche models.

Data Availability Statement
All model input and output files as well as the preprocessing and postprocessing scripts used to generate model input and the figures in the study are available online (https://doi.org/10.5880/PIK.2020.008). The source code for the model used in this study is archived at the Potsdam Institute for Climate Impact Research and is made available upon request.