On the Sensitivity of the Devonian Climate to Continental Configuration, Vegetation Cover, Orbital Configuration, CO 2 Concentration, and Insolation

During the Devonian (419 to 359 million years ago), life on Earth witnessed decisive evolutionary breakthroughs, most prominently the colonization of land by vascular plants and vertebrates. However, it was also a period of major marine extinctions coinciding with marked changes in climate. The cause of these changes remains unknown, and it is therefore instructive to explore systematically how the Devonian climate responds to changes in boundary conditions. Here we use coupled climate model simulations to investigate separately the influence of changes in continental configuration, vegetation cover, carbon dioxide (CO2) concentrations, the solar constant, and orbital parameters on the Devonian climate. The biogeophysical effect of changes in vegetation cover is small, and the cooling due to continental drift is offset by the increasing solar constant. Variations of orbital parameters affect the Devonian climate, with the warmest climate states at high obliquity and high eccentricity. The prevailing mode of decadal to centennial climate variability relates to temperature fluctuations in high northern latitudes which are mediated by coupled oscillations involving sea ice cover, ocean convection, and a regional overturning circulation. The temperature evolution during the Devonian is dominated by the strong decrease in atmospheric CO2. Albedo changes due to increasing vegetation cover cannot explain the temperature rise found in Late Devonian proxy data. Finally, simulated temperatures are significantly lower than estimates based on oxygen isotope ratios, suggesting a lower δ18O ratio of Devonian seawater.


Introduction
The Devonian (419 to 359 Ma, 1 Ma = 1 million years ago) is a key period in Earth's history characterized by fundamental changes in the atmosphere, the ocean, and the biosphere. With respect to the evolution of life, the Devonian is best known for the diversification of fish as well as the colonization of the continents by vascular plants and vertebrates. Although microbial mats and nonvascular plants could be found on land even before the Devonian (Boyce & Lee, 2017), the appearance of vascular land plants beginning in the Early and Middle Devonian was certainly an important first step toward modern land ecosystems. By the Late Devonian, vascular plants had vastly diversified, going hand in hand with the evolution of more advanced leaves and root systems (Algeo & Scheckler, 1998;Algeo et al., 1995). In the ocean, coral stromatoporoid reefs reached their largest extent during the Phanerozoic in the Middle Devonian (Copper & Scotese, 2003) and fish evolved into rich diversity (Dahl et al., 2010). Finally, in the Late Devonian, the first tetrapods moved from ocean to land (Brezinski et al., 2009;Clack, 2007).
While the Devonian is best known for these evolutionary breakthroughs, it is also a period of species mass extinctions. The extinction rate during the Devonian is marked by three distinct peaks (Bambach, 2006), with the highest pulse (the Frasnian-Famennian mass extinction, 378-375 Ma) ranking among the "Big Five," that is, the five most severe mass extinctions in Earth's history (Bambach, 2006). The cause of these extinctions, which mostly took place in the ocean (Bambach, 2006), is still under debate, as it is challenging to explain the episodic nature and duration of the extinctions (Algeo et al., 1995). Discussed causes for the Frasnian-Famennian extinction are a bolide impact (McGhee, 1996;McGhee et al., 1984), volcanic activity (Ma et al., 2016), changes in sea level (Bond & Wignall, 2008;Ma et al., 2016), rapid temperature variations (McGhee, 1996;Ma et al., 2016), and the development of ocean anoxic waters (Bond & Wignall, 2008;Ma et al., 2016).
These changes in atmospheric composition, and in particular the drop in atmospheric CO 2 concentrations, also resulted in climatic changes, which in turn affected the biosphere. Indeed, 18 O oxygen isotope data from marine microfossils (a proxy for seawater temperature shown in Figure 1b) indicate a greenhouse climate in the Early Devonian and much cooler temperatures in the Middle Devonian (Joachimski et al., 2009). For the Late Devonian, proxy studies (Joachimski et al., 2009;van Geldern et al., 2006) indicate rising temperatures again which are still challenging to reconcile with the decreasing CO 2 concentrations. During the late Famennian, Earth's climate cooled again (Brezinski et al., 2009;Joachimski et al., 2009), with some studies even indicating glaciations (Brezinski et al., 2008(Brezinski et al., , 2009Caputo, 1985;Caputo et al., 2008). Furthermore, the Late Devonian was a period of fast sea level variations (Brezinski et al., 2009;Haq & Schutter, 2008;Johnson et al., 1985;Ma et al., 2016;Sandberg et al., 2002), with periods of sea level transgression being linked to the development of ocean anoxic conditions (Bond & Wignall, 2008;Johnson et al., 1985).
As any other period in Earth's history, the 60 million years spanning the Devonian is characterized not only by long-term changes in temperature but also by large fluctuations in temperature around these longer-term trends (see Figure 1b). Furthermore, there are several studies stressing the importance of astronomical forcing for the climate during the Devonian. The geologic record of the Devonian shows cyclic structures (De Vleeschouwer et al., 2013, 2017 which can be interpreted as the result of astronomical cycles according to Milankovitch theory (Milankovitch, 1941). The configuration of Earth's orbit and rotational axis determines the total amount as well as the spatial and temporal distribution of solar radiation and therefore impacts climate. In addition, identifying astronomical cycles in the geologic record can help assign a timescale to cyclic features observed in the geologic record and thus a timescale for palaeoclimatic changes (De Vleeschouwer et al., 2013, 2014, 2017. In an effort to link the changes in the various components of the Devonian Earth system, there are several studies investigating potential connections between land plant evolution, climate change, and the oceanic mass extinction (Algeo & Scheckler, 1998;Algeo et al., 1995;Berner, 1994;Goddéris & Joachimski, 2004;Le Hir et al., 2011;Simon et al., 2007). It is suggested, for example, that the increase of weathering due to the spreading of plants on land could have been a cause of the decrease in atmospheric carbon dioxide (Algeo & Scheckler, 1998;Algeo et al., 1995;Berner, 2006). Additionally, the increased weathering rates could have lead to a higher transport of phosphorus to the ocean, promoting eutrophication with its negative consequences for life in the ocean (Algeo & Scheckler, 1998;Algeo et al., 1995).
Given the multitude of changes in the Earth system during the Devonian and the intricate coupling of atmosphere, ocean, and biosphere, it is challenging to disentangle causes and effects and to determine which forcings are most important for Devonian climate change. In this study, we test the sensitivity of the Devonian climate to different forcings in order to quantify their relevance. We therefore set up simulations with a coupled ocean-atmosphere model considering three different continental configurations representing the Early, Middle, and Late Devonian. In addition, changes in the solar constant, in atmospheric carbon dioxide concentrations, in orbital parameters, and in the spatial distribution of vegetation are systematically tested. Concerning the influence of land plants, this study focuses on their impact on the climate via biogeophysical effects, in particular changes in albedo and evapotranspiration. This paper is organized as follows. Section 2 describes the coupled climate model used in this study, the various sensitivity experiments which serve to investigate the influence of changes in different parameters throughout the Devonian as well as the boundary conditions for three simulations which are representative for climate states of the Early, Middle, and Late Devonian. Results from the sensitivity simulations are presented and discussed in section 3. Section 4 highlights an ocean-sea ice mechanism resulting in climate oscillations at high northern latitudes. Section 5 discusses the climate evolution through the Devonian based on simulations for three time slices and compares the results to other modeling as well as proxy studies. Finally, Section 6 summarizes our key findings and discusses the broader implications of our results.

Model Description
To be able to assess the sensitivity of the Devonian climate to a wide range of parameters, we use CLIMBER-3 , a relatively fast coupled Earth system model of intermediate complexity (Montoya et al., 2005). Its main component consists of a modified version of the ocean general circulation model MOM3 (Hofmann & Morales Maqueda, 2006;Pacanowski & Griffies, 1999) used at a horizontal resolution of 3.75 • × 3.75 • with 24 vertical levels. The coupled model further includes a dynamic/thermodynamic sea ice model (Fichefet & Morales Maqueda, 1997) operated at the resolution of the ocean model and a fast statistical-dynamical atmosphere model (Petoukhov et al., 2000) with a coarse spatial resolution of 22.5 • in longitude and 7.5 • in latitude. Atmosphere-surface interactions (including the biogeophysical effect of terrestrial vegetation) are modeled by a land surface model distinguishing six prescribed surface types (open ocean, sea ice, forest, grass, bare soil, and glaciers) characterized by albedo, roughness, and evapotranspiration. CLIMBER-3 compares well with other models in model intercomparison projects (e.g., Eby et al., 2013;Zickfeld et al., 2013). For the simulations of the Devonian climate presented in this paper, we employ a model version with improved lapse rate parametrization and more realistic ice and snow albedo values Note. Settings for the continental configuration, the solar constant S, the atmospheric CO 2 concentration, the vegetation distribution, the obliquity of Earth's axis, the eccentricity e of its orbit, and the precession angle are listed for all simulations. (Feulner & Kienert, 2014) already used in studies of a wide variety of deep-time palaeoclimate problems (e.g., Brugger et al., 2017;Feulner, 2017;Feulner et al., 2015). The surface model of CLIMBER-3 is also used in the coarser-resolution model CLIMBER-2 and compares well with more sophisticated models when simulating the biogeophysical effects of vegetation on climate (Ganopolski et al., 2001).

Boundary Conditions and Design of the Numerical Experiments
The main aim of this study is to simulate the climate evolution through the Devonian by designing three "best guess simulations" representative of the Early, Middle, and Late Devonian. These simulations use the most likely boundary conditions in terms of the distribution of continents, vegetation, the solar constant, and the atmospheric carbon dioxide level. In addition, we systematically investigate the sensitivity of the Devonian climate to changes in continental configuration, vegetation distribution, solar constant, atmospheric carbon dioxide levels, and orbital parameters by varying each of these boundary conditions separately while keeping all others fixed. The insights gained from these sensitivity experiments help to interpret the results from the best guess simulations and thus to identify the dominant drivers to climate changes during the Devonian. Apart from the sensitivity experiments investigating the influence of changes in orbital parameters (section 2.2.5), we use orbital parameters of an idealized (or median) orbit with obliquity = 23.5 • and eccentricity e = 0. This obliquity is similar to the present-day value; using a circular orbit makes this configuration independent of precession; that is, it minimizes seasonality for a given obliquity. All model simulations are initialized from an ice-free state with constant present-day ocean salinity and a smoothed sea surface temperature (SST) profile which is symmetrical about the equator and approximates modern observations. The simulations are integrated for 4,000 model years until climate equilibrium is approached. The design of the sensitivity experiments (sections 2.2.1 to 2.2.4) as well as of the best guess simulations (section 2.2.6) will be described in more detail in the following.

Continental Configurations
To test the influence of the continental configuration on the Devonian climate, we set up three model runs using continental configurations representing three time slices during the Early Devonian (415 Ma), Middle Devonian (380 Ma), and Late Devonian (360 Ma), respectively, based on reconstructions (Scotese, 2014). The three different Devonian continental configurations, shown in Figure 2, are characterized by an arrangement of the continents' largest fraction in the Southern Hemisphere and no land northward of 60 • N. Throughout the Devonian, the most significant changes affect the distribution of land and ocean areas around the South Pole. For our sensitivity experiments with respect to these three Devonian continental Paleoceanography 10.1029/2019PA003562 configurations, it is assumed that there are no land plants; that is, the land surface type is set to bare land throughout, and all other boundary conditions remain fixed; see Table 1.

Solar Constant
The effect of changes in the solar constant is investigated by adjusting the solar constant for each Devonian timeslice based on the best estimate (Kopp & Lean, 2011) for the present-day solar constant of 1,361 W/m 2 and a scaling according to a standard solar evolution model (Bahcall et al., 2001), yielding values of 1,315.3, 1,319.1, and 1,321.3 W/m 2 for the Early, Middle, and Late Devonian, respectively. All other parameters are fixed: We use continental configuration and a CO 2 concentration (1,500 ppm) representative for the Middle Devonian (380 Ma) orbital parameters for the median orbit and set vegetation to bare land.

CO 2 Concentration
We test the influence of changes in the atmospheric carbon dioxide level on the Devonian climate using CO 2 concentrations of 2,000, 1,500, and 1,000 ppm. These particular values have been chosen because they are representative for the values reconstructed for the Early (415 Ma), Middle (380 Ma), and Late Devonian (360 Ma) periods, each lying within one standard deviation of the pCO 2 probability maximum given in Foster et al. (2017); see Figure 1a. Note that we chose to use approximate values rather than reading off the fitted values for the three time slices because of the sparse temporal sampling and large error bars of the CO 2 reconstruction. All other parameters are fixed: We use the continental configuration and solar constant (S = 1, 319.1W/m 2 ) representative for the Middle Devonian (380 Ma), idealized orbital parameters, and set vegetation to bare land.

Vegetation Distributions
The biogeophysical influence of the spread of vascular land plants on the Devonian climate is investigated by different idealized vegetation distributions. These are roughly representative of the situation during the Early, Middle, and Late Devonian and will thus be also used in the simulations representing the three Devonian time slices described in section 2.2.6. The sensitivity experiments with respect to the effect of vegetation cover alone, however, are based on the continental configuration for 380 Ma to isolate the biogeophysical effect of vegetation cover from the changes due to continental drift.
The design of the vegetation distributions is based on palaeontological evidence for the evolution of vascular plants during the Devonian (Boyce & Lee, 2017). In particular, the evolving ability of land plants to spread upcountry is taken into account, which we quantify by imposing limits in terms of distance from the coast and elevation at the original resolution of the continental reconstruction. Furthermore, the evolution of plant communities from bryophyte-dominated coastal vegetation toward the first forests is considered by varying the mix of plant types and their biogeophysical properties albedo, roughness length, and evapotranspiration; see Table 2. Finally, vegetation distribution is constrained by climate. The derivation of the vegetation distributions for the different Devonian time slices is explained in detail in the following and shown in Figure 3.
Nonvascular land plants arose as early as in the Mid-Ordovician (Rubinstein et al., 2010) and may have spread widely already during the Ordovician (Porada et al., 2016). The first vascular plants evolved during the Late Ordovician (Steemans et al., 2009). During the Late Silurian and Early Devonian, the first vascular land plants were distributed across all latitudes (Boyce & Lee, 2017), but due to shallow roots and no water conducting cells, they were restricted to wet lowland areas (Algeo & Scheckler, 1998;Algeo et al., 1995;Boyce & Lee, 2017). For the very Early Devonian, we therefore assume that land plants are restricted to areas closer than 500 km to the coast and up to 500 m in altitude. Areas closer than 300 km to the coast and/or with an altitude lower than 300 m are covered by 60% vegetation. Areas with a coastal distance between 300 and 500 km and/or an altitude between 300 and 500 m are covered by 10% vegetation. To represent the dominance of nonvascular vegetation, we assume the plant cover to consist of two-thirds bryophytic vegetation mixed with one-third grass-like vegetation for the Early Devonian (see Table 2). Note that grass has evolved much later in Earth's history and that the grass-like surface type in our model approximates the basic biogeophysical properties of any shrub-like vegetation. For the remaining surface area the land surface type in the model is set to bare land.
By the Middle Devonian, vascular land plants gained in importance and the first arborescent vegetation evolved (Algeo & Scheckler, 1998). Thus, for our shrub-like vegetation, we mix four fifths of the vegetation of shrub-like characteristics with one-fifth vegetation of bryophyte characteristics. We add a second vegetation type representing an early type of trees (see Table 2). We assume 70% Middle Devonian shrub-like vegetation Note. Parameters of bare soil, shrubs, and trees are the standard values in our model. Albedo and roughness values for the early trees are based on the woodland parameters in Dickinson et al. (1993); the evapotranspiration value is an estimate based on the values used for shrubs and trees. The bryophytes' albedo is an average value of the moss types in Stoy et al. (2012); roughness length and LAI are taken from Dickinson et al. (1993) for short grass. Evapotranspiration is proportional to LAI; see equation (45) in Dickinson et al. (1986). LAI = Leaf Area Index.
cover and 10% early tree vegetation cover in areas lower than 500 m in altitude and closer than 500 km to the coast and set the remaining surface area to bare land.
In the Late Devonian, trees grew larger and diversified (Algeo & Scheckler, 1998) and vegetation was more widespread but did not extend in upland areas (Boyce & Lee, 2017). Therefore, for this model setup, we use the same shrub-like vegetation mix as in the Middle Devonian but a tree type with characteristics of an evergreen forest after Dickinson et al. (1993). Coastal areas lower than 500 m in altitude and closer than 500 km to the coast are covered by 40% of shrub-like vegetation and 40% of trees. On areas from 500 to 1,000 m in altitude and/or between 500 and 1,000 km from the coast we assume 30% shrub-like vegetation and 30% trees.
As vegetation and climate do not interact dynamically in our model, we implement temperature and precipitation limits to avoid vegetation on land areas with unsuitable climate. Trees do not cover cells for which the average temperature of the coldest month is below −25 • C, whereas the growth of shrub-like vegetation is not restricted by temperature (Schaphoff et al., 2018; supporting information Table S4). The value for the temperature limit for trees is an estimate based on an average of different types of trees. Additionally, we restrict all vegetation to areas with annual precipitation above 700 mm. This value is based on Lieth (1973, Figure 3) taking into account that Devonian plants can store less water compared to present-day plants. The areas unsuitable for vegetation according to these limits are determined using a simulation with all land cover set to bare land and are indicated in Figure 3. There is a marked increase in the fraction of land covered by vegetation from 15.1% during the Early Devonian to 28.5% during the Middle Devonian and 37.3% during the Late Devonian.
The resulting vegetation distributions at 1 • spatial resolution are then used to calculate vegetation fractions on the coarser atmospheric grid of the model. The bias due to the coarse model resolution should be lim- ited, however, because the atmosphere model takes subgrid fractions into account when calculating surface fluxes.
In order to assess the influence of the different vegetation types separately, we add three simulations representing extreme scenarios in which all land is covered by bare soil, shrub-like vegetation, and trees, respectively, irrespective of coastal distance, elevation, temperature, and precipitation.
All sensitivity simulations for the different vegetation distributions are run for a fixed carbon dioxide concentration of 1,500 ppm, a solar constant of S = 1, 319.1 W/m 2 and orbital parameters of a median orbit; see also

Influence of the Continental Configuration
First, we investigate the sensitivity of the Devonian climate to changes in continental configuration while keeping all other boundary conditions fixed (CO 2 concentration 1,500 ppm, solar constant 1,319.1 W/m 2 , bare land, and median orbital parameters; see Table 1 and section 2.2.1). We compare three simulations using the continental configurations shown in Figure 2 and find that the change in the distribution of continents throughout the Devonian leads to a decrease of global annual mean surface air temperature over time. However, the global mean temperature differences between the climate states for the three continental configurations are small: for the Early Devonian, global annual mean surface air temperature is 21. Although the decreasing trend in surface air temperature for the three Devonian continental configurations is small, regional effects caused by changes in continental drift are significant: The largest temperature differences are found at high southern latitudes and can mostly be attributed to the differences in the distribution of land and ocean areas (see Figure S1). Furthermore, differences in orography between the three continental configurations result in regional temperature differences. Finally, part of the differences in temperature can be attributed to changes in ocean circulation between the different continental configuration as indicated by the ocean surface velocities and temperatures shown in Figure S2. Note that changes in continental configuration could in principle effect the climate-carbon feedback which is not represented in our model, possibly inducing stronger climatic changes than modeled here.

Influence of the Solar Constant
Annual global mean surface air temperatures increase from 20.9 • C to 21.2 • C and 21.4 • C with increasing solar constant (1,315.3 W/m 2 in the Early, 1,319.1 W/m 2 in the Middle, and 1,321.3 W/m 2 in the Late Devonian) while keeping all other boundary conditions fixed (see Table 1 and section 2.2.2). These temperature changes are consistent with what one would expect from the change in solar forcing. The effect of an increasing solar constant is largest in the high northern latitudes (see Figure S3), in particular during the winter months. This is closely linked to differences in Arctic sea ice triggered by the differences in global mean temperature. Note that, on a global scale, the warming due to the brightening Sun is very similar in magnitude to the cooling due to continental drift over the Devonian described in section 3.1 above. Both effects will therefore roughly cancel each other. Regionally, however, the cooling due to continental drift is stronger at high southern latitudes, whereas the warming due to the increasing solar constant is more pronounced at high northern latitudes.

Influence of CO 2 Concentrations
Next, we evaluate the impact of the decrease in atmospheric CO 2 concentrations throughout the Devonian (2,000 ppm in the Early, 1,500 ppm in the Middle, 1,000 ppm in the Late Devonian; see section 2.2.3) on the Devonian climate while keeping all other boundary conditions fixed (see Table 1). Global annual mean surface air temperature decreases from 22.4 • C for 2,000 ppm to 21.2 • C for 1,500 ppm and 19.4 • C for 1,000 ppm. The temperature difference of 2.9 • C between the climate states with 2,000 and 1,000 ppm of atmospheric carbon dioxide is, of course, equivalent to the model's equilibrium climate sensitivity to CO 2 doubling for Devonian boundary conditions. This value is in good agreement with the range for the present-day climate sensitivity given in IPCC (2013IPCC ( , p.1110. For the relatively warm Devonian climate state, climate sensitivity is smaller than for cooler configurations analyzed with the same model (Feulner & Rahmstorf, 2010;Feulner & Kienert, 2014), predominantly due to the smaller ice-albedo feedback.
With respect to regional climate change, differences in surface air temperature between different CO 2 concentrations are significantly larger for higher latitudes, especially in the Northern Hemisphere (see Figure  S4). This is similar to the phenomenon of Arctic amplification under present-day global warming (Serreze & Barry, 2011) which is mainly caused by the ice-albedo feedback effect. We will encounter the same pattern when analyzing the best guess simulations in section 5.
Looking at the results of the sensitivity studies so far, it is important to note that the magnitude of the global temperature change due to CO 2 is considerably larger than the effects of continental drift and changes in the solar constant.

Influence of the Vegetation Distribution
To assess the potential climate impact of the colonization of land by vascular plants during the Devonian, we investigate the sensitivity of the climate to the biogeophysical effects of different vegetation distributions: We use the three vegetation distributions constructed to represent the three Devonian periods on the 380-Ma continental configuration (see section 2.2.4) but first compare the three extreme scenarios of bare land, shrubs, and trees covering the entire continents of the 380-Ma configuration.
Global mean surface air temperatures are 21.2 • C for the bare land case, 20.7 • C for the shrub-covered continents, and 21.6 • C for the tree-covered continents. The continental temperature differences are much more significant with 16.2 • C for the bare land case, 13.4 • C for the shrub-covered continents, and 15.0 • C for the tree-covered case.
Maps of differences in surface air temperature, evaporation, surface albedo, and atmospheric water vapor content for the shrub-covered case minus the bare land case and the tree-covered case minus the bare land case (see Figure 4, first two rows) help to understand the origin of temperature differences. Comparing patterns of evaporation and surface air temperature differences over continents, it is obvious that changes in evaporation determine continental surface air temperature differences for the shrub-covered case. For the tree-covered case, this is only true for latitudes up to 60 • . Stronger evaporation results in lower temperatures due to evaporational cooling and dominates over the warming effect of a decreasing surface albedo. This trend is in agreement with the climatic effect of modern tropical forests (Betts, 2000;Bonan, 2008). For modern boreal forests, in contrast, the negative climate forcing of evaporation is counteracted by the low albedo (Betts, 2000;Bonan, 2008). (The effects of modern temperate forests are very uncertain.) Furthermore, simulated total cloud cover is higher over continents in case of vegetation cover, with only small differences between the shrub-covered and the tree-covered case, leading to additional cooling ( Figure S5).
Interestingly, surface air temperatures over the ocean are higher for the tree-covered case than for the bare land case. Trees increase the atmospheric water vapor content which leads to an increased greenhouse effect, resulting in higher surface air temperatures. Although this is a global trend, it dominates surface air temperature differences only for oceanic regions (see Figure 4e). For continental regions, this effect is counteracted by the evaporative cooling described above. The warmer temperatures over the ocean for the tree-covered case lead to less snow in the high southern latitudes, decreasing the albedo and hence increasing continental temperatures. In contrast, for the high southern latitudes of the shrub-covered case there are a small increase in albedo as the continental snow cover increases due to lower temperatures. In the Arctic temperature and albedo changes are closely linked to changes in sea ice cover for both the shrub-covered and the tree-covered case. Still, the regional temperature differences on land are more pronounced (Figure 4, third to fifth row) and can largely be understood with the explanations derived from the more extreme idealized simulations discussed above. Plant evolution during the Devonian is accompanied by an increase in evaporation, leading to cooling over the continents. This cooling effect dominates changes in continental temperature except for the high southern latitudes where the climatic limits (see Figure 2) preclude vegetation. Here, temperature differences are negligible for the Early and Middle Devonian vegetation distribution compared to the bare land case. For the Late Devonian the higher atmospheric water vapor content caused by the higher fraction of trees increases the atmospheric water vapor content globally which dominates surface air temperature patterns over the ocean and over the continents at high southern latitudes. Changes in albedo are very small over continental areas, and therefore, their effect is insignificant in our simulations. Note that land surface albedo may not have changed too much during the Devonian due to the early presence of microbial mats (Boyce & Lee, 2017). As observed for the extreme scenarios, albedo changes in the Arctic region are linked to changes in temperature and sea ice cover for all three Devonian vegetation distributions.

10.1029/2019PA003562
In the following we compare our results of the influence of the Devonian vegetation changes with an earlier modeling study. Le Hir et al. (2011) investigate the Late Devonian climate using a coupled climate-carbon-vegetation model focusing on the impact of changing surface properties on climate. Similar to our approach, they use three representative settings for the Early, Middle, and Late Devonian. As forcing, they use the solid Earth degassing rate, continental configurations, and an increasing solar constant. Although they simulate a strong decrease in CO 2 (5,144 ppm for the Early Devonian, 4,843 ppm for the Middle Devonian, and 2,125 ppm for the Late Devonian), they find an almost unchanged temperature of 23.4 and 23.5 • C for the Early and Middle Devonian, and considering the strong CO 2 decrease, the temperature decrease to 22.2 • C for the Late Devonian is comparatively small. Note that in Le Hir et al. (2011) CO 2 decreases by more than a factor of 2 from the Early to the Late Devonian; the expected cooling for a halving of the CO 2 concentration is about 3 • C (Section 3.3). Le Hir et al. (2011) conclude that the decrease in albedo accompanying the evolution of land plants is counteracting the carbon dioxide decrease. In contrast to our simulations, however, they use relatively high values for the bare soil albedo ( To investigate the impact of the choice of albedo for the different surface types, we repeat our bare-land simulation and our tree-covered simulation with the albedo values used in Le Hir et al. (2011). The comparison with our original vegetation simulation results illustrates how strongly the modeled climate impact of differences in vegetation is determined by the choice of albedo values: Due to the albedo increase for bare land global mean surface air temperature decreases by 1.7 to 19.5 • C and mean continental temperature by 4.2 to 12.3 • C. The strong warming over the continents is shown in Figure S6a. For the tree-covered case, the small differences in albedo have only negligible effects on global mean temperatures, but small regional differences can be seen in Figure S6b.
We note, however, that the soil albedo value used in Le Hir et al. (2011) is very high; given the early presence of microbial mats (Boyce & Lee, 2017) which have a much lower albedo (Sanromá et al., 2013) than the one used in Le Hir et al. (2011), our low value for the soil albedo is better suited to represent rocks covered with microbial mats and therefore appears to be more plausible.
It is important to bear in mind that, in addition to this biogeophysical effects, changes in the weathering process due to the evolution of deeper roots strongly influence the carbon cycle and therefore climate. Hence, these results only give a rough idea of the influence of the land plant evolution on climate. The combined effect of vegetation distribution and decreased CO 2 concentration is discussed in the context of the best guess simulations in section 5.

Influence of the Orbital Configuration
As a further sensitivity test, we investigate the influence of orbital configuration on the Devonian climate using simulations for a range of different orbital parameters while keeping all other boundary conditions fixed (continental configuration for 380 Ma, carbon dioxide concentration 1,500 ppm, solar constant 1,319.1 W/m 2 , and bare land). In these simulations, we systematically test different values for the obliquity ( = 22.0 • , 23.5 • , and 24.5 • ), the eccentricity (e = 0, 0.03 and 0.069), and precession ( from 0 • to 315 • ); see section 2.2.5. Concerning the influence of the shape of Earth's elliptical orbit on the Devonian climate, global annual mean surface air temperatures increase with eccentricity due to the higher annual mean insolation for more eccentric orbits. Surface air temperature differences for different precession angles are negligible in our sensitivity experiments.
Lowest global annual mean surface air temperatures are reached for the configurations with = 22.0 • and e = 0 (21.0 • C), whereas the maximum is reached for the configurations with = 24.5 • and e = 0.069 (21.6 • C). Hence, for our median orbit simulations the uncertainty arising due to the unknown exact orbital configuration is +0.4 −0.2 We note that the differences in climate between orbital configurations in the Devonian might be larger in reality than in our idealized simulations at fixed carbon dioxide levels due to feedbacks involving the carbon cycle.
We can compare our results with an earlier modeling study on the impact of orbital parameters on the Devonian climate (De Vleeschouwer et al., 2014). This study finds a strong dependence of the Devonian climate on orbital parameters and supports a link between ocean anoxic events during the Devonian and orbital forcing. Except for a small eccentricity offset of the coldest state and a slight dependence on precession, De Vleeschouwer et al. (2014) also find that global mean temperature increases with obliquity and eccentricity. However, they report significantly larger temperature variations with orbital parameters. Indeed, their coldest state for low obliquity and eccentricity has a global mean temperature of 19.5 • C (1.5 • C lower than our value for corresponding orbital parameters), while their warmest state for high obliquity and eccentricity has a temperature of 27 • C values (more than 5 • C higher than the warmest state in our simulations).
Using selected model output kindly provided by David De Vleeschouwer, we are able to investigate the discrepancies between De Vleeschouwer et al. (2014) and our model results in more depth. We find generally good agreement for the median orbit simulations of the two models and for patterns over continents. Differences primarily arise for other orbital configurations and for ocean areas, in particular over the Arctic ocean where we find significant differences in the sea ice distribution between the two studies. For example, De Vleeschouwer et al. (2014) report no sea ice at all in their = 24.5 • , e = 0 simulation, whereas our model simulates sea ice, leading to temperature differences at high latitudes. Although we are not able to fully explain these discrepancies based on the available data and although we cannot rule out a contribution from our simplified atmosphere model, the observed differences in the sea ice distribution make it reasonable to assume that the differences arise from differences in ocean circulation, ocean heat flux, and ocean heat transport and their interplay with sea ice formation and dynamics. Note that, in contrast to our study, the slab ocean model used in De Vleeschouwer et al. (2014) does not explicitly model ocean dynamics and uses the same prescribed ocean heat transport for all orbital configurations. Furthermore, they use different boundary conditions with a higher CO 2 concentration (2,180 vs. 1,500 ppm) and solar constant (1,324 vs. 1,319.1 W/m 2 ) as well as a Late Devonian continental configuration and vegetation cover, resulting in a warmer climate state, leading to differences in sea ice cover, for instance. Test simulations with our model using similar boundary conditions show, however, that the amplitude of temperatures changes between different orbital configurations remains small and that sea ice is present at high latitudes in all cases.

Devonian Climate Variability: Flips Between Two Climate States in the Arctic
While studying the sensitivity of the Devonian climate, we have discovered that the dominant mode of variability on centennial timescales in our model simulations of the Devonian climate relates to quasi-periodic variations in the global mean surface temperature on the order of 0.2 • C with an average periodicity of about 300 years. As an example, Figure 6a shows these variations for a representative time period in one of our sensitivity simulations. These fluctuations in global mean temperature originate from the Arctic region where they are driven by the pace of North Polar SST changes. Indeed, the values of annual mean surface air temperatures averaged over the Arctic between 60 and 90 • N oscillate between about 3.5 and 5.0 • C (Figure 6b). The large regional temperature fluctuations in the Arctic and the waxing and waning of winter sea ice cover with a periodicity of several centuries have the potential to impact Arctic marine ecosystems considerably. For example, Harada (2016) reports significant changes in the biogeochemical cycling due to the anthropogenic climate change in the western Arctic Ocean at all trophic levels. As marine life in the Devonian underwent considerable changes (Copper & Scotese, 2003;Dahl et al., 2010;McGhee, 1996), it is important to further investigate the described mechanism, although verification from proxy data might prove difficult Figure 7. Mechanism describing the flip between the two climate states in the Arctic: In the warm mode, the North Polar water column is destabilized by the northbound transport of warmer and saltier waters to the ocean surface. This inhibits the growth of sea ice and is connected to the establishment of a North Polar deep water convection. The accompanying Arctic meridional overturning circulation transports cold and fresh water northward, leading to the flip into the cold mode. The formation of sea ice stops the deep water convection and the Arctic meridional overturning circulation. The resulting meandering density and ocean surface gradients lead to a northbound salt transport, increasing surface salinity at the North Pole and thereby triggering the flip back into the warm mode.
due to the short timescale and the regional focus of the effect. Furthermore, it is crucial for a meaningful interpretation of the data of this study to take into account these temperature variations. Therefore, in the analysis of our sensitivity studies (section 3) and the best guess simulations (section 5) we evaluate long-term mean values (averaged over 1,000 years) in order to represent an average climate state rather than one of the extreme cases.
The flips between cold and warm climate states in the Arctic can be attributed to coupled changes in sea ice cover (shown in Figure 6b), ocean convection, and ocean overturning (Figure 6c) in the region around the North Pole. A schematic overview of the mechanism is shown in Figure 7 and described in more detail in the following. During the discussion we will often refer to Figure 8 where salinity, northern sea ice cover, ocean surface velocities, and the mixed layer depth are shown for the warm mode, the cold mode, and the two transition states.
During the warm mode (Figures 8a-8c) the situation is characterized by the presence of high-salinity surface water at the North Pole, causing a destabilization of the water column and thus open ocean convection. As a consequence, warmer and saltier water masses are injected into the surface ocean layer while inhibiting the growth of sea ice during the warm mode. Deep ocean convection at the North Pole is, however, intimately connected with the establishment of an Arctic meridional overturning circulation (ArcMOC); see Figures 6c and 9a. The water masses from the surface of the North Pol sink toward the bottom while spreading southward and resurface between 60 • N and 70 • N.
While flowing back at sea surface toward the North Pole, the water masses transported by the ArcMOC deliver relatively fresh waters from the latitudinal band between 70 • N and 80 • N, which promotes the recurrence of sea ice formation and thus the transition from the warm into the cold mode (Figures 8d-8f). As seasonal sea ice grows, deep convection at the North Pole weakens and ultimately ceases because the layer of sea ice inhibits heat losses from sea surface. As a consequence, the ArcMOC collapses (Figures 6c and 9b) and the system lapses into the cold mode.
Within the time interval of the cold mode (Figures 8g-8i), the formation of seasonal sea ice in the Arctic leads to a lowering of the sea surface salinity around the North Pole caused by the impacts of brine rejection during sea ice formation and a freshening of sea surface waters during its melting period (Bouttes et al., 2010). Consequently, the cold mode is characterized by a shallow mixed layer and strong stratification of Polar water masses where cold and relatively fresh water lies above warmer and salty water masses. The surface velocities follow almost a cyclonic flow pattern south of 60 • N and a weaker anticyclonic pattern north of 70 • N, combined with a circularly symmetric shaped dynamical free surface .
However, due to deviations from perfect circular symmetry of water mass properties around the North Pole, the cold mode is unstable, and the system begins the transition back into the warm mode (Figures 8j-8l). As sea ice forms and melts over the years, the relatively fresh surface water layer accumulating around the North Pole spreads slightly asymmetrically toward the south. Therefore, the concomitant formation of a meandering density front within the latitudinal band between 60 • N and 80 • N reshapes the dynamic free surface according to the geostrophical balance (Marshall & Plumb, 2008). The resulting deflection of the Oscillations like the one described here could in principle also result from numerical instabilities caused by a long simulation time step (12 hr by default), the splitting of tracer, and dynamic time steps employed in our ocean model (Montoya et al., 2005) or by rounding errors due to optimization during code compilation. We have performed a large set of test experiments including simulations with short time steps (as short as 1 hr) and without splitting of ocean time steps, experiments on different platforms and with different compiler versions, and optimization settings to verify that the oscillations described above are a real phenomenon rather than a numerical artefact.
Our results suggest that similar modes of climate variability could be a more general feature of climate states with open water at one pole. Indeed, Poulsen and Zhou (2013) describe a climate state-dependent Arctic variability pattern in simulations of the mid-Cretaceous climate. Investigating the influence of different atmospheric CO 2 concentrations, they find a link between temperature variability and changes in sea ice cover for low (preindustrial) CO 2 but additionally observe temperature variability for warmer, ice-free climate states (10× preindustrial CO 2 concentration). For both cases, they describe an intensification of the Northern Hemisphere meridional overturning circulation during the warm mode, similar to the maximum ArcMOC in our simulations, connected to northward transport of lower-salinity and lower-density water and coupled to a strong feedback involving low-level clouds.

Climate States of the Early, Middle, and Late Devonian
After having investigated the sensitivity of the Devonian climate system to a range of parameters, we performed three simulations with the aim to represent the most likely state for the Early, Middle, and Late Devonian by choosing the most appropriate combination of parameter values (Table 1), as explained in section 2.2.6.
The simulated global temperature decreases significantly throughout the Devonian: During the Early Devonian, global annual mean surface air temperature is 22.0 • C, decreasing to 21.1 • C for the Middle Devonian The contributions of the different forcings to the global cooling over the Devonian seen in the set of best guess simulations can be quantified and compared with the help of the sensitivity experiments described in section 3; see the summary shown in Figure 11. The cooling is dominantly driven by the strong decrease in CO 2 concentration (section 3.3) from 2,000 ppm for the Early Devonian to 1,000 ppm for the Late Devonian. The temperature decrease due to changes in the continental configuration (section 3.1) and the temperature rise caused by the increasing solar constant (section 3.2) approximately cancel each other. Finally, the influence of differences in vegetation cover on global mean surface air temperature, as described in section 3.4, can be neglected compared to the temperature changes induced by the other forcings.
While the impact of continental drift on the global climate during the Devonian is relatively small, the effect on continental and regional temperature distributions and ocean surface circulation is significant ( Figure 10). A comparison of temperature patterns and land-ocean masks for the three continental configurations shows that the temperature distribution in the southern high latitudes is largely due to the differences in continental configuration. It is difficult to disentangle the influence of the different parameters on land temperature, but the sensitivity study of the vegetation's influence showed that vegetation has a cooling effect on the continental temperature north of 60 • S which therefore likely contributes to the continental temperature decrease and the continental temperature patterns of the best guess simulations. Finally, the stronger cooling effect of the CO 2 decrease at high northern latitudes observed for the CO 2 sensitivity Figure 11. Simulated temperature change between the Early and the Late Devonian caused by changes in individual boundary conditions described in section 3 and by the combined forcing used in the best guess simulations.
simulations (section 2.2.3) is visible in the temperature maps of the best guess simulations as well and dominates over the warming expected in this region due to the increasing solar constant (section 2.2.2).
Note that we have fixed orbital parameters for the best guess simulations to a median orbit as the main aim of these simulations is the investigation of longer-term trends. This results in an uncertainty of ±0.3 • C which can be derived from the sensitivity experiments with respect to changes in orbital configuration (see section 3.5).
In the following, we want to compare our results with other proxy and modeling studies considering the climate throughout the Devonian. In their modeling study investigating the impact of orbital forcing on the Devonian climate (see the discussion in section 3.5) De Vleeschouwer et al. (2014) provide maps of December, January, and February and June, July, and August surface temperatures for a Late Devonian median orbit ( = 23.5 • and e = 0) climate simulation. Keeping in mind that they use higher values for the CO 2 concentration and the solar constant, the seasonal temperature patterns agree well with our Late Devonian simulation. Furthermore, the patterns for the difference of evaporation and precipitation in our simulations (see Figure S7) match the model results and lithic indicators of palaeoclimate presented in De Vleeschouwer et al. (2014) reasonably well. We can therefore conclude that our model exhibits a good latitudinal representation of the major climate zones during the Devonian. Simon et al. (2007) investigate the evolution of the atmospheric CO 2 concentration during the Devonian using a global biogeochemical model coupled to an energy balance climate model, taking into account the changes in continental configuration, sea level changes, and the increase of vegetation covered by land plants. The modeled CO 2 concentration of 3,000 ppm for the Early Devonian is very high (see Figure 1) and shows a general decreasing trend throughout the Devonian going along with the evolution of land plants, reaching low concentrations of 1,000 ppm already around 390 Ma and then fluctuating around this value during the following 30 million years. For the three Devonian periods covered by our best guess simulations, the low-latitude temperatures of Simon et al. (2007) compare well with their equivalent in our simulations although both studies are based on very different modeling approaches.
Considering proxy data, several studies determine tropical SSTs using 18 O. Van Geldern et al. 2006 use 18 O from brachiopod shells to derive SSTs ranging from 24 to 27 • C in tropical to subtropical latitudes for the Early and Middle Devonian. For the Late Devonian, they discuss that their temperature estimates of 31 to 41 • C are significantly too warm for providing good living conditions for the diverse marine life found in their investigated successions, possibly suggesting an increase in 18 O of seawater (Jaffrés et al., 2007) from the Devonian to present-day. Joachimski et al. (2009) point out difficulties in the interpretation of 18 O brachiopod shell signals and provide a 18 O record derived from conodonts (see Figure 1b) which they argue to be a more reliable palaeotemperature record. In Figure 12 we compare their tropical SST estimates for different calibration standards and temperature equations with the SSTs from our best guess simulations. We averaged our SSTs between 10 • S to 30 • S since all the sections used for the temperature reconstruction in Joachimski et al. (2009) are located at palaeolatitudes in this range. Our tropical SSTs decrease from 28.9 • C for the Early Devonian to 28.3 • C for the Middle Devonian to 27.1 • C for the Late Devonian. For the Early Devonian, the new SST curve (M. Joachimski, private communication, January, 30th, 2018) calculated using the new calibration standard and temperature curve (Lécuyer et al., 2003(Lécuyer et al., , 2013 shows warm SSTs of 35 to 37 • C, with a cooling trend, continuing in the Middle Devonian, to temperatures of 23 to 27 • C. This is followed by an increase to 36 • C in the early Late Devonian and a slight decrease to 30 • C toward the end of the Devonian. We note that the large spread of the data as well as the model uncertainty make the comparison of the model results with the temperatures from 18 O difficult. Nevertheless, except for the Middle Devonian, we can state that the modeled temperatures from our best guess simulations shown in Figure 12 generally compare better with the lower SSTs values given in Joachimski et al. (2009) rather than the SSTs calculated using the new calibration standard and temperature curve. Together with the fact that these very high SSTs are  (Lécuyer et al., 2003) rather than 22.6 (Vennemann et al., 2002) and assume a different temperature equation (Lécuyer et al., 2013). The solid red and dashed blue lines show local regression fits to the data points. The filled black circles are modeled SSTs from our best guess simulations averaged from 10 • S to 30 • S since all the sections used for the temperature reconstruction in Joachimski et al. (2009) are located at palaeolatitudes in this range. The error bar indicates the temperature range from 10 • S to 30 • S for each time slice. Note that the additional error due to uncertainties in the boundary conditions is much smaller than this temperature range; in the sensitivity experiments, SSTs at 10 • S and 30 • S vary by ±1 • C at most for the range of CO 2 concentrations, orbital parameters, and vegetation distributions tested. Additional seasonal temperature variations (not shown) would mostly affect the lower end of the indicated temperature ranges; typical seasonal variations are roughly ±2 • C at 30 • S and ±0.5 • C at 10 • S. SST = sea surface temperature.
beyond the mortality limit for diverse ecosystems of higher marine life (van Geldern et al., 2006), this could indicate lower 18 O levels of Devonian ocean water as compared to the present-day ocean (van Geldern et al., 2006;Jaffrés et al., 2007).
Considering the temperature trend, our model results, like other modeling studies (Simon et al., 2007), cannot reproduce the higher temperatures of the late Frasnian and early Famennian compared to the Middle Devonian, suggested by SST proxy data (Joachimski et al., 2009;van Geldern et al., 2006) and palaeontological studies (Streel et al., 2000). This could suggest missing forcings in the models, for example, the forcing induced by a change in O 2 concentration or an increase in the CO 2 concentration not resolved in the CO 2 proxy compilation.
Despite important differences, the general agreement of our results with the modeling studies of De Vleeschouwer et al. (2014) and Simon et al. (2007) is encouraging: We find similar temperature patterns for the median orbit simulations of De Vleeschouwer et al. (2014) and the low-latitude temperatures of Simon et al. (2007) compare well with ours. It is important, however, to keep the limitations of our model in mind. First, our it does not model vegetation dynamically. For this reason, climatic limits for the Devonian vegetation distributions had to be imposed from a simulation without vegetation cover (see section 2.2.4). We would argue, however, that our modeling approach is a reasonable approximation and complementary to studies with dynamic vegetation models since the limited knowledge of Devonian plant and soil properties could make the accurate representation of Devonian vegetation dynamics in these models difficult. Second, the increase in weathering due to the spread of land plants and the development of deeper roots are not modeled explicitly but only by imposing a reduction of CO 2 during the Devonian. Hence, the climate weathering feedback is not taken into account. The main effect, however, is captured by the decreasing atmospheric CO 2 concentration throughout the Devonian; hence, this seems to be a reasonable approach. Third, we employ a coupled climate model with a simplified, coarse-resolution atmosphere model. However, there is good agreement between our model and models with a more sophisticated atmosphere in a number of studies for a wide range of different time periods, for example, for the Neoproterozoic (Liu et al., 2013;Feulner & Kienert, 2014), for the last millennium (Feulner, 2011;Schurer et al., 2014), or for the 21st century (Feulner & Rahmstorf, 2010;Meehl et al., 2013). Furthermore, Ganopolski et al. (2001) investigate the effect of vegetation surface cover modeled by the coarser-resolution model CLIMBER-2, but with the same surface scheme as CLIMBER-3 , and find good agreement with results from more sophisticated models. Therefore, we do not consider the atmosphere model to be a serious limitation of this study.

Conclusions
In an attempt to disentangle the various factors influencing Earth's climate during the Devonian, we have systematically performed and analyzed a large set of sensitivity simulations with a coupled climate model, focussing in particular on the impact of changes in continental configuration, the biogeophysical effects of increasing vegetation cover, and variations related to Earth's orbital parameters. The key findings of our paper can be summarized as follows: 1. The impact of changes in orbital parameters (at a fixed CO 2 concentration) on the Devonian climate is in line with earlier findings in the sense that annual global mean temperature increases with obliquity and eccentricity. However, the differences in global mean temperature are much smaller in our simulations using an ocean general circulation model (rather than a simple slab ocean model with fixed heat transport) because of changes in ocean heat transport counteracting the insolation changes.
2. The most important mode of Devonian climate variability on centennial timescales relates to coupled changes in Arctic meridional ocean overturning, northern sea ice cover and deep ocean convection around the North Pole. Our results and a similar variability pattern described for the mid-Cretaceous (Poulsen & Zhou, 2013) suggest that this mode of climate variability is a more general feature of climate states with open water at one pole. 3. Best guess simulations for three time slices approximately representing the Early, Middle, and Late Devonian based on estimates of continental configuration, vegetation cover, solar constant, and atmospheric carbon dioxide concentration show a general climatic cooling trend in accordance with proxy estimates. This cooling is predominantly driven by decreasing levels of atmospheric CO 2 as derived from reconstructions. Furthermore, the biogeophysical impact of an increasing vegetation cover on the Devonian climate is small. In particular, the associated albedo changes cannot completely compensate for the falling CO 2 concentrations even under the assumption of a relatively high bare-soil albedo. The increase in temperatures around 390 Ma observed in reconstructions therefore remains difficult to reconcile with reconstructions indicating falling CO 2 levels. 4. Finally, simulated tropical SSTs for our Devonian time slices are generally significantly lower than the estimates based on 18 O and the latest temperature calibration. This discrepancy and the fact that reconstructed temperatures are at times beyond the lethal limit for higher life (van Geldern et al., 2006) could indicate a general difference in the oxygen isotope ratio between the Devonian and modern oceans.
Our results on the Devonian climate highlight that much remains to be done to improve our understanding of this crucial period in Earth's history. Open issues include a better quantification of the biogeochemical impact of the spread of vascular land plants, the puzzling warming at the beginning of the Late Devonian despite falling atmospheric carbon dioxide levels, the climate response of an increase in the atmospheric oxygen concentration during the Devonian, and the interplay between terrestrial changes and the marine extinction events.

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 paper are available online (http://www.pik-potsdam.de/data/doi/10. 5880/PIK.2019.002/). 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.