Low atmospheric CO2 levels before the rise of forested ecosystems

The emergence of forests on Earth (~385 million years ago, Ma)1 has been linked to an order-of-magnitude decline in atmospheric CO2 levels and global climatic cooling by altering continental weathering processes, but observational constraints on atmospheric CO2 before the rise of forests carry large, often unbound, uncertainties. Here, we calibrate a mechanistic model for gas exchange in modern lycophytes and constrain atmospheric CO2 levels 410–380 Ma from related fossilized plants with bound uncertainties of approximately ±100 ppm (1 sd). We find that the atmosphere contained ~525–715 ppm CO2 before continents were afforested, and that Earth was partially glaciated according to a palaeoclimate model. A process-driven biogeochemical model (COPSE) shows the appearance of trees with deep roots did not dramatically enhance atmospheric CO2 removal. Rather, shallow-rooted vascular ecosystems could have simultaneously caused abrupt atmospheric oxygenation and climatic cooling long before the rise of forests, although earlier CO2 levels are still unknown.

Supplement S2.2 (p.11, l 207 ff.): there is no need to solve equations 1 to 3 numerically; equating equations 1 and 2 leads to a quadratic equation in c_a which can easily be solved by paper and pencil methods.
Reviewer #2 (Remarks to the Author): Tais et al propose that atmospheric CO2 in the Devonian between 410 and 380 Mya ranged between 460 and 660 ppm. These new paleo-CO2 estimates are between 4 to nearly 10 times lower that the prevailing estimates based on multiple proxies. Based on the revised paleo CO2 estimates Tais et al go on to overturn the prevailing paradigm that the emergence of forested ecosystems in the Devonian and the acquisition of evolutionary adaptations for life on a terrestrial planet ( leaves, roots, arborescence etc etc) enhanced chemical weathering of silicate rocks thus leading to enhanced sequestration of carbon as carbonate in the ocean system. In order to support this paradigm shift the authors use the COPSE model output to reinforce the idea that the evolution of deep rooted terrestrial ecosystems did not enhance silicate weathering but they offer no mechanistic basis or empirical data set for their contrary hypothesis. Science always builds on previous scientific endeavors and sometimes it smashes the ideas and studies that have gone before it however if a study, such as that presented by Tais et al is proposing a complete paradigm shift in understanding on biosphere -atmosphere interaction in the Paleozoic then it must hold up to scrutiny and have a firm empirical/observational/ logical/modelling foundation. My appraisal of this manuscript indicate that the provocative and potentially exciting and paradigm shifting conclusion of this study are not sufficiently supported by the data, analysis or model revisions presented. 1) From my understanding of the Supplementary files, the new paleo-CO2 estimates are based on the analysis of 4 rock specimens, on which 25 fragmented fossil plant specimens were analyzed for stable carbon isotopic composition. Although the authors argue for the importance of sourcing stomatal data and carbon isotopic composition from the same plants to parameterize the mechanistic stomatal proxy model used in this study -this was not undertaken in the study. Instead, stomatal data, including stomatal density, pore length and depth were taken from different specimens of the same species growing in different locations and from my reading of the supplementary files (although it is not at all clear) from the literature. In addition the preservation of the fossil specimens indicates some diagenesis. I acknowledge that the authors have taken many steps to argue for the pristine nature of the carbon isotopic signal from the 25 fossil fragments however it would have been good to provide further supporting data from charcoalified or purely coalified preserved specimens to bolster their arguments. Parameterization of the mechanistic model with stomatal and isotopic data is therefore unsatisfactory.
2) In addition to the initial fossil stomatal and isotopic dataset being unsatisfactory, the authors have completely re-calibrated the original mechanistic stomatal proxy model of Frank et al. using two tropical epiphytic lycophytes. The choice of tropical epiphytes, which grow on other plants, usually woody trees rather than terrestrial lycophytes to recalibrate the mechanistic proxy model is not sufficiently justified by the authors. They argue throughout the manuscript that a focus on lycophytes rather than extinct taxa is preferable to historical studies which took an analogous trait (plants from similar environments and with similar morphology likely function similarly) rather than homologous trait approach (plants which are related likely function similarly). Of course when dealing with nearly 400 million years of evolution separating the fossils from the modern taxa used to calibrate them neither approach is 100% satisfactory and many unknowns remain. In order for a complete recalibration of a published model to stand up to scrutiny the authors should have undertaken a verification in both ambient AND elevated CO2 conditions. The others only present a test of their recalibration of the Frank model from a few plants grown in a greenhouse under ambient CO2. An elevated CO2 calibration was not undertaken and it is a flaw in the study design. The re-calibration of an existing model is not therefore sufficiently supported with the test examples currently presented.
3) There are strong hints that the new epiphyte lycophyte calibration of the Franks model does not work on lines 221 -224 of the manuscript. The revised model estimates assimilation rates for a number of charismatic early Devonian fossils plants that suggest they are barely autotrophic without any discussion or citation of existing literature arguing to the contrary. Are the authors suggesting that these fossils are almost non-photosynthetic non independent sporophytes living on the gametophytes with such low A0 values or could it be that the recalibration of the model requires further work using a lycophyte which is rooted in the ground or at least not an epiphyte. Have the authors considered using additional lycophytes like Selaginalla in their re-calibration of the mechanistic stomatal model? Selaginalle have double the observed stomatal densities of Huperzia. How would this influence the paleoCO2 estimates of the Devonian lycophytes (I believe it would double their estimates in line with most current CO2 reconstructions for this interval). This demonstrates the scale of uncertainty in estimating paleoCO2 from fossil stomata for the earliest land plants for extinct and extant lineages. The choice of a living taxon to recalibrate the Franks model will drastically alter paleo-CO2 estimates from fossils. The choice of taxon must be robustly supported. 4) Throughout the manuscript the authors are bombastically confident about both the precision and accuracy of the new paleo-CO2 estimates however they have not considered unknown errors that cannot be propagated or the fact that other studies show non -linear errors under different CO2 concentrations. See for example tests of the mechanistic model in recent papers by Porter which demonstrate that relatively small errors under ambient CO2 but very significant underestimates and thus large errors ( 500 to 1000 ppm) under elevated CO2 conditions. Incidentally, the manuscript does not consider experimental findings or suggested correction factors for lycophytes and ferns in the Porter papers. Why were lycophyte/ phylogenetic correction factors to the mechanistic model not considered? 5) In general too much material and methods is buried in the supplementary files which are a little chaotic. They jump between materials and methods, personal notes remain in some of the figure legends and the referencing section for ref 1 to 10 do not seem to align with the supplementary text. It is particularly difficult to work out who collected the stomatal data, whether the data is a meta-analysis from the literature. Inclusion of scaled micrographs with examples of the pore width and length measurements would be a good addition. 6) In parts of the manuscript the authors selectively cite the literature and facts to support their arguments and this has resulted in a number of contradictions in the logic and framing. For example the authors argue that because the stomata of fossil and modern lycophytes look similar then the fossils must have functioned the same as their modern relatives. Yet the authors reject similar arguments made in relation to extinct Devonian taxa-they suggest that even though the have similar morphology to some living species because they are unrelated thet they cannot share ecophysiological likeness or functioned the same. Generally, the authors present in black and white and miss much of the nuance, weaknesses and strengths of past studies. As an example, there are multiple references to the productivity of the paleoecosystems throughout sections of the paper without any explicit data or tests on paleo productivity or citations. These in turn weaken the arguments on paleo weathering.
In conclusion, the authors present very interesting data and analysis that is preliminary in nature only. The new dataset and estimates of paleo-CO2 suggest that our understanding of paleozoic atmospheric evolution may be incomplete or event incorrect, however the preliminary nature of the data sets and remodeling presented by Tais et al do not provide sufficient confidence to completely overturn all existing concepts, hypotheses and previous paleoCO2 estimates in relation to plant-atmospheres and weathering etc etc in the Paleozoic.
Reviewer #3 (Remarks to the Author): In this paper, Dahl and colleagues refine a mechanistic gas exchange model for modern lycophytes and apply this information to closely-related fossilised lycophytes in order to constrain atmospheric CO2 concentrations from 410 -380 Ma. They then use this information in a paleoclimate model to derive surface temperatures, snow coverage and precipitation for a given CO2 concentration and various orbital configurations. Finally, they use the biogeochemical COPSE model to explore possible reasons behind changing CO2 concentrations across the Ordovician to early Carboniferous, and try to match their data from the mechanistic model for the early to mid-Devonian.
It is clear that a huge amount of work has gone into this paper, which is greatly appreciated. There are a distinct lack of geochemical proxies for CO2 in the early to mid Paleozoic, so the work here on developing the mechanistic model and using relationships between modern and ancient fossilised lycophytes to generate CO2 values is invaluable. The paleoclimate modelling is also pretty neat. In general I agree with what the authors conclude with regards to shallow rooted vascular systems causing significant changes to the Earth system. However, I would note that this has been predicted before. For example COPSE Reloaded predicts a first decline in CO2 and rise in O2 in the Ordovician, as do Krause et al. (2018) for O2 well before the first forests and even before the first fossil evidence for vascular plants.
For the paper as a whole, I have a few questions and comments which should be easily taken care of (see below), but I have some serious concerns with regards to the COPSE modelling section. I have tried my best to outline these clearly below.

COPSE section (main text and SI)
First of all, it would be useful to describe in a bit more detail the differences between COPSE Reloaded (CR) and COPSE 2016 (C16). Why is it important that the authors change forcings in C16 and show what they do to CO2 and O2 predictions, as well as in CR? For people not particularly familiar with the model they are going to question the utility of doing this because they do not know the scale of differences between the iterations of the model. However, from looking at this figure, I'm wondering, are the authors combining the results from their revisions to the two different COPSE iterations? Or did they make further changes to CR and these ended up in Figure 3, but not in S18b? Because looking at the SI Figures, S17 and S18, the pCO2 from the COPSE 2016 (C16) revision (S17b) is above 5000 ppm at 480 Ma, whereas the revised CR is 4000 ppm, yet in Figure 3 the main pCO2 evolution curve is less than 4000 ppm at 480 Ma. The overall trend of the curve in Figure 3 looks slightly different to those in S17b and S18b and doesn't come as close to the data. The pO2 evolution curves are also very different. The revised C16 shows ~5% at 480 Ma, the revised CR at ~10% but Figure 3 shows it to be closer to 15% atm at that time. C16 is consistently above 20% atm within the fire window, while CR is barely within it, yet Figure 3 predicts pO2 to generally be around 20% atm within the window.
Furthermore, the pO2 curve from the original COPSE Reloaded model that is shown in Figure  3 does not look right at all here. For a start it should be around 5% atm from ~460 to 480 Ma. It should dip below the wildfire threshold during the Devonian (at ~406 Ma) rising above it again at ~382 Ma. The curve in Figure 3 does not match that in Figure S18a, and that curve in S18a doesn't look right compared to the output from the original model either (see figure below, which for O2 is figure 13d in Lenton et al. (2018)). The original COPSE 2016 model pO2 curve (S17a) also does not look right when compared to any of the modelled scenarios presented in the Lenton et al. (2016) paper.
And while there is no such discrepancy between Figure 3 and Figure S18a for the CO2 curve, again when compared to the output from the original CR model (below and/or figure 13a in Lenton et al. (2018)) it is very different. The original model predicts CO2 in line with the authors' final data point (at ~380 Ma) though admittedly it is higher than the data prior to that point, but not as high as suggested by the authors' figures.
These points are important because if the authors ran the original CR model in order to plot a curve for the paper, it would seem that they may have run a version in a mode with parameters set to explore some underlying uncertainties but not the mode which produces the best guess (i.e. the main black lines in Fig It is not surprising then, that the CO2 and O2 curves for the original C16 and CR are not the same as those in the respective papers. The authors are proposing changes to forcings such as degassing and plant evolution to enable a decline in CO2. But as their version of the original CR has relatively high CO2 at 480 Ma and only gradually declines through the timeframe of interest (through to 320 Masee Fig. S18a) then they probably needed to make more dramatic changes to the forcings than is necessarily needed, considering that the actual CR forcings are different and thus CO2 output has a 2-step decline in CO2. It also means that the decline in CO2 initiated by the first vascular plants is not as dramatic (~2000 to ~500 ppm, instead of ~3800 to ~500 ppm). Indeed, other models such as GEOCARBSULF (Berner, 2008;Royer et al., 2014) predict CO2 levels of ~2800 ppm at 440 Ma.
This leads me on to more questions about changes to the forcings. The plant enhancement of weathering (W) is now assumed to scale with changes to mudrock proportions, based on McMahon and Davies (2018), but there's no mention of how that is done in the discussion, methods, or SI sections. In McMahon and Davies, in the main text at least, they only present data to the end of the Carboniferous, so is it assumed that the end Carboniferous value is the one by which normalisation occurs?
For degassing, it reads as if the authors have decreased degassing from 1.5x to 0.9x Present so that the CO2 prediction then better fits their data, and then cite it being plausible due to atmospheric CO2 input being primarily sourced from continental arc volcanism. Is this the case? Or have the authors come up with a new degassing forcing by somehow converting the zircon abundance/ages in McKenzie et al. (2016) into a normalised forcing? While I agree that degassing should not be solely attributed to ocean spreading rates as has been the case in CR (2018) and other such models, I am unsure as to whether a reduction in continental arc volcanism would result in such a big decline in degassing rates. For the last 200 million years at least, the range in CO2 emitted from arc volcanoes has equalled that from mid-ocean ridges and are dwarfed by continental rifting (see Wong et al. (2019) in Frontiers in Earth Science) and this may have been the case further back in time. So even with a decline in arc activity, this would need to be mixed in with a modestly declining ocean spreading rate, which would probably result in a smaller overall decline than 1.5x to 0.9x.
The same with the other forcings that are changed, this bit needs to be fleshed out in the Methods or SI. And again, this leads me on to another line of thought: although it's shown in Figure 3 what a different degassing rate does to the CO2 predictions, which of the other forcings that have been changed (W, E, F, CP) are the main driver in producing the solid blue line in Figure 3? Is it changes to plant evolution, or are they all contributing equally? This is unknown but should be explored here.
In the text for the SI (lines 623-627) the authors state that selective P weathering (forcing: F) is omitted for simplicity and because it is not necessary, however, they mention in the main text (line 406) that this forcing is changed, and it is plotted in S18c and d implying that it is still in use in the modelso which is it? In fact, it is plotted in S17c and d as well, but called epsilon instead of F in the legend.
Why are the degassing and uplift forcings not plotted in S17c and d? And why is plant evolution called EVO in S17 but E in S18? Why are the normalised weathering fluxes dashed lines in S18c and d but the legend has them as solid lines?
For Figures 3, S17 and S18 what reference does the red circle for pO2 come from? As far as I can remember the inertinite data from Glasspool and Scott (2010) and Glasspool et al. (2015) only extends back to 400 Ma. Scott and Glasspool (2006) do note the first evidence for fire towards the end of the Pridoli Epoch (~420 Ma) but don't ascribe a pO2 value to that, but the fire window is based off of that first evidence in combination with the experimental work of Belcher and McElwain (2008). All this is to say, have I missed some new data about this, or is it a combination of the above? Some references are needed here.
Also, if you want some other pO2 proxies to plot in Figures 3, S17 and S18 then Sønderholm and Bjerrum (2021) have some estimates from fossil roots that can be plotted. . Is there any reason to have the CO2 outputs grey instead of blue as in main text Figure 3, and also why are the CO2 proxies from paired stomata and isotope data in S17 and S18 not coloured the same as Figure 3?
Figures S17 and S18 are quite low resolution compared to the equivalent figure in the main text, this could be improved because the main text mentions that degassing is decreased from 1.5x to 0.9x but in Figure S18 it looks like a drop from 1.5x to 1x with a minor rise again at ~400 Ma to ~1.1x.

Other comments and questions
Was the lighting in the greenhouse purely natural sunlight, or also aided by artificial light? Were the lighting conditions optimal for the species in question? There is some evidence (though I am unsure how robust) that light spectral quality can affect photosynthetic activity and carbon isotope fractionation (e.g. Tarakanov et al. (2022)).
The stomatal density used in section 2.5 in the SI is 28 +/-5 which appears to be based primarily on H. phlegmaria at an RH of ~80%? But for H. squarrosa it appears that changes in relative humidity can affect stomatal density. Might this also be the same for H. phlegmaria? How does this then affect the carbon-based results if CLIMBER predicts a RH of ~65% for the Baragwanathia flora at the Yea fossil site in Victoria, Australia?
How do the average surface temperatures from the climate model (Table S4) compare to sea surface temperatures from oxygen isotopes?
Lines 274-275: I think this sentence could be made tighter by saying it's post vascular colonisation, because the unforested Earth could mean anytime in the previous ~4 billion years.
Line 173: I'm not sure if underpins is the correct word here. The new data don't really underpin the model predictions, but are used as a comparison, aren't they? You're not using the data to force the model. lineage of lycophytes (Order: Selaginales) that evolved later, and further S. Kraussiana is anatomically much smaller than the Early Devonian lycophytes Drepanophycus, Asteroxylon and Baragwanathia.
We agree with the reviewer that the choice of taxon for calibration is critical. We have chosen the most logical and most appropriate set of extant taxa possible.
We thank the reviewer for pointing attention to the experiments with S. Kraussiana grown in controlled Conviron chambers at high pCO2 (1899±60 ppm) for which there is stomatal and isotope data at various O2 and CO2 levels (Porter et al. 2017(Porter et al. , 2019. In these experiments, the stomatal density and size of fresh microphyll from S. Kraussiana (and many other plant lineages) show no significant difference at elevated pCO2 to that of equivalent plants grown at 400 ppm. Yet, their carbon isotope composition did change, and applying the Franks proxy shows that they adapted and operated at a faster CO2 assimilation rate, which is qualitatively correct. Yet, the response of S. Kraussiana deviated quantitatively from that calculated by the Franks equations using the observed stomatal size and density.
There are several reasons why this does not affect our interpretation.
1) The two species of Huperzias that we used for the proxy calibration are physiologically more similar to Drepanophycus, Asteroxylon and Baragwanathia. Huperzias belong to the Order Lycopodiales, where Selaginella belong to a sister group (Order: Selaginales) that evolved later (Late Devonian). Thus, the Early Devonian lycophytes belong to the Drepanophycales that might be closer related to the Lycopodiales. ). This suggest that Selaginella plants studied to date must have other ways to cope with gas exchange. 3) For example, in contrast to other lycophytes, Selaginella has ligules and their function is not accounted for in the Franks model. Likely, ligules prevent water loss from the microphyll and can allow for a higher stomata density. Therefore, Selaginella are likely not passively adjusting their stomata to ambient pCO2 level in the simple manner that more basal lycophytes do, but gas exchange is also affected by ligule function. 4) Although this is beside the point of our paper, we do raise concerns regarding the interpretation of experiments under Conviron conditions. In Franks model, it is critical that stomata are either fully open or fully closed, and that the measured pore size on a harvested plant corresponds to that of stomata in the living leaves. Alternatively, one could adopt a "stomata size correction term". We find that if Selaginella had a factor of 2-4 larger effective stomata size in the Conviron chambers than in the harvested leaves, then that would be sufficient to correct for pCO2-proxy overprediction in the experiment (see Table below).

Table. Predicted pCO2 in controlled experiments with the lycophyte Selaginella kraussiana grown at 1899±60 ppm in Conviron chambers if a stomata size correction factor of 2-4 was applied.
stomata size correction factor Finally, -and most importantly -our paper demonstrates that Devonian lycophytes have low SD, small stomata size, and small carbon isotope fractionation similar to that of equivalent modern plants. Experiments with lycophytes show that they respond qualitatively to high pCO2, but the magnitude of the change deviates for S. kraussiana some experiments performed in Conviron settings. Regardless of how the Drepanophycales responded quantitatively if grown under high ambient pCO2, the observed stomata density, size and carbon isotope fractionation is compatible only with low pCO2 similar to that of our calibration. Therefore, the uncertainty of the calibration at high pCO2 does not affect our main conclusion, and we must still reject the current paradigm and rule out 4-10-fold higher pCO2 levels in the Early Devonian compared to today.

Figure S22. COPSE reloaded simulations (as in Fig. S19) with revised forcing of Weathering (W), Plant evolution (E), selective weathering (F), and volcanic outgassing (D) scaled linearly to Continental Volcanic Arc (CVA) emission, which in turn is scaled in proportion to Young vs. Old grains in sedimentary deposits (McKenzie et al. 2019). The baseline CR model is shown in dotted curves in all panels.
Previous version (forcings are not updated to CR baseline)

Redacted
Editorial Note: Parts of this Peer Review File have been redacted as indicated to remove third-party material where no permission to publish could be obtained.

Our palaeoclimate model predicts Sea Surface Temperatures (SST) in the Devonian oceans that can be compared to proxy data obtained from the oxygen isotope composition of marine brachiopods and conodonts, provided we know the 18 O composition of the Devonian oceans in which the shells formed.
The tabulated data is summarized in Data S2.
We compiled 18 Ocalcite data from well-preserved (e.g. non-luminescent shells)  Overall, the manuscript presented by Dahl et al is much improved in all aspects of the analysis, discussion and presentation. The authors have undertaken a huge amount of work and the results are paradigm changing and will be of huge interest to the community. I remain concerned with the estimated CO2 values in relation to full error propagation and treatment of some of the scaling factors used within the CO2 model (as outlined below). I want to reiterate that I support this work but as a reviewer want to be thorough and make sure that the authors work stands up to scrutiny and reanalysis. If the authors could provide a simple table for all of their input values (as in Table 1 below) in the model following that specified by Franks for each species/sample investigated so that their entire analysis could be repeated by others then I think this research and manuscript is worthy of publication. As it stands there remain a number of more minor concerns (listed below) all of which I believe could be addressed and revised as shown below.
Tabel 1 Dab, eDab  Dad  eDad  GCLab  eGCLab  GCLad  eGCLad  GCWab  eGCWab  GCWad  eGCWad  d13Cp  ed13Cp  d13Ca  ed13Ca  CO2_0  A0  eA0  gb  egb  s1  es1  s2  es2  s3  es3  s4  es4 s5 es5 In Supplementary: (1) Here, the operational conductance is to a good approximation a fixed proportion ( conductance (full daylight) and no conductance (night time), g c(op) = ·g c(max) . -Currently is unknown for Lycophytes. The authors have chosen to select 20%. What this means in plain English is that the authors assume that ancient lycophytes operated at 20% of their maximum theoretical stomatal conductance which is a function of the stomatal pore size and geometry. Full propagation of errors would include a wider range of values or at least more fully justify their choice -(2) we found A 0 values of 2.17 and 3.5 for H. squarrosa and H. phlegmaria, respectively -How were the photosynthetic rates measured for the greenhouse specimens? Were they measured under ambient CO2 of the greenhouse or 400 ppm CO2.
(3) Full error propagation of errors would include modelling multiple values as shown below eg modelling a circle ( = 1) as well as varying proportions of a circle ( = .8, .6 etc) for the geometry of the stomatal pore.
From reading the supplementary materials I think the authors used a value of 0.6 but they state they have used a value of 0.2. This may be my error however please could the authors clarify which value they used in the model and provide full justification From running the Franks model wit the inputs above and an An of 3 to test the authors results it looks like the model is very sensitive to values. It would be extremely helpful if the authors could report the input paameters in a table format following the same method as Franks. This would allow others to run the model with all possible iterations of the input variables and repeat the authors work. This is currently not possible. ) and stomatal pore length (p = 20±3 µm).
(7) Figure S8 for observational data (H leaf , SD, p; upper row) to calculate probability distributions for model parameters Please provide sample number for original observational data in table format as for Baragwanthia above. It is still not possible to work out how many observations of fossil stomatal density and pore length and depth were made on fossil Asteroxylon. Were the measurements made by this author team or are they from the literature? If they are from the literature please provide citation

(8) S2.5 Sensitivity analyses
It looks like no sensitivity analysis was conducted on pore length measurments or modelled pore area usin different geometries for stomatal opening or for different values -see previous comment (9) operational conductance efficiency the use of this term throughout the manuscript is incorrect it is operational conductance (10) Four parameters can produce significantly (>10%) higher pCO 2 estimates, including lower SD (and higher p), wer operational conductance efficiency (~ < 0.20), and higher CO 2 assimilation rate at reference CO 2 concentration (A 0 > 3.5 qmol m -2 s -1 ) in the ancient flora compared to modern lycophytes.
Lower SD and lower, not higher P would be expected to produce higher CO2 estimates? but in the table below this paragraph (s3) the alternate values include a substantially higher SD and unsurprisingly this leads to a lower CO2. Could the authors clarify how they did the sensitivity analysis and what the justification for more than doubling SD as an alternate scenario? All other alternate factors used seem reasonable and clear.
Most Devonian fossils with carbon isotopic data do not simultaneously preserve stomatal data, but the stomata of Drepanophycus spinaeformis from Emsian, Givetian and Frasnian deposits carry stomata with similar size and density as modern plants, including SD = 18±4 mm -2 and p = 17±3 µm. provide citation if data is from meta analysis or specimen numbers if observed in this study In Table S3 and alternate value of 89.5±13.5 is used for SD but as reported above alternate SD values for other lyvophyte like taxa indicate SD's of = 18±4 and as low as 11.5±1.5 mm -2

In main manuscript
Errors I roughly recalculated CO2 using Franks with A0 value of 3, and using the median values for all other inputs reported in Fig 2 for Drepanophycus and estimated the following CO2 with 16% and 85% errors. I acknowledge that this is roughly done calculation but what strikes me is the very narrow error range reported by the authors as quartiles compared to the error range in the output from the Franks model below. It would be good if the authors could clarify why and how their errors are so small compared to all other published error output from the Franks model to date. I appreciate the extra work undertaken to generate the temperature records from oxygen isotopes with which the Climber model results are compared to, this strengthens the paper considerably.
Revising the paper to just use COPSE Reloaded (thus omitting COPSE 2016) tightens up the paper and improves its clarity, and showing the effect of the iterative changes on the model predictions is very nice. However, I do still have some concerns with the COPSE section. The CO2 and O2 predictions for the baseline of COPSE Reloaded (CR2018) as plotted in Figure  Here, CO2 and O2 levels start to change at around 460 Ma. However in Figure 3, the CR2018 baselines (dashed) as plotted do not begin to change until just before 440 Ma:

Redacted
Editorial Note: Parts of this Peer Review File have been redacted as indicated to remove third-party material where no permission to publish could be obtained.
This paper has shown a massive transformation since initial submission. It is superb and will be paradigm shifting. Dahl et al present robust and groundbreaking results on Devonian pCO2 estimates using a mechanistic stomatal proxy model. The authors report much lower pCO2 values than previously estimated and published by any proxy method. In support of these new low paleo CO2 estimates the author team have provided multiple lines of independent evidence that now strongly bolsters their bold claims. I am thoroughly convinced that the newly reported CO2 estimates and the implications for earth system processes (in particular biosphere-atmosphere interactions) will hold up to scrutiny . The paper now reads exceptionally well. My congratulations to the author team. I have only two minor comments -I agree that the commentary on the life cycle implications of this work for Aglaophyton and Rhynia etc should be kept in the paper as currently written.
I have one very minor suggestion -line 141 -I would add ' presumed' or assumed in front of 'physiologically similar. I strongly support acceptance in its current form.
Reviewer #3 (Remarks to the Author): I would like to thank the authors for their comprehensive response to my concerns and those of the other reviewers. I have a few additional minor concerns which need to be addressed and are listed below. Once these have been completed, I am happy to recommend it for publication and look forward to reading the final version.
-Alex Krause L38-41: I think the wording of this sentence is a bit confusing at the beginning. I'm not sure if you mean something like, 'Enhanced chemical weathering is suggested to have caused a decline in atmospheric CO2 pressure…'? The way it is written at the moment it seems like you're saying silicate weathering levels were both enhanced and declined during the Devonian-Carboniferous transition.
L218-219: Needs a reference. Figure 3: The wildfire minimum line is missing. The previous version of this figure showed the goethite and phytane estimates. I thought this was useful to see, but is not necessary, so I will leave it to the authors to decide as to whether to include them again when revising the figure to add the wildfire line.
In the SI: L964: C16 needs to be defined.
L968: The CO2 range has formatted to time and date. Figure S20 (and Fig. S21) and L993-1002: Is the revised E forcing plotted correctly? It appears that you're ramping up the forcing later (~445 Ma) than in COPSE Reloaded (CR, ~465 Ma) but state in the text the opposite occurs. Is it the case that you're ramping up at the same time as CR but to 0.5 instead of 0.25? Or are you starting to ramp up to 0.5 at 485 Ma?