Abstract
Ecosystem process models provide unique insight into terrestrial ecosystems by employing a modern understanding of ecophysiological processes within a dynamic environmental framework. We apply this framework to deep-time ecosystems made up of extinct plants by constructing plant functional types using fossil remains and simulating—as close as possible—the in vivo response of extinct taxa to their paleoclimatic environment. To accomplish this, foliar characteristics including maximum stomatal conductance, distance from leaf vein to stomata, and cuticular carbon and nitrogen were input as model parameters derived from measurements of well-preserved Pennsylvanian-age fossil leaves. With these inputs, we modeled a terrestrial tropical forest ecosystem dominated by “iconic” plant types of the Pennsylvanian (∼323–299 Ma) including arborescent lycopsids, medullosans, cordaitaleans, and tree ferns using a modified version of the process model BIOME-BGC, which we refer to as Paleo-BGC. Modeled carbon and water—and, for the first time, nitrogen—budgets of a tropical ecosystem from Euramerica driven by daily meteorology are simulated using the Global Circulation Model GENESIS 3.0. Key findings are:
lycopsids have the lowest daily leaf water potential, soil water content, surface runoff, and degree of nitrogen leaching indicating an intensive water use strategy compared to medullosans, cordaitaleans, and tree ferns that have increasingly lower simulated water use, greater surface, and nitrogen loss in this order;
modeled vegetation response to aridification, which was caused by reduced precipitation and intensified through the close of the Carboniferous and into the Permian shows that lycopsids and medullosans have the lowest tolerance for precipitation decrease compared to cordaitaleans and tree ferns, consistent with the paleobotanical record of occurrence of floral turnovers through the Middle Pennsylvanian through earliest Permian;
elevated atmospheric pO2, hypothesized as characteristic for the latter half of the Pennsylvanian and early Permian (∼299–272 Ma), caused higher atmospheric pressure reducing plant transpiration, higher surface water runoff rates, and increased nitrogen export for all plant types simulated, manifested most strongly in the lycopsid dominated ecosystems—with overall only a small reduction in net daily assimilation (≈1 μmol CO2 m−2 s−1).
Both aridification and elevated atmospheric oxygen reduced transpiration, increased water retention in soils, with higher surface runoff. With more discharge, enhanced and higher short-term surface soil loss and silicate weathering would have been possible in broad regions of the paleotropics during the late Carboniferous and early Permian. These results are only obtainable by integrating multiple, fossil-derived measurements into the simulation framework of an ecosystem model that utilizes daily meteorology.
INTRODUCTION
The evolution of terrestrial plants transformed Earth's surface over the Phanerozoic Eon, affecting albedo, geochemical flux rates, and establishing a new, non-aquatic organic carbon pool in the global carbon cycle. Understanding the influence of terrestrial ecosystems on Earth's atmosphere is critical for investigating the dynamics of greenhouse gas-induced radiative forcing on both short and long-time scales. Expanded computational abilities increasingly permit modeling the response of ecosystem functions to past climate change. Such studies, however, raise new challenges. For example, Earth system models of paleo-ecosystems are typically based on modern angiosperms and conifers as physiological surrogates for extinct plants (Kutzbach and Ziegler, 1993; Beerling and others, 1998; Rees and others, 2002; Poulsen and others, 2007; Heavens and others, 2015). Although appropriate for simulations of the Cenozoic (Harrison and Prentice, 2003), the approach is less valid in deep-time because plants of these paleo-ecosystems are increasingly distinct, morphologically and functionally, and evolutionarily distant from extant plants. Recent attempts to include time-appropriate vegetation (Le Hir and others, 2011; Lenton and others, 2016) have been limited to employing basic types (for example, lichens, woody plants) rather than using the traits of extinct plants to infer physiology. As interest increases in modeling vegetation-climate feedbacks during past intervals of climate change, so does the need to reevaluate how paleo-vegetation is represented in these models. Given the wealth of paleobotanical records for many deep-time intervals of Earth history—particularly those records that provide constraints on leaf traits and morphologies (Boyce and Knoll, 2002; Taylor and Taylor, 2009; Peppe and others, 2011; Peppe and others, 2014) and, in some cases, cuticular chemistry (Mösle and others, 1998; Gupta and others, 2007; Montañez and others, 2016; Royer and Hren, 2017)—an opportunity exists for assessing plant function in deep time particularly where models are mechanistic and process-based.
Understanding plant function in deep time is based on inference from a) nearest living relatives (NLR), b) nearest living equivalent (NLE), based primarily on morphology, and/or c) modeling processes from structural characteristics based on limited fossil evidence. It is tempting to use modern plant analogs in deep time. However, plants of the past evolved in unique climates, atmospheres, and community compositions, resulting in potentially novel individual combinations of structural and biochemical configurations that impact their basic physiological function (Wilson and others, 2017). Using mechanistic models to assess the physiological capacity of extinct plants avoids evolutionary progressionist assumptions that infer a persistent, linear trend towards presumed optimization over evolutionary scales of time (Gould, 1994; Romero, 2001; You and others, 2017). Without direct observation of extinct plants, understanding the physiological capacity and ecosystem function of ancient plant communities is best approached through process-based modeling founded on biophysical principles.
Why are mechanistic, process-based models important for simulating the function of paleo-plants and their interaction with climate processes? Simulation at increasingly finer time scales comes closest to bridging the gap between an in silico representation of plant physiology, (that is, the use of computer models to represent biological activity; Peck, 2004) and the unfeasible, in vivo investigation of how ancient flora functioned. Living plant metabolism operates at rapid timescales, from seconds to days. Plant biochemical signatures, archived in leaves, stems, and roots, integrate the dynamic homeostatic mechanisms of the plant's biological activity as it responds to environmental variability over months to decades. In well-preserved plant fossil material, physical structure and chemical composition can be analyzed successfully to reveal both rapid and integrated physiological signals (DeNiro and Cooper, 1989; McElwain and Chaloner, 1995; Bai and others, 2009; Peppe and others, 2011; Franks and others, 2014; Diefendorf and Freimuth, 2017; Soh and others, 2017). Plants, however, filter environmental signals in varying degrees depending upon their evolutionary/taxonomic relationships and functional characteristics (for example, individual plant physiology). Consequently, the accuracy and precision at which environmental drivers can be inferred from plant fossil material is contingent not only on their preservation quality (that is taphonomy) but also on the degree to which species responded to environmental conditions at different timescales while living (Seibt and others, 2008). In the absence of “living fossils”, mechanistic models allow investigation of how paleo-plant taxa may have functioned in terms of daily water use, carbon uptake, and nutrient cycling, among other ecosystem processes.
Here, we utilize an ecosystem model to simulate the ecophysiological functioning of major plant groups (arborescent lycopsids, medullosans, cordaitaleans, and tree ferns) that dominated the late Paleozoic paleotropical forests during Earth's penultimate icehouse (that is, Late Paleozoic Ice Age [LPIA]; Carboniferous Period through early Permian Period: 323 – 290 Ma) and to assess their carbon, water, and nitrogen budgets. The LPIA is a significant event in Earth's history with massive sequestration of biologically fixed carbon, much of it as coal generated from forests, peaked during the overall cool icehouse conditions. The LPIA is also a historical bookend to the current climate warming of which much of this previously fixed carbon is now being released back into the atmosphere. We modified BIOME-BGC, a mechanistic modeling system with ∼45 years of development and testing across modern ecosystems (Waring and Running, 1976; Running, 1984; Running and Coughlan, 1988; Running and Gower, 1991; Hunt and Running, 1992; Kimball and others, 1997; Churkina and Running, 1998; White and others, 2000; Nemani and others, 2002; Thornton and others, 2002; Turner and others, 2007; Bond-Lamberty and others, 2009; Ooba and others, 2010; Wang and others, 2011; Raj and others, 2014; Reeves and others, 2014; Schwalm and others, 2015; Tian and others, 2017; Neumann and others, 2018) to represent paleo-plant characteristics and accommodate paleo-atmospheric conditions and paleo-meteorology (at a daily time-step) for a tropical forest at the end of the Carboniferous. Using this model, we assess three question:
How different were plants possibly of the Late Pennsylvanian in terms of water use and productivity from each other and from modern plant types based on limited information derived from plant fossils?
How might elevated atmospheric oxygen levels that existed during the Late Pennsylvanian directly and indirectly affect water and carbon exchange?
Based on the evidence of aridification during the Carboniferous-Permian boundary, were some Late Pennsylvanian plants better adapted physiologically for this globally significant climate change?
METHODS
Rationale for Model Selection
Paleo-vegetation is not currently represented in fully coupled biosphere-atmosphere climate models (for example, Community Earth System Model; Blackmon and others, 2001). Rather, modern vegetation types have has been used as a surrogate for deep-time climate simulations (Zhou and others, 2012; Heavens and others, 2012; Heavens and others, 2015; Wade and others, 2019). Because plants of the past have different characteristics, especially those associated with water transport and use (Wilson and others, 2017), further consideration is necessary on how to best represent extinct plants in GCMs when assessing the biosphere's potential feedback to climate. For this reason, we utilize BIOME-BGC as a modeling starting point for three primary reasons. First, BIOME-BGC shares a historical connection to the Community Land Model (CLM), the widely-utilized vegetation sub-model of the Earth system model CESM (Bonan and others, 2013) particularly the photosynthesis and soil carbon/nitrogen cycling routines. Second, the stomatal conductance routine of CLM is based on a modified Ball-Woodrow-Berry scheme (Ball and others, 1987), although very accurate for modern plants, requires gas-exchange data for accurate parameterization—values that are impossible to infer from the paleobotanical record. In BIOME-BGC, stomatal conductance is estimated according to Jarvis (1976) where a maximum stomatal conductance (gsmax) value is reduced by scalars representing environmental constraints (for example, light, temperature, vapor pressure deficit, CO2). Where the Jarvis model may be less sensitive to internal leaf CO2, use of gsmax is well-suited for paleo-plant simulations as this parameter can be directly measured from leaf fossils (for example, Franks and Beerling, 2009; Montañez and others, 2016). Third, use of BIOME-BGC permits “off-line” modification of code to evaluate physiological representation of extinct plants as a step towards subsequent integration into fully coupled models.
Given the long-term increase in atmospheric pO2 through the late Carboniferous and early Permian (Berner, 2009; Glasspool and others, 2015; Krause and others, 2018; Lenton and others, 2018; Schachat and others, 2018) to values ∼43 percent higher (maximum pO2 ranges between studies from ∼26–30%) than present-day oxygen concentrations, we assess the effect of higher pO2 on the paleo-plants due to photorespiration (Igamberdiev and Lea, 2006) and increased atmospheric pressure on evaporation (Poulsen and others, 2015). Evidence for increased atmospheric pressure associated with higher oxygen is primarily indirect from insect and bird flight evolution biomechanics (Graham and others, 1995; Dudley, 1998), however remains as an important question given the physical contraints on the atmosphere with flucuating oxygen levels (Berner, 2006). During this time, pCO2 varied (<200–650 ppm), through 100,000 year, eccentricity glacial-interglacial cycles, and through million-year scale cycles (Montañez and others, 2016). With oscillating pCO2 levels, a trend of pantropical aridification occurred, which intensified during the demise of the late Paleozoic icehouse glaciation (Tabor and others, 2008). Repeated plant community turnover in tropical vegetation occurred during this ∼20 million-year period, including the extinction of most kinds of arborescent lycopsids throughout Euramerica (∼304 Ma) (DiMichele, 2014). By using a fossil-parameterized ecosystem-process model, we critically evaluate the sensitivity of plants and ecosystem responses to atmospheric and climatic changes that drove these community turnovers.
Modeled Taxonomic Groups
Throughout the Pennsylvanian, lowland, paleo-tropical forests contained diverse assemblages of plant life, though with persistent dominant taxa (323 – 299 Ma) (DiMichele, 2014; Montañez, 2016). To represent these forests, we chose four plant groups for our simulations commonly found in fossil assemblages for this time period. The first group was the arborescent lycopsids, a distant tree-sized relative of the extant lycopsid Isoetes. Key lycopsid taxa include the stem “form genera” (fossil genera) Lepidophloios and Sigillaria, whereas roots are assigned to the “form genus” Stigmaria and leaves have been assigned to a variety of genera, including Lepidophyllum, which is employed in this study (DiMichele, 1979; DiMichele, 1980). Another group was the medullosan seed ferns, a diverse group of seed plants with a wide range of ecologies, ranging from shade to sun tolerant. Medullosan stems contained high-conductivity xylem cells arranged in anastomosing vascular segments (Stewart and Delevoryas, 1952; Delevoryas, 1955; Basinger and others, 1974; Wilson and others, 2008; Wilson, 2013; Wilson and others, 2015) attached to meter- or multi-meter-long compound leaves (Wnuk and Pfefferkorn, 1984; Laveine, 1986; Laveine and Belhis, 2007). Next, the cordaitaleans were selected as an early coniferophyte lineage, with a wide range of growth habits, from scrambling forms to small and large trees with trunks that contained dense, pycnoxylic wood and bearing strap-shaped leaves that could approach a meter in length (Florin, 1951). Cordaitalean trees were abundant within the central Pangaean wetland flora at certain times and places (Phillips and Peppers, 1984; Raymond and others, 2010; Montañez, 2016), where larger forms dominated many seasonally dry areas (Falcon-Lang and Bashforth, 2004). The last group were marattialean tree ferns, represented by the well-known stem “form” genus Psaronius and the common leaf form genus used here, Pecopteris (although see Cleal, 2015, regarding nomenclature of fossil marattialean foliage), which rose to importance in the late Middle Pennsylvanian and then to dominance in wetland tropical forests after the Middle-Late Pennsylvanian boundary (∼306 Ma; Phillips and others, 1974; DiMichele and Phillips, 2002, Montañez, 2016).
These taxonomic groups were chosen for simulation, in part, due to their importance in major biome shifts that occurred in the latter half of the Pennsylvanian. Montañez (2016) summarized these changes beginning with the transition from lycopsid dominance in the wetland biome to woody cordaitalean–lycopsid forests coincident with an episode of intensification of glaciation (∼311 Ma). Subsequently, two biome turnovers occurred that were concurrent with two brief but deep glacials (∼308 and 306 Ma), resulting in extirpation of cordaitaleans and lycopsid taxa from the wetland biome. The climate warming that followed was associated with increased abundance and dominance of marattialean tree ferns in the wetland habitat. Finally, as drier, seasonal conditions became more frequent, cordaitalean and conifer woodlands became a permanent biome that persisted across the Carboniferous–Permian boundary (299 Ma). Lastly, use of these four taxa captures the functionally important components of the Late Pennsylvanian age forests and permits assessment of the importance of anatomical differences that existed among these plants in the context of a common environment.
BIOME-BGC
Model description and development.
The Biogeochemical Cycles terrestrial ecosystem process model (version 4.1.2; http://www.ntsg.umt.edu/project/biome-bgc.php; Thornton and Running, 2002) was designed to simulate carbon, nitrogen, and water dynamics for the vegetation and soil of a terrestrial ecosystem (Schimel and others, 2000; Thornton and Rosenbloom, 2005). It has been extensively applied to research on modern climate change and disturbance in different ecosystems (Kimball and others, 1997; White and others, 2000; Wang and others, 2005; Lagergren and others, 2006; Chiesi and others, 2007; Ye and others, 2016). The BIOME-BGC model simulates stands of general plant types based on defined ecophysiological attributes of component plant functional types (PFT) and assumes no successional dynamics.
The model simulates primary production using the Farquhar-von Caemmerer-Berry (FvCB) model (Farquhar and others, 1980); in both the sunlit and shaded leaf C pools, assimilation rate is constrained by stomatal conductance based on leaf temperature, CO2 and O2 concentration, leaf C to N ratio (C:Nleaf; kg C kg N−1), and incoming solar radiation. Evaporation and transpiration are estimated using the Penman-Monteith equation for the soil and sunlit/shaded portions of the canopy (Penman, 1948; Monteith, 1965; Monteith and Unsworth, 2008). Autotrophic respiration includes apportionment for maintenance (metabolic) and growth (biosynthesis), where maintenance respiration is based on a Q10 relationship calculated from tissue mass, N concentration, and temperature (Ryan, 1995) and growth respiration is modeled based on static C requirements for different tissues (De Vries, 1975). Fixed C is allocated to plant structures such as leaves, stems, and coarse and fine roots, based on the pool of available N using fixed fractions defined by the user. Projected leaf area is calculated by dividing leaf C by input values of specific leaf area (SLA; leaf m2 kg C−1). Differences in SLA values among plant types are an important, intrinsic variable that mediates canopy radiation absorption, water interception, photosynthesis, transpiration, and litter inputs to detrital pools. Plant biomass C is transferred to litter that eventually becomes part of the soil organic matter. Litter and soil C and N are each divided into four pools reflecting differing rates of input-biomass chemical degradation associated with lignin and hemicellulose concentrations assigned by the user for different plant parts. Nitrogen immobilization and mineralization from litter and soil pools are based on threshold C:N ratios reflecting the different requirements of decomposer biota within those pools. As C is allocated to different parts of the plant (that is leaves, stems, roots), mineralized N is taken up from the litter and soil pools as part of the N cycle based on the input C:N of specified for different plant parts representing tissue sink strength.
Paleo-BGC.
For this study, we created a modified version of the BIOME-BGC model, referred to as Paleo-BGC, to include paleo-environmental and fossil plant structural attributes, the latter of which can be assessed and quantified from the fossil record (fig. 1). Inputs included: A) leaf hydraulic conductivity based on measured leaf venation and stomatal density, B) estimated mesophyll conductance, C) foliar C:N values and D) input of pCO2 and pO2 with related changes in atmospheric pressure, the density of air, and specific heat of air. While geologic evidence low atmospheric pressure has been estimated from rain droplet imprints in Archean sediments (Som and others, 2016) and modern lava basaltic vesicles related to altitude (Sahagian and Maus, 1994), evidence for high pressure is limited to inference from biological evolution of flight and respiratory requirements of insects (Graham and others, 1995; Dudley, 1998).
Simplified flowchart of the Paleo-BGC model including inputs from fossil leaf data, GENESIS-derived daily meteorology, and input atmospheric gas concentrations of CO2 and O2.
Leaf hydraulic conductivity.
Vascular tissue in leaves affects leaf water and, particularly, as changes in guard cell turgor pressure cause opening and closing of stomata, stomatal conductance. Because of non-analog morphologies in extinct taxa, it is critical to consider the coordinated evolution of features important in leaf water supply and distribution (conductivity) and the control of stomatal aperture size (conductance; McElwain and others, 2016a). Leaf hydraulic conductivity (KL; mol H2O s−1 m−2 MPa−1) constrains leaf water potential through resupply of water lost via transpiration to leaf cells. McKown and others (2010) and others have shown that KL has two components, leaf xylem and outside xylem conductivity (Kx and Kox, respectively; mol H2O s−1 m−2 MPa−1) where KL can be rewritten as:
We considered two options for solving KL including 1) as a function of the modified Hagen–Poiseuille equation by summing up the conductance of each individual xylem cell in all of the veins of the leaf (McCulloh and others, 2009) or 2) empirically, as a function of mesophyll path length (Dm). Mechanistically, inclusion of conductivity specific to an individual vein seems favorable; however, as pointed out in Scoffoni and others (2011), the derivation of conductance from the application of first principles based on measurements of vein size, vein dimensions, xylem cell dimensions, et cetera results in a large error predicting Kx when applied to a single leaf. This error propagates when extrapolated to the canopy scale. Therefore, when faced with two tradeoffs (that is mechanistic vs. empirical and specificity vs. generality), we chose the empirical and generalizable, as this form is particularly appropriate for simulation of paleo-flora with limited available anatomical information. The empirical formulation of total leaf conductance (including both Kx and Kox) by Brodribb and others (2007) based on empirically-derived Dm values provide a robust analysis for all plant evolutionary groups, where:
However, this value is static, whereas, operationally, KL varies as a function of stem water potential (Sack and Holbrook, 2006). Based on the curves from figure 2 of Blackman and others (2009), we define the change in KL calculated from the maximum leaf hydraulic conductivity (KLmax) based on the shortest mesophyll path length (Dmmin):
and the minimum conductivity value (KLmin), which is calculated from the longest mesophyll path length (Dmmax):
The change between maximum and minimum KL is assumed to co-occur with maximum and minimum stem-water potentials (ΨSmax, ΨSmin), which were set at values of −0.5 and −2.5 MPa, respectively. The value of −2.5 MPa is a general water potential value for loss of plant cell turgor pressure, also associated with leaf hydraulic restrictions (Bartlett and others, 2016). Plant water transpired is replaced by the shortest pathlength that becomes increasingly depleted from sources closest to the stomata, such as veins. Approaching the wilting point (ΨSmin), only the longest paths from source to stomata are likely to continue to rewet cells. To calculate the daily value of KL, we first define a linear desiccation response function based on the path length from the plant's anatomical features related to plant water availability. This is determined by:
where m is the slope. The intercept (b) is defined by:
with the daily value of KL calculated by:
where ΨS is the daily stem-water potential, assumed to be equivalent to the soil-water potential from the previous day. Soil-water potential is calculated as a function of soil textural characteristics, soil depth, and the net-water balance of the soil, which includes precipitation, runoff, evapotranspiration, and vertical drainage.
A schematic illustration of the feedback of stem water potential on leaf hydraulic conductivity. Upper panel (A): maximum leaf hydraulic conductivity (KLmax) is derived from the minimum mesophyll hydraulic pathway (Dmmin) and minimum leaf hydraulic conductivity (KLmin) is derived from the longest mesophyll hydraulic pathway (Dmmax). KLmax will be reached when stem water potential is highest and KLmin will be reached when stem water potential is lowest. In Paleo-BGC, leaf conductivity decreases linearly with a slope of m as stem water potential decreases. Lower panel (B): three possible paleo-plant leaf hydraulic conductivity response to stem water potential. For Plant 1 (black squares), short minimum and long maximum values of Dm, leaf hydraulic conductivity decrease from high values to low values with stem water. Plant 2 (gray triangles), with moderate but consistent lengths of Dm, have lower maximum leaf hydraulic conductivity but higher minimum KL values. Plant 3 (open circles), with longer mesophyll pathways, will have the lowest KLmin and KLmax values. See methods for details.
This linearized model approximates leaf behavior during evaporation, mimicking the “desiccation curve” response defined by experimental leaf hydraulics (fig. 2A). In leaves with short paths from veins to stomata, KLmax will be high and vice versa. The distribution of path sizes plays a role in defining leaf behavior: plants with consistently short paths from veins to stomata will be able to maintain a high KL as the leaf desiccates, whereas in plants with a large range between the shortest and longest pathway, KL decreases at a faster rate. Finally, these leaf pathways are preserved in the fossil record and can be measured directly from fossilized material.
Leaf-water potential (ΨL) is modeled as a function of leaf conductivity and balance between the internal leaf water supply and transpiration. Daily, ΨL is calculated as:
where gsmax is the maximum stomatal conductance of water vapor (mol H2O m−2 s−1) supplied as an input variable for each plant type. The rationale for using gsmax is that this sets the upper bounds of transpiration and assumes that plants, in the morning, may open their stomata fully in response to blue light illumination (Roelfsema and Hedrich, 2005). Note that Paleo-BGC does not include the effect of plant height on gravitational water potential (a minor component of plant water potential in most trees), which will be added to a forthcoming version incorporating stem hydraulic characteristics to limit leaf water supply.
In Paleo-BGC, operational stomatal conductance of water vapor (gs) is partially based on a water-potential scalar (mΨL), ranging from 0 to 1:
where ΨLmin and ΨLmax are the minimum and maximum leaf-water potential values at the start and completion of stomatal closure, set for these simulations at −0.5 and −2.5 MPa, respectively (water potential is expressed in negative units of pressure, such that the maximum value is the most negative, rather than the most positive, number). In addition, other scalars representing the effect of light, VPD, and temperature are all multiplied to reduce gsmax to daily gs; these are retained from the original formulation of BIOME-BGC (Jarvis, 1976; Thornton and Running, 2002).
Leaf conductance.
Leaf conductance is constrained by several factors, including characteristics of the leaf boundary layer, conductance of stomata, cuticle thickness, and mesophyll tissue. Originally, BIOME-BGC only included conductance of the boundary layer, stomata, and cuticles with mesophyll conductance (gm) assumed to be included as part of gsmax. Mesophyll tissue, however, may affect gas exchange independently of gsmax especially for leaves with dense internal tissue filled with cells (Flexas and others, 2012) and may have varied across geologic time (Yiotis and McElwain, 2019). Direct measurements of mesophyll conductance are not available for paleo-plants; therefore, the gas exchange pathway must be derived from measurements on or analysis of material preserved in the fossil record. Using the 1-D diffusion model of Niinemets and Reichstein (2003), we estimate gm from the fraction of air space in the leaf, mesophyll tissue thickness, mean mesophyll cell width, and cell-wall thickness determined from published cross-sections of fossil leaves of representative Late Pennsylvanian taxa (See Appendix B). From these data, a mean value of 0.273 mol H2O m−2 s−1 was derived and used for these simulations. Using all this information, total leaf conductance (gL; including gm) is calculated in Paleo-BGC as:
where gb (mol H2O m−2 s−1) is the boundary layer conductance and gc (mol H2O m−2 s−1) is the cuticular conductance. The values of gb were specified for each plant type to account for leaf size effects on gas exchange, and based on the assumption of forced convection transfer (Campbell and Norman, 2012):
where d is 0.72 times the mean width of a leaf, leaflet, or pinnule (for fern fronds), and u is the wind speed, assumed to be a constant value of 2.0 m s−1. Leaf widths for each plant type were derived from published sources (Somner, 1989; Hernandez-Castillo and others, 2001; Montañez and others, 2016). The value of gc was set to a value of 4 × 10−4 mol H2O m−2 s−1 (Jones, 1992) for all plant types as a conservative estimate in the absence of measured values.
Photosynthesis.
The FvCB photosynthesis model used in BIOME-BGC predicts CO2 assimilation based on limitation by either the rate of carboxylation or by substrate concentrations (Farquhar and others, 1980). This model is useful for paleobotanical studies because it includes the Michaelis-Menten kinetics of enzyme competition of ribulose-bisphosphate carboxylase/oxygenase (Rubisco), which is affected by both internal leaf CO2 (Ci; Pa) and O2 (Oi; Pa) concentrations. For Paleo-BGC, we updated the temperature-dependent enzyme kinetic sub-models based on the versions included in CLM 4.5 (Oleson and others, 2013).
Calculation of assimilation begins with assessing the maximum rate of carboxylation (Vcmax; μmol CO2 m−2 s−1) and the maximum potential rate of electron transport associated with the regeneration of the reactant substrate ribulose-bisphosphate (RUBP) (Jmax; μmol CO2 m−2 s−1). The value of Vcmax at a 25°C (Vcmax25; μmol CO2 m−2 s−1) is approximated as the product of leaf N concentration on an area basis (Na; g N m−2 leaf area), the fraction of leaf N partitioned by the plant into the Rubisco enzyme (fLNR; g N in Rubisco per g leaf N), the fraction of the total molecular mass of Rubisco that is N (fNR; g total mass of Rubisco per g N in Rubisco), and the mass-specific activity rate of Rubisco at 25 °C (aR25; μmol CO2 g−1 Rubisco s−1):
where (1) daily values of Na; are derived from growth, (2) N allocation to foliage is constrained by plant-specific SLA; values, (3) fLNR; is a static input variable, and (4) fLNR; and aR25 are constants with values of 7.16 g Rubisco g−1 N in Rubisco and 60.0 μmol CO2 g−1 Rubisco s−1, respectively. Leaf N is calculated from daily net photosynthetic allocation and constrained using input foliar carbon to nitrogen ratios (C:Nleaf). This ratio, set by the user, limits maximum carbon allocated for leaves based on the plant-available N, which is derived from soil N mineralization and is required to maintain a minimum leaf N. Unconstrained allocated carbon is used to calculate a new daily leaf N concentration on a leaf mass basis (Nm; g N g C−1) with Na = Nm/SLA (eq 13). Specific Leaf Area (SLA), as a species-intrinsic variable that specifies plant's foliar area relative to its mass, is related to C:Nleaf (White and Scott, 2006), cuticle thickness (Soh and others, 2017), adaxial stomatal density (Haworth and Raschi, 2014), petiole width (Royer and others, 2005), and petiole size (Li and others, 2008). The value used here for fLNR was set for all plant types to 0.06 g N in Rubisco g−1 N, an intermediate value of those used in CLM 4.5 for all modern forest types.
Estimated Vcmax25 is then used to calculate a realized value of Vcmax based on the change in entropy, activation, and deactivation energy of the carboxylase reaction (Harley and others, 1992; Oleson and others, 2013). Similarly, Jmax was calculated to include temperature kinetic effects by assuming that Jmax scales with Vcmax by a factor of approximately 1.97.
Atmospheric oxygen affects baseline photosynthetic activity of plants through the CO2 compensation point (Γ*; Pa):
where Kc and Ko are the Michaelis-Menten half-saturation coefficients for carboxylase/oxygenase that are adjusted for temperature, and Vomax is the maximum rate of Rubisco oxygenase.
Daily average net assimilation (An; μmol CO2 m−2 s−1) is calculated based on the minimum of the limiting rates of carboxylation (Ac) and RUBP regeneration (Aj). Formulation of these rates are:
where Rd is the leaf respiration rate (μmol CO2 m−2 s−1), calculated as a function of nitrogen concentration and temperature, and:
where J is the actual rate of electron transport calculated from Jmax as a function of absorbed photosynthetically active radiation and the quantum efficiency of electron transport. Because both Ac and Aj include Ci, we solved both quadratically, assuming that assimilation is proportional to the CO2 gradient between plant and atmosphere, and stomatal conductance to CO2 [A = gCO2s(Ca − Ci); eq 17]. Thus we allow for external factors to dominate gCO2s (gCO2s = 
; mol CO2 m−2 s−1; eq 18), which is typical of the “Jarvis” form for leaf gas-exchange (Jarvis, 1976).
Canopy precipitation interception.
The BIOME-BGC model derives plant canopy water interception from input precipitation based on the calculation of a maximum canopy storage capacity. Throughfall (qthru; kg H2O m−2 s−1) is calculated based on precipitation exceeding this maximum capacity. In Paleo-BGC, qthru is calculated using the function described by Oleson and others (2013) as:
where prcp is site precipitation derived from the meteorological data, α is a scale conversion factor set to a value of 0.25, and LAI is the leaf area index (m2 m−2) derived from the carbon budget of Paleo-BGC.
Paleo-floral Plant Functional Types and Physiological Characterization
The four broad taxonomic groups representing dominant Late Pennsylvanian paleo-floral PFTs include lycopsids (for example, Lepidodendron, Lepidophloios, Ulodendron), pteridosperms (specifically Macroneuropteris scheuchzeri), cordaitaleans (for example, Cordaites, Cordaadaxicutis, Cordaabaxicutis), and marattialean tree ferns (for example, Pecopteris, Psaronius). To simulate these plant types in Paleo-BGC, general parameters were derived by modifying the evergreen broad-leaved forest PFT, a default modern plant type distributed with BIOME-BGC model. Paleo-floral PFTs are defined by modifying four specific leaf parameters: 1) leaf conductance, involving both gsmax and gb, 2) C:Nleaf, 3) SLA, and 4) Dm.
Median values of gsmax and gb (and minimum/maximum values) were estimated from plant leaf measurements [Somner, 1989; Hernandez-Castillo and others, 2001; Montañez and others, 2016 (table 1)]. For paleo-plants, initial C:Nleaf values are derived from C and N and isotopic analyses of a set well-preserved leaf fossil cuticles (n=120) representative of the four taxonomic groups simulated in this study as well as several other taxonomic groups characteristic of late Carboniferous tropical forests (Montañez unpublished data). The C:N ratios of the plant fossils cuticles range between 24 and 71, a similar range as that of modern vascular plants (Meyers and others, 1995) including deciduous needle- (median C:N of 26; range of 17–37) and broad-leaf (24; range of 16–36) and evergreen needle-leaf forests (41; range of 23–70) and shrubs (33; range of 17–57) (White and others, 2000). Moreover, the plant fossil C:N ratios define subpopulations by taxonomic group with the inferred longer leaf life span of lycopsids (median of 41; range of 35–47) and cordaitaleans (median of 42; range of 39–44), compared to marattialean tree ferns (avg. of 36; range of 29–43) and Macroneuropteris (median 30; range of 28 – 32). As these plant fossil C:N ratios overlap with values from modern vascular plants, have high-quality morphology, and have reasonable δ13C (−26.5 to −23‰) and δ15N (−0.5 to +6.5‰) values, this supports the excellent preservation of the geochemical compositions of the fossil cuticles (Mösle and others, 1998; Tu and others, 1999; Gupta and others, 2007). That said, the C:N ratios of leaf fossils can be shifted during litter decomposition and subsequent burial, and additionally may differ from overall leaf C:N both in vivo and taphonomically. Decomposition by soil fauna and/or under anaerobic conditions has been shown to decrease C:N ratios (Rice and Tenore, 1981; Meyers and others, 1995), although experimental studies have largely documented that litter decomposition typically increases C:N ratios by 1.5- to 2-fold (Rice and Tenore, 1981; Benner and others, 1990; Lanuza and others, 2019). A synthesis of litter C:N ratios associated with modern evergreen needle- (ENF) and deciduous needle-(DNF) and broad-leaf (DBF) forests as well as shrubs (White and others, 2000) indicates a greater than doubling of C:Nlitter ratios over living leaf values (that is, median of 93 vs. 42 for ENF; 120 vs. 27 for DNF; 55 vs. 25 for DBF; and 75 vs. 35 for shrubs). These modern litter C:N ratios are well above the average values of their fossil taxonomic group counterparts (for example, median of 49 for evergreen needle-leaf voltzian and walchian conifers vs. 93 for the C:Nlitter ratio of the modern ENF analog). Differential increase in cuticle C:N ratios likely reflect differences in macromolecule composition and resistance to decomposition of the original epicuticular waxes of the plant fossils (compare Rice and Tenore, 1981; Mösle and others, 1998; Gupta and others, 2007).
Median (± minimum and maximum values) stomatal conductance (gsmax, mol H2O m−2 s−1) calculated for each plant type based on stomatal measurements from fossil cuticles based on the methods described in Franks and Farquhar (2001)
In Paleo-BGC, fossil-derived C:Nleaf values were also used to initialize the C:N values for other biomass pools of leaf litter (C:Nlitter), fine roots (C:Nfineroot), live wood (C:Nlivewood), and dead wood (C:Ndeadwood). These pools are important to consider given that living biomass C:N drives sink strength for soil uptake in the modeled growth, and C:N of non-living biomass drives environmental N cycling, including immobilization and mineralization rates, in the simulated soil pools. In the absence of tissue-specific samples, we used the ratio of C:N values of the different biomass pools [that is, (C:Nleaf):(C:Nlitter):(C:Nfineroot):(C:Nlivewood):(C:Ndeadwood)] for the modern evergreen broad-leaved forest PFT in BIOME-BGC (1: 1.17: 1: 1.19:7.14) to initialize the C:N values for each paleo PFT based on measured leaf fossil cuticle C:Nleaf (Montañez unpublished data). For example, the fossil C:Nleaf of the lycopsid group is 41.0, which using the above ratio produced C:Nlitter, C:Nfineroot, C:Nlivewood, and C:Ndeadwood values of 47.8, 41.0, 48.8, and 292.9, respectively (table 2).The SLA values for the paleoplant types were derived from the C:Nleaf values based on relationships between SLA and C:Nleaf developed for modern New Zealand podocarps and tree ferns (White and Scott, 2006; Montañez and others, 2016) (table 2; See Appendix B).
Median carbon to nitrogen ratios for leaves derived from fossil cuticles for paleoflora
To calculate Dm, we followed Brodribb and others (2007) and measured flow path distance from vein to stomata (μm) from micrographs published in literature sources for relevant taxa using ImageJ (https://imagej.nih.gov/ij/) (table 3). This method specifies that for a given measurement of cell width (Cx) and height (Cy) where Cx > Cy, the mean horizontal path length (X) is evaluated as X = ν/Cx[(Cx − Cy) + πCy/2] (eq 20) where v is the average longest length from veins to stomata and the vertical path length (Y) is Y = t/Cy(πCy/2) (eq 22) where t is the leaf thickness from vein to stomata. In cases where Cx < Cy, the dimensions are switched in these functions. Finally, Dm is derived from Dm = 
(eq 23). When no leaf cross-sectional data exists, Dm was calculated as the mean intervein distance measured from photographs multiplied by 2.0 to account for path tortuosity.
Mesophyll pathlength (Dm; μm) for paleo- and modern flora derived from literature sources
Paleometeorology
Daily inputs required by Paleo-BGC were derived from a climate simulation made using the GCM GENESIS 3.0, coupled to a land-surface model, including topography, vegetation, soil, and snow (Peyser and Poulsen, 2008). Briefly, the atmospheric portion GENESIS has a resolution of 3.75° × 3.75° with 18 vertical levels to 100 km above the Earth's surface with the land-surface grid having a spatial resolution of 2° × 2°. Sea surface temperatures were computed based on a diffusive heat flux using a 50-m slab oceanic layer (Pollard and Thompson, 1995; Pollard and Thompson, 1997). Late Paleozoic simulations were constructed, following Horton and others (2010), using an early Permian (Sakmarian, about 290 Ma) continental configuration and topography (after Ziegler and others, 1997). Solar luminosity was reduced by 2.5 percent relative to modern values in accordance with solar evolution models (Gough, 1981). The climate simulation represent an interglacial scenario with no land ice and pCO2 of 600 ppm and pO2 of 0.21 (following Montañez and others, 2016), a circular orbit with an obliquity of 23.5°. Within GENESIS, a “place-holder” vegetation cover is required to represent land surface feedbacks to the climate system that was prescribed as a uniform needle-leaved evergreen forest type to generate the daily meteorological data. From the GENESIS simulation, ten years of daily minimum/maximum temperature (°C), precipitation (m s−1 - converted to cm day−1), and solar radiation (W m−2) were generated for input into Paleo-BGC.
Point Modeling
For required meteorological data, values were extracted from the GENESIS output for a single paleo-location within the present-day Illinois Basin (central Euramerica) where tropical wetland forests, likely refugia, containing our plant species existed during interglacial periods of the Carboniferous (fig. 3). In addition to the meteorological inputs from GENESIS, vapor pressure deficit (VPD, in Pa) and day length (in seconds) were generated from the MT-CLIM model (Running and others, 1987; https://www.ntsg.umt.edu/project/mt-clim.php) using the GENESIS-derived daily minimum/maximum temperature, precipitation, and paleolatitude as input data.
Annual (A) temperature (°C), (C) precipitation (cm yr−1), (E) net shortwave radiation (GJ yr−1), and (G) vapor pressure deficit (kPa) are derived from the global GENESIS GCM for 300 Mya for the interglacial scenario with pCO2 = 600 ppm and pO2 = 0.21 mol mol−1. In addition, (B) temperature (°C), (D) precipitation (cm d−1), (F) radiation (W m−2), and (G) vapor pressure deficit (kPa) are presented as daily averages from the 10-years of GENESIS simulation extracted for location “A” representative of the tropical moist forest of the Illinois Basin flora. The shaded areas on the daily average graphs denoting ± one standard deviation.
Forest stands were simulated using Paleo-BGC for each plant type based on 50-year model runs, in which the 10-year GENESIS-derived meteorological data were cycled five times. The length of these simulations was set to allow the model to reach a steady-state for the carbon, water, and nitrogen pools. The Paleo-BGC model initialization routine, referred to as a spin-up, allows ecosystem carbon and nitrogen pools to reach equilibrium, which was assigned a maximum value of 10,000 years (Thornton and Rosenbloom, 2005).
Simulation Sets
All simulation sets included utilizing a single set of meteorological data extracted from the GENESIS simulation as input into Paleo-BGC (fig. 4). With the fossil-derived parameters to define the four paleo plant types, three simulated sets were run. The objective of modeling Late Pennsylvanian tropical wetland forests during late glacials (DiMichele, 2014; Montañez and others, 2016) under interglacial atmospheric compositions and meteorology was to evaluate the sensitivity and response of the wetland plants to changing environmental conditions during turnovers from glacial-to-interglacial conditions and longer-term climate shifts. An initial, baseline set of simulations were run for each of the four paleo-floral PFTs (that is lycopsids, the medullosan - Macroneuropteris scheuchzeri, cordaitaleans, and tree ferns) to assess model function including outputs of carbon, water, and nitrogen budgets for each of the plant groups with the emphasis on understanding differences between the plants based on the limited fossil-derived characteristics input into Paleo-BGC. These baseline simulations assumed atmospheric pCO2 of 600 ppm, representing interglacial values (Montañez and others, 2016), and pO2 of 0.21 mol mol−1. The high CO2 and modern oxygen (pO2 = 0.21 mol mol−1) setting is chosen as a demonstration of the model for interglacial, Early Pennsylvanian conditions.
Flowchart showing the basic model interaction and data requirements for each simulation set conducted. All simulations were run assuming an atmospheric pCO2 = 600 ppm.
As pO2 increased linearly through this period from near present-day levels (0.21) to between 0.26 and 0.28 mol mol−1 during the Late Pennsylvanian (Krause and others, 2018), the second set of models was run for the four paleo-floral PFTs using the same meteorology but with an input pO2 of 0.28 mol mol−1. Because oxygen is the second-largest component of Earth's atmosphere, the molar mass of dry air (Ma), the specific heat of air (cp), and atmospheric pressure (P; Graham and others, 1995) were adjusted accordingly following the approach described by Poulsen and others (2015). As a result, an atmosphere with a pO2 value of 0.28 has an estimated dry air mass of 29.03 g mol−1 and a sea-level atmospheric pressure of 113.57 kPa, compared to the modern values of 28.96 g mol−1 and 101.33 kPa, respectively. Whereas these differences in values may seem subtle, they are important for calculating evaporation processes in Paleo-BGC. The mass of dry air is directly related to the specific heat of air, where cp = (7/2)(R/Ma) (eq 23), which decreased to 932.7 J kg−1 K−1 for pO2 = 0.28 mol mol−1 from 1010.0 J kg−1 K−1 for pO2 = 0.21 mol mol−1. Both Ma and atmospheric pressure are important for calculating evaporation processes in Paleo-BGC, and cp affects calculation of the psychometric “constant” (γ; K−1):
where λ is the temperature-dependent latent heat of vaporization, and Mw is the molar mass of water (18.01 g mol−1). The value of γ is an input in the Penman-Monteith equation for estimating evapotranspiration based on net radiation, temperature, and vapor pressure deficit (Monteith, 1965). Higher pO2, through cp and P, results in increased atmospheric pressure and reduced rates of evaporation and transpiration.
The third set of models were used to assess physiological limitations specific to each paleo-plant PFT due to the effect of long-term intensifying pantropical aridification during the Pennsylvanian through the early Permian. For this, precipitation data, initially derived from GENESIS as input into Paleo-BGC, were sequentially reduced by 80, 60, 40, 20, and 10 percent within the Paleo-BGC initialization file to produce five in separate simulations that ranged from 142.9 to 17.8 cm yr−1 for each paleo-plant type and atmospheric compositions of pCO2 = 600 ppm and pO2 = 0.28 mol mol−1. Whereas the Pennsylvanian to earliest Permian component of the long-term aridification occurred within the context of declining pCO2, from 600 to < 200 ppm, this modeling is intended to isolate the effects of water availability on the different plant groups for a given atmospheric CO2 concentration, as would be the case for spatial-environmental heterogeneity at any given time.
In addition, two modern vegetation types, evergreen broad-leaved and needle-leaved forests, were also simulated for comparison. The input parameter values of C:N, SLA, et cetera for the modern vegetation types were taken from the original ecophysiological constant files in the BIOME-BGC code (http://www.ntsg.umt.edu/project/biome-bgc.php) that were derived from the data summary in White and others (2000). These modern vegetation types were selected because, as evergreen vegetation, they are phenologically closest to the Late Pennsylvanian flora and have been used in previous studies as paleoflora proxies for Paleozoic modeling (Poulsen and others, 2007; Horton and others, 2010; Montañez and others, 2016; Wade and others, 2019).
RESULTS
Modeling Set #1: pCO2 = 600 ppm and pO2 = 0.21 mol mol−1 Scenario
The GENESIS 3.0-derived meteorology for paleo-location of the simulation site had an average daily temperature of 27.6 °C and mean annual precipitation of 178.6 cm (figs. 3B and 3D). Spatial variability in the GENESIS meteorological data around the simulation site based on an analysis of a 6° × 6° grid, or approximately 450, 000 km2 area, showed mean annual temperature ranged from 18.5 to 29.2 °C. Annual precipitation values ranged from 22.6 to 235.6 cm per year. Precipitation was seasonal with rainfall primarily between simulated days 75 (early spring) and 325 (early winter). Based on the rainfall having exceeded simulated evapotranspiration for 7 months per year for this location, the climate was classified as sub-humid (Cecil and others, 2004).
The simulated seasonal precipitation affected plant-water use by the paleo-floral PFTs, with lycopsids having the lowest and tree ferns the highest based on our analysis of leaf water potential and soil water content. Given that lycopsids have the highest maximum conductance (gsmax; table 1), this resulted in lower average daily leaf water potential (ΨL, −1.67 MPa) for a given stem water potential (ΨS) (fig. 5). The high water use by lycopsids and consequent soil water depletion reduced the simulated days per year that the lycopsid ΨL values were above the wilting point (−2.5 MPa) (fig. 6), thus decreasing the growing season by approximately two months when compared to the other paleo PFTs. Despite a shortened growing season, the lycopsid PFT had the lowest average annual volumetric soil-water content (0.56) of all PFTs (fig. 7), having transpired a significant amount of soil water. As a consequence of the low soil water content, the lycopsid PFT also had the lowest amount of predicted surface runoff (184.2 L m−2 yr−1) from its simulated ecosystem.
Simulated soil and leaf water potentials (MPa) for the paleo-floral plant functional types by Paleo-BGC for each plant type for 50 years of simulations with pCO2 = 600 pm and pO2 = 0.21 mol mol−1. The plant wilting point (−2.5 MPa) is shown as a dotted line for reference.
Average daily leaf water potential (MPa) predicted by Paleo-BGC for each plant type for 50 years of simulations with pCO2 = 600 pm and pO2 = 0.21 mol mol−1. The plant wilting point (−2.5 MPa) is shown as a dotted line for reference. The shaded areas represents ± one standard deviation from the mean based on the 50 years of simulation for each plant type.
Average daily fraction of available soil predicted by Paleo-BGC for each plant type for 50 years of simulations with pCO2 = 600 pm and pO2 = 0.21 mol mol−1. The shaded areas represent ± one standard deviation from the mean based on the 50 years of simulation for each plant type.
For non-lycopsid PFTs, lower gsmax produced higher ΨL values as a result of less water lost through transpiration. In particular, the tree fern PFT had the highest average ΨL (−0.94 MPa). Less depleted soil water resulted in higher average soil water content (0.67) and higher surface runoff (740.1 L m−2 yr−1). For all paleo-floral PFTs, variability in surface runoff mirrors seasonal precipitation, with peak flows occurring during mid-summer and complete cessation of flow during the late autumn (fig. 8).
Average daily surface runoff (L H2O m−2 yr−1) predicted by Paleo-BGC for each plant type for 50 years of simulations with pCO2 = 600 pm and pO2 = 0.21 mol mol−1. The shaded areas represent ± one standard deviation from the mean based on the 50 years of simulation for each plant type.
Differences in simulated N leaching among the paleo-floral PFTs showed the lycopsid PFT had the lowest annual N loss (46.1 g N ha−1yr−1), compared to the tree fern PFT with the highest (135.5 g N ha−1 yr−1; table 4). In contrast, the modern PFTs had higher N leached for the same meteorological conditions (191.6 and 163.7 g N ha−1 yr−1 for the evergreen broad- and needle-leaved forest PFTs, respectively). Annual N mineralization rates for all paleo-floral PFTs were similar with an average value of 40.9 kg N ha−1 yr−1 (table 4). Mineralization for the modern evergreen broad-leaved forest PFT were comparable to values derived from the paleo-floral PFTs (42.3 kg N ha−1 yr−1), whereas the average mineralization rate for the modern evergreen needle-leaved forest PFT (29.5 kg N ha−1 yr−1) was lower than for the paleo-PFTs.
Average net inorganic nitrogen mineralization (kg N ha−2 yr−1) and leaching (g N ha−2 yr−1) simulated from Paleo-BGC for paleo- and modern PFTs for different pO2 levels with pCO2 = 600 ppm
For carbon, average net assimilation (An) ranged from 8 to 10 μmol CO2 m−2 s−1 for all paleo-floral PFTs (fig. 9). With the same meteorology, An) for the modern PFTs was higher ranging from 10 to 12 μmol CO2 m−2 s−1, with the evergreen broad-leaved forest PFT having the highest rate. Relative to photosynthesis, large differences in biomass were found. Biomass predicted from the end of the simulations (50 years) was used to represent mature stands, with the medullosan, M. scheuchzeri, having the lowest biomass of 85.3 kg C m−2 compared to the cordaitalean PFT with the highest value of 171.4 kg C m−2 (table 5). Overall, the average biomass values of the paleo-PFTs were lower than those of the evergreen broad- and needle-leaved forests (179.5 and 156.2 kg C m−2, respectively). Leaf area index (LAI), a canopy variable related to canopy biomass and cover, varied among PFTs, ranging in value from a low of 3.1 m2 m−2 for the M. scheuchzeri PFT to a high of 5.4 for the tree fern PFT. For comparison, LAI values for modern broad- and needle-leaved evergreen forests were 2.1 and 4.8 m2 m−2, respectively (table 5).
Daily net assimilation (μmol CO2 m−2 s−1) predicted by Paleo-BGC for each plant type for 50 years of simulations with pCO2 = 600 pm and pO2 = 0.21 mol mol−1, excluding zeros. Boxes represent the 25th and 75th percentiles with average shown as a line. Whisker lines represent the 95% confidence interval with outliers shown as individual points. The evergreen broad- and needle-leaved forest simulations are based on modern flora using the same GENESIS meteorology shown for reference.
End-of-simulation leaf area index (LAI; m2 m−2) and vegetation biomass values (kg C m−2) simulated from Paleo-BGC for paleo- and modern PFTs for different pO2 levels with pCO2 = 600 ppm
Combining the water and carbon budgets, average intrinsic water-use efficiency values (WUEi = An/gs; μmol CO2 mol H2O−1) greatly differed among the PFTs. Lycopsids and M. scheuchzeri had the lowest WUEi values (7.3–13.7 μmol CO2 mol H2O−1), whereas the cordaitaleans and tree fern PFTs had higher average water-use efficiencies (41.8–61.5 μmol CO2 mol H2O−1) (fig. 10). In these simulations, then, lower transpiration from cordaitaleans and tree ferns resulted in more assimilation per volume of water transpired.
Daily average intrinsic water-use efficiency values (WUEi; μmol CO2 mol H2O−1) predicted by Paleo-BGC for each plant type for 50 years of simulations with pO2 = 0.21 and pO2 = 0.28 mol mol−1. All simulations were run with pCO2 = 600 pm and whisker lines represent the 95% confidence interval.
Modeling Set #2: pCO2 = 600 ppm and pO2 = 0.28 mol mol−1 Scenario
With higher atmospheric oxygen, plants used less water and had higher ΨL, resulting in more water remaining in the soil and greater surface runoff for all PFTs. The lycopsid PFT, like the previous scenario with modern pO2, used the most water of all PFTs—though with slightly higher ΨL (−1.58 MPa), higher average soil-water content (0.60), and a doubling of surface runoff value of 321.2 L m−2 yr−1 (as compared to simulations with modern pO2). Conversely, the tree fern PFT had the lowest water use and the highest average ΨL (−0.86 MPa), with slightly more soil water content (0.66) and runoff (828.5 L m−2 yr−1) than when atmospheric oxygen concentrations were similar to the present.
Increased soil water content, a consequence of higher O2 levels, resulted in higher annual N leaching rates. The lowest annual leaching rate among PFTs in the high O2 scenario (72.2 g N ha−1 yr−1 with lycopsid) was more than 1.5× the leaching rate under the present-day pO2 scenario (table 4). Net N mineralization rates do not change appreciably between the two pO2 scenarios.
Average daily net assimilation of CO2 in the high O2 scenario for all paleo-floral PFTs was reduced by ∼ 1 μmol CO2 m−2 s−1, that is, 93 percent of An simulations with modern pO2 levels (fig. 11). Despite lower mean assimilation, paleo-floral leaf area and biomass values were equivalent to those for the modern pO2 simulations, both in terms of value and range (table 5). Average WUEi with high pO2, however, was lower for the cordaitalean and tree fern PFTs by ∼ 5–6 μmol CO2 mol H2O−1 compared to the modern pO2 simulations (fig. 10).
Daily assimilation (μmol CO2 m−2 s−1) predicted by Paleo-BGC for each plant type for 50 years of simulations with pO2 = 0.21 and pO2 = 0.28 mol mol−1. All simulations were run with pCO2 = 600 pm and whisker lines represent the 95% confidence interval.
Modeling Set #3: Aridification and Precipitation Reduction
Lower precipitation resulted in lower simulated biomass and LAI overall, but in varying degrees among PFTs. The lycopsid and M. scheuchzeri PFTs exhibited the greatest degree of sensitivity to reduced precipitation (fig. 12) with reduced growth beginning at around 80 percent of the initial precipitation. Both PFTs become non-viable, meaning no biomass was simulated, at a precipitation reduction to 20 percent of the initial input (35.7 cm yr−1). In contrast, the cordaitalean and tree fern PFTs continued to have stimulated growth with precipitation as low as 10 percent of original levels (17.8 cm yr−1). Furthermore, the cordaitalean and tree fern PFTs maintained a constantly high LAI until precipitation was reduced to 60 percent of the initial input level (107.2 cm yr−1).
End-of-simulation values of (A) total vegetation carbon (kg C m−2) and (B) leaf area index (LAI; m2 m−2) for Paleo-BGC-predicted values for each plant type for 50 years of simulations with pCO2 = 600 ppm and pO2 = 0.28 mol mol−1 for individual scenarios with reduced input precipitation (%) from the starting original value of 178.6 cm yr−1. The shaded areas in each graph represent ±one standard deviation from the mean based on the 50 years of simulation for each plant type.
DISCUSSION
Baseline Pennsylvanian Ecosystem Simulations
Water budget differences among paleoplants.
The range of water use by simulated Pennsylvanian plants reflects a range of particular ecophysiological habits, from those that lose water rapidly and require hydraulic resupply to those with low water loss and high resistance to drought. As a general trend, maximum stomatal conductance was inversely proportional to residual annual soil water and the length of the growing season for stands of the different plant types. However, if gsmax was the only important variable driving water use and gas-exchange, we would expect the medullosan PFT to show similar water limitations to lycopsid PFT (1.53 mol H2O m−2 s−1; table 1), as their gsmax values are of the same order of magnitude. Instead, the medullosan PFT behaved more like cordaitaleans and tree ferns in terms of water usage. This showed the influence of leaf hydraulic conductivity (KL) feedback on modeled operational stomatal conductance.
Differences in modeled paleo-plant KL reflects the magnitude and range in average mesophyll path distances measured among the paleo plants types (table 3) accentuated by the non-linear relationship with these distances and KL (for example, equations (3) and (4); table 6). As a result, the slope value (m), calculated from the range of KL divided by the difference in stem water potentials for open and closed stomata (eq 5), showed that medullosans had the lowest (that is most negative) value (−27.43; table 6). From a modeling perspective, this slope means that medullosans, with their large range in mesophyll path distances, were the most sensitive to changes in stem water potential of all plant types simulated (for example, Plant type 1, lower graph, fig. 2B). Despite a higher maximum stomatal conductance, fluctuating stem water content affected calculated daily leaf water potential (ΨL) for medullosans, reducing operational stomatal conductance, depleting soil water, and extending the growing season by maintaining potential values above the wilting point. Tree ferns also had a low m value (−14.99; table 6) indicating greater sensitivity to stem water; however, this PFT also had the lowest gsmax, which limits maximum daily water loss, imparting a different adaptation to environmental water limitation.
Values of maximum and minimum leaf hydraulic conductivity
Both lycopsids and cordaitaleans had higher m values (−2.25 and −1.42 respectively), indicating a leaf water supply strategy comparable to plant types 2 (plants with moderate but similar minimum and maximum mesophyll pathlengths) and 3 (plants with long but similar minimum and maximum mesophyll pathlengths) illustrated in figure 2. For lycopsids, high gsmax with low sensitivity to stem water supply supports rapid water loss. For survival, however, this strategy requires a continual supply of soil water, and perhaps explains (or is explained by) their abundance and ecological dominance in Carboniferous swamp habitats. For cordaitaleans, the combination of low gsmax and low sensitivity to stem water potential suggests an adaptation to drought resistance: low water loss from leaves and internally conservative water resupply. This strategy is consistent with the increased prevalence of cordaitaleans in lowland basins of Euramerica during the seasonally dry interglacials, which occurred on eccentricity timescales and within long-term aridification from the Late Pennsylvanian to the early Permian (Tabor and others, 2008).
Vegetation-specific water-use differences in these simulations indicate that vegetation, in addition to climate, could influence regional environmental conditions. The difference in the surface-water discharge values among the four paleo-PFTs highlights the potential importance of watershed hydrology of both ecophysiology and spatial distribution of paleo-vegetation in the Late Pennsylvanian—early Permian Pangaean landscape (325–290 Ma). By modern example, the African Congo River has an average drainage area of 4 × 106 km2, an annual precipitation of 150 to 200 cm yr−1 (Munzimi and others, 2015), and an average flow near the river outlet to the Atlantic Ocean of 1457 km3 yr−1 (Runge, 2008). Using the average surface runoff per unit area predicted from Paleo-BGC, if the Congo Basin were dominated solely by our simulated lycopsids, predicted flow at the mouth of the river would be less than half of the current flow (737 km3 yr−1) for a comparable amount of precipitation. In contrast, was the Basin covered entirely by tree ferns, the estimated simulated flow would be 4-fold higher (2963 km3 yr−1). However, we did find that significant variability in precipitation from the GENESIS simulations particularly to the north of the selected paleolocation (fig. 3C) in which mean annual precipitation was approximately 13 percent of that used for the Paleo-BGC simulations. Future, spatially explicit simulations will allow us to investigate validity of this basinal hydrology projection.
Modeling Pennsylvanian forests as monoculture stands represent end-member values from among the range of actual possibilities. However, ecosystems of the Euramerican tropical basins, reconstructed using a rich paleobotanical record (Pfefferkorn and Thomson, 1982; DiMichele and Phillips, 2002; Cleal and Thomas, 2005; DiMichele and others, 2009; DiMichele, 2014; Bashforth and others, 2016b; DiMichele and others, 2017) though dominated by either lycopsids or tree ferns, were composed of mixed plant communities (Willard and others, 1995). When the freshwater surface runoff values with all PFTs present together are extrapolated to an extent of 2.0957 × 106 km2 (the estimated average drainage area of the Central Appalachian Basin during the Middle Pennsylvanian and when all the modeled PFTs were present together) (Archer and Greb, 1995) the resulting values of 671 to 1732 km3 yr−1 are similar to that estimated for precipitation-driven discharge volume of 800 to 1500 km3 yr−1 into the Late Pennsylvanian Midcontinent Sea of central tropical Euramerica (Algeo and Heckel, 2008). Notably, simulation results show how the replacement of one dominant taxon would have substantially impacted the surface hydrology of tropical Euramerica.
The major vegetation turnover across the Middle to Late Pennsylvanian boundary (∼306 Ma), which included the loss of most lycopsids and rise to dominance of tree ferns (Phillips and Peppers, 1984; DiMichele and others, 2009; Montañez, 2016), based on our model results would have increased freshwater surface runoff potentially across broad regions of the Euramerican paleo-tropics. Physiological differences between the dominant PFT may, therefore, have led to higher stream flow and during more seasonal periods of the glacial-interglacial cycles. At times, this would have resulted in increased sediment yield and stronger silicate weathering (for example, increased weatherability, discussed later). This major vegetation turnover was coincident with intensifying pantropical aridification, which began in the mid-Pennsylvanian and continued through the early Permian (Tabor and others, 2008; DiMichele and others, 2009; Tabor and others, 2013). For example, vegetation-driven decrease in transpiration in response to the shift to lower ecosystem-scale stomatal conductance would have led to reduced recycling of water vapor to the continental interior possibly intensifying the aridification (compare Ibarra and others, 2019). Given the substantial differences in water use, maximum stomatal conductance, and water-use efficiency between the paleo-PFTs, it is likely that vegetation feedbacks to the climate (that is, decreased precipitation and increased seasonality) occurred. These hypothesized vegetation-climate effects and feedback need to be further tested through integrated climate and fossil data-ecosystem coupled modeling.
Nitrogen cycles in Pennsylvanian vegetation.
The nitrogen budgets of Pennsylvanian forests are a relatively unexplored topic; however, they are of particular importance in terms of linking terrestrial and aquatic (for example, marine) systems through eutrophication potential and ocean productivity. Nitrogen mineralization rates, an indicator of nitrogen recycling, were lower in these simulations than observed rates for modern tropical rain forests (91.3–186.9 kg N ha−1 yr−1; Neill and others, 1995; Wolf and others, 2011) for both paleo- and modern evergreen flora (including broad and needle-leaved). Rather, simulated paleo-flora N mineralization rates were more similar to those observed from existing temperate and boreal conifer forests (12.5–120.0 kg N ha−1 yr−1; Reich and others, 1997; Aber and others, 1998; Mattsson and others, 2003; Kielland and others, 2006). While no study has verified nitrogen cycle predictions from BIOME-BGC for tropical forests, we found that estimated net nitrogen mineralization rates were underestimated by BIOME-BGC for modern dry oak forests in Spain. We attribute our lower values to the seasonal soil water availability as a combined effect of seasonal precipitation and high water use of the paleoplants that led to moisture limitations on modeled soil respiration. Between the two modern taxa simulated, the N mineralization rates of the paleo-PFTs were found to be similar comparable in magnitude to those of the evergreen broad-leaved forest PFT (≈40.0 kg N ha−1 yr−1) rather than the evergreen needle-leaved forest PFT (≈30.0 kg N ha−1 yr−1; table 4). This reflects that measured C:Nleaf of the paleoflora (table 2) are more comparable to current Amazonian forests (32 kg C kg N−1; Luizão and others, 2004) than temperate coniferous forests (51 kg C kg N−1; McGroddy and others, 2004). We acknowledge that in modern tropical forests, phosphorus (P) is also a potentially limiting nutrient (Turner and others, 2018) but is not included in Paleo-BGC due to the difficulty in accurately representing the geologic distribution of P-bearing rocks.
Ecosystem nitrogen lost through leaching was also lower than expected for a tropical forest (320.0–2400.0 g N ha−1 yr−1; McDowell and Asbury, 1994; Lewis and others, 1999; Taylor and others, 2015) and was more similar to modern temperate and boreal conifer forests (86.0–300.0 g N ha−1 yr−1; Hart and Firestone, 1989; Strader and others, 1989; Aber and others, 1998; Pérez and others, 1998). Part of this difference may be attributed to the GENESIS-simulated precipitation; higher rainfall ultimately leads to higher N losses. Yasin, Kroeze, and Mayorga (2010) estimated an average dissolved inorganic N flux from the Congo basin of 805 g N ha−1 yr−1 with a range of total basin load of 244 to 400 Pg yr−1, although this includes some fraction of anthropogenic inputs of waste and fertilizer. Again, using a modern example, the calculated loading from a lycopsid-dominated Congo Basin would be <10 percent of this amount or 18 Pg yr−1. If the basin were dominated by marattialean tree ferns, the N loading would be on average 56 Pg yr−1. Such low nitrogen flux from late Paleozoic terrestrial systems, as compared to the modern, would have negatively impacted marine productivity particularly in tropical marginal or epicontinental seas. This is consistent with the inferred low productivity and the documented allochthonous character of most organic matter in the Late Pennsylvanian North American Midcontinent Sea (Algeo and Heckel, 2008). However, the 3-fold increase in simulated N leaching rates between lycopsids and tree ferns suggests that the loss of lycopsids and rise to dominance of tree ferns in the Euramerican tropics across the Middle-Late Pennsylvanian boundary (DiMichele and others, 2009; Montañez and others, 2016) had the potential to increase not only surface water discharge to the Midcontinent Sea, but also to terrigenous N export. These results suggest that vegetation physiology had the potential to influence increased marine productivity and widespread benthic anoxia that has long been recognized in late Paleozoic epicontinental seas, but is not well understood (Algeo and Heckel, 2008; Montañez and others, 2018).
Productivity of Pennsylvanian vegetation tied to water and nitrogen availability.
Maximum stomatal conductance and CO2 assimilation rates across modern global vegetation have been found to broadly correspond across very different plant types and climates (Schulze and others, 1994). However, despite an order of magnitude difference among paleo-PFT input gsmax values (table 1), we found daily average assimilation rates to be quite similar, reflecting a within-climate variability of different vegetation. For example, gsmax ranged from a high of 3.66 mol H2O m−2 s−1 for the lycopsid to a low 0.34 mol H2O m−2 s−1 for the tree fern, yet daily average assimilation values ranged from 9.41 and 7.72 μmol CO2 m−2 s−1, for these same plant types, respectively. This demonstrates that realized assimilation is manifested within an ecosystem through interactions of plant water use, soil water availability, nutrient recycling, and leaf nutrition. Notably, such organism-ecosystem feedbacks suggest that inferring ecosystem processes from a single plant trait variable leads to biased results. Due to lycopsids' high gsmax and low leaf hydraulic conductivity, our simulated arborescent lycopsids would have been considered to be low productivity plants were they to have grown on terra firma soils that experienced only limited moisture recharge. In such habitats, lycopsid trees literally might “drink themselves dry”. However, occurrence data suggest that nearly all arboreal lycopsid, with the exception of specialized sigillarians (Pfefferkorn and Wang, 2009), were confined to swamp habitats with nearly constant water availability (Cleal and Thomas, 2005; DiMichele, 2014).
Gas-exchange is only part of the story as leaf nutrients also affects photosynthetic capacity, highlighting the importance of utilizing measured C:N for our different PFTs. As an example, M. scheuchzeri had the highest daily average assimilation rate of 10.08 μmol CO2 m−2 s−1, though only the second highest gsmax (1.53 mol H2O m−2 s−1). The key to the high assimilation of M. scheuchzeri is having the lowest input C:Nleaf value (29.6 kg C kg N−1), leading to the highest calculated maximum theoretical photosynthetic rate (Vcmax25) of all PFTs (31.66 μmol CO2 m−2 s−1).
However, having a high average assimilation rate does not necessarily translate to higher simulated biomass as the modeled results demonstrate. Rather, initial C:Nleaf for each plant type (table 2) is proportional to total vegetation production (table 5) indicating that the plants with the nitrogen-rich leaves tend to produce less biomass. This seems counter-intuitive but is explained by considering the processes of carbon metabolism and allocation as they relate to nitrogen. Plant tissue nitrogen concentration is positively correlated with maintenance respiration rates as modeled in Paleo-BGC (Ryan, 1995). Following the previously considered example, as a result of their low input C:Nleaf M. scheuchzeri had the highest tissue respiration rate (1.61 μmol CO2 m−2 s−1 kg biomass C−1 ; including both maintenance and growth) as compared to tree ferns (1.39 μmol CO2 m−2 s−1 kg biomass C−1), lycopsids (1.31 μmol CO2 m−2 s−1 kg biomass C−1), and cordaitaleans (0.92 μmol CO2 m−2 s−1 kg biomass C−1). Despite having the highest Vcmax25, M. scheuchzeri had higher carbon loss and, therefore, less carbon to allocate to biosynthesis of leaves, stems, and roots.
Biomass allocation is also tied to C:Nleaf based on tissue sink-strength. For total biomass, the stem proportion is the largest component; therefore, if less model carbon is allocated to stems, it is either due to lower net carbon availability or high tissue-specific nitrogen requirement, which, under low nitrogen availability, leads to overall lower total biomass. For the four paleo-PFTs, we did not specify a specific leaf, root, and stem allometry but rather allow the model to apportion these based solely on the prescribed C:Nleaf. We note that the medullosans were an ecologically diverse group, some growing in the understory, others forming low canopies in open, sunnier microhabitats. Of these diverse forms, M. scheuchzeri trees were considered to be comparatively short, with leaves a meter or somewhat more in length (Laveine and Belhis, 2007) suggesting that they were understory members of tropical wetland forest stands (Wnuk and Pfefferkorn, 1984; Wnuk and Pfefferkorn, 1987; Falcon-Lang, 2009) with larger lycopsids and cordaitaleans in the upper canopy. These paleobotanically-inferred architectural characteristics agree with our biomass simulations, wherein, for M. scheuchzeri, with its low input C:Nleaf, required more nitrogen for biomass construction than in the other PFTs. In contrast, lycopsids and cordaitaleans had the highest simulated stem carbon (table 4) as higher input C:Nleaf for these plant types with less nitrogen required for carbon allocation allowing for larger biomass relative to the other PFTs. The low nitrogen requirement by lycopsids is also is compatible with the interpreted oligotrophic swamp habitat of lycopsids inferred from sedimentological analysis (Cleal and Thomas, 2005; DiMichele, 2014).
The Effect of Elevated pO2
Higher atmospheric pO2 has been suggested to reduce primary productivity based on the effect of increased photorespiration, the competitive oxygen fixation associated with the photosynthetic enzyme Rubisco (Beerling and Woodward, 1997). This hypothesis, however, does not consider the possible mitigating effect of increased atmospheric pressure associated with higher pO2 on evaporation physics. As the second-largest gas component of the Phanerozoic atmosphere, oxygen could not be added or removed without changing the pressure and other properties of Earth's troposphere (for example, Poulsen and others, 2015; Wade and others, 2019). Thus, because high atmospheric oxygen of the Pennsylvanian likely increased atmospheric pressure, we included the physics of changing pressure on the mechanisms associated with evaporation, as the most likely effect of elevated pO2. For these sets of simulations, higher pO2 did reduce average daily net assimilation of CO2 (An) slightly, though also suppressed transpiration by plants thereby increasing soil water throughout the growing season.
Although elevated oxygen (0.28 mol mol−1) reduced plant productivity, the effect was substantially less than expected. The negative photorespiratory effects of elevated O2 were counter-balanced by the increased number of days that were physiologically conducive for stomatal opening (62% for pO2 = 0.21 mol mol−1 compared to 64% for pO2 = 0.28 mol mol−1) due to pressure-induced reduction in transpiration and increase in soil water. In this way, elevated O2 lengthened the growing year—an important environmental change, particularly for lycopsids. Despite overall lower daily assimilation (fig. 11), in fact, the lycopsid PFT had slightly higher biomass under elevated pO2 (table 4).
Despite modeled physical suppression of water loss due to elevated oxygen, WUEi values for all paleo-PFTs are lower for simulations with elevated pO2 compared to those with present-day pO2, with the greatest modeled change for cordaitaleans and tree ferns (fig. 10). This indicates that under elevated pO2, photorespiratory effects on CO2 assimilation (An), while balanced, are not completely compensated for by reduced transpiration.
In addition to productivity difference, reduced evaporation allowed for higher soil water contents and surface water discharge rates for all paleo-PFTs and subsequent increased N leaching rates. Surface water discharge increased by 11 percent for tree ferns and nearly doubled for lycopsids in elevated pO2 simulations, corresponding to proportional increases in N loss through leaching.
Modeled Drought and Vegetation Response
The mechanistic, process-based modeling results presented here are consistent with previously proposed turnovers in Late Pennsylvanian community compositions driven by changes in water availability inferred from paleobotanical and taphonomic records (for example, Falcon-Lang, 2004; DiMichele, 2014; Willis and McElwain, 2014). Specifically, due to combined effects of high transpiration and low leaf sensitivity to water loss our simulations indicate that small changes in site water budgets would have had a larger impact on the arborescent lycopsids, which dominated large parts of Early and Middle Pennsylvanian wetland forests, than the other three modeled taxa. The simulation of lycopsids under decreasing precipitation further points to specific, and perhaps novel mechanisms for their drought intolerance.
The high gsmax of lycopsids, even under elevated pO2, indicates they could have used large amounts of water. High gsmax coupled with long mesophyll path length (Dm), however, were disadvantageous foliar characteristics for lycopsids. Lycopsids would have experienced more negative leaf-water potential than other taxa due to leaf water transport and distribution constraints. In order for lycopsids to maintain sufficient leaf-water potential and stomatal conductance for requisite CO2 assimilation to obtain ‘viable’ biomass (fig. 12), landscapes with persistent water sources like saturated lowlands were likely key to their survival. Furthermore, if vegetation turnover occurred in adjacent sites, then this would have reinforced these wet refugia with added freshwater runoff. Although this environmental scenario has been long proposed based on the paleobotanical record (for example, Pfefferkorn and Thomson, 1982; Gastaldo, 1987; DiMichele and others, 2002; Cleal and others, 2012; DiMichele, 2014; Bashforth and others, 2016b; DiMichele and others, 2017) this is, to our knowledge, the first modeled test of lycopsid growth conditions that provides mechanism-based quantitative constraints on their viability under a range of site water conditions.
The simulated changes in precipitation evaluated in this study indicate that lycopsids and medullosans (represented by M. scheuchzeri) required higher mean annual precipitation than cordaitaleans and tree ferns to maintain growth (fig. 12). However, because we simulate relatively low precipitation thresholds (35.7 cm yr−1) for the lycopsids and medullosans as compared to cordaitaleans and tree ferns (17.8 cm yr−1), we regard these thresholds as only the absolute, not realized, physiological limit for these plants. Such low precipitation amounts define modern arid ecosystems and are unlikely to have supported the Late Pennsylvanian vegetation we modeled. Notably, these results, (that is, cordaitaleans and tree ferns were more drought tolerant than lycopsids and medullosans) are consistent with a plethora of paleobotanical evidence. Aridification, which intensified through the Late Pennsylvanian, is inferred to have impacted global hydrologic cycling given geologic and modeling evidence for a shift during this time from less seasonal and overall higher effective moisture during late glacials and earliest interglacials (referred to here as ‘glacials’) to more highly seasonal climate during interglacials and early glacials (referred to here as ‘interglacials’) (DiMichele and others, 2009; Horton and others, 2012; Rosenau and others, 2013; Cecil and others, 2014). The lycopsids and medullosans, which dominated the lowland tropical landscape during glacials, were repeatedly locally or regionally replaced during interglacials by a suite of plants that included drought-tolerant cordaitaleans, conifers, taeniopterids, and other taxa (Bashforth and others, 2014; DiMichele, 2014; Bashforth and others, 2016a). Moreover, certain wetland vegetation members were able to locate and thrive in fragmented microhabitats with high soil moisture, such as tree ferns and calamitaleans (DiMichele and others, 2006; Falcon-Lang and Dimichele, 2010). The disadvantageous foliar characteristics (high gsmax and high Dm) of the lycopsids likely led to the retreat of these taxa during interglacials to fragmented microhabitats within lowland basins, where adjacent topographic impoundment areas allowed adequate site water supply. Consequently, the return of the lycopsid-dominated forests represents a reordering of these ecosystems in concert with the return of increased rainfall and of widespread humid climatic conditions in regions of the tropics, where there were adequate lycopsids surviving in refugia to provide reproductive stock.
We further hypothesize that a precipitation threshold, or thresholds, led to a tipping point that began a cascade of regional and local environmental changes that drove the repeated vegetation turnovers during the latest Mississippian, Pennsylvanian, and earliest Permian. Restructuring of the paleo-tropical terrestrial ecosystems of Euramerica occurred during long-term pantropical aridification that intensified through the Pennsylvanian and into the early Permian. A precipitation threshold hypothesis is particularly relevant for the major vegetation turnover across the Middle-Late Pennsylvanian boundary (∼306 Ma), which involved the permanent loss of most arborescent lycopsids from paleo-tropical Euramerican wetlands and the rise to dominance of tree ferns in those habitats. This has been widely hypothesized to have occurred as a result of intensifying aridification (Phillips and Peppers, 1984; DiMichele and others, 2009; Falcon-Lang and others, 2018). During this event, the lycopsids and, perhaps many of the medullosans, would have become physiologically non-viable from a plant water balance perspective as precipitation thresholds were reached.
Our ecosystem simulations are primarily based on the foliar characteristics of the now-extinct relatives of early-diverging vascular pteridophytes and seed plants, which dominated late Paleozoic tropical terrestrial ecosystems. We fully appreciate, however, that simulated whole-plant characteristics are critical for robust estimation of the ecophysiological functioning and biogeographical distributions of these early vascular plants (compare Wilson and others, 2017). A modeled representation of the complete hydraulic pathway of these plants is required to better delineate the precipitation threshold amounts required to avoid critical hydraulic failure in stems and thereby tree mortality. That is, canopy water supply is constrained by the risk associated with stem and root xylem embolism occurring at a daily time-scale (Sperry and others, 2016). Further model development will include aspects of stem and root hydraulic conductivity, particularly those structurally important characteristics that can be derived easily from fossil specimens (Wilson, 2013; Wilson and others, 2015; Wilson and others, 2017). Secondly, there are undoubtedly differences in drought resistance within the taxonomic groups in this study that remain unaddressed. Increased taxonomic resolution and additional fossil analysis will improve precision: median gsmax values were derived for each PFT based on many fossils that collectively represent a few genera. Therefore, detailed site water fidelity—such as inhabitation of freshwater and brackish wetlands among lycopsids, as described by Thomas and Dimitrova (2017), were not addressed. Lastly, we acknowledge that the Paleo-BGC model does not represent other important biological mechanisms of plant survival and mortality, including reproduction, defense, and other factors such as tolerance to waterlogged soil and/or salinity stress.
Effects of Late Paleozoic Vegetation on Silicate Weathering, and CO2 Consumption
The evolution of vascular plants in the mid-Paleozoic has long been hypothesized to have impacted continental silicate weathering flux, and, in turn, atmospheric CO2 through increased CO2 consumption (Berner, 1991; Berner, 1992; Algeo and Scheckler, 1998; Beerling and Berner, 2005; Boyce and Lee, 2011; McMahon and Davies, 2018). This impact of vascular plants on silicate weathering has been estimated as a 7 to possibly 10-fold amplification relative to that of non-vascular plants, which covered earlier landscapes. Recent empirical and modeling studies have argued that the efficiency of silicate weathering and the regulation of CO2 consumption are primarily linked to the intensity of hydrologic processes, of which surface runoff is a key component (Dessert and others, 2001; Godsey and others, 2009; Maher, 2010; Maher, 2011; Maher and Chamberlain, 2014; Ibarra and others, 2016). In the context of this linkage, the 105-yr-scale shifts in dominant plant types across the landscape during glacial-interglacial cycles, and on the longer-term (106-yr) through the Pennsylvanian, had the potential to affect silicate weathering and thus CO2 consumption. This reflects two important characteristics of the paleotropical forests. The first is the aerial extent (up to 2400 × 103 km2) of the paleo-forests (Cleal and Thomas, 2005) and their vast expanse throughout tropical Euramerica and the North China Block. The second is the notable differences in modeled water-use and soil-water amounts between the extinct vascular seed plants (medullosans) and lycopsids (high water use), which dominated the landscape during wet periods, and the pteridophytes and cordaitaleans (and conifers) (low water use), which dominated during drier and more seasonal intervals, and the consequent up-to four-fold differences in simulated surface-water discharge rates revealed by this study. Low water-use plants (that is, dryland plants) that dominated the tropical landscape during drier and more seasonal periods would have lived in soils with relatively high volumetric soil water content. With up to ∼4 times higher surface-water discharge rates associated with dryland plants coupled with less dense vegetation coverage during these drier periods (DiMichele and others, 2009; DiMichele, 2014), whether on the 105- or 106 yr-scale, could have induced greater local- to regional-scale soil loss due to increased runoff. This potentially would have exposed more bedrock and fresh minerals to weathering (Drew, 1983; Istanbulluoglu and Bras, 2005). This increased surface water discharge, which undoubtedly varied seasonally, however, could have been sustained through the duration of the drier intervals—over 10s of thousands to millions of years thus impacting silicate weathering rates and global CO2 consumption.
This potential vegetation-driven influence on silicate weathering would have been amplified during the aforementioned ‘threshold’ vegetation turnovers given the associated large-scale changes in vegetational composition and architecture (DiMichele and others, 2009; DiMichele, 2014; Falcon-Lang and others, 2018). In contrast, overall surface water discharge rates would have been low during intervals when the paleo-tropical forests were dominated by the wetland communities rich in lycopsids and medullosans (for example, eccentricity scale glacials; the first half of the Pennsylvanian). Moreover, silicate weathering in standing water bodies (peatlands), which may have persisted for 100s of years or more, would have been temporally limited given the short fluid transport times (reactive flow path length normalized for porosity/runoff amount) and long residence times relative to the time required for silicate mineral weathering reactions to reach equilibrium (compare Maher, 2011; Maher and Chamberlain, 2014). Overall, the water-use and soil water attributes of the plant dominants during drier and more seasonal intervals of the Pennsylvanian had the potential to make silicate weathering more efficient in the tropics (increased ‘weatherability’) and thus increase chemical denudation rates (compare Maffre and others, 2018; Maher, 2010). In turn, more efficient CO2 consumption by increased silicate weathering in the warm tropical climate belt would have promoted a tighter coupling between atmospheric CO2 and climate, thus increasing the strength of the silicate weathering negative feedback during eccentricity to million year-scale periods of drier and more seasonal conditions (Ibarra and others, 2016).
Undoubtedly, the Central Pangaean Mountains, which spanned 1000s of km of the paleotropics, would have been a primary control on continental silicate weathering flux and CO2 consumption during the late Paleozoic given the effects of relief on erosion, soil thickness and fluid residence time and length of fluid transport (Jacobson and Blum, 2003; Jacobson and Blum, 2003; Riebe and others, 2004; Torres and others, 2005; West and others, 2005; Maher and Chamberlain, 2014; Caves and others, 2016; Goddéris and others, 2017). Importantly, our simulations indicate that the presence or absence of different dominant plant types on the Pennsylvanian—early Permian (325–290 Ma) landscapes had the potential to greatly alter surface water discharge rates in broad regions of the lowland paleotropics as well. This is particularly important for the supercontinent Pangaea given that expansive continental interiors are typically characterized as having low effective moisture and thus low runoff rates and lowered weatherability. Furthermore, it is feasible that the high water-use and transpiration rates of wetland plant dominants, the lycopsids and medullosans, modeled in this study could have expanded the aerial extent of the Pangaean tropics accessible to silicate weathering through increased downwind transport of water vapor to continental interior (compare Boyce and Lee, 2011; Ibarra and others, 2019). In turn, this influence would have been dampened during times of predominance of cordaitaleans, tree ferns, and conifers in the tropical lowlands (drier intervals) and overall, in the drier climate of western tropical Pangaea. At times when dryland PFTs dominated the Pangaean lowlands, recycling of water vapor to Pangaean interior regions would have been reduced due to their higher water-use efficiency and reduced transpiration rates.
Our elevated pO2 (0.28 mol mol−1) simulations indicate additional mechanisms for increasing surface water discharge rates and potentially the silicate weathering flux through the Pennsylvanian and early Permian. All four of the floral dominants of the paleo-tropical forests that we modeled would have had higher soil-water contents with elevated atmospheric pO2 due to the atmospheric pressure-induced reduction in transpiration. In turn, decreased transpiration rates under elevated pO2 translate to higher surface-water discharge rates for all four PFTs in comparison to under present-day (21%) pO2 conditions. The change in surface-water discharge rate for a pO2 of 28 percent ranges from an 11 percent increase for tree ferns to a doubling of discharge for lycopsids, all with concurrently higher modeled soil leaching rates. Overall, the intensifying aridification through the Pennsylvanian and early Permian, the aforementioned repeated vegetation shifts, and the long-term increase in atmospheric pO2 (from ∼ present-day in the Early Pennsylvanian to 26–30% by the end Pennsylvanian-early Permian), may have counteracted the dampening effects of supercontinentality on surface hydrology and silicate weathering, given that our simulations indicate that increased pO2, along with the effect of changing vegetation type, had the potential to increase soil water contents and higher surface water discharge rates.
SUMMARY
This study found evidence for significant differences in physiological function of four dominant Late Pennsylvanian plant types, lycopsids, medullosans, cordaitaleans, and tree ferns, based on process-based ecosystem modeling with fossil-based input parameters. Rather than using modern proxies for the fossil plants, measurements of the dominant leaf characteristics were made directly from fossilized remains. These measurements include stomatal conductance, C:N, leaf size, average mesophyll path distances. Major ecophysiological differences were observed among the four PFTs. Lycopsids were more water-use intensive, produced lower nitrogen leaching, and had moderate growth when compared the tree ferns, which exhibited low water use, high nitrogen loss, and lower growth rates.
The two different atmospheric compositions modeled represent different parts of the late Paleozoic, one with present-day levels of atmospheric oxygen, and one with elevated oxygen levels hypothesized to characterize the latter part of the Pennsylvanian. With elevated oxygen levels, the modeled physiological characteristics of the PFTs and their effects on the Earth System were modified. All PFTs transpired more water with elevated oxygen because of higher atmospheric pressure; however, elevated oxygen had a limited effect on overall plant productivity. This modeled result is compatible with the independently inferred ecological distributions of Pennsylvanian plants (Cleal and Thomas, 2005; Bashforth and others, 2014; DiMichele and others, 2017; Thomas and Dimitrova, 2017), and with the repeated major turnovers of wetland vegetation thought to characterize the Pennsylvanian (∼315–300 Mya) tropics. These turnovers resulted in the permanent establishment of seasonally dry vegetation as the sole biome in tropical lowlands by the early Permian. This biome establishment was accompanied by the diminished prevalence of certain wetland floral lineages, which were relict to isolated high-moisture microhabitats (Montañez, 2016). Consequences of these turnovers may have included reduced water use from plant communities, higher silicate weathering, and transport of nitrogen from land to sea, the latter of which could have stimulated marine productivity.
Our simulations demonstrate that the ecological differences among plants are not due to any single attribute, but rather are associated with multiple traits in their anatomy and ecophysiological responses that lead to feedbacks in water use, nutrition, and assimilation within a dynamic environment. Modeling extinct plant ecosystems necessitates the use of fossil data. Our approach presents a novel perspective on paleobotanical data as a previously untapped wealth of information about the late Paleozoic and other Eras, made available by process-based ecosystem modeling. Paleo-ecosystem models that use modern plants as analogs are not truly representative of late Paleozoic ecosystems. Given the results of recent studies, different combinations of paleo-plants traits affect their physiological function and render them distinct from extant groups (Wilson and others, 2015; Wilson and others, 2017; Wilson and others, 2020). Representing these differences in Earth System models is important to understanding biogeographical distributions and erosional potential of land area over the course of geologic time.
SOURCE CODE
The original BIOME-BGC model version 4.2 is available from the Numerical Terradynamics Simulation Group at the University of Montana at the following address: https://www.ntsg.umt.edu/project/biome-bgc.php. The source code for the Paleo-BGC model presented here, in addition to the climate and ecophysiological files for each vegetation type can be found at: https://github.com/josephdwhite/paleo-bgc.
ACKNOWLEDGMENTS
The authors are grateful for the foundation, development, and availability of the BIOME-BGC model specifically provided by Steve Running, Peter Thornton, Ramakrishna Nemani, John Kimball, Mike White, other past and current members of the Numerical Terradynamic Simulation Group at the University of Montana, USA. We also thank the insightful and helpful reviews of an earlier version of this manuscript by Prof. Danny Rye, Yale University, and Dr. Daniel Ibarra, University of California Berkeley, and the anonymous reviewer. The data and findings presented in this paper were supported by NSF grants EAR-1338247 (J.D.W.), EAR-1338281 (I.P.M.), EAR-1338200 (C.J.P.), and EAR-1338256 (M.T.H.), and ERC-2011-StG and 279962-OXYEVOL to J.C.M.
APPENDIX A
Variables used as part of this study related to simulations utilizing the Paleo-BGC model
APPENDIX B
ESTIMATING MESOPHYLL PROPERTIES FROM FOSSIL LEAVES
Mesophyll Conductance
This section describes the basics of estimating mesophyll conductance (gm; m s−1) from anatomical features. Derivation of, gm is from a 1-D diffusion model (Niinemets and Reichstein, 2003) where gm is composed of conductance affecting diffusion of CO2 through leaf inner air spaces (gias; m s−1) and cellular liquid (gliq; m s−1):
where H is the equilibrium gas/liquid-phase distribution coefficient for CO2 (Henry's law constant, 2948.6 Pa m3 mol−1 at 298 K), R is the gas constant (8.314 J K−1 mol−1), and Tk is the temperature (Kelvin). The inner air space, gas-phase conductance (gias) is:
where Da is the diffusion coefficient for CO2 in the gas phase (1.51 × 10−5 m2 s−1 at 25 °C), fias is the fraction of leaf air space, ΔLias is the average gas-phase thickness assumed to be half of the mesophyll thickness in m (Tomás and others, 2013), and ς is the diffusion path tortuosity (m m−1). The path tortuosity has been estimated at 1.57 m m−1 for Macademia sp. (Syvertsen and others, 1995). The value of gliq (m s−1) is estimated by:
where Sc is the surface area of chloroplast, S is the leaf surface area, and the subscripts of resistance (r) of cw, pl, cyt, en, and st are for cell wall, plasmalemma, cytosol, chloroplast envelope, and stroma, respectively. For rpl and ren, a value 281.7153 s m−1 has been estimated for the phospholipid bi-layer membrane (Evans and others, 1994; Tosens and others, 2012b). The value of Sc/S is likely impossible to be derived from fossilized leaf sections. However, Tosens and others (2012b) found a conservative relationship between Sc/S and Sm/S (fig. B1) where Sm is the surface area of mesophyll cells, a feature that can be identified and measured in some well-preserved specimens.
Calculated mesophyll conductance values to both carbon dioxide and water vapor from morphological values derived from various literature sources
Average mesophyll surface area exposed to intercellular air space per unit leaf area (Sm/S, m2 m−2), chloroplast surface area exposed to intercellular air space per unit leaf area (Sc/S, m2 m−2) for leaves of Populus tremula exposed to varying levels of water and light (Tosens and others, 2012a).
The values of rcw, rcyt, and rst are evaluated from the following general equation for liquid phase resistances (ri; s m−1):
where ΔLi is the pathlength (m), Dw, is the aqueous phase diffusion of CO2 taken (1. 79 × 10−9 m2 s−1 at 25°C) rf,I is the relative decrease in Dw relative to free diffusion in water, and pi is the effective porosity (m3 m−3). For rcyt, and rst, the value of pi can be assumed to be equal to 1.0 and the value of ri,f set to 0.3. For rst, the value of ΔLst is approximated as half the thickness of the chloroplast with an average value of 1.65 μm (Oguchi and others, 2003). For rcyt, the value of ΔLcyt may be derived from leaf specimens as the 0.5 of the mean distance within mesophyll cells. For rcw, rcyt,f may be assumed to be equal to 1.0. The value pi can vary with cell thickness based on an empirical formula developed by Tosens and others (2012b) where pi = 0.1775 – 0.3274 (Tcw), where Tcw is the cell wall thickness (μm), a value derived from microscopy of fossil leaf section cross-sections. Final values of gm can be converted into molar units:
where TL is the leaf temperature (°C), P is atmospheric pressure at the location (Pa), and Pmsl is the pressure at sea level (Pa).
Ideally, these models, using assumed values presented, may be populated with fossil measurements of four key attributes including: fias, the fraction of leaf air space, ΔLias, the average half of the mesophyll thickness, ΔLcyt, the mean mesophyll cell width, and Tcw, the cell wall thickness. An example of these measurements is shown from a cross-section of an Alethopteris sp. pinnule from the upper Pennsylvanian (Mickle and Rothwell, 1982) (fig. B2).
Alethopteris sp. cross-section copied from Mickle and Rothwell (1982), Plate 1, #9, ×30. Shown are potential measurement opportunities including fias, the fraction of leaf air space, ΔLias, the average half of the mesophyll thickness, ΔLcyt, the mean mesophyll cell width, and Tcw, the cell wall thickness.
Mesophyll conductance (gm) values were calculated based on morphological variables from images of leaf transverse sections from the literature for Carboniferous species. The average gm was 0.436 mol CO2/m2/s. Using reported values from Niinemets and others (2009) average values for modern angiosperm and gymnosperm species was 0.192 mol CO2/m2/s. Distribution of Niinemets and others (2009) data and those calculated for extinct species are shown figure B3.
Distribution of mesophyll conductance and modern plants derived from Niinemets and others (2009) and paleoplants.
Estimating Specific Leaf Area from Fossil Leaves
Utilizing data from a study on native New Zealand tree species collected within canopies of intact primary forests (White and Scott, 2006), linear regression models were constructed for podocarp and tree ferns, separately with leaf carbon to nitrogen ratios (C:N) as the independent variable and specific leaf area (expressed as m2/kg C) as the dependent variable (fig. B4).
Data from White and Scott (2006) are shown with linear regression models of C:N and SLA for both podocarps and tree ferns from that data set. Also shown are the other plants samples (angiosperms) for reference.
