A dynamic local-scale vegetation model for lycopsids (LYCOm v1.0)

Lycopsids represent a distinct lineage of vascular plants with a long evolutionary history including numerous extant and extinct species that started out as tiny herbaceous plants and later went on to grow into forests with tree-like structures. They enriched the soil carbon pool through newly developed root-like structures and promoted soil microbial activity by providing organic matter. These plants enhanced soil carbon dioxide (CO₂) via root respiration and also modified soil hydrology. These effects potentially promoted the dissolution of silicate minerals, thus intensifying silicate weathering. The weathering of silicate rocks is considered one of the most significant geochemical regulators of atmospheric CO₂ over a long (hundreds of thousands to millions of years) timescale. The motivation for this study is to achieve an increased understanding of the realized impacts of vascular plants, represented by modern relatives of the most basal plants with vascular tissues and shallow root systems, on silicate weathering and past climate. To this end, it is necessary to quantify physiological characteristics, spatial distribution, carbon balance, and the hydrological impacts of early lycopsids. These properties, however, cannot be easily derived from proxies such as fossil records. Hence, as a first step, a process-based model is developed to estimate net carbon uptake by these organisms at the local scale. The model includes key features such as the distribution of biomass above and below ground, along with a plausible root distribution in the soil affecting water uptake by plants. The stomatal regulation of water loss and its immediate implications for photosynthesis are considered. Moreover, root respiration plays a crucial role in the model by affecting soil carbon dioxide and weathering rates. The model features ranges of key physiological traits of lycopsids to predict the emerging characteristics of the Lycopsida class community under any given climate by implicitly simulating the process of selection. In this way, extinct plant communities can also be represented. In addition to physiological properties, the model also simulates weathering rates using a simple limit-based approach and estimates the biotic enhancement of weathering by these plants. We run the Lycopsid model, called LYCOm, at seven sites encompassing various climate zones under today's climatic conditions. LYCOm can simulate realistic properties of lycopsid communities at the respective locations and estimates values of net primary production (NPP) ranging from 126 to 245 g carbon m⁻² yr⁻¹. Our limit-based weathering model predicts a mean chemical weathering rate ranging from 5.3 to 45.1 cm ka⁻¹ of rock with lycopsids varying between different sites, as opposed to 0.6–8.3 cm ka⁻¹ of rock without these plants, thereby highlighting the potential importance of such vegetation at the local scale for enhancing chemical weathering. Our modeling study establishes a basis for assessing the biotic enhancement of weathering by lycopsids at the global scale and also for the geological past. Although our method is associated with limitations and uncertainties, it represents a novel, complementary approach towards estimating the impacts of lycopsids on biogeochemistry and climate.

1 Introduction

The increasing carbon dioxide (CO₂) content of the atmosphere is one of the most discussed topics in climate research and pertains to rising concerns over the aggravating impacts of its greenhouse effect. Since societal well-being is closely connected to the level of CO₂ in the atmosphere, it is important to look at past developments and understand the key events which ensued changes in the level of CO₂ in the atmosphere. The fluctuations in atmospheric CO₂ content of the Earth through history have been significant (Berner, 1998; Foster et al., 2017). A rapid drawdown of carbon dioxide is observed at around 420 Ma, around a similar timeline, the advent of the earliest lycopsids recorded (Garratt et al., 1984). Very few attempts have been made to explain this dip, and silicate weathering has often been cited as an important factor contributing to this rapid fall in CO₂ over geologic timescales (Hilley and Porder, 2008; Colbourn et al., 2015) but the influence of these evolving flora has not been explored in detail.

Lycopsids are regarded as one of the earliest forms of extant vascular plants (tracheophytes) with a long evolutionary history (Kenrick et al., 2012; Morris, 2018; Dahl and Arens, 2020). This ancient group of vascular plants (Kenrick and Crane, 1997; Qiu et al., 2006; Wickett et al., 2014) dates back to the Silurian (Garratt et al., 1984; Kotyk et al., 2002) and were initially relatively simple, but later evolved into tree stature in the Middle Devonian (Stein et al., 2012) before surpassing a height of 50 m during the Carboniferous Period (Thomas and Watson, 1976; Taylor et al., 2009). This long transformative history along with a rich fossil record and extant diversity have made studies of lycopsids central to plant evolution. In any case, lycopsids were not only abundant during the Carboniferous Period but became ecologically significant much earlier in the late Silurian–Lower Devonian (Wellman et al., 2013) and later dominated the flora. The dominant role of these plants as terrestrial vegetation in those periods could be attributed to the appearance of root-like structures (Matsunaga and Tomescu, 2016) coupled with stomatal control (Brodribb and McAdam, 2011) for controlling water loss.

Root evolution has been a gradual process and Fig. 1 depicts the transformation from its early elementary stages of appearance in lycopsids to a sophisticated structure comparable to the roots of modern plants during the Upper Devonian (Hetherington and Dolan, 2017; Matsunaga and Tomescu, 2016). It is without a doubt that the transformation took place over a considerable period. In fact, this evolution led to the foundation of modern plants and trees which are characterized by complex tissue structure allowing transport of water and nutrients, thus making high rates of net primary production (NPP) possible. Root-like structures have been observed in fossil records dating back to as early as 413 Ma (Hao et al., 2010) in the Lower Devonian (Gensel et al., 2001). This innovation might have been central to the dominance of flora on the planet and enabled plants to grow in size and overcome the physiological drawback of their dependence on moist habitats (Bateman et al., 1998).

Figure 1 A conceptual diagram of evolution of roots from its earliest stage in the Silurian Period to a sophisticated form during the Devonian.

Roots provided the lycopsids with a key interface for interaction with the terrestrial environment for their nutrient and water supply, as they could utilize the underground reserves which had remained previously untouched. Roots provided an anchorage, which may have been a precondition for the development of the first trees. This effect also increased soil stability, and hence it potentially reduced the erosion rate (Vannoppen et al., 2017). Roots influenced hydrological processes through transpiration due to the evolution of deep roots, which accessed groundwater that would otherwise remain shielded from uptake and/or subject to loss through evapotranspiration from plants according to atmospheric and physiological demand. This potentially could have enhanced weathering (Lee and Boyce, 2010). The other notable consequence of root development is believed to be an enhancement of chemical weathering (Algeo and Scheckler, 1998, 2010), which could have had strong impacts on atmospheric CO₂ content and the global carbon budget (Le Hir et al., 2011).

In general, the weathering regime is often influenced by several factors and plants, and the associated microorganisms have been suggested to have had a substantial impact on weathering, through respiration in the soil, hence increasing the partial pressure of soil CO₂ to several-fold greater than the atmospheric CO₂ level (Kelly et al., 1998; Montanez, 2013). In addition to plant-derived organic matter respired by these microorganisms, root respiration also adds to soil CO₂ content, thus increasing soil acidity and influencing weathering rates (Moulton et al., 2000; Andrews and Schlesinger, 2001; Berner et al., 2003; Ibarra et al., 2019).

Roots were not the only suite of innovations from these kinds of plants: they further developed stomatal control (Kenrick and Strullu-Derrien, 2014; Brodribb and McAdam, 2011) of water loss to be able to photosynthesize under adverse climatic conditions, which might have given lycopsids yet another advantage over their contemporaries. The evolution of stomata followed the innovation of roots closely, at least 420 Ma (Chater et al., 2011; Ruszala et al., 2011). Water loss and NPP are regulated by the opening and closing of stomata (McAdam and Brodribb, 2013). By changing the aperture of the stomata, plants inhibit water loss to prevent desiccation, at the cost of reduced photosynthetic productivity. The combination of adjustable stomata with an internal water transport system, i.e., vascularity, was a turning point in plant evolution that might have enabled vascular plants to invade most terrestrial environments, tolerating water stress and exploiting favorable conditions (Raven, 2002). The development of such traits in lycopsids was gradual and was initially relatively primitive. Brodribb and McAdam (2011) suggested that lycophyte (a broader classification encompassing lycopsids) and fern stomata lacked key responses to abscisic acid and epidermal cell turgor, making stomatal behavior highly predictable. The stomatal regulation of water loss was determined solely by the surrounding environmental conditions due to a poorly developed hydraulic system (Boyce and Lee, 2010), unlike today's plants which have a more complex mechanism driving stomatal control.

Stomatal development primarily impacted and regulated water loss through transpiration, which may have had a significant role in altering the global water cycle, in particular by increasing the land surface area affected by precipitation through the recycling of water vapor (Berner, 1992; Drever, 1994; Ibarra et al., 2019). This may have then indirectly contributed to the alteration of weathering rates via an enhancement of local and global water recycling, as well as increasing the area available for weathering in the interiors of continents (Shukla and Mintz, 1982; Boyce and Lee, 2010, 2017). Moisture recycling likely also increased plant productivity, thus increasing soil CO₂ content across continents. Consequently, the presence of land plants increased the quantity of water in continental interiors available for weathering and increased primary productivity, further contributing to an increase in global weathering.

Silicate weathering is considered one of the most significant geochemical regulators of atmospheric CO₂ over a long (hundreds of thousands to millions of years) timescale (Berner, 1990). Silicate weathering is primarily mediated by the hydrological cycle (Walker et al., 1981; Bluth and Kump, 1994; Maher, 2011; Maher and Chamberlain, 2014; Green et al., 2017; Ibarra et al., 2019) and it serves as a dominant carbon sink for Earth in the long run (Berner, 1998; Bergman et al., 2004; Gibling et al., 2014). A precondition for weathering reactions is the production of weak carbonic acid from rainwater and soil CO₂, which dissociates into bicarbonate (HCO${}_{\mathrm{3}}^{-}$) ions during the dissolution of silicate minerals. The ions subsequently enter streams via runoff and eventually the oceans, where they are used by various organisms to build shells and skeletons. These then serve as a long-term carbon sink after sedimentation. Hence the weathering of silicate rocks corresponds closely to the flux of carbon from the atmosphere into seafloor sediments (Algeo and Scheckler, 1998; Elick et al., 1998).

This silicate weathering feedback might have been significantly influenced over various periods of Earth's history following the further evolution of plants and the development of roots and stomata. The increased soil carbon pool had the potential to promote the dissolution of silicate minerals (Walker et al., 1981). The study of Matsunaga and Tomescu (2016) demonstrates substantial plant–substrate interaction that was already underway in the Devonian, causing weathering by root penetration (Algeo and Scheckler, 1998), and thereby confirming the significance of such an advent of vascular vegetation. Further studies have already claimed an enhancement of chemical weathering by vascular plants to have exceeded six fold as compared to their earlier counterpart, i.e., non-vascular vegetation (Cochran and Berner, 1993), along with improved pedogenesis (Quirk et al., 2015).

Although lycopsids were abundant in the Upper and Middle Devonian, these plants thrived even in the Lower Devonian, but their effects on weathering rates remain unclear. Detailed information on physiology, productivity, spatial distribution, and hydrological properties of these early plants, which is crucial to determine their effects on weathering rates in the root zone, is lacking. This makes it difficult to estimate the impacts of early lycopsids on the climate of the past. Lycopsids today are primarily represented by Lycopodiaceae (club mosses), Selaginellaceae (spike mosses), and Isoetaceae (quillworts) and are typically found as undergrowth beneath the forest canopy (Spencer et al., 2021). They are morphologically quite different from those of the past. Also, fossil records and proxy studies, which provide key information on morphological properties of early lycopsids and their approximate geographical distribution, are usually not sufficient to derive the extent of biogeochemical impacts of early vegetation in a quantitative way. It is thus necessary to develop an approach that can estimate vegetation properties for the late Silurian–Lower Devonian, based on available information, such as climate reconstructions and known physiological mechanisms.

Here, we introduce a process-based model, called LYCOm, which simulates lycopsid properties and is able to simulate diverse extinct plant communities based on given climatic conditions. This is complemented by a simple limit-based weathering model (Arens and Kleidon, 2011), which is utilized to determine the biotic enhancement of weathering by these plants. We test if the model can predict realistic lycopsid properties for a range of prescribed current climatic conditions, and assess whether the organisms have a significant impact on weathering rates. The current study focuses on the impacts of roots and stomata of Lycopsida on hydrological processes and the weathering of rocks of today, as a first step towards a more quantitative understanding of the impacts of early vascular vegetation on global biogeochemical cycles and climate of the past.

2 Methods

2.1 Lycopsid model (LYCOm)

The process-based vegetation model used in the current study aims at a general representation of lycopsids in order to estimate the productivity and physiological properties of these organisms under a broad climatic range. LYCOm follows other dynamic global vegetation models (DGVMs) closely regarding the representation of plant properties and functioning (Randerson et al., 2009). Vegetation is represented as a set of carbon pools, namely below and above ground, corresponding to roots and shoots with leaves. The balance of these pools depends on NPP, which is simulated as a function of climate forcing.

The current version of LYCOm differs from other vegetation models in some aspects. Stomatal conductance, for instance, is described as a function of soil water content and potential evaporation, instead of using the vapor pressure deficit (VPD) as the controlling factor for the calculation of photosynthesis (Sellers et al., 1996; Lawson et al., 2010). The use of VPD fails to accurately capture the fact that a plant’s stomata may also react to water loss due to transpiration into saturated air, which is warmed up by radiation and thus may take up water while staying close to saturation during this process. Hence, VPD does not capture the possibility of plants' water stress due to transpiration of water into saturated air (at low VPD), which is driven by the energy balance, or water stress due to low soil moisture. Hence, efficient handling of a broader range of environmental conditions is possible. Since we ultimately want to apply LYCOm to the Devonian Period, we prefer a general approach that considers all relevant factors affecting stomatal conductance over a specific parameterization, which may be more accurate for the present day, but might not hold for the geological past.

Another difference between LYCOm and other DGVMs is the representation of leaves and stems. To account for the accurate morphology of lycopsids, the aboveground biomass is distributed in the form of inclined cylinders, representing the stems, with small-sized leaves attached to the sides (see Fig. 2) instead of the detailed and complicated branching structure occurring in most plants. The vertical distribution of roots simulated in LYCOm is similar to other DGVMs with an exponential decrease in biomass distribution with depth.

Figure 2 The root profile and general soil scheme comprising the five layers of top soil up to 15 cm, followed by a bucket extending up to 0.65 m, are represented along with the hydrological profile and carbon pathways.

In contrast to most global vegetation models, which use plant functional types, LYCOm explicitly represents physiological variation concerning several key characteristics of plants. This means that the ranges of observed physiological parameters derived from the literature are sampled by a Monte Carlo algorithm to generate hypothetical “species”, which will be referred to as strategies in the article. These do not correspond directly to any particular lycopsids in the real world but instead represent the diversity of physiological strategies for lycopsids.

In this way, the response of a plant community of unknown composition with regard to a large range of environmental conditions can be simulated, and the effects of the adapted community on biogeochemical processes can be estimated. This modeling approach is more flexible compared to simulating a low number of “average” functional types with fixed properties, and it has already been successfully implemented in JeDi-DGVM (Pavlick et al., 2013), and was further developed in LiBry-DGVM for lichens and bryophytes (Porada et al., 2013). LYCOm is based in parts on the LiBry model. Another advantage of explicitly simulating diversity is the possibility of representing extinct species in the model, for which direct current analogs are lacking. The LiBry model, for instance, has already been applied to simulate the biogeochemical effects of early non-vascular vegetation in the Ordovician Period (Porada et al., 2016).

Variation in the physiological features of different lycopsid strategies in LYCOm results in differences regarding their carbon balance. To acquire new carbon, lycopsids have to invest their carbon resources into structures for carbon uptake. Since carbon can be allocated to different parts of the plant (leaves, roots), this represents a trade-off. The various strategies simulated by the model differ in their carbon allocation and are thus adapted to different climatic conditions. Hence, through analyzing the relative performance of each strategy, we can identify the key ecophysiological properties which determine the distribution of lycopsids under any given recent or past climate.

2.2 Model description

To determine both the spatial distribution of lycopsids and their impacts on weathering rates, it is necessary to quantify the carbon balance of the organisms and their dependence on climatic conditions. LYCOm includes three carbon pools for lycopsids: belowground (root) biomass and aboveground biomass, which is subdivided into a pool for leaves and another for stems. Gains in biomass result from NPP, while losses are due to mortality, which includes both aboveground litterfall and root turnover. NPP is calculated as the difference between photosynthesis and respiration (partitioned into leaf and root respiration), which are computed in LYCOm as functions of environmental climatic conditions, such as available water, light, and temperature. Thereby, root water uptake and subsequent transpiration through leaves connect the carbon to the water balance in the model. In the following sections, we describe the various equations aimed to describe biochemical and physical processes in LYCOm.

2.2.1 Soil scheme and hydrology

To simulate the water balance at the land surface and its connection to vegetation processes, we developed a hybrid soil scheme, consisting of a root zone, which is subdivided into five layers, and a “bucket” below the root zone. The layers have a thickness of 3 cm each and the bucket has a depth of 65 cm, resulting in a total depth of 0.8 m. This scheme aims to take into account the effects of the vertical distribution of soil water on the success of simulated strategies of class Lycopsida that differ in their root profiles. Since extant lycopsids, however, can only reach a maximum rooting depth of approximately 15 cm (Matsunaga and Tomescu, 2016), we resolve the soil water profile and its effects on the root profile up to this depth (Fig. 2). The inclusion of a bucket as water storage below the layers is necessary to simulate a realistic overall partitioning of rainfall at the land surface into evapotranspiration and runoff, i.e., avoid overestimation of runoff due to frequent saturation of the soil.

Rainfall serves as the primary water source for the soil and enters the first layer of the root zone. From there, water percolates down into the next layer, as a function of the relative water content of the respective layer, limited to the actual amount of water contained in the layer (Eqs. 1–3). This sequence is repeated until the bucket is reached. In each layer, surplus water, which results from the exceeded water storage capacity of a layer, directly enters the bucket. This is justified by the spatial resolution of the model, which is designed to ultimately run at the regional or global scale. If the surface soil is saturated locally, overflowing water usually infiltrates the surrounding soil, and does not immediately contribute to surface runoff. Only in the case the soil is saturated at a large scale, which corresponds to the bucket being near saturation, is surface runoff generated in the model, which then contributes to total runoff.

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathrm{Qin}=\mathrm{Rainfall}+\text{snow melt}(..\text{layer 1})\text{(2)}& {\displaystyle }\mathrm{Wu}=\mathrm{WO}+(\mathrm{Qin}-\mathrm{Qp})\cdot p_{\mathrm{dt}}\text{(3)}& {\displaystyle }\mathrm{Qp}=min(\mathrm{Qp}\mathrm{0}\cdot S,\mathrm{Wu}/p_{\mathrm{dt}})\end{array}$$

Here, Qin = incoming water into layers (rain + snow melt − interception by surrounding vegetation (for first layer)), p_dt = time step, Qp = percolation (which becomes the Qin for the following layer), QP0 = ($\mathrm{0.5}\times {\mathrm{10}}^{-\mathrm{7}}$ [m s⁻¹]) Percolation under water saturated condition, S = relative water saturation level of the layer, Wu = updated layer water content, and WO = initial layer water content.

In addition to infiltration and percolation, root water uptake changes the water content of the layered soil in the model. Root water uptake is calculated for each layer as a function of potential evapotranspiration (Monteith, 1981), which is then modified depending on the share of the root biomass contained in the respective layer on the total root biomass (Eq. 4). Additionally, root uptake is limited to the available water in each layer.

$$\begin{array}{}\text{(4)}& \mathrm{QR}=min(\mathrm{QRpot}\cdot \mathrm{Br}/\sum (\mathrm{Br}),\mathrm{Wu}/p_{\mathrm{dt}})\end{array}$$

Here, QR = root uptake from layer, QRpot = potential evapotranspiration [m s⁻¹], Br = layer biomass content, $\sum \left(Br\right)=$ total root biomass, and Wu = water content of layer.

$$\begin{array}{}\text{(5)}& \mathrm{Qb}=min(\mathrm{Qb}\mathrm{0}\cdot S\mathrm{2},\mathrm{Wx}/p_{\mathrm{dt}})\end{array}$$

Here, Qb = base flow, Qb0 = base runoff rate [m s⁻¹], S2 = average saturation of the bucket soil part, and Wx = water content of bucket.

The water balance of the bucket is then calculated according to the following equation:

$$\begin{array}{}\text{(6)}& \mathrm{Wx}=\mathrm{WIx}+(\mathrm{Qp}+\mathrm{Qs}-\mathrm{Qb}-\mathrm{Qo})\cdot p_{\mathrm{dt}},\end{array}$$

where Wx = updated bucket water content, WIx = initial bucket water content, Qp = percolation from upper layers, Qs = overflow from top layer, Qb = actual base flow, and Qo = large-scale runoff comprising the excess water from the saturated bucket.

2.2.2 Stomatal conductance

In LYCOm, stomatal conductance depends on two main factors: firstly, it is controlled by the evapotranspiration ratio (ETratio), which means the ratio of the potential evapotranspiration at a given time of the day and the average potential evapotranspiration for the overall timeframe (Eq. 8). Through the ETratio (Eq. 7), a normalized estimate of the atmospheric water demand (and, thus, potential water stress) is defined, which is independent of empirical parameters and considers both the VPD and net radiation as drivers of evaporation. The motivation for this approach is the potential dependence of empirical parameterizations on recent climatic conditions and also on physiological properties of the current vegetation community, which is likely adapted to the climate. Consequently, when applied to the geological past, these parameterizations may no longer be valid. Our approach is flexible enough to be applied to paleo-climatic scenarios. This thereby depends on the assumption that the plant community has adapted to the prevailing average conditions, and thus using our normalized measure for atmospheric water demand will lead to a realistic response of stomatal conductance to short-term variations in the climate forcing.

The average atmospheric evaporation demand for the full time span of the run is determined in the pre-processing. The run-time atmospheric demand is then evaluated using the following (Eq. 7):

$$\begin{array}{}\text{(7)}& \mathrm{ETratio}=\mathrm{ETdata}/\mathrm{ETdavg},\end{array}$$

where ETratio = relative atmospheric demand for the particular run-time, ETdata = run-time atmospheric water demand, and ETdavg = average atmospheric water demand for the whole runtime.

The stomatal conductance is then computed as follows:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}\mathrm{gsleaf}& =\mathrm{gs}\mathrm{0},\text{when},\mathrm{ETratio}\le \mathrm{1},\\ & \text{else},\mathrm{gsleaf}=\left(\mathrm{gs}\mathrm{0}\right)/(\mathrm{ETratio}{)}^{\log p_{\mathrm{gs}\mathrm{1}}/\log \mathrm{ETrmax}},\end{array}\end{array}$$

where gs0 = maximum stomatal conductance (varies between maximum derived from Soni et al. (2012) and minimum derived from (Tosens et al., 2016); 20 % more than the actual value stipulated for Pteridium aquilinum), ETratio = potential evapotranspiration at the given time step / average evapotranspiration, p_gs1 = shape parameter (1000 used here), ETrmax = maximum ETratio for the entire simulation, and gsleaf = stomatal conductance due to atmospheric demand.

Secondly, stomatal conductance is limited by maximum conductance (dependent on species) and soil water availability in our model. To this end, the conductivity for water flow at the soil root interface is computed, which is proportional to the soil water content (Eq. 10). The realized stomatal conductance is then limited to root conductance. In this way, also the impacts of soil water stress on stomatal conductance can be taken into account, in addition to the atmospheric drivers.

$$\begin{array}{}\text{(9)}& {\displaystyle }\mathrm{gsroot}=\mathrm{gs}\mathrm{0}\cdot S\mathrm{1}\text{(10)}& {\displaystyle }\mathrm{gsfinal}=min(\mathrm{gsleaf},\mathrm{gsroot})\end{array}$$

Here, gsroot = stomatal conductance determined by average layer water content, gsfinal = stomatal conductance limited both by soil water availability and atmospheric demand, and S1 = average saturation of the layered soil part.

2.2.3 Above- and belowground biomass

In LYCOm, plant growth corresponds to NPP, which is based on the simulated net photosynthesis averaged over 1 month. Thereby, net photosynthesis is computed as gross photosynthesis minus respiration according to Farquhar and Von Caemmerer (1982).

$$\begin{array}{}\text{(11)}& {\displaystyle }{f}_{\mathrm{GPP},\mathrm{L}}={\displaystyle \frac{J\cdot ({\mathrm{CO}}_{\mathrm{2}}-P)}{\mathrm{4.5}\cdot {\mathrm{CO}}_{\mathrm{2}}+\mathrm{10.5}\cdot P}}\text{(12)}& {\displaystyle }{f}_{\mathrm{GPP},{\mathrm{CO}}_{\mathrm{2}}}=\mathrm{vcm}\cdot {\displaystyle \frac{({\mathrm{CO}}_{\mathrm{2}}-P)}{{\mathrm{CO}}_{\mathrm{2}}+\mathrm{Kc}\cdot \frac{\mathrm{1.0}+{\mathrm{O}}_{\mathrm{2}}}{\mathrm{Ko}}}}\end{array}$$

Here, f_GPP,L = light-limited gross primary production, $f_{\mathrm{GPP},{\mathrm{CO}}_{\mathrm{2}}}$ = carbon dioxide-limited gross primary production, vcm = maximum specific carboxylation rate, J = actual rate of electron transport, CO₂ = CO₂ concentration in the chloroplast, O₂ = oxygen concentration in the chloroplast, P = CO₂ compensation point, and Kc, Ko = Michaelis–Menten constants of carboxylation and oxygenation reactions.

Gross photosynthesis is simulated as the minimum of a light-limited and a CO₂-limited rate. The light-limited rate is an increasing function of the absorption of light and is constrained by the maximum rate of the electron flux (J_max). The CO₂-limited rate is an increasing function of the CO₂ concentration in the chloroplasts and saturates at high values of CO₂ concentration, the maximum rate being VC_max. The detailed formulation can be found in Porada et al. (2013).

$$\begin{array}{}\text{(13)}& \mathrm{GPP}=min(f_{\mathrm{GPP},\mathrm{L}},f_{\mathrm{GPP},{\mathrm{CO}}_{\mathrm{2}}})\end{array}$$

Here, GPP = gross primary production.

Photosynthesis peaks around an optimum surface temperature estimated for the individual strategies (June et al., 2004), and the estimation of respiration is done using a common empirical temperature response framework with respect to Q₁₀ (Eq. 15) (Vanderwel et al., 2015). GPP and NPP are calculated as follows:

$$\begin{array}{}\text{(14)}& {\displaystyle }\mathrm{NPP}=\mathrm{GPP}-R,\text{(15)}& {\displaystyle }R=R\mathrm{0}\cdot Q_{\mathrm{10}}^{\left(\right(T-\mathrm{Topt})/\mathrm{10.0}))},\end{array}$$

where NPP = net primary production, R = Respiration, R0 = Reference respiration rate at 10 ^∘C (Tosens et al., 2016; Soni et al., 2012; Carriquí et al., 2019), Q₁₀ = Q₁₀ value of respiration (Atkin et al., 2005; Proctor, 2014), T = temperature on the surface of the plant, Topt = optimum temperature.

The pore space CO₂ concentration inside the leaves is a function of the stomatal conductance which is in turn controlled by the water availability in the soil (Eq. 10). Thereafter, we considered a steady-state assumption of pore space CO₂ between diffusion and net photosynthesis which is in turn dependent on respiration. This drives the light- and CO₂-limited Farquhar and Von Caemmerer (1982) assimilation model of photosynthesis. The rest of the scheme is analogous to the LiBry model (Porada et al., 2013).

Plant growth is often limited by water availability, which causes the organisms to invest carbon from NPP into root structures to increase access to soil water. Since plant growth is also limited by light, carbon needs to be invested into stems and leaves, too. This leads to a trade-off regarding the allocation of the assimilated carbon to either above- or belowground biomass. To address this key physiological trade-off, LYCOm includes a flexible allocation scheme that is adaptive to the prevailing atmospheric conditions, i.e., whether the plant photosynthesizes under light-limited or CO₂-limited (i.e., water-limited) conditions. The light-limited photosynthesis enhances the allocation of NPP into shoots and leaves while the CO₂ limitation leads to the accumulation of root biomass (Eqs. 16, 17 and 18). The distribution and amount of root biomass and aboveground biomass are key to the quantitative estimation of impacts of vegetation on multiple biogeochemical processes, such as microbial respiration of soil organic matter being closely linked to chemical weathering.

Equations (16) and (17) determine the amount of assimilated carbon (NPP) assigned to root biomass with depth depending on the relative amount of limiting factors for photosynthesis. The available root biomass B_root is then distributed layer-wise ($B_{{\mathrm{root}}_i}$​​​​​​​ in Eq. 18), with decreasing fractions at increasing depth, which is determined by the fracR₀ factor (Fan et al., 2016). This ensures a decreasing cumulative root distribution profile with depth.

$$\begin{array}{}\text{(16)}& {\displaystyle }{B}_{{\mathrm{root}}_{\mathrm{in}}}=\mathrm{NPP}\cdot (\mathrm{1.0}-\mathrm{fracL})..\left(\text{layer 1}\right)\text{(17)}& {\displaystyle }{B}_{{\mathrm{root}}_{\mathrm{in}}}=(\mathrm{1.0}-{\mathrm{fracR}}_{\mathrm{0}})\cdot B_{\mathrm{root}}..\left(\text{layers below}\right)\end{array}$$

Here, NPP = net primary productivity and fracL = fraction of time where light limitation dominates during photosynthesis. In this case, we allocate more biomass into leaves in order to capture more light for photosynthesis. fracR₀ = shape constant (0.25) for an exponential decrease of root biomass allocation with depth and B_root = available root biomass that is still available for allocation in this layer.

$$\begin{array}{}\text{(18)}& B_{{\mathrm{root}}_i}={\mathrm{IB}}_{{\mathrm{root}}_i}+{\mathrm{fracR}}_{\mathrm{0}}\cdot B_{{\mathrm{root}}_{\mathrm{in}}}-{\mathrm{Mortality}}_i\end{array}$$

Here, $B_{{\mathrm{root}}_i}$ = root biomass of the layer i, IB${}_{{\mathrm{root}}_i}$ = Root biomass contained in the ith layer, $B_{{\mathrm{root}}_{\mathrm{in}}}$ = root biomass available, and Mortality_i = a term based on baseline mortality and the biomass allocation in soil layers.

The assigning of root biomass to litterfall (Eq. 19) is an important aspect to focus on as it determines the silicate weathering in the soil. When the soil layer reaches its maximum root biomass holding capacity, the excess is directly assigned to litterfall (Eq. 19). This adds to the mortality calculated in Eq. (18), which is in turn utilized to determine the soil carbon pool and later the weathering regime.

LYCOm considers the root respiration as well, which potentially could further affect the soil CO₂ content (Eq. 20). Root respiration (Eq. 20) is basically comprised of root maintenance respiration, which is used for keeping roots alive, root growth respiration for growing new roots as well as new root tissues, along with microbial respiration.

$$\begin{array}{}\text{(19)}& \mathrm{RtoM}=max(\mathrm{0.0},B_{{\mathrm{root}}_i}-{\mathrm{Brmax}}_i)\end{array}$$

RtoM = root to direct litterfall, $B_{{\mathrm{root}}_i}$ = root biomass of layer, and Brmax_i = maximum layer root biomass.

$$\begin{array}{}\text{(20)}& R_{\mathrm{resp}}=R\mathrm{0}\cdot Q_{\mathrm{10}}^{\left(\right(T-\mathrm{Topt})/\mathrm{10.0}))}\cdot (B_{{\mathrm{root}}_i}/\sum (B_{\mathrm{root}}\left)\right)\end{array}$$

Here, R_resp = root respiration, R0 = respiration determined by strategy-specific parameters, Q₁₀ = Q10 value, T = soil temperature, Topt = optimum temperature, $B_{{\mathrm{root}}_i}$ = root biomass of layer, and $\sum \left(B_{\mathrm{root}}\right)=$ total root biomass.

Besides roots, the hypothetical plants in the model also consist of shoots and leaves partitioned at a 3:7 ratio, which constitute the total available aboveground biomass, represented as a cylinder to emulate the vertical structure of the given plant strategies, as well as the leaf area index (LAI) calculated in Eq. (22). The leaf-to-shoot ratio is analogous to the horizontal to vertical biomass ratio value derived from Zier et al. (2015). Lycopsids exhibit root-bearing axes and true roots, leafy shoot axes, along with scale-like leaves called microphylls covering the aboveground part of the plant (Hetherington et al., 2021). In contrast to trees or shrubs, the scale-like nature of leaves is hard to distinguish from the stem in the model and hence we assume that photosynthesizing parts of the plant in the model comprise approximately 70 % of the aboveground biomass.

This is a distinct feature of the model which departs from the traditional simple LAI approach while representing the aboveground share of biomass.

The leaf biomass (Eq. 21) is used thereafter to estimate the LAI which is computed using Eq. (22).

$$\begin{array}{}\text{(21)}& \mathrm{Bl}=\mathrm{Bli}+\mathrm{NPP}\cdot \mathrm{fracL}-{\mathrm{Mortality}}_{\mathrm{shoot}}\end{array}$$

Here, Bl = leaf biomass, Bli = initial leaf biomass, and Mortality_shoot = a term based on baseline mortality and the biomass production and allocation into aboveground biomass.

$$\begin{array}{}\text{(22)}& \mathrm{LAI}=(\text{Total Leaf Area}/\mathrm{GrA})\end{array}$$

Total leaf area = Aleaf ⋅ (Bl$/$Drw), Aleaf = leaf area using measurements from Valdespino (2015a) and Drw = dry leaf weight as measured in Leopold et al. (1981), GrA = ground area.

The projection of the ground area covered by plants used for the calculation of the LAI in the model assumes a random value ranging from 30 to 50 cm². The range is close to the reconstruction of Zosterophyllum shengfengense (Hao et al., 2010) but we assume a maximum coverage of a lateral area equivalent of a square of 70 mm edge.

2.2.4 Weathering and soil carbon dioxide

To estimate chemical weathering rates and their biotic enhancement, we apply a simple limit-based approach (Arens and Kleidon, 2011), which has been used in other studies to quantify large-scale weathering, both for the current climate and also for the geological past (Porada et al., 2016). This approach assumes that chemical weathering can be derived from the minimum of a supply-limited rate, which corresponds to the provision of primary minerals into the dynamic soil profile via uplift, and a transport-limited rate, which describes the removal of dissolved weathering products from the soil via runoff, named the eco-hydrological limit. Without runoff, the soil solution would remain in a state of chemical saturation, and the weathering reactions would stop. Consequently, the chemical weathering rate cannot exceed the rate of weathering products exported via runoff. Furthermore, the rate of chemical weathering cannot exceed the rate of exposure of unweathered rock material at the surface, which is equal in the steady state to the erosion rate (see Table 3).

While the supply limit is not affected by biotic processes, unless the impact of vegetation on erosion and uplift are considered, the transport limit may be altered through the effect of CO₂ in the soil solution on the equilibrium concentration of weathering products, and also through the effects of vegetation on runoff. The erosion is simply a function of the elevation of the sites (see Table 2), which is derived from the database of the National Center for Atmospheric Research called “TerrainBase” (TBASE).

Plants affect soil CO₂ mainly via two processes: firstly, they directly release CO₂ into the soil as a result of root respiration, and, secondly, they enhance microbial respiration due to the decomposition of plant biomass originating from root turnover and litter input into the soil. In LYCOm, we assume a steady state between the incoming soil CO₂ via mortality (Eqs. 16–18) and root respiration (Eq. 19), and the outgassing of the CO₂ from the soil. The overall carbon going into the soil as litter is subject to diffusion, and this is the primary carbon pool responsible for facilitating weathering. Subsequently, the concentration of CO₂ in the soil can be determined from the known atmospheric CO₂ concentration and a soil diffusivity equation (Jabro et al., 2012) as shown in Eq. (23).

$$\begin{array}{}\text{(23)}& {\mathrm{Soil}}_{\mathrm{C}}=\left(B_{\mathrm{r}}\right)\cdot S_{\mathrm{d}}/\left(D_{\mathrm{c}}\right)\cdot p_{\mathrm{a}}+\left({\mathrm{CO}}_{{\mathrm{2}}_{\mathrm{p}}}\right)\cdot S_{\mathrm{d}}/\left(D_{\mathrm{c}}\right)\cdot p_{\mathrm{a}}+{\mathrm{CO}}_{\mathrm{2}}\end{array}$$

Here, B_r = root turnover and litter input to the soil, D_c = $T\cdot \mathit{ϵ}\cdot D\mathrm{0}\cdot \mathrm{hmon}$, T=0.66 (tortuosity constant), D0=0.05, ϵ = air-filled soil porosity, hmon = hours in month, CO${}_{{\mathrm{2}}_{\mathrm{p}}}$ = soil carbon dioxide concentration from the previous time step, p_a = 0.51 (L g⁻¹), CO₂ = atmospheric CO₂ content, and S_d = maximum root depth.

2.3 Physiological strategies

The variation of physiological properties in the model is summarized in Table 1. Some photosynthetic parameters such as the stomatal conductance responsible for the diffusion of atmospheric CO₂ into the leaves is limited as shown in the aforementioned table. The upper limit of the optimal temperature is limited to between 0 and 50 ^∘C along with a restricted Q₁₀ value between 1.5 and 2.3. The RuBisCO kinetic limits are adjusted from Galmes et al. (2014). This is done keeping in mind that the current model serves as a basis for paleoclimate (Devonian) simulations when the lycopsids were abundant and the mean surface temperature was comparatively higher than today (Song et al., 2019). Furthermore, the physical realism of the lycopsid is asserted by constraining certain parameter ranges such as the leaf area (Valdespino, 2015b) or the dry weight of leaves (Leopold et al., 1981) based on literature that includes various taxonomic details of lycopsids. The rest of the scheme is analogous to the model described in Porada et al. (2013).

Zhang et al. (2014)Hikosaka and Hirose (1997)Galmes et al. (2014)Savir et al. (2010)Galmes et al. (2014)Savir et al. (2010)Medlyn et al. (2002)Medlyn et al. (2002)Kattge and Knorr (2007)Kattge and Knorr (2007)

Table 1 Parameters used to define the lycopsid-specific model.

2.4 Weighting scheme

The weighting scheme used to calculate the resulting productivity depends on the allocation of root and aboveground biomass (Eqs. 24, 25, 26). To obtain a representative value, we scale the productivity of all surviving strategies by their weights. This yields an averaged representative NPP (Table 4) and soil CO₂, which would be difficult once the dominant species is solely chosen for the purpose. To eliminate overestimation, we consider the root as well as leaf biomass for the weighing scheme. In Table A1 a comparative chart of the productivity by weighted mean, mean, and that of the dominant species is summarized.

$$\begin{array}{}\text{(24)}& {\displaystyle }\text{weighted NPP}=\sum ({\mathrm{NPP}}_{ji}\cdot W_i)\text{(25)}& {\displaystyle }{\mathrm{NPP}}_j={\displaystyle \frac{\mathrm{1}}{\mathrm{30}}}\sum _{n=\mathrm{1}}^{\mathrm{30}}\left({\mathrm{NPP}}_{nj}\right),j=\mathrm{1},\mathrm{2}\mathrm{\dots }\mathrm{12}\text{(26)}& {\displaystyle }{W}_i={\displaystyle \frac{(M_i/\stackrel{\mathrm{‾}}{M})+({\mathrm{LAI}}_i/\stackrel{\mathrm{‾}}{\mathrm{LAI}})}{{\sum }_i\left(\right(M_i/\stackrel{\mathrm{‾}}{M})+({\mathrm{LAI}}_i/\stackrel{\mathrm{‾}}{\mathrm{LAI}}\left)\right)}}\end{array}$$

Here, NPP_ji = monthly mean NPP of species i, W_i = weight calculated for species i, NPP_j = monthly NPP averaged over 30 years, M_i = root biomass accumulated over 30 simulation years for species i, $\stackrel{\mathrm{‾}}{M}$ = mean root biomass for all surviving species, LAI_i = leaf area index at the last time step for species i, and $\stackrel{\mathrm{‾}}{\mathrm{LAI}}$ = mean leaf area index of the surviving species at the final time step.

3 Model setup and validation

The performance of LYCOm is assessed using gas exchange measurements of photosynthesis by lycopsids under varying ambient conditions, such as temperature, carbon dioxide concentration, and photosynthetic photon flux density (PPFD) (Soni et al., 2012). The data used for the purpose is obtained for Selaginella bryopteris (see Fig. 5 dashed lines). This aims to confirm the ability of LYCOm to replicate real species and affirms the realism of our approach. Moreover, certain parameter ranges, for example, the habitable temperature range (Topt in Table 1) of lycopsids, are determined using the data derived from Soni et al. (2012).

3.1 Model setup

We run LYCOm locally for seven locations (listed in Table 2) with 2.8125^∘ spatial resolution data for 30 years between 1958 and 1989 with an hourly time step, simulating 100 different physiological strategies one at a time, and evaluate their survival rates. The meteorological forcing data set used for the purpose is derived from WATCH Forcing Data (WFD) by making use of the ERA-Interim reanalysis data (Weedon et al., 2011). In addition, the model setup is initialized with a close-to-present atmospheric composition of CO₂ concentration of 400 ppm and oxygen content of 210 000 ppm. The soil porosity is constrained to a relative value of 0.45 with a percolation rate of 0.5 × 10⁻⁷ m s⁻¹ between the layers and an average baseflow of 1.0 × 10⁻⁸ m s⁻¹ from the bucket.

Table 2 Various characteristics of the sites where the model is run, including the local relief, and a comparison of the runoffs from the model and historical runoff data between 1958 and 1988.

3.2 Model validation

The chosen study sites encompass various climate zones to assess whether the model can capture the productivity and viability of the local lycopsid communities. The chosen sites are characterized by a high occurrence of lycopsids, as determined based on the Global Biodiversity Information Facility (GBIF, https://www.gbif.org, last access: 10 June 2021).

The study of Ghiggi et al. (2019) provides data that is suitable to evaluate the hydrological cycle achieved via the hybrid soil hydrological scheme in LYCOm.

Although the physiological evolution of lycopsids is well documented, there is a significant lack of literature concerning biochemical properties and the characteristic productivity range of lycopsids under various climate regimes. Hence, our validation of the model is limited. However, the on-site measurements of productivities of Lycopodium annotinum and Lycopodium clavatum in Estonia by Tosens et al. (2016) in May and early June provides an order-of-magnitude estimate in this regard, which improves our ability to evaluate LYCOm. Both Lycopodium annotinum and Lycopodium clavatum have a widespread distribution across several continents today and are the most common species in Estonia (https://www.gbif.org, last access: 10 June 2021). A similar study, by Campany et al. (2019) of Costa Rica, affirms the robust model performance.

3.3 Parameters of LYCOm

As described in Sect. 2, LYCOm represents the physiological diversity of lycopsids via generating multiple strategies, which are closely analogous to thriving lycopsid species, and the intra-specific diversity of lycopsids is accounted for in the physiological processes implemented in the LYCOm model. The distinguishing features of lycopsid strategies are represented by 14 parameters and their corresponding ranges. To generate the lycopsid strategies, these 14 characteristic parameters are assigned through randomly sampling ranges of possible values (Table 1). Assignment of parameter values is performed in two steps: (a) for each strategy, a set of 14 random numbers uniformly distributed between 0 and 1 is sampled. The random numbers are generated by a Latin hypercube algorithm (McKay et al., 2000). (b) These random numbers are used to map values of the parameters from permissible ranges particular to lycopsids derived from the literature. The particulars of the equations utilized in the model are briefly described in Table 1.

4 Results

The modeled lycopsids cope well at all the sites except in India, where the plants die due to water stress. The aforementioned run is undertaken to incorporate the microclimate of Selaginella bryopteris from where the sample was collected by Soni et al. (2012). The study forms the basis of the light, CO₂, and temperature dependence of productivity, for model calibration (Fig. 4). Out of 100 strategies, a maximum of 46 species survive in New Zealand followed by the USA which features 41 surviving strategies. The climate of Sweden suits 24 varieties of the generated species, and in Costa Rica the number is 30. The prevailing weather conditions in Japan and Peru support 32 and 24 strategies, respectively. In India, none of the strategies survive 30 simulation years as shown in Table 2, and Estonia features the least number of surviving strategies at 18.

The measured productivity of Lycopodium clavatum and Lycopodium annotinum in Estonia by Tosens et al. (2016) range from 1.4 to 1.7 µmol CO₂ m⁻² s⁻¹ for the respective species. The productivity range in LYCOm model varies from 0.5 to 1.23 µmol CO₂ m⁻² s⁻¹ in May and early June for five dominant species possessing the maximum weights (Eq. 26). The ecophysiological study of Selaginella, an early lineage vascular plant group from a tropical forest understory in Costa Rica (Campany et al., 2019) in June exhibits a range of net photosynthesis from 1.8 to 5.8 µmol CO₂ m⁻² s⁻¹. The mean productivity of the top 10 strategies (by weights) over the 30-year simulation period for June varies between 1.07 and 2.42 µmol CO₂ m⁻² s⁻¹ in the model. The higher productivity range in the study results from the prevailing atmospheric temperature, which ranges from 24 to 29 ^∘C during the on-site measurement, while in the model the temperature range in June varied from 13 to 25 ^∘C between 1958 and 1988 with a mean temperature of 21 ^∘C which is significantly lower than the aforementioned study.

In general, the LYCOm simulates high weighted yearly productivity (see Table 4) at the tropical site in Costa Rica of 20.45 mol CO₂ m⁻² (245.37 g C m⁻²), followed closely by the area in New Zealand estimating around 14.86 mol CO₂ m⁻² (178.35 g C m⁻²), which shows a reverse seasonality compared to other locations (Fig. 4) situated in the Northern Hemisphere. The mainland around the Great Lakes in the USA shows annual productivity of 10.68 mol CO₂ m⁻² (128 g C m⁻²) and values in Sweden and Estonia of around 12.22 mol CO₂ m⁻² (146.7 g C m⁻²) and 10.5 mol C m⁻² (126.03 g C m⁻²), respectively, which is comparatively low as a result of the local climate that restricts high productivity to the summer and renders poor productivity in other times of the year. Peru, with a higher average year-round temperature, provides a reasonably better climate for lycopsids, resulting in productivity of 14.17 mol CO₂ m⁻² (170.09 g C m⁻²) per year.

We averaged the weighted monthly productivity over a 30-year seasonal cycle to obtain an average measure of annual plant growth (Table 4). We use a weighting scheme (see methods) to generate weighted NPP to obtain a mean seasonal cycle over 30 years (Fig. 4). LYCOm can capture the seasonality with the highest productivity over the summer months. The productivity is solely a function of the local climate and is clear from Fig. 4.

The assessment of the net assimilation rate is undertaken to evaluate the model performance, as mentioned in the section on model setup and validation (Sect. 3.2), to verify its capability to reproduce the lycopsid-specific lab measurements with our newly imposed parameter ranges (Table 2). The model performs reasonably, and increased photosynthetic productivity is observed with an increase in photosynthetic photon flux density (PPFD), reaching saturation at 600 µmol m⁻² s⁻¹ for all strategies with a steeper slope, as opposed to Soni et al. (2012), where they observed a saturation around 800 µmol m⁻² s⁻¹. The temperature response studies in Soni et al. (2012) showed a gradual increase until 40 ^∘C and it dropped marginally beyond that, a behavior well captured in our calibration (see Fig. 3) with a slightly lower maximum (around 7 µmol m⁻² s⁻¹) in all strategies. Photosynthetic response to carbon dioxide levels in close consistency with that of Selaginella bryopteris for strategy 53, while the rest of the strategies exhibit saturation at higher maxima. We use four strategies to calibrate since the generated strategies are not a direct representation of a single species of lycopsid. This gives a relatively good idea of the model behavior under various conditions for unique strategies.

Figure 3 The locations used for the simulations are found in the following regions: United States of America (USA) [47^∘ N 78^∘ W], Costa Rica (CR) [10^∘ N, 85^∘ W], Sweden (Swe) [58^∘ N, 13^∘ E], Estonia (Est) [58.75^∘ N, 24^∘ E], Japan (Jap) [33^∘ N, 133^∘ E], New Zealand (NZ) [40^∘ S, 176^∘ E], Peru [2^∘ S, 74^∘ W].

Runoff is generated in LYCOm when water input by rainfall or snowmelt exceeds the water-storage capacity of the soil, with consideration of water loss via plants according to the climatic demand and stomatal regulation. Hence climate and vegetation are pivotal to the calculation of the runoff of the model. Table 2 confirms the realistic nature of our hydrological approximations and we obtain values of runoff (mean yearly runoff) from LYCOm that are comparable to historical runoff data between 1958 and 1988 (Ghiggi et al., 2019)​​​​​​​. In general, the model runoff is overestimated when compared with the literature (Ghiggi et al., 2019) except for the location in New Zealand where LYCOm slightly underestimates the annual runoff per unit area (m²). The model renders an average runoff of 2605 mm yr⁻¹ for Peru and 2134 mm yr⁻¹ for Costa Rica, followed by 823 mm yr⁻¹ for Japan (see Table 2). This overestimation can be attributed to the fact that only the layered topsoil is subject to transpirational loss as opposed to the conditions in real forests or vegetated regions since the roots of modern plants penetrate deeper into the soil. The loss by transpiration is several times higher than we capture in the current model, which simulates only one kind of vegetation, i.e., lycopsids. In addition, the soil water saturation determines the amount of CO₂ diffusing out of the soil. The water in the soil acts as an inhibitor to the diffusion of CO₂ through pathways of available pore space. The more the soil saturates, the lower the CO₂ emissions. The concentration of soil CO₂ (see Fig. 4) is thus coherently linked with the rainy season and productivity (monthly productivity shown in Fig. 4). This effect is prominent in some locations, namely, in Costa Rica, when both the productive season and most rainfall coincide, leading to a high content of soil CO₂ in June. The effect of both productivity and rainfall on carbon content in soil can be seen in some cases as well, such as in Sweden and Estonia, where the most productive months for lycopsids are May and June but soil CO₂ peaks slightly in August in the rainy months.

The resulting soil CO₂ concentration after diffusion forms the basis of the weathering estimate in the post-processing of the model.

Figure 4 Annual cycle of weighted NPP (solid lines) and carbon dioxide concentration (dashed) in the soil.

Soil CO₂ (see Table 4) concentration, which is primarily a function of productivity and soil water saturation, is seen to reach the maximum level in Costa Rica on a yearly average scale of 77 330 ppm over the 30 simulation years, followed closely by the location in Peru, where the average carbon dioxide content in soil is averaged at 53 640 ppm on a yearly scale (Table 4). Even though the seasonality of the soil CO₂ shows a close relationship with productivity (Fig. 4) on a seasonal scale, we identify a broad range of soil CO₂ values over the study sites, with the overall minimum at the study site in the USA, where NPP is low as well. The values in Table 4 are the yearly aggregation of the weighted monthly soil carbon dioxide content over 30 years. Figure 4 shows the monthly fluctuating nature of the monthly soil CO₂ content in the soil averaged over the full period of simulation (in ppm).

The supply-limited calculation of weathering rates (Arens and Kleidon, 2011) is a unique approach, where chemical weathering of silicate minerals in the soil is limited either by erosion rate or by runoff, the latter being the “eco-hydrological limit”. Furthermore, LYCOm accounts for the effect of soil CO₂ on the maximum saturation state of the soil solution (see Sect. 2.2.4). At the respective sites, the chemical weathering is evaluated using the soil CO₂ and runoff from LYCOm along with erosion. The eco-hydrological limits play a crucial role in determining the chemical weathering rates at all sites except for Peru and Costa Rica, where the weathering is limited by erosion. Areas with shallow relief structures show consequently low erosion, such as lowland Peru, where weathering by lycopsids becomes limited by a layer of highly weathered material blocking access to the unweathered rock substrate. We represent this potential limitation of chemical weathering implicitly by calculating the erosion rate, which is essentially a function of elevation of the location. Maximum ecohydrologically limited weathering rates of 0.418 and 0.451 mm of rock per year for every square meter is found in Costa Rica and Peru, which is representative of the high levels of productivity. Weathering in these areas is strictly restricted by erosion (see Table 3). The area of study in Japan exhibits a weathering rate of 0.086 mm rock m⁻² per year, which is comparatively higher than the simulated areas in the USA, Sweden, and New Zealand, which showcase 0.045, 0.055, and 0.055 mm rock m⁻² weathering per year, respectively.

Table 3 Weathering regimes post processed using Arens and Kleidon (2011) approach from LYCOm products.

Table 4 LYCOm outputs: survival rates, weighted leaf area index (LAI), weighted NPP and weighted soil CO₂ content.

A regional influence of ecohydrologically limited weathering is conceived in the sensitivity runs undertaken to obtain insight into the influence of lycopsids on weathering. In addition to the runs with (a) vegetation (lycopsids), we run the model in (b) a vegetation-free scenario. The model run with no lycopsids is undertaken to determine the large-scale runoff, and the soil CO₂ concentration is assumed to be in equilibrium with that of the atmosphere for computation of the weathering when the location is devoid of vegetation (scenario (b) in Table 5). In general, an elevated level (see Table 5) of weathering is observed at all sites with vegetation as compared to the case without vegetation. The weathering was observed to be 6 times higher in Costa Rica and almost 5 times in Peru with vegetation. The location-specific enhancement due to vegetation is a result of different levels of productivity and varying hydrology. The least enhancement is recorded around the great lakes in the USA, followed by Sweden and Japan.

Table 5 Comparison of supply limited eco-hydrological-controlled weathering regimes with vegetation (a) and without vegetation (b). In addition, the percolation rate is doubled in the with vegetation (a) scenario, which incorporates the impact of roots on the soil (making it more porous).

5 Discussion

LYCOm represents a community of lycopsids at several local sites. The respective locations provide suitable climatic conditions for the growth of lycopsids (GBIF), and the model shows relatively high survival rates there, as pointed out in Table 2. The tropical climate of Costa Rica and Peru ensure continuous productivity throughout the year, while the temperate climate at the rest of the sites results in a strong seasonal pattern (see Fig. 4). Increased productivity in warmer months as compared to colder months varies over latitudes as well. Indirectly, differences in climatic conditions between the locations affect productivity levels of the plants also through soil water availability. In addition, a prominent impact of soil water can be noticed in the soil CO₂ concentration. The microclimate at the site chosen for emulating the lab sample in India is likely to be substantially different from the large-scale climate derived from the global dataset. The climate of the region is markedly drier compared to the other chosen locations in the study. Therefore, the rainfall is insufficient for the sustenance of the generated lycopsid strategies. As a result, no strategies survive through the whole simulation period at the location.

The runoff in the LYCOm model is comparable to Ghiggi et al. (2019) and the lycopsid productivity in the model is comparable to Tosens et al. (2016) and Campany et al. (2019), even though the generated strategies in the model are not analogous to the species utilized for observation. This establishes LYCOm as a robust process-based vegetation model with reasonable eco-hydrological representation (Table 2) given the complexity of the model. Weathering is classified in LYCOm according to the dominant limiting factor, which is either the supply of parent material or the eco-hydrological conditions. In the sensitivity study, it is apparent that the influence of lycopsids on chemical weathering rates is substantial at the local scale (Table 3). It is in fact, significantly higher than lichens and bryophytes, which colonized terrestrial ecosystems earlier than lycopsids (Wellman et al., 2003; Wellman and Gray, 2000; Chang and Graham, 2011; Yin-Long, 2008). Aghamiri and Schwartzman (2002) explore weathering by lichens and bryophytes on rock surfaces and estimate a range of 0.0004 to 0.01 mm rock yr⁻¹ weathering, comparatively lower than that of lycopsids, which exhibit a greater weathering potential of 0.026 to 0.418 mm rock yr⁻¹ (see Table 3). The change in weathering potential emphasizes the impacts of lycopsids and their corresponding effects on the biogeochemical cycle.

The approach of simulating strategies encompassing various physiological properties and trade-offs adds flexibility to the model, which is advantageous in representing already extinct or evolving species. In Fig. 5, four species generated in LYCOm are plotted under varying ambient conditions, against assimilation rates measured in the laboratory by Soni et al. (2012). This figure depicts the potential (test species in figure, shown in dashed lines) of the model to closely replicate real-world species, which can be achieved with an appropriate selection of parameters within the defined ranges of LYCOm.

Figure 5 CO₂, light, and temperature responses of four lycophyte strategies simulated (colored solid lines) in LYCOm against the four best test strategies (dashed lines) generated from LYCOm parameterization ranges in order to replicate the lab data (in black solid line; Soni et al. (2012)) as closely as possible.

LYCOm emulates well the behavior of the organisms under local climate at the local scale. The chosen locations, suggested to be suitable for such plant communities, do not , however, allow for direct extrapolation of the effects of lycopsids at the global scale. For further assessment of biogeochemical impacts of lycopsids, it is necessary to extend the model to a large scale.

Even though we cannot account for all biochemical processes existing in nature, the model accounts for the most relevant properties of lycopsids and represents the local interactions sufficiently to a certain degree. It is crucial to ensure that our modeling approach remains consistent when we upscale it in the future. LYCOm, paired with the weathering model, forms a crucial component towards understanding the long-term impacts of such plants on Earth's climate, as influenced by the CO₂ content of the atmosphere. Although our model can provide subsequent insights into the impacts of such plants, the study is restricted by a lack of data on the properties of lycopsids. Data concerning plant physiology and measurements involving productivity will contribute towards better optimization and validation of the model. The light, CO₂, and temperature curves of lycopsid species other than Selaginella bryopteris would not only strengthen the model performance further, but they would also help us constrain the parameter ranges with better precision.

One option for improvement would be the use of high-resolution meteorological data to reduce the uncertainty of the model estimates. This would not only provide an enhanced possibility to evaluate the model performance but also resolve the biogeochemical processes on a finer scale. Local enhancement of weathering rates could be captured better with such a setup since runoff and soil CO₂ concentration are dependent on productivity and hydrological balance at the respective sites. Table 6 summarizes the sensitivity of weathering with increases and decreases in these components. We do not observe any disproportional enhancement factor in the weathering and thus assume that our weathering estimates are robust. Although we consider the light and water interception by the surrounding, the influence of the large forests on the water uptake from soil and humidity of the surrounding needs to be accounted for explicitly when upscaling our estimates. This may be achieved with a higher spatial resolution of the model input data representing environmental conditions.

Table 6 Sensitivity of eco-hydrological-controlled weathering regime at respective locations with (+) 20 % enhancement and (−) 20 % decrease of runoff and soil CO₂ concentration relative to the original enhancement due to vegetation from Table 5.

6 Conclusions

The LYCOm model draws in information from fossil records as it helps in setting the physiological limitations of the extinct species. Modern lycophytes, such as terrestrial lineages of Lycopodiaceae, maintain some degrees of functional and morphological conservatism concerning Lower Devonian herbaceous lycopsids and their forerunners. On this basis, modern lycophytes can be used as a suitable analog system for widespread elements of Lower Devonian plant communities. Lycopsids survive under a variety of climatic conditions (Petersen and Burd, 2018) and have been represented well by the local version of the model. The productivity from the model shows strong agreement with the on-site measurements for the contemporary lycopsids. LYCOm, therefore, has the potential to be utilized beyond the range of current conditions. The study aims to lay the foundation for a global lycopsid model which could also estimate the variability of atmospheric CO₂ levels and climate on geological timescales on coupling with a global mass balance model of carbon pools. The current study hints at a potential enhancement of the weathering by the lycopsids considerably higher than that of lichens but requires further exploration, especially keeping in mind the nutrient limitation that has not been incorporated in LYCOm. This model therefore can be applied toward investigations targeting weathering dynamics in early terrestrial ecosystems. Hence, lycopsids could turn out to be a crucial contributor to the enhancement of silicate weathering on a global scale, especially, during the periods in which the species dominated the flora. Hence, this implication of the model for the fossil record is imminent and we, therefore, try to explore the past using information from such records. Further studies delving into the geochemical changes resulting from the advent of lycopsids could give us detailed insights into understanding the gradual evolution of the composition of the atmosphere today.

The study is restricted by two primary factors. Firstly, the current coarse resolution of the dataset in use poses the risk of averaging the real behavior of the plant species and therefore might mask a higher potential impact of the lycopsids on weathering rates. Secondly, poorly sampled physiological properties of lycopsids make it difficult to evaluate the model. Further enhancement and extension of the lycopsid datasets are thus required, which would require a collective effort of the research community. The current work attempts to map out a relatively new field of research that has shown to hold much potential based on this preliminary work. The work aims to draw scientific interest in the field and try to incorporate these details in works such as the Paleoclimate Modelling Intercomparison Project (Otto-Bliesner et al., 2021).

Appendix A

Table A1 Comparison of NPP: (i) mean, which refers to the average of all surviving species; (ii) weighted, according to Eq. (24), averaged by root and leaf biomass; (iii) dominant species at the location which accumulates the maximum biomass during the simulation period in gram carbon per square meter per year.