A dynamic local scale vegetation model for lycophytes (LYCOm)

Lycophytes (club mosses) represent a distinct lineage of vascular plants with a long evolutionary history including numerous extant and extinct species which started out as herbaceous plants and later evolved into woody plants. 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 (CO2) via root respiration and also modified soil hydrology. These effects had the potential to promote the dissolution of silicate minerals, thus intensifying silicate weathering. The weathering of 5 silicate rocks is considered one of the most significant geochemical regulators of atmospheric CO2 on 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 to the most basal plants with vascular tissue and shallow root system, 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 lycophytes. These properties, however, cannot be easily 10 derived from proxies such as fossil records, for instance. Hence, as a first step, a process-based model is developed here to estimate net carbon uptake by these organisms at the local scale, considering key features such as biomass distribution above and below ground, root distribution in soil regulating water uptake by plants besides, stomatal regulation of water loss and photosynthesis, and not withholding respiration in roots. The model features ranges of key physiological traits of lycophytes to predict the emerging characteristics of the lycophyte community under any given climate by implicitly simulating the process 15 of selection. In this way, also extinct plant communities can 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 lycophytes. We run the Lycophyte model, called LYCOm, at seven sites encompassing various climate zones under today’s climatic conditions. LYCOm is able to simulate realistic properties of lycophyte communities at the respective locations and estimates values of Net Primary Production (NPP) ranging from 126 g carbon m−2 year−1 to 245 g carbon m−2 year−1. Our 20 limit-based weathering model predicts a mean chemical weathering rate ranging from 5.3 to 45.1 cm ka−1 rock with lycophytes varying between different sites, as opposed to 0.6 8.3 cm rock ka−1 without lycophytes, thereby highlighting the potential importance of such vegetation at the local scale for enhancing chemical weathering. Our modeling study establishes a basis for 1 https://doi.org/10.5194/gmd-2021-202 Preprint. Discussion started: 21 July 2021 c © Author(s) 2021. CC BY 4.0 License.

. A conceptual diagram of evolution of roots from its earliest stage in the Silurian era to a sophisticated form during the Devonian notable consequence of root development is believed to be an enhancement of chemical weathering, which could have had strong impacts on atmospheric CO 2 content and the Global carbon budget.
In general, the weathering regime is often influenced by several factors, and plants, and their associated microorganisms have been suggested to have a substantial impact on weathering, through respiration in the soil, hence increasing the partial 60 pressure of soil CO 2 to several-fold greater than the atmospheric CO 2 level (Kelly et al., 1998;Montanez, 2013). In addition to plant-derived organic matter respired by these microorganisms, the root respiration also adds to the soil CO 2 content, thus increasing soil acidity and influencing weathering rates (Moulton et al., 2000;Andrews and Schlesinger, 2001;Berner et al., 2003).
Roots were not the only suite of innovations from these kinds of plants, they further developed stomatal control of water loss to be able to photosynthesize under adverse climatic conditions which might have given lycophytes yet another advantage over their contemporaries. The evolution of stomata followed the innovation of roots closely, around 400 million years before the present. Water loss and NPP are regulated by the opening and closing of stomata. By changing the aperture of the have enabled vascular plants to invade most terrestrial environments, tolerating water stress and exploiting favorable conditions (Raven, 2002). The development of such traits in lycophytes was gradual and was initially, relatively primitive. Brodribb and McAdam (2011) suggested that lycophyte 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 75 a more complex mechanism driving the stomatal control.
The 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 80 weathering in the interiors of continents (Shukla and Mintz, 1982;Lee, 2010, 2017). The moisture recycling likely also increased plant productivity, thus increasing soil CO 2 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 as one of the most significant geochemical regulators of atmospheric CO 2 on a long (hun-85 dreds 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 the weathering reactions is the production of weak carbonic acid from rainwater and soil CO 2 , which dissociates into bicarbonate (HCO 3 − ) ions during the dissolution of silicate minerals. The ions subsequently enter streams via 90 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.
This so-called 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 95 potential to promote the dissolution of silicate minerals (Walker et al., 1981). The study by Matsunaga and Tomescu (2016) demonstrates a substantial plant-substrate interaction that was already underway in the Devonian, causing weathering by root penetration, 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-folds as compared to their earlier counterpart, i.e. non-vascular vegetation (Cochran and Berner, 1993), along with an improved pedogenesis (Quirk et al., 2015).

100
Although lycophytes were abundant in the Silurian, 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 lycophytes Here, we introduce a process-based model, called LYCOm, which simulates lycophyte properties and is able to simulate 110 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 lycophytes. We test if the model can predict realistic lycophyte 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 lycophytes on hydrological processes as well as weathering of rocks, as the first step to a more quantitative understanding of 115 the impacts of early vascular vegetation on global biogeochemical cycles and climate of the past.

Lycophyte Model (LYCOm)
The process-based vegetation model used in the current study aims at a general representation of lycophytes for estimating the productivity and physiological properties of these organisms under a broad climatic range. LYCOm follows other dynamic 120 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 the roots and shoots with leaves. The balance of these pools depends on Net Primary Productivity (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 125 Deficit (VPD) as the controlling factor for the calculation of photosynthesis (Sellers et al., 1996;Lawson et al., 2010). The main reason for this choice is that VPD does not capture the possibility of water stress of plants due to transpiration of water into saturated air (at low VPD), which is driven by the energy balance, and also water stress due to low soil moisture. Since we ultimately want to apply LYCOm to the Silurian 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 130 not hold for the geological past. Another difference between LYCOm and other DGVMs is the representation of leaves and stems. To account for the particular morphology of lycophytes, the above-ground 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 a 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.  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 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 lycophytes in the real world but represent the diversity of physiological strategies for lycophytes instead.   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 the JeDi-DGVM (Pavlick et al., 2013), and was further developed in the LiBry-DGVM for lichens and bryophytes (Porada et al., 2013). LYCOm is based in parts on the LiBry model. Another  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 lycophytes under any given recent or past climate.

155
To determine both the spatial distribution of lycophytes as well as their impacts on weathering rates, it is necessary to quantify the carbon balance of the organisms and its dependence on climatic conditions. LYCOm includes three carbon pools for lycophytes: below-ground (root) biomass and above-ground biomass, which is subdivided into a pool for leaves and another for stems. Gains in biomass result from Net primary production (NPP), while losses are due to mortality, which includes both above-ground litterfall and root turnover. NPP is calculated as the difference between photosynthesis and respiration 160 (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.

165
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 5 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 lycophytes that differ in their root profiles. Since extant lycophytes, however, can only reach a maximum rooting depth of approximately 15 cm (Matsunaga and 170 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 175 contained in the layer (equation 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 into the surrounding soil, and does not immediately contribute to surface runoff. Only in case, the soil which then contributes to total runoff.
where, Qin= Incoming water into layers (rain + slow melt -Interception by surrounding vegetation (for first layer)),p dt In addition to infiltration and percolation, root water uptake changes the water content of the layered soil in the model.
Uptake of water by roots 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 (equation 4). Additionally, root uptake is limited to the available water in each layer. The water balance of the bucket is then calculated according to the following equation: where, Wx = Updated bucket water content, WIx= Initial bucket water content, Qp= Percolation from the upper layers, Qs= Overflow from top layer, Qb=Actual base flow , Qo= Large scale runoff comprising of the excess water from the saturated bucket 205

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 vapor pressure 210 deficit 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 not be valid anymore. Our approach is flexible enough to be applied to palaeo-climatic scenarios. It thereby depends on the assumption that the plant community is adapted to the prevailing average conditions, and thus using our normalized measure 215 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 equation 7.
ETratio=Relative atmospheric demand to the whole run-time, ETdata= run-time atmospheric water demand, ETdavg=average atmospheric water demand for the whole runtime The stomatal conductance is then computed as: gsleaf = gs0, when, ET ratio <= 1, where, gs0 = maximum stomatal conductance, ETratio= (Potential evapotranspiration at the given timestep/ Average evapo-225 transpiration), p gs1 =shape parameter , ETrmax= Maximum ETratio for the entire simulation, 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 230 soil water stress on stomatal conductance can be taken into account, in addition to the atmospheric drivers. gsroot = gs0 · S1 (9) where, gsroot= stomatal conductance determined by average layer water content, gsfinal= stomatal conductance limited both 235 by soil water availability and atmospheric demand, S1 = average saturation of the layered soil part.

Above and below ground biomass
In LYCOm, plant growth corresponds to Net Primary Production (NPP), which is based on the simulated net photosynthesis averaged over one month. Thereby, net photosynthesis is computed as gross photosynthesis minus respiration according to Farquhar and Von Caemmerer (1982).
where, f GP P,L =Light limited gross primary production , f GP P,CO2 =Carbon dioxide limited gross primary production, vcm= maximum specific carboxylation rate, J = actual rate of electron transport, CO2=CO2 concentration , O 2 = Oxygen concen- Gross photosynthesis is simulated as the minimum of a light-limited and a CO 2 -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 2 -limited rate is an increasing function of the CO 2 concentration in the chloroplasts and saturates at high values of CO 2 concentration, the maximum rate being VC max . The detailed formulation can be found in Porada et al. (2013).
where, 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 10 (eq. 15) (Vanderwel et al., 2015). The Gross and Net Primary Production (GPP and NPP) are calculated: where, NPP=Net primary production, R= Respiration, R0= Reference respiration rate at 10 • C, Q 10 = Q 10 value of respiration, T= Temperature on the surface of the plant, Topt= Optimum temperature 260 The pore space CO 2 concentration inside the leaves is a function of the stomatal conductance which is in turn controlled by the water availability in soil (eq. 10). Thereafter, we have considered a steady-state assumption of pore space CO 2 between diffusion and the net photosynthesis which is in turn dependent on the respiration. This drives the light and CO 2 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 below-ground 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 2 -limited (i.e. water-limited) conditions. The light- where NPP = Net primary productivity, 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 0 = shape constant [0.25] for a exponential decrease of root biomass allocation with depth,B root = available root biomass which is still available for allocation in this layer where, B rooti = root biomass of the layer i, IB rooti = initial root biomass of the i th layer, B rootin = Rootbiomassavailable, M ortality i = a term based on baseline mortality and the biomass allocation in soil layers The assigning of root biomass to litterfall (equation 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 (Eqn. 19). This adds to the mortality calculated in equation 18 which is in turn utilized to determine the soil carbon 290 pool and later the weathering regime.
LYCOm considers the root respiration as well, which potentially could further affect the soil CO 2 content (equation 20).
Root respiration (equation 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.
R resp = Root respiration, R0= Respiration determined by strategy-specific parameters, Q 10 = Q10 value, T= Soil temperature, Topt= Optimum temperature, B rooti =Root biomass of layer, sum(B root )=Total root biomass Besides roots, the hypothetical plants in the model also consist of a shoot and leaf partitioned in three to seven ratio of the 300 total available above-ground biomass which is represented as a cylinder to emulate the vertical structure of the given plant strategies along with Leaf Area Index (LAI). The leaf to shoot ratio is analogous to the horizontal to vertical biomass ratio value derived from Zier et al. (2015). Lycophytes exhibit true roots, stems, and also scale-like leaves called microphylls which cover the above-ground part of the plant. 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 approximately comprise 70 305 % of the above-ground 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 (Eqn. 21) is used thereafter to estimate the LAI which is computed using equation 22.

310
where Bl = leaf biomass , Bli= initial leaf biomass, Mortality s hoot =a term based on baseline mortality and the biomass production and allocation into above ground biomass Total Leaf Area= Aleaf · (Bl / Drw), Aleaf= Leaf area using measurement from Valdespino (2015a) and Drw=Dry leaf weight as measured in Leopold et al. (1981), GrA=Ground area (below)

315
The projection of the ground area covered by plants used for the calculation of the LAI used in the model assumes a random value ranging between 30 to 50 square centimeters. 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.

Weathering and soil carbon dioxide
To estimate chemical weathering rates and their biotic enhancement, we apply a simple limit-based approach Arens and Kleidon 320 (2011), which has already 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 325 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 4).
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 2 in the soil solution on the equilibrium concentration of weathering 330 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 2 mainly via two processes: Firstly, they directly release CO 2 into the soil as a result of root respiration, and, secondly, they enhance microbial respiration due to 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 2 via mortality (equation 16 -335 18) and root respiration (equation 19), and the outgassing of the CO 2 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 the weathering. Subsequently, the concentration of CO 2 in the soil can be determined from the known atmospheric CO 2 -concentration and a soil diffusivity equation (Jabro et al., 2012) as shown in Eqn. 23.

Physiological strategies
The variation of physiological properties in the model is summarised in table 1. Some photosynthetic parameters such as the 345 stomatal conductance responsible for the diffusion of atmospheric CO 2 into the leaves is limited as shown in the aforementioned table. The upper limit of the optimal temperature is limited between 0 and 50 • C along with a restricted Q 10 value between 1.5 -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 (Ordovician) simulations when the lycophytes appeared and the mean surface temperature was comparatively higher than today (Cocks and Torsvik, 2020). Furthermore, the physical realism of the lycophytes is asserted 350 by constraining certain parameters 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 lycophytes. The rest of the scheme is analogous to the model described in (Porada et al., 2013).

Weighting scheme
The weighting scheme used to calculate the resulting productivity depends on the allocation of root and above-ground biomass This yields an averaged representative NPP (table 3 ) and soil CO 2 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 summarised.

375
We run LYCOm locally for seven locations (listed in table 2) with a 2.812-degree spatial resolution data for 30 years between 1958 through 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., 2018). Besides, the model setup is initialized with close to present atmospheric composition of CO 2 concentration of 400 ppm and an oxygen content of 210000 ppm. The soil porosity 380 is constrained to a relative value of 0.45 with a percolation rate of 0.5x10 −7 m s −1 between the layers and an average baseflow of 1.0x10 −8 m s −1 from the bucket.

Model validation
The chosen study sites encompass various climate zones, to assess if the model can capture the productivity and viability of the local lycophyte communities. The chosen sites are characterized by a high occurrence of lycophytes, as determined based 385 on the Global Biodiversity Information Facility (GBIF, www.gbif.org).
The study by 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 lycophytes is well documented there is a significant lack of literature concerning biochemical properties and also the characteristic productivity range of lycophytes under various climate regimes. Hence,

390
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 (gbif.org). A Similar study by Campany et al. (2019) in Costa Rica affirms a robust model performance.  Table 1.

Results
The modeled lycophytes 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 2 , and temperature dependence of productivity, for model calibration ( Figure   5). Out of 100 strategies a maximum of 46 species survive in New Zealand followed by the USA which features 41 surviving 410 strategies. The climate of Sweden suits 24 varieties of the generated species, and that in Costa Rica the number fixates to 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  We averaged the weighted monthly productivity over a 30-year seasonal cycle to obtain an average measure of annual plant growth (table 3). We use a weighting scheme (see methods) to generate weighted NPP to obtain a mean seasonal cycle over 30 years (Figure 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 figure 4.
The assessment of the net assimilation rate is undertaken to evaluate the model performance as mentioned in the sec-435 tion model calibration, to verify its capability to reproduce the lycophyte-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) and reached saturation at 600 µmol m −2 s −1 for all strategies with a steeper slope as opposed to Soni et al. (2012) where they observed a saturation around 800 µmol m −2 s −1 . The temperature response studies in Soni et al. (2012) showed a gradual increase until 40°C and dropped marginally beyond that, a behavior 440 well captured in our calibration (see Fig. 3) with a slightly lower maximum (around 7 µmol m −2 s −1 ) in all strategies. Photosynthetic response to carbon dioxide levels is 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 lycophyte. This gives a relatively good idea of the model behavior under various conditions for unique strategies.

445
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  i.e. lycophytes. In addition, the soil water saturation determines the amount of CO 2 diffusing out of the soil. The water in the soil acts as an inhibitor to the diffusion of CO 2 through pathways of available pore space. The more the soil saturates, the lower are the CO 2 emissions. The concentration of soil CO 2 (see fig 4) is thus coherently linked with the rainy season and productivity (see 4). This effect can be prominently noticed at some locations, namely, in Costa Rica, when both the productive season and most rainfall coincide leading to a high content of soil CO 2 in June. The effect of both productivity and rainfall on 460 carbon content in soil can be seen in some cases as well, such as in Sweden and Estonia, where the most productive months for lycophytes are May and June but soil CO 2 peaks slightly in August in the rainy months.
The resulting soil CO 2 concentration after diffusion forms the basis of the weathering estimate in the post-processing of the model.  LYCOm represents a community of lycophytes at several local sites. The respective locations provide suitable climatic conditions for the growth of lycophytes (GBIF), and the model shows relatively high survival rates there, as pointed out in table 2.

495
The tropical climate of Costa Rica and Peru ensure continuous productivity throughout the year, while the temperate climate in 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. Besides, a prominent impact of soil water can be noticed in the soil CO 2 concentration. The microclimate at the site chosen for emulating the lab sample in India is likely to be substantially 500 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 lycophyte 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) (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 Lycophytes on chemical weathering rates is substantial at the local scale (table 4). It is in fact, significantly higher than their predecessors namely, early forms of lichens and bryophytes. Aghamiri and 510 Schwartzman (2002) explores weathering by lichens and bryophytes on rock surfaces and estimates a ranges of 0.0004 to 0.01 mm rock yr −1 weathering and is comparatively lower than that of lycophytes which exhibit a greater weathering potential of 0.026 to 0.418 mm rock yr −1 (see table 4).
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 figure 5 four species generated in LYCOm 515 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 (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.   restricted by a lack of data on the properties of lycophytes. Data concerning the plant physiology as well as measurements involving productivity will contribute towards a better optimization and validation of the model. Light, CO 2 , and temperature curves of lycophyte species other than Selaginella bryopteris would not only strengthen the model performance further but help us constrain parameter ranges with better precision.

530
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 2 concentration are dependent on productivity as well as hydrological balance at the respective sites. Table 6 summarises the sensitivity of weathering with increase and decrease in these components. We do not observe any 535 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.

540
The process-based lycophyte model LYCOm performs reasonably at the local scale. It is possible for lycophyte traits to survive in various climatic conditions and these can be represented well by the model. The productivity from the model shows strong agreement with the on-site measurements for the contemporary lycophytes. LYCOm, therefore, has the potential to be utilized beyond the range of current conditions. The study aims to lay the foundation for a global lycophyte model which could also estimate the variability of atmospheric CO 2 levels and climate on geological time scales on coupling with a global mass balance 545 model of carbon pools. The current study hints at a potential enhancement of the weathering by the lycophytes considerably higher than that of lichens but requires further exploration especially keeping in mind the nutrient limitation that has not been incorporated in LYCOm. Hence, lycophytes could 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. Further studies delving into the geochemical changes as a result of the advent of lycophytes could give us detailed insights into understanding the gradual 550 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 a risk of averaging the real behavior of the plant species and therefore might mask out a higher potential impact of the lycophytes on weathering rates. Secondly, poorly sampled physiological properties of lycophytes make it difficult to evaluate the model.
Further enhancement and extension of the lycophyte datasets are thus required, which would require a collective effort of the 555 research community. The current work attempts to map out a relatively new field of research that has shown to hold a lot of potentials 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).