Discerning the thermodynamic feasibility of the spontaneous coexistence of multiple functional vegetation groups

Can the Second Law of Thermodynamics explain why ecosystems naturally organize into a complex structure composed of multiple vegetation species and functional groups? Ecosystem structure, which refers to the number and type of plant functional groups, is the result of self-organization, or the spontaneous emergence of order from random fluctuations. By considering ecosystems as open thermodynamic systems, we model and study these fluctuations of throughput signatures on short timescales to determine the drivers and characteristics of ecosystem structure. This diagnostic approach allows us to use fluxes of energy and entropy to calculate an ecosystem’s estimated work and understand the thermodynamic behavior of the system. We use a multi-layer canopy-root-soil model to calculate the energy and entropy fluxes of different scenarios for field sites across various climates. At each site, scenarios comprised of native individual plant functional groups and a coexisting multi-group composition scenario including all functional groups observed at the site are compared. Ecosystem-scale calculations demonstrate that entropy fluxes and work efficiency—the work performed for the amount of radiation entering the ecosystem—are greatest in the multi-group scenario when its leaf area is significantly larger than each of its individual functional groups. Thus, we conclude that ecosystems self-organize towards the vegetation structure with the greatest outgoing entropy flux and work efficiency, resulting in the coexistence of multiple functional groups and performing the maximum amount of work within the constraints of locally available energy, water, and nutrients.


Scientific Reports
| (2020) 10:18321 | https://doi.org/10.1038/s41598-020-75050-4 www.nature.com/scientificreports/ sites are modeled with two functional groups based on the composition of their dominant vegetation, and the TAP vegetation is divided into four functional groups (based on the details in Domingues et al. 35 ). For simplicity, the functional groups are abbreviated as: understory (UN), mid-canopy trees (MT), overstory trees (OT), and lianas (L). The scenario with multiple functional groups is abbreviated as MG. Site-specific classifications can be found in Table S1 of the Supplementary Information. At each location we compare the existing scenario of vegetation consisting of multiple functional groups with hypothetical scenarios of each one of the individual functional groups present. For all scenarios, we calculate the entropy flux ( J eco ; see Eq. 7), work (W; see Eq. 11), and work efficiency (WE; see Eq. 13). We define thermodynamic advantage as the production of larger entropy fluxes as well as greater work efficiency by the ecosystem as a whole. We hypothesize that the multiple-functionalgroup systems are more thermodynamically advantageous than or similar to their respective single-group scenarios.

Results
The distributions of entropy fluxes at each timestep and the distributions of daily work efficiency for the 2-year study period are shown in Fig. 1. The three sites have different ranges of entropy fluxes and work efficiencies due to distinctions in the local availabilities of water, energy, and nutrients. WCR, SRM, and TAP have relatively energy-limited, water-limited, and nutrient-limited environments, respectively. Thus, the entropy and work efficiencies should not be directly compared across sites. However, within each site the ranges gradually change as we look across the different functional group scenarios. SRM has the least variability amongst its functional groups, and TAP has the most. Considering the entropy fluxes and work efficiencies, the multi-group scenarios all appear to have distributions similar to or consisting of larger values than the other scenarios for each site.  Table S1). The distributions of entropy fluxes are developed based on each time step of simulation for the 2-year study period, and the distributions of work efficiency are calculated based on daily energy fluxes (see Eq. 13).
Means are shown as black diamonds. www.nature.com/scientificreports/ Using the Miller Jacknife and Kolmogorov-Smirnov tests, the distributions of entropy fluxes and work efficiency for the multiple functional group scenarios at each site are compared for statistically significant differences with each of their individual functional groups. Overall, the results indicate that multiple-functional-groups either have a thermodynamic advantage over single-groups or they are not at a disadvantage due to greater or similar values of entropy flux and work efficiency. There was a statistically significant difference in the distributions of entropy and work efficiency between the MG scenarios and each of the individual functional groups except for WCR-OT. This case is unique since WCR-UN only contributes to 3% of the total leaf area of WCR-MG, meaning that the WCR-OT and WCR-MG are very similar in vegetation composition. For all other cases, the values of entropy and work efficiency were significantly larger in the MG scenarios, indicating thermodynamic advantage. The "Statistical analysis" section in the "Methods" lays out the statistical tests and related analyses for comparing these entropy flux and work efficiency distributions.
Factors impacting entropy flux. The behavior of an ecosystem's entropy is determined by the combined variation of its individual energy fluxes leaving the system, such as: shortwave radiation (SW), longwave radiation (LW), latent heat (LE), and sensible heat (H). For each scenario at each site, an entropy per unit energy (EUE) value is computed for all energy fluxes (see Eq. 8) and reported as an average over the simulation period in Fig. 2a. Each ecosystem has its own partitioning of energy among these categories, leading to differences in total entropy per unit energy leaving the ecosystem ( EUE eco ), corresponding to "Total Out" in Fig. 2a.
The EUE for each energy category can be explained by the temperature of its source. SW originates from the sun, so its EUE is based on the temperature of the sun. However, outgoing longwave radiation ( LW out ), LE, and H originate from the leaves as well as the soil surface. The total outgoing energy of the ecosystem is a resultant of the energy leaving each of the canopy layers including the soil surface within MLCan, each with its own temperature. Therefore, we calculate temperature equivalences ( T eq ) for each of these energy categories based on the weighted average of the temperature of each leaf and soil layer contributing to the overall energy flux of that category (Fig. 2b, Eq. 5). Entropy is calculated directly from temperature (Table 1), leading to an important inverse relationship between EUE and T eq . Figure 2b shows that for each of the sites, across almost all energy categories the UN scenario has the highest T eq and the multi-group scenario has the lowest. Due to the inverse relationship between T eq and EUE, all of the EUE averages for each category in Fig. 2a are smallest for the UN scenarios and largest in the multi-group scenarios, though some of these differences are marginal. From this pattern, one would expect that this would lead to a clearly greater overall EUE eco for all multi-group scenarios ("Total Out" in Fig. 2a). Yet, this is not the case for all sites, as TAP is the only one in which the average EUE eco varies considerably among functional group scenarios.
Although EUE averages tend to increase with leaf area index (LAI) and when multiple functional groups coexist, each energy category has different relative values of EUE. Further, EUE can be interpreted as an indicator of how degraded a particular form of energy is and the ability of that energy to perform additional work. SW has the least EUE across all sites and scenarios (Fig. 2a), so it has the greatest capacity for work to be done; plants are able to use this energy to perform work. Yet, SW is radiative energy, so it does not perform work itself. On the other end, LW out has the highest EUE of the energy categories studied here, meaning it is the most degraded with little capacity for additional work to be performed from it. In the middle, LE and H are more degraded than SW, but they are still able to perform physical work in the ecosystem through convection and conduction and contribute to the redistribution of heat throughout the vertical profile of the ecosystem. Since we are only considering the fluxes that enter and leave the ecosystem control volume in this analysis and all irreversible work releases heat, LE and H are also proxies for the work performed within the ecosystem internally (i.e. by the vegetation itself). Overall, when an ecosystem has low thermodynamic efficiency and high LW out , it degrades the incoming SW quickly without performing much work within the ecosystem and constitutes wasted energy. However, higher values of LE and H leaving the ecosystem mean that more work has been performed.
Since each scenario partitions the incoming energy differently throughout the ecosystem, the considerably higher overall EUE for LW out has important implications for overall EUE eco and entropy fluxes, yielding larger values when more outgoing energy is allocated towards LW out . Figure 3a displays the average partitioning of incoming radiation among the various energy categories at TAP, and Fig. 3b indicates the corresponding entropy fluxes presented as percentages of the incoming entropy flux and disaggregated into different energy categories. The energy and associated entropy entering the ecosystem are the same for all functional groups, but the outgoing fluxes vary significantly among them. Due to conservation of energy, the outgoing energy for all scenarios corresponds to 100% of the incoming radiation. Alternatively for entropy, this percentage is greater than 100%, indicating entropy production (shaded in grey). At TAP and similarly at the other two sites, the multi-group scenario produces more total entropy on average than the other two scenarios.
However, the contribution of each energy category towards the total entropy differs among functional groups. As LAI increases from UN to OT and with the addition of multiple functional groups, the proportion of LE and H, or work, increases while the outgoing radiation decreases (Fig. 3). The UN scenario partitions more energy towards LW out than the OT and multi-group scenarios. Since LW has the largest EUE of all the energy categories (Fig. 2a), its percentages in Fig. 3b are larger than in Fig. 3a. This observation, consistent across all sites, indicates that UN scenarios are able to make up for lower performing EUE values by partitioning more energy towards the higher entropy-producing energy category, LW out . Thus, even though the UN scenario has the lowest EUE for all energy categories, its total EUE eco -and overall ability to degrade the incoming SW-and entropy production are similar to the other scenarios since it has more outgoing energy partitioned towards the largest EUE category, LW out . However, EUE eco is an indication of energy degradation, not work performed. Thus, the additional assessment of work and work efficiency is necessary for the interpretation of thermodynamic advantage. www.nature.com/scientificreports/ Work as an indicator of self-organization. As discussed in the previous section, partitioning of energy and entropy fluxes are important for understanding the overall thermodynamic behavior of ecosystems. Energy fluxes with large EUE values result in greater entropy production for an ecosystem but do not always yield more work performed. Work-estimated as the sum of the ecosystem's latent and sensible heat fluxes (W; Eq. 11)represents the ability of an ecosystem to diminish the temperature gradient through the ecosystem. Despite having the largest EUE value, LW out is not a component of work since it is a passive response to the temperature state. Figure 4 displays the relationships of work versus temperature gradient between the atmosphere and land surface ( �T/�z ; Eq. 12) for the functional group scenarios at each site. Work performed within an ecosystem has a positive nonlinear relationship with temperature gradient across all scenarios. Further, each functional group scenario is fitted to a power function: W = a(�T/�z) b . At each site a power law is observed; higher powers (b) correspond to functional groups with larger LAI, and except for the SRM-UN scenario, the highest power at each site corresponds to the multiple-functional-group scenario. This means that the work performed by the ecosystems with multiple functional groups has an exponentially greater response to marginal changes in temperature gradient as it increases.  www.nature.com/scientificreports/ To better understand the relationship between work and temperature gradients, we refer to the 1994 paper in which Schneider and Kay 10 theoretically explored Silveston's 36 Bénard cell experiments of heating an enclosed fluid from below. They demonstrated that without self-organization, a system's work performed from conduction alone has a linear relationship with the temperature gradient. However, when self-organization in the form of convection occurs at a critical point, Bénard cells form, the relationship becomes nonlinear, and more work is performed for each additional unit increase in gradient 10 . In this experiment, analysis of data in which no self-organization occurred at all demonstrated a linear relationship between work and gradient. Alternatively, when Bénard cells formed (i.e. self-organization occurred) power law relationships between work and gradient emerged with higher-degree powers corresponding to greater dissipation rates due to increased convection. We give this example not to study the emergence of the phenomena, but as a means of comparing the behavior of possible end states: with and without self-organization. If we take the functional group scenarios as possible end states with various levels of self-organization, the interpretation for each of our sites is consistent with those of the Bénard cells: with self-organization, there is a nonlinear relationship between work and gradient with larger exponents in the power law relationship corresponding to more advanced levels of self-organization. This yields exponentially more work performed in the more highly organized multiple-functional-group ecosystem scenarios. From the combined results of Schneider and Kay and Fig. 4, we infer that self-organization is the leading driver of the nonlinearity shown in work-gradient plots. This supports the proposition that the existence of multiple functional groups reflects a higher degree of self-organization that results in nonlinear increases of work performed in response to marginal increases in temperature gradient, reflecting thermodynamic advantage.

Discussion
Through the concepts of entropy and work efficiency, this paper establishes a framework for identifying thermodynamic advantage for the spontaneous self-organization of ecosystems towards a vegetation structure that includes multiple functional groups. We identify decreased canopy temperature, increased LAI, and greater partitioning of energy towards LW out as important factors amplifying the entropy production of an ecosystem. From these factors, one can deduce the relative changes in entropy flux and, thus, changes in thermodynamic behavior of an ecosystem. Entropy provides insights into the total disorder of a system. The second law of thermodynamics requires that closed systems yield increases in entropy over time. Although, ecosystems are open thermodynamic systems, scientists have correlated greater entropy production as a driver of self-organization 11,37 . However, not all energy and hence entropy is productive in terms of thermodynamic advantage. Thus, we provide work efficiency as another important component of thermodynamic advantage and the directionality of self-organization.
Work efficiency captures the ability of an ecosystem to perform work based on the energy throughflow of the system. Work without context of incoming radiation doesn't tell much about the performance of an ecosystem relative to others. Since all functional group scenarios at each site receive the same incoming radiation in this study, either work or work efficiency can be used as a metric of thermodynamic advantage. However, work efficiency is a more attractive metric for further study as it normalizes an ecosystem based on its local availability of energy, allowing for comparison of ecosystems across multiple climates. Thus, work efficiency provides promise for future research to compare ecosystems with varying energy availabilities and external environments directly.
Additionally, work efficiency measures an ecosystem's ability to rapidly convert incoming radiation into alternate forms of energy that disperse throughout the ecosystem control volume and diminish the imposed temperature gradient. Work efficiency helps us understand the reorganization of available energy entering an ecosystem towards thermodynamically-productive uses-meaning depleting the imposed temperature gradient. According to thermodynamic theory, all systems work to decrease gradients of their state variables, which in turn drive the movement of energy and mass from high concentrations to low concentrations (i.e. high to low temperatures) 13 . In this study, work exhibits a nonlinear power law relationship with temperature gradient. This means that exponentially more work is performed to combat the greater temperature differences from the Earth surface to the atmosphere above the canopy. For the sites studied here, scenarios with multiple functional groups exhibit the highest power law, meaning that the MG structure is more efficient at depleting the driving temperature gradient. This is a demonstration of high work efficiency. Since structures that perform more work for a given temperature gradient have a thermodynamic advantage over those with lower efficiencies, ecosystems have a natural tendency to self-organize to this MG structure.
Because ecosystems are formed and evolve through a process of random fluctuations, there is a nonzero statistical probability for the existence of any possible vegetation structure or state. The state exhibiting thermodynamic advantage identifies the state with the highest probability of occurrence. This does not mean that the advantageous structure will always result. In all of our sites, the highest work efficiency corresponds to the existing MG scenario. Thus, the ecosystems exist in the highest probability thermodynamically-advantageous state.
The outcomes of this work provide valuable insight into the self-organization of natural ecosystems. Thus far, we have identified that when multiple functional groups coexist this structure exhibits a thermodynamic advantage over other possible individual functional group scenarios. Thus, ecosystems will have a higher probability of self-organizing towards this greater work efficiency state. Additionally, this work highlights new areas for further study. The framework of thermodynamic advantage through greater entropy production and work efficiency could be applied to other ecosystem structures, such as the existence of individual functional groups in nature. Further, this framework could help scientists understand how human-induced perturbations could impact thermodynamic behavior and alter the most advantageous state. Thus, we propose the concepts of entropy and work efficiency as valuable contributions to the basic understanding of the existence of a particular vegetation structure and present thermodynamic advantage as a tool for future use in understanding and studying the stability and behavior of ecosystem self-organization.  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 humidity [32][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 liana 35 . See Table S1 of the Supplementary Information for functional group abbreviations.
The LAI data for all sites are taken from MODIS 39 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 considered 30,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. where dQ is change in heat and T is temperature 49 . Thus, the flux of entropy for a given energy flux (E) across a boundary is:

Site descriptions. The SRM site is located on the Santa Rita Experimental Range in southern Arizona
However, thermal radiation (SW and LW) cannot be treated this simply. The entropy flux for blackbody radiation is: where σ is the Stefan-Boltzmann constant, and E BR is the blackbody radiation flux defined as σ T 4 from the Stefan-Boltzmann Law 48,49 .
(3) www.nature.com/scientificreports/ 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(ǫ) to account for the entropy produced during the 'diluted emission' of radiation given by an object's emissivity, ǫ . This factor is defined as 45,46 : Although SW diffuse is still a blackbody radiation, it has been demonstrated 47 that the entropy flux due to SW diffuse can be treated similarly to non-blackbody radiation with a new variable, ξ , in place of emissivity. ξ 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 space 47 . Since diluted blackbody radiation ( SW diffuse ) is mathematically equivalent to non-blackbody radiation (LW) when the dilution factor is equal to the emissivity, ξ can also be plugged into Eq. 4 to solve for the amplifying factor of entropy production due to scattering, X(ξ ) 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: where T eq,j is the equivalent temperature of energy category j such that j ∈ {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 ω j,k is the weight of energy category j coming from layer k given by: 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: 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: 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: 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 www.nature.com/scientificreports/ 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, Q: where Q is significantly less than LE and H and can be ignored. So work can be simplified to: 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: 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: 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 test 50,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 < 0.05 ) for all scenarios except for the WCR-OT scenario. This indicates that for these ecosystems the distributions of entropy fluxes consist of larger values when multiple functional groups are present. www.nature.com/scientificreports/ Using the KS test to compare work efficiency distributions for each site, we test the null hypothesis, H 0 , that the multiple-functional-group measures of work efficiency and those for each of its single-functional-groups are from the same continuous distribution, or population. The alternate hypothesis, H A3 , states that the CDFs of the entropy flux from the multi-group scenario are smaller than those from the corresponding single-groups, meaning that the multi-group scenarios consist of values that are larger than their associated single-group scenarios. The results from this test, shown in Table 2, indicate that H 0 is rejected in favor of H A3 at the 5% level ( p < 0.05 ) for all scenarios except for the WCR-OT scenario. This indicates that the distributions for work efficiency are indeed larger when multiple functional groups are present in these ecosystems, as indicated by a smaller CDF (Fig. S3 of the Supplementary Information).
However, the tests of comparison for the WCR-OT scenario for both work efficiency and entropy flux distributions have p values larger than 0.05 (i.e. H 0 cannot be rejected at the 5% level). This means that we cannot say that the WCR multi-group entropy flux distribution has a variance larger than the OT single-group distribution or the multi-group work efficiency scenario comes from a larger distribution than the OT single-group scenario. This is not entirely surprising, as there is very little difference in LAI between these two scenarios; the maximum difference in LAI is about 0.2 (Fig. S4 in the Supplementary Information), only 3% compared to the total WCR-MG LAI. This small increase in LAI from the single to the multi-group scenario provides less opportunity for increased energy dissipation and hence entropy production due to the smaller understory. Thus for completeness, we also perform the MJ and KS tests in the opposite direction for the WCR-OT scenario with the following alternative hypotheses. For the MJ test on entropy flux variances, H A2 states that the distribution of entropy fluxes from the WCR multiple-functional-group has a smaller variance than that of the OT singlefunctional-group, meaning that the two populations do not belong to the same distribution and the multi-group scenario consists of smaller values than the single-group scenario. For the KS test on work efficiency, H A4 states that the CDF of the entropy fluxes from the WCR multiple-functional-group is larger than the CDF from the OT single-functional-group, meaning that the multi-group scenario consists of values that are smaller than the single-group scenario. The results from both tests, shown in Table 2, indicate that we again cannot reject H 0 at the 5% significance level for WCR-OT. Thus, although the WCR multi-group scenario compared to the OT scenario does not have a significantly greater entropy flux variance or a greater work efficiency distribution, it also does not have less variance or a smaller distribution. Overall, the test results indicate that multiple-functionalgroups have either greater or similar values of entropy flux and work efficiency than the modeled scenarios of their individual functional groups.