Multi‐hypothesis comparison of Farquhar and Collatz photosynthesis models reveals the unexpected influence of empirical assumptions at leaf and global scales

Abstract Mechanistic photosynthesis models are at the heart of terrestrial biosphere models (TBMs) simulating the daily, monthly, annual and decadal rhythms of carbon assimilation (A). These models are founded on robust mathematical hypotheses that describe how A responds to changes in light and atmospheric CO2 concentration. Two predominant photosynthesis models are in common usage: Farquhar (FvCB) and Collatz (CBGB). However, a detailed quantitative comparison of these two models has never been undertaken. In this study, we unify the FvCB and CBGB models to a common parameter set and use novel multi‐hypothesis methods (that account for both hypothesis and parameter variability) for process‐level sensitivity analysis. These models represent three key biological processes: carboxylation, electron transport, triose phosphate use (TPU) and an additional model process: limiting‐rate selection. Each of the four processes comprises 1–3 alternative hypotheses giving 12 possible individual models with a total of 14 parameters. To broaden inference, TBM simulations were run and novel, high‐resolution photosynthesis measurements were made. We show that parameters associated with carboxylation are the most influential parameters but also reveal the surprising and marked dominance of the limiting‐rate selection process (accounting for 57% of the variation in A vs. 22% for carboxylation). The limiting‐rate selection assumption proposed by CBGB smooths the transition between limiting rates and always reduces A below the minimum of all potentially limiting rates, by up to 25%, effectively imposing a fourth limitation on A. Evaluation of the CBGB smoothing function in three TBMs demonstrated a reduction in global A by 4%–10%, equivalent to 50%–160% of current annual fossil fuel emissions. This analysis reveals a surprising and previously unquantified influence of a process that has been integral to many TBMs for decades, highlighting the value of multi‐hypothesis methods.


| INTRODUC TI ON
As the gateway for carbon entry into terrestrial ecosystems, photosynthesis plays the keystone role in the biosphere of transferring atmospheric CO 2 into terrestrial ecosystems. Since its inception 40 years ago the Farquhar et al. (1980;hereafter FvCB) model of C3 photosynthesis has revolutionized photosynthesis research (>5,000 citations, at time of writing). The FvCB model describes photosynthetic carbon assimilation (A) using a set of mathematically described hypotheses that represent the enzymatic subprocesses of photosynthesis and their integration, including: light-stimulated electron transport, CO 2 fixation in the Calvin-Benson cycle and photorespiration. The FvCB model is an integrated set of mathematically described hypotheses, a system hypothesis, that yields quantitative predictions to accurately describe the dynamics of A in response to incident radiation (I), carbon dioxide concentration (C a ) and temperature. Observation, experiment and model-based photosynthesis research has seen substantive advances due to the availability of this mathematically rigorous hypothesis. However, variants of the FvCB model exist, chief among them is that proposed by Collatz et al. (1991;hereafter CBGB). Differing hypotheses for three key subprocesses distinguish the models: (a) electron transport, (b) limiting process selection, and (c) triose phosphate use (TPU).
Terrestrial biosphere models (TBMs)-which simulate global land ecosystems and their role in the Earth System-rely on the FvCB and CBGB models (as well as various hybrids and additions to these core models) to simulate leaf-scale photosynthesis and its response to global change, in particular, increasing C a . The CBGB model and hybrids with the FvCB model are employed by several prominent TBMs (Table 1; e.g. IBIS, JULES, CLM and ELM; Clark et al., 2011;Foley et al., 1996;Oleson et al., 2013). Yet despite the keystone role of photosynthesis in the biosphere and the wide variation in TBM simulations of photosynthesis (Anav et al., 2015;Rogers, Medlyn, et al., 2017), a rigorous, quantitative comparison of the FvCB and CBGB models has not been undertaken. In part this is because rigorous methods to compare and evaluate competing sets of hypotheses have not been available until recently.
Sensitivity analysis (SA) is used to determine the sensitivity of model output to variability in individual model components.
Variability in model output can be introduced through a number of K E Y W O R D S carbon assimilation, high-resolution A-C i curve, multi-hypothesis modelling, photosynthesis, process sensitivity analysis, terrestrial biosphere model

Model
Carboxylation TPU

M1111
1. Equations (1), (4) and (5) 1. Absent 1. Equation While all of these models assume TPU limitation, they all assume α tpu = 0 so that TPU limited A = 3TPU. c IBIS uses a unique TPU formulation (Foley et al., 1996). TA B L E 1 Factorial list of possible models from the given hypotheses National Science Foundation, Grant/Award Number: 1852977 and 1552329 sources (Beven, 2016), two key sources are parameter choice and differences in mathematical representations of the multiple processes that a model simulates, e.g. photosynthetic electron transport. SA methods to assess model sensitivity to various parameters in TMBs are common (e.g. Dietze et al., 2014;Gupta & Razavi, 2018;Koven et al., 2019;Ricciuto et al., 2018;Zaehle et al., 2005) while SA methods to assess model output sensitivity to alternative process representations are rare. Parameter SA methods can be applied in the context of multiple models and sensitivity indexes averaged to get an overall assessment of parameter influence under model uncertainty (Dai & Ye, 2015). However, these methods miss a key element of model output variability-the difference in the means among models, or between-model variation. Parameter SA can only account for within-model variation, necessitating a process SA that is designed to account for both within (parameter) and between (process representation) model variation (e.g. Dai et al., 2017).
A further obstacle to rigorous process SA is that the majority of commonly used modelling codes are not sufficiently flexible to switch between all of the various hypotheses for all of the various processes under investigation. The Multi-Assumption Architecture and Testbed (MAAT) has been designed to overcome this issue (Walker et al., 2018). MAAT is a hyper-modular, multi-hypothesis modelling framework designed to easily pose multiple alternative models by combining alternative mathematically described hypotheses for the processes that form the building blocks of a model, or system hypothesis (Walker et al., 2018). Through hyper-modularity MAAT allows a factorial combination of each hypothesis across all processes and subprocesses, exploring the full range of possible models and ensuring that each representation of each process is evaluated against all other representations for all other processes, i.e. it is fully comprehensive. MAAT's ability to combine models at the scale of individual process hypotheses enables the application of rigorous process (SA) algorithms, such as that of Dai et al. (2017) which was designed to work with modelling codes like MAAT.
In this study we use MAAT to formally compare the leaf-scale enzyme-kinetic models of C3 photosynthesis by comparing the FvCB and CBGB model formulations. The choice of electron transport model, limiting-rate selection, TPU limitation and parametric variability are examined. We ask the questions: (a) which processes are most influential for simulating carbon assimilation at various levels of atmospheric CO 2 concentration and incident radiation, (b) which parameters are most influential, and (c) are the influential process and parameters different when considering absolute assimilation or the response of assimilation to a change in CO 2 ? We further evaluate the outcome of this SA using global TBM simulations and measurements of leaf-scale photosynthesis.

| Comparison of the FvCB and CBGB models of photosynthesis
Enzyme-kinetic models of photosynthesis (Farquhar et al., 1980) simulate net CO 2 assimilation (A-µmol CO 2 m −2 s −1 ) in response to CO 2 concentrations in the intercellular airspace of the leaf (C i -Pa) and incident photosynthetically active radiation (I-µmol photons m −2 s −1 ). The model scales the gross carbon assimilation rate (A g -µmol CO 2 m −2 s −1 ) to account for photorespiration, minus 'dark' respiration (R d -µmol CO 2 m −2 s −1 ): where Γ * is the photorespiratory CO 2 compensation point (Pa), the C i at which A g is equal to the rate of CO 2 release from oxygenation. A g is simulated as a change point model (Gu et al., 2010) where one of two (FvCB) or three (CBGB) potentially limiting processes (A c,g , A j,g , and A p,g -µmol CO 2 m −2 s −1 ), described in detail below, are selected. FvCB simply identifies the minimum of the potentially limiting rates: resulting in discontinuities in the derivative of the A-C i or A-I curves at the change points where A c,g = A j,g and A j,g = A p,g . In order to 'introduce a more realistic, gradual transition from one limitation to another, and to allow for some co-limitation', CBGB proposed non-rectangular hyperbolic (quadratic) smoothing between the three potentially limiting rates: where A cj,g is a latent variable resulting from smoothing between A c,g and A j,g . Parameters θ cj and θ cjp (θ and β in CBGB's original notation) are curvature parameters that take a value 0-1 with lower values leading to greater smoothing. The FvCB method is a special case of the CBGB method where both θ cj and θ cjp take the value 1, while if θ cj and θ cjp take the value 0 smoothing becomes rectangular hyperbolic (Johnson & Thornley, 1984).
A c,g , A j,g and A j,g are modelled as Michaelis-Menten functions of C i . For A c,g , V cmax (µmol CO 2 m −2 s −1 ) determines the asymptote: where O i is the internal O 2 partial pressure (kPa), and K c (Pa) and K o (kPa) are the Michaelis-Menten half-saturation constants of the RuBisCO enzyme for CO 2 and for O 2 . For A j,g , the asymptote is proportional to the electron transport rate (J-µmol e m −2 s −1 ) where: A number of hypotheses to represent J exist, most commonly used are the following three. Based on Smith (1937), two representations of J saturate at a maximum electron transport rate (J max ), (a) Farquhar and Wong (1984) used non-rectangular smoothing (Equation 6a), (b) Harley et al. (1992) use an alternative non-rectangular hyperbola ), while (c) CBGB proposed a linear model that has no maximum (Equation 6c) respectively: where a is the leaf absorptance (the fraction of I absorbed by the leaf, unitless), α i is the intrinsic quantum efficiency of electron transport (the number of electrons transported per absorbed photon, e photon −1 ) and θ j is a non-rectangular hyperbola smoothing parameter. aα i is the apparent quantum efficiency of electron transport (electrons transported per incident photon). Following FvCB, in this study we define α i as 0.5(1 − f), where f is the fraction of photons absorbed by the leaf but not absorbed by the light harvesting complexes, and 0.5 represents the requirement of two photons for full linear transport of a single electron.
Subsequent to the development of the FvCB model, a third potential limitation was identified under high C a and high irradiance (I). TPU in sucrose and starch synthesis releases phosphate needed for the regeneration of RuBP, thus low rates of sucrose and starch synthesis can limit RuBP regeneration and therefore A (Sharkey, 1985). CBGB proposed this additional rate-limiting cycle in their model (A p,g ), which was refined (von Caemmerer, 2000) to account for reversed sensitivities of A to C i and O i in the TPU limiting state (Harley & Sharkey, 1991): where α T is the fraction of glycolate exported but not returned to the chloroplast during photorespiration. CBGB assumed a closed photorespiratory cycle (α T = 0) such that multiplication by the first term in Equation (1) yields 3TPU.

| ME THODS
As described above, the MAAT (Walker et al., 2018; https://github. com/walke ranth onyp/MAAT) is a hyper-modular, multi-hypothesis modelling framework. MAAT is written in R (R Core Team, 2020) and provides a general framework and code structure for building models that allows for new processes to be added easily. Higher-level 'system models' integrate multiple processes into a coherent representation of a given system. A number of these system models have been coded into MAAT and in this study we use the leaf-scale enzyme-kinetic model of photosynthesis (Walker et al., 2018). MAAT also encodes process and parameter SA algorithms. In this study we used MAAT (tag: v1.2.1_walkeretal2020_GCB, commit hash: 09b1479). Note that in the original studies, FvCB and CBGB describe the enzyme-kinetic models using different units and slightly different parameter definitions. In this study we use the unification of (6a) 0 = j J 2 + a i IJ max J + a i IJ max ,  Walker et al. (2018) that are built upon Gu et al. (2010) and predominantly follow FvCB. Building on the model details presented in the introduction, Γ * can be simulated as a function of K c and K o :

| Additional model details
where k o and k c are the turnover numbers for the oxygenase and carboxylase functions of RuBisCO. At 25°C k o is 0.21 times k c , and their activation energies are the same so this ratio is preserved across a range of temperatures (Farquhar et al., 1980). J max can be represented in a number of ways but many studies have demonstrated the tight correlation between J max and V cmax (Leuning et al., 1995;Walker, Beckerman, et al., 2014;Wullschleger, 1993) and for the basis of this study we use the commonly employed linear relationship based on the relationship proposed by (Wullschleger, 1993): where a jv and b jv are calculated from linear regression of J max,25 on V cmax,25 (which are J max and V cmax at a reference temperature of 25°C).
We use a similar relationship to calculate TPU: We focus on the core components of the FvCB and CBGB models and do not consider choices of CO 2 -diffusion resistance models (i.e. leaf boundary layer, stomata and internal; see Collatz et al., 1991;Walker et al., 2018) nor temperature response models, which would substantially increase scope. Further, while stomatal models do influence C i , the response to I and C a of many stomatal models used by TBMs is similar, i.e. g s = f(A/C a ). In order to preserve realism in the calculation of A we used the unified stomatal model (Medlyn et al., 2011) with a g 1 of 4.3, the global C3 mean (Lin et al., 2015), and a typical value of 0.01 mol m −2 s −1 for g 0 . We assume that leaf internal resistance and leaf boundary layer resistance are 0, vapour pressure deficit is 1 kPa, and no soil water limitation. Leaf temperature was set to a standard of 25°C (which assumes that all temperature sensitive parameters were at 25°C reference values).
To generate C a and I response curves, the two models were run from a C a of 50-1,500 µmol/mol in 50 µmol/mol increments at I of 960 μmol m −2 s −1 , and from I of 10-1,960 μmol m −2 s −1 in increments of 50 μmol m −2 s −1 at C a of 400 µmol/mol.

| Sensitivity analysis
We use statistical, variance-based process SA and parameter SA which both rely on an ensemble of model simulations to calculate model output variance and ascribe this variance to variation in processes and parameters. The photosynthesis model process SA was broken into four processes: limiting-rate selection (two representations, Equations 2 or 3), electron transport (three representations, Equations 6a-c), TPU limitation (included or not included, Equation 7) and carboxylation (one representation, Equation   4). Factorial combination of the alternative process representations gives a total of 12 individual models. For both the process SA and parameter SA, 14 parameters were varied ±10% from a central, commonly used value (Table 2). Thus the model ensemble comprises 12 individual models and variation of 14 parameters.
Both methods calculate sensitivity indexes that represent the proportion of variance in model-ensemble output (both between and within-model variance for the process SA, and just within-model variance for the parameter SA) caused by variance in either a specific process or a specific parameter. Process SA provides a quantitative assessment of the influence of processes on model output and includes both variation caused by alternative hypotheses for the mechanics of a given process, and variation caused by parameters within a given process, i.e. between-model and within-model variability (Dai et al., 2017). Parameter SA can assess the influence of In this analysis we calculate and focus on first order sensitivity indexes, analogous to main effects in ANOVA. We do not calculate two-way or higher order interactions among processes or parameters. While some interactions are likely to be interesting, the sum of first order sensitivity indexes sum to around 0.95 in many cases indicating that over 95% of variance in model output was explained by first order effects.
For both analyses, we investigated sensitivities of simulated carbon assimilation (A, Equation 1) across a number of environmental scenarios that were a factorial combination of atmospheric CO 2 (C a ; 280, 400 and 600 µmol/mol) and incident light (I; 200, 500, and 1,000 µmol photons m −2 s −1 ) conditions. The sensitivity of the absolute assimilation response (ΔA) to changes in C a (from preindustrial to present-day, 280-400 µmol/mol, and from present-day to projected future concentrations, 400-600 µmol/mol) were also calculated under the three incident light conditions. Variance-weighted means of the sensitivity indexes across different environment combinations allow us to quantify the general influence of a process or parameters across environment combinations. For parameters this can also be done across model combinations, but again still only account for within model variance.
Parameter samples were drawn from uniform distributions and Other values have since been used by TBMs; e.g. 0.9 in IBIS (Foley et al., 1996) and 0.83 in JULES (Clark et al., 2011). The possible values for these smoothing parameters are from 0 ≤ θ ≤ 1 and so to maintain values within this range and preserve the ±10% variation in all parameters, we use a central value of 0.9 for θ cj and θ cjp . These data are publically available (Walker, Lu, et al., 2020).
To assess convergence in sensitivity index calculations, preliminary SAs were run with an n of 1,000 for the process SA and 1,000,000 for the parameter SA. Subsampling and bootstrapping indicated that an n of 300 for the process SA and 300,000 for the parameter SA were ample to achieve convergence (standard deviations of the sensitivity indexes were less than 0.001). The results of the SA shown in this study were generated using these smaller values of n, resulting in a total of 4,320,000 executions for the process SA and 50,400,000 for the parameter SA. These total number of executions are larger than n as n is the base number of samples and the full SAs require multiple sets of iterations that are a function of the number of model combinations and parameters investigated (see Walker et al., 2018).

| Estimation of θ cj from high-resolution A-C i curves
High-resolution A-C i curves (Anderson et al., 2020)  Preliminary measurements identified saturating I and C a where photosynthesis transitioned from RuBP saturated (A c,g ) to RuBP limited (A j,g ) photosynthesis. These preliminary measurements informed a commonly used A-C i response protocol , developed to include a high density of measurements around the transition point (when A c,g = A j,g ). Leaves were first acclimated to chamber conditions (I = 2,000 μmol m −2 s −1 , C a = 400 μmol/mol, flow rate = 600 μmol/s, relative humidity = 70%-75%, leaf temperature = 30°C) and measurements began once steady-state gas exchange was achieved. C a was taken from 400 to 50 µmol/mol then returned to a conservative estimate of the start of the transition zone (305 µmol/mol) and raised progressively in 5 µmol/mol increments to 1,000 µmol/mol (a value comfortably higher than the end of the transition zone). C a was then raised in larger increments to capture the full extent of a standard A-C i curve.

TA B L E 2 Comparison of the parameters used in the original papers by Farquhar and Collatz unified to common units
Abbreviation: TPU, triose phosphate use. a Calculated from Collatz CO 2 :O 2 specificity ratio, τ in their notation, of 2,600 where k o : b Parameters were not originally specified in Farquhar but values in parentheses featured in Farquhar and Wong (1984).
Bayesian machine-learning, Markov chain Monte Carlo (MCMC), algorithms were used to numerically approximate the posterior distribution of the smoothing parameter, θ cj (Equation 3a), and V cmax and J max , from these high-resolution A-C i curves. Numerical approximation is achieved by randomly sampling from a specified prior distribution, stochastically generating a proposal for the parameters to these were then thinned to 1% to remove auto-correlation.

| TBM simulations
We use three TBMs to test the impact of quadratic smoothing on Decadal average annual GPP from the two simulations were then compared to determine the impact of the smoothing parameters.
Model results were re-gridded to a common 0.5° × 0.5° spatial grid, using bilinear interpolation where necessary (ELM). A 0.5° land mask was then applied to constrain model output to a common areal extent on which to base annual calculations and maps of global GPP (Walker, Fisher, et al., 2020).

| RE SULTS
In their original parameterizations, the C a response of CBGB is smoother and more sensitive than FvCB at both low and high C a ( Figure 2a). With unified parameters the models are similar at low to intermediate C a but at high C a the CBGB model is again more sensitive to C a and the difference in A approaches 10 μmol m −2 s −1 at 1,500 µmol/mol (Figure 2b). Comparison of A implied by each of the two or three potentially limiting rates (i.e. calculating A from A c,g , A j,g , and A p,g in Equation 1) explains these responses (Figure 2c,d).    Across all environmental conditions, limiting-rate selection was responsible for 57% of the variation in A, carboxylation responsible for 22%, electron transport 10% and TPU 2% (Table 3). The strong influence of limiting-rate selection was borne out across the majority of environmental conditions (Figure 3b). However, at saturating I (1,000 μmol m −2 s −1 ) the process of carboxylation was most influential at preindustrial C a (280 µmol/mol), and at present-day C a (400 µmol/mol) the influence of carboxylation was about equal to limiting-rate selection. This pattern was similar at close to saturating I (500 μmol m −2 s −1 ), but limiting-rate selection had a generally higher influence. At low values of I (200 μmol m −2 s −1 ), the process of electron transport was more influential than carboxylation with the sensitivity of A to electron transport increasing as C a becomes less limiting. Bimodality was most apparent when limiting-rate selection was most influential (the alternative modes corresponding with the two alternative representations of limiting-rate selection).
When all 12 models and all nine environmental scenarios were combined, variation in A of more than 5% was caused by only five of the 14 parameters ( Figure 3c; Table 4). V cmax was the most influential parameter responsible for 35% of the variation in A, followed by K c with 22%, θ cj with 19%, K o with 7%, and a with 7%.
The influence of these parameters somewhat reflects the influence of the processes to which they belong. However, if we were to only consider variability in A caused by parameter variation, the influence of carboxylation would be over-estimated.
V cmax , K o and K c are all parameters in the process of carboxylation and together they were responsible for 64% of the variation in A in the parameter SA (total variance: 0.94). Together θ cj and θ cjp , the parameters of limiting-rate selection, were responsible for only 23% of variance in the parameter SA. On the other hand, the process SA suggested that carboxylation was responsible for a more modest 22% of the variation in A (total variance: 2.59), while limiting-rate selection was responsible for 57%. The results F I G U R E 3 Sensitivity of carbon assimilation (A, μmol m −2 s −1 ) to variability in processes and parameters across various C a and I environmental conditions. (a) Semi violin plots showing distributions of A against C a (μmol/mol, x-axis) and I (μmol m −2 s −1 , panels) boxes represent the interquartile range and median, whiskers the full range. (b) First order sensitivity index of A to variability in the four processes against C a (μmol/mol, x-axis) and I (μmol m −2 s −1 , panels). (c) First order sensitivity index of A to variability in the 14 parameters, indexes integrated across the 12 models and nine environmental conditions. (d) First order sensitivity index of A to variability in the 14 parameters, indexes integrated across the 12 models and colour coded for each of the nine environmental conditions. (e) First order sensitivity index of A to variability in the 14 parameters, indexes integrated across the nine environmental conditions and colour coded for each of the 12 models a a jv Sensitivities for individual models (integrated across environmental scenarios) showed that, given equal variation (±10%), V cmax and θ cj shared similar maximum sensitivities: 48% and 52% respectively ( Figure 3e; Table 4). θ cjp had a maximum sensitivity of 21% in the model with smoothing, TPU and no J max term in electron transport (M1223). Leaf light absorption, a, featured in all models of electron transport, and therefore all models, and sensitivity to a varied between 9% and 13% or 3%-6% depending on minimum or smoothing limiting-rate selection respectively.
Similarly, V cmax , K o and K c were all more influential in the models that used the minimum for limiting-rate selection. Sensitivities for individual environmental conditions (integrated across models) showed that as expected a was influential under low-light conditions while V cmax , K o and K c were influential under high-light conditions (Figure 3d). Sensitivities for individual model and environmental condition combinations showed that for some cases some of the previously unmentioned parameters were influential (e.g. θ j , f), while others remained with very little influence (e.g. Figure 4a shows the distribution of ΔA in response to changes in C a (from preindustrial to present-day, 280-400 µmol/mol, and from present-day to projected future concentrations, 400-600 µmol/mol) at three levels of I. Variation of ΔA was greatest at intermediate light levels and going from present to future C a (range about 2 to 5 μmol m −2 s −1 ). Variation was similar at high light. For both C a changes and at both high and intermediate light, the distribution of ΔA was highly bimodal, with stronger bimodality going from preset to future C a . As for A, the alternative modes were associated with the alternative representations of limiting-rate selection.

| SA of the assimilation response (ΔA) to changes in C a
When all environmental scenarios were combined, limiting-rate selection was responsible for 65% of the variation in ΔA (total variance: 0.50), carboxylation responsible for 5%, electron transport 13% and TPU 3% (Table 3). The strong influence of limiting-rate selection was borne out across the majority of environmental conditions ( Figure 4b). For both C a changes and at high and intermediate I, sensitivity of ΔA to limiting-rate selection ranged from 64% to 76%, with the higher sensitivities at intermediate I. At low I electron transport was the most influential process, accounting for 50% of the variation in ΔA at the lower C a change and 32% at the higher C a change. At low I sensitivity of ΔA to limiting-rate selection increased from 7% at the lower C a to 37% at the higher C a .
In contrast with the sensitivity of A, seven parameters were responsible for over 5% variation in ΔA (total variance-0.09) when models and environmental scenarios were combined. Of these seven

TA B L E 4 First order sensitivity to parameters
parameters, four were in common with A (V cmax with 26%, K c with 9%, θ cj with 20% and a with 9%; Figure 4c; Table 4), K o did not feature, and θ cjp with 11%, b rv with 7% and θ j with 5% were also influential.

| Consequences of limiting-rate selection assumptions
Given the sensitivity of A and ΔA to the processes of limitingrate selection, we now investigate these models in more detail.

Mathematical analysis shows that FvCB sets the upper limit for
A while CBGB smoothing always reduces A below that of the minimum (see Supporting Information). The greatest reduction in A caused by smoothing is when all three limiting rates-A c,g , A j,g and A p,g -are equal, and yields Equation S3 (see Supporting Information). With θ cj = 0.95, θ cjp = 0.98, Equation S3 shows that the smoothing scalar on A g is 0.77, i.e. when A c,g = A j,g = A p,g quadratic smoothing reduces A g by 23%. When only A c,g and A j,g are equal and substantially lower than A p,g (so A p,g effectively has no influence on smoothing), A g is reduced by 18%. Figure 5a shows that A g is reduced below the minimum rate across a wide range of A c,g and A j,g values and that the reduction in A g approaches 0 monotonically as the difference between the minimum rate and the larger rate increases.  year, or 4.4%-9.7%).

F I G U R E 5
Relative reduction in calculated A g of GPP (%) using quadratic smoothing compared to the minimum of the limiting rates. (a) Relative reduction in A g (%) against the relative difference in A c,g and A j,g (% increase relative to the minimum of the two rates) when triose phosphate use is not simulated or simulated, colours represent relative difference in A p,g and the minimum of A c,g and A j,g with orange representing the lowest difference and therefore the largest reduction in A g . (b) The relative reduction in A g (%) as a function of both A c,g and A j,g when A p,g is 20 μmol m −2 s −1 (see colour scale for d). (c) Global GPP (g C/m 2 ) simulated by the three terrestrial biosphere models (TBMs): ELM, FATES and SDGVM. (d) Relative reduction in GPP (%) caused by non-rectangular hyperbolic smoothing in the three TBMs. (e) High-resolution A-C a curves used to estimate θ cj . (f) MCMC posterior distributions for θ cj estimated from the curves in (e). The vertical grey line in (e) represents the C a cutoff for the single high-resolution curve that showed a drop in A at high C a . Where values of GPP in (c) were less than 250 g C/m 2 , values in (d) were screened to avoid over-emphasizing high relative changes on small absolute rates of GPP that do not contribute substantially to the global carbon cycle

| Discriminating among hypotheses for limitingrate selection
High-resolution A-C i curves were taken on three individuals of Populus canadensis (Figure 5e) in order to estimate the θ cj parameter,  Figure 5f). Thus the only data-driven estimate of θ cj made to-date is not significantly different from 1.0, providing support for the FvCB method of limiting-rate selection.

| D ISCUSS I ON
A novel, mathematically rigorous, process SA that accounts for both hypothesis (process representation) and parameter variability in common models of photosynthesis has shown that limiting-rate selection was the most influential process, accounting for 57% of variation in A and 65% of variation in ΔA in response to a change in CO 2 . For simulating A, carboxylation was the next most influential process (which was all due to parameter variability) followed by electron transport.
When simulating ΔA, electron transport was the next most influential

| The influence of limiting-rate selection
Finding the smaller roots of the quadratics described by Equations (3a,b) is intended to smooth the abrupt transition between A c,g , A j,g and A p,g , and has been described as representing co-limitation between limiting rates. In so doing, smoothing also imposes a reduction of modelled A (Figure 5) ). This scenario, when all potentially limiting rates are equal, results in the greatest smoothing-related reduction in A g and A, but it is not an extreme physiological scenario. Potential assimilation rates A c,g and A j,g are often observed to be close to co-limiting in saturating light (e.g. Ainsworth et al., 2003;Bernacchi et al., 2005) and the co-ordination hypothesis assumes that A c,g and A j,g are equal under mean environmental conditions (Maire et al., 2012). Co-ordination hypotheses are used in a number of optimization schemes (Smith et al., 2019;Wang et al., 2017), which maintain photosynthesis close to the transition given a changing environment. What process is smoothing intended to represent? CBGB introduce smoothing to represent a more realistic transition in the light response of A, and this could represent imperfect coupling among the cycles of electron transport, carboxylation and photorespiration (Farquhar et al., 1980). It is also possible that smoothing is accounting for different limiting states of individual chloroplasts within a leaf (Buckley et al., 2017;Kull & Kruijt, 1998

| Influence of other processes
Electron transport, the process most commonly thought of as the key difference among FvCB and CBGB, was not a strongly influential process (10% for A and 13% for ΔA when integrated across environmental scenarios). The influence of electron transport as a process was not greater than the sum of the influence of its parameters: a, a jv , b jv , f and θ j (sum of sensitivities indexes 11% for A and 16% for ΔA). That is to say, the alternative representations of electron transport did not result in appreciable between-model variability under the environmental and other conditions of the SA.
This result suggests that under the conditions of this SA, a linear electron transport rate or a saturating rate with J max simulated as a linear function of V cmax had very little effect on simulated assimilation rates. This inference is supported by the small sensitivity indexes of a jv and b jv (0% and 0% integrated across models and scenarios, for A).
Our setup, based on a commonly used relationship in TBMs (Wullschleger, 1993)   .
Nevertheless, the central values of a jv and b jv that we used are commonly employed by TBMs and a V cmax of 50 μmol m −2 s −1 is a fairly representative value, so the parameter space of our SA is likely representative of a substantial proportion of parameter space across multiple TBMs.
Electron transport in the FvCB models also uses empirical smoothing between J max and an electron transport rate that is linear with I. In models that employ non-rectangular smoothing in electron transport the parameter θ j has a sensitivity index of 4%-8% when integrating across environmental scenarios. Smoothing of J in response to light has less influence than smoothing for limiting-rate selection because it only smooths the light response, thus can only influence the light-limited rate, while CBGB smoothing is applied to every single calculation of A. Furthermore the empirical smoothing used to represent the response to I is the best representation of that relationship currently available and lacks the mechanistic understanding and representation possible for A c,g using Michaelis- Menten kinetic theory (von Caemmerer, 2000).
With the addition of TPU as a limiting rate, smoothing further decreases A g compared with using the minimum (Figure 5a). The SA suggested that the inclusion of TPU limitation as a process was not strongly influential (2% for A and 3% for ΔA when integrated across environmental scenarios). Due to the co-ordination of photosynthetic apparatus, the rate of TPU export is usually simulated as a proportion of V cmax (Equation 9). Collatz et al. (1991) used a value of b tv equivalent to 0.167, which is commonly used and is the central value used in this SA. Lombardozzi et al. (2018) pointed out that this value of b jv may be too high based on Wullschleger (1993) and demonstrated a 9 Pg C (about 9%) smaller increase in global terrestrial ecosystem C between 1850 and 2100 in CLM4.5 under RCP8.5 when a value of 0.083 was used for b jv . However, 0.083 is >1 standard error lower than the Wullschleger (1993) mean and lower than the 95% CI at 25°C from a recent synthesis (Kumarathunge, Medlyn, Drake, Rogers, & Tjoelker, 2019). Ellsworth et al. (2015) showed that TPU can be limiting under high-light and high-CO 2 conditions, concluding with a general recommendation that modellers interested in simulating A should consider TPU limitation (as formulated by Equation 7). However, their results demonstrate that TPU limitation is primarily influential under conditions of low O 2 (2%) or saturating C i (>100 Pa). Given that these conditions are rather extreme, the low ratio of TPU to V cmax chosen by Lombardozzi et al. (2018), and the results of our SA, we suggest that calculating TPU at the photosynthetic core of TBMs is probably an unnecessary computational cost.
Despite the relative lack of influence of the process of carboxylation, V cmax was still the most influential parameter (i.e. accounted for more of the within-model variance than any other parameter) when models and environment were combined. This discrepancy highlights the importance of considering the variability in model process representation when conducting model SA, as illustrated by the different variances calculated by the two SA types (for A, parameter SA variance = 0.94, while process SA variance = 2.59, despite equal means).

| CON CLUS IONS
At the heart of TBMs lies the surprising dominance of the nonmechanistic, limiting-rate selection process. While empirical smoothing among limiting photosynthetic rates may account for a number of mechanistic processes at various scales, it is unsatisfying that empirical functions have such influence in a model that is intended to be highly mechanistic. Indeed the FvCB model is at the core of many TBMs specifically because of its mechanistic simulation of the primary response of the terrestrial biosphere to rising CO 2 concentration, a principal driver of global change. In this SA, limiting-rate selection accounts for 65% of the variance in the CO 2 response of A.
That this empirically driven variation lies within what is assumed to be a highly mechanistic process representation at the core of TBMs, it is perhaps not surprising that there is such a vast range of disagreement in Earth System model projections of the future terrestrial carbon sink. While FvCB limiting-rate selection represents selection of A at its upper bound, we suggest that this is a more defensible assumption than a highly influential non-mechanistic function with an essentially arbitrary choice of parameter values that are not supported by data.
To increase confidence in our understanding and future projections of the carbon cycle, and thus climate, we need to understand how the process representations and parameters used by TBMs drive variation in TBM simulations . Previous methods to evaluate process representations have relied on model inter-comparison projects (MIPs) either with multiple models (e.g. Anav et al., 2015;Arora et al., 2019;Walker, Hanson, et al., 2014) or comparison of alternative representations of submodels or processes within a single higher-level system model (e.g. Burrows et al., 2020). Both of these methods sample only an extremely small fraction of possible model combinations (Abramowitz & Bishop, 2014;Fisher & Koven, 2020). By allowing a fully factorial combination of models, MAAT and the process-level SA (that includes hypothesis and parameter variability) in this study represents a new frontier for model analysis and development. We investigated a small but influential component of TBMs, finding a surprising leaf-scale sensitivity that has global-scale implications. Yet the analysis presented here is just the beginning of what is possible. The quantitative multi-hypothesis tools provided by MAAT, and by other multi-hypothesis modelling groups, will help to provide rigorous advances in process-level understanding of the dynamics of complex ecosystems.