Discerning the thermodynamic feasibility of the spontaneous coexistence of multiple functional vegetation groups
Experimental design
A multi-layer canopy-root-soil model (MLCan)24,26,27 is used to calculate the energy and entropy fluxes for three climatologically-different ecosystems containing multiple functional groups: water-limited Santa Rita Mesquite (SRM), energy-limited Willow Creek (WCR), and nutrient-limited Tapajos National Forest (TAP)38.
MLCan takes site-specific parameters and weather forcing data and computes the energy and entropy fluxes and temperatures for each of the ecosystem layers. Entropy calculations are based on both the energy fluxes and temperature of soil, air, and leaves (see Entropy Calculations). The model is run for a simulation period of 2 years (2004–2005) at a half-hourly timescale for SRM and WCR and an hourly timescale for TAP due to data availability. Weather forcing data was downloaded from FLUXNET2015: air temperature, air pressure, global radiation, precipitation, wind speed, friction velocity, and relative humidity32,33,34. Additional model input parameters can be found in Table S2 of the Supplementary Information.
The initial soil moisture and temperature profiles for each of the sites—and snow properties for WCR—were produced from a spin-up of the model. The WCR and TAP sites used 2004 LAI with 2003 forcing data for a spin-up of 2 years to provide the initial conditions for the beginning of the 2004 simulation. For the SRM site, the FLUXNET2015 data was not available for 2003, so 2004 data was used instead.
At each site, the model splits up the vegetation into plant functional groups. Domingues et al.35 demonstrates the importance of modeling ecosystems based on functional groups. For WCR and SRM, the vegetation is represented by understory herbaceous species and overstory trees. For TAP, a high biodiversity ecosystem in Amazonia, the vegetation is further divided and represented by four groups: understory tree, mid-canopy tree, upper-canopy tree, and upper-canopy liana35. See Table S1 of the Supplementary Information for functional group abbreviations.
The LAI data for all sites are taken from MODIS39 and calibrated based on site documentation (Fig. S4 of the Supplementary Information). The LAI is then partitioned into two or four components based on the number of functional groups at each site. Additional LAI information can be found in the Supplementary Information.
MLCan has been previously validated for each of the sites considered30,40. Since entropy cannot be directly measured, we provide a comparison of the model outputted latent heat fluxes with the observed fluxes at each site in Fig. S5 of the Supplementary Information for additional validation.
Site descriptions
The SRM site is located on the Santa Rita Experimental Range in southern Arizona ((31.8214^{circ }hbox{N}), (110.8661^{circ }hbox{W})). SRM has a hot semi-arid climate and consists of woody savannas with mesquite trees (Prosopis velutina Woot.) and C4 grasses and subshrubs40,41.
The WCR (Willow Creek) site is located within the Chequamegon-Nicolet National Forest in northern Wisconsin ((45.8059^{circ }hbox{N}), (90.0799^{circ }hbox{W})) with a northern continental climate. It is a deciduous broadleaf forest dominated by sugar maple (Acer saccharum Marsh.) with understory shrubs, including bracken ferns (Pteridium aquilinum), and overstory seedlings and saplings42,43,44.
The TAP (Tapajos National Forest) site data is taken from the Santarem Km 67 Primary Forest site located in Belterra, Pará, Brazil ((2.8567^{circ }hbox{S}), (54.9589^{circ }hbox{W})). This evergreen broadleaf forest in Amazonian Brazil has a tropical monsoon climate with vegetation consisting of dozens of known tree species and lianas30,35.
Entropy calculations
Entropy calculations are based on model-simulated temperature and energy at each of the 20 canopy layers and the soil-surface layer, and results are scaled up to the ecosystem level. No lateral exchange of fluxes are considered. The net sum of energy fluxes from all layers of the ecosystem is equivalent to the total flux of energy across the boundary of the control volume (Fig. S1 of the Supplementary Information). These energy fluxes include shortwave radiation (SW), longwave radiation (LW), latent heat (LE), and sensible heat (H). All results are categorized as the flux of energy at the boundary entering ((SW_{in}), (LW_{in})) or leaving ((SW_{out}), (LW_{out}), LE, H) the ecosystem. Because the total energy flux across the ecosystem boundary is equal to the sum across the canopy layers in the model, the total entropy flux across the boundary can also be taken as the cumulative sum of the entropy fluxes from all layers of the ecosystem.
Entropy flux calculations are summarized in Table 1. All energy variables have units of (hbox{W/m}^2), entropy variables are in (hbox{W/m}^2hbox{K}), and temperatures are in K.
Table 1 Entropy calculations
Full size table
Entropy for LE and H calculations are based on simple heat transfer. The change in entropy is:
$$begin{aligned} dS=frac{dQ}{T} end{aligned}$$
(1)
where dQ is change in heat and T is temperature49. Thus, the flux of entropy for a given energy flux (E) across a boundary is:
$$begin{aligned} J=frac{E}{T} . end{aligned}$$
(2)
However, thermal radiation (SW and LW) cannot be treated this simply. The entropy flux for blackbody radiation is:
$$begin{aligned} J_{BR} = frac{4}{3} sigma T^3 = frac{4}{3} frac{E_{BR}}{T} end{aligned}$$
(3)
where (sigma) is the Stefan–Boltzmann constant, and (E_{BR}) is the blackbody radiation flux defined as (sigma T^4) from the Stefan–Boltzmann Law48,49.
SW is considered blackbody radiation, and entropy fluxes for direct shortwave radiation ((J_{SW,direct})) can be obtained by Eq. 3. However, LW is considered non-blackbody radiation, also called ‘diluted blackbody radiation’, which must include an additional factor (X(epsilon )) to account for the entropy produced during the ‘diluted emission’ of radiation given by an object’s emissivity, (epsilon). This factor is defined as45,46:
$$begin{aligned} X(epsilon ) = 1-Big [frac{45}{4pi ^4}ln {(epsilon )}(2.336-0.26epsilon )Big ]. end{aligned}$$
(4)
Although (SW_{diffuse}) is still a blackbody radiation, it has been demonstrated47 that the entropy flux due to (SW_{diffuse}) can be treated similarly to non-blackbody radiation with a new variable, (xi), in place of emissivity. (xi) is the ‘dilution factor’ of radiation due to scattering, meaning it is the ratio of diffuse solar radiance on Earth’s surface to solar radiance in extraterrestrial space47. Since diluted blackbody radiation ((SW_{diffuse})) is mathematically equivalent to non-blackbody radiation (LW) when the dilution factor is equal to the emissivity, (xi) can also be plugged into Eq. 4 to solve for the amplifying factor of entropy production due to scattering, (X(xi ))37,45,46,48.
Each of the entropy calculations in Table 1 have a temperature value corresponding to the temperature of the energy’s source. For instance, shortwave radiation originates from the sun, so the source temperature in its entropy equations is (T_{sun}). Likewise, longwave radiation is assumed to originate from the atmosphere, leading to a corresponding temperature of (T_{atm}). However, (LW_{out}), LE, and H do not have a single source location, so we must calculate an equivalent temperature ((T_{eq})) for each energy category based on the modeled temperatures and weighted contribution of each layer to the total energy flux at the ecosystem boundary. The equivalent temperatures for these three energy categories are calculated as follows:
$$begin{aligned} T_{eq,j} = sum _{k=1}^{21}[T_{k} times omega _{j, k}] end{aligned}$$
(5)
where (T_{eq,j}) is the equivalent temperature of energy category j such that (j in {LW_{out}, LE, H}). k refers to the layer in the ecosystem such that layers 1-20 are the canopy layers, and layer 21 refers to the ground surface. (T_k) is the temperature of layer k, and (omega _{j,k}) is the weight of energy category j coming from layer k given by:
$$begin{aligned} omega _{j,k} = frac{E_{j,k}}{E_{j,eco}} end{aligned}$$
(6)
where (E_{j,k}) is the energy j leaving layer k, and (E_{j,eco}) is the total energy j leaving the ecosystem.
The total entropy flux of the ecosystem ((J_{eco})) is calculated by summing the energy categories:
$$begin{aligned} J_{eco} = sum J_j + J_{SWout} end{aligned}$$
(7)
where (J_{SWout}) is the entropy flux of diffuse shortwave radiation leaving the ecosystem. The entropy flux per unit energy (EUE) is another way to view the thermodynamic state of ecosystem vegetation. EUE is calculated as:
$$begin{aligned} EUE_{j} = frac{J_{j}}{E_{j}} end{aligned}$$
(8)
where (EUE_{j}) is the entropy per unit energy in 1/K of energy category j. It follows that the corresponding (EUE_{SWout} = J_{SWout}/E_{SWout}), and the total ecosystem EUE is:
$$begin{aligned} EUE_{eco} = frac{sum J_j + J_{SWout}}{sum E_j + E_{SWout}}. end{aligned}$$
(9)
Work calculations
Work in an ecosystem represents the energy required to directly perform motion in the form of heat, effectively decreasing the temperature gradient within the ecosystem. We assume that LE and H are the primary regulators of temperature within a natural ecosystem, and (LW_{out}) is wasted energy. Additionally, we assume that the bottom of the control volume is sufficiently deep such that the temperature at the boundary is consistent and there is no loss of heat (i.e. ground heat flux is ignored). Thus, work is estimated and calculated directly from LE, H, and change in internal energy due to photosynthesis, (Delta Q):
$$begin{aligned} W = LE+H+Delta Q end{aligned}$$
(10)
where (Delta Q) is significantly less than LE and H and can be ignored. So work can be simplified to:
$$begin{aligned} W = LE+H. end{aligned}$$
(11)
Since work represents the ability of an ecosystem’s vegetation to deplete the driving temperature gradient imposed upon the ecosystem, our analysis compares work with temperature gradient. We define temperature gradient as:
$$begin{aligned} frac{Delta T}{Delta z} = frac{T_{surf}-T_{air}}{h_e} end{aligned}$$
(12)
where (T_{surf}) is the temperature of the soil surface, (T_{air}) is the temperature of the air in the top layer of the ecosystem, and (h_e) is the ecosystem height (see Table S2 in the Supplementary Information).
Work efficiency is the work performed for the amount of radiation entering the ecosystem defined as:
$$begin{aligned} WE = frac{LE+H}{E_{SWin} + E_{LWin}} = frac{W}{E_{in}}. end{aligned}$$
(13)
Since each vegetation functional group partitions energy differently among the energy categories, work efficiency is a good way to compare thermodynamic behavior across model scenarios at each site in a normalized way.
Statistical analysis
To determine if the differences of entropy flux and work efficiency among scenarios at each site are statistically significant, we perform two separate tests for entropy flux and work efficiency. Since entropy flux distributions are positively skewed (Fig. 1a), we use the variance as an indicator of the difference between them. To this end we use the distribution-free Miller Jackknife (MJ) significance test50,51 for variance that does not assume that the distributions come from populations with the same median. However, the work efficiency distributions exhibit no such pattern (Fig. 1b), and, therefore, we use the two-sample Kolmogorov–Smirnov (KS) test, which measures the maximum absolute difference between two empirical cumulative distribution functions (CDF)52,53,54.
First, the entropy flux variances are compared with the MJ test. Because functional group scenarios at each site are bounded on the lower end by similar values, if a distribution has a larger variance than another, then the two populations cannot be considered as coming from the same continuous distribution, and the distribution with a larger variance generally consists of larger values. For each site we test the null hypothesis, (H_0), that the distribution of multiple-functional-group entropy fluxes and the distribution for each of its single-functional-groups have the same variance. This is done with each functional group present at each site (Table S1). The alternate hypothesis, (H_{A1}), states that the distribution of entropy fluxes from the multiple-functional-group has a larger variance than that of the corresponding single-functional-group, meaning that the two populations do not belong to the same distribution and the multi-group scenario consists of larger values than the single-group scenario. The results from this test, shown in Table 2, indicate that (H_0) is rejected in favor of (H_{A1}) at the 5% level ((p More
