Estimating canopy gross primary production by combining phloem stable isotopes with canopy and mesophyll conductances.

Gross primary production (GPP) is a key component of the forest carbon cycle. However, our knowledge of GPP at the stand scale remains uncertain because estimates derived from eddy covariance (EC) rely on semi-empirical modeling and the assumptions of the EC technique are sometimes not fully met. We propose using the sap flux/isotope method as an alternative way to estimate canopy GPP, termed GPPiso/SF , at the stand scale and at daily resolution. It is based on canopy conductance inferred from sap flux and intrinsic water-use efficiency estimated from the stable carbon isotope composition of phloem contents. The GPPiso/SF estimate was further corrected for seasonal variations in photosynthetic capacity and mesophyll conductance. We compared our estimate of GPPiso/SF to the GPP derived from PRELES, a model parameterised with EC data. The comparisons were performed in a highly instrumented, boreal Scots pine forest in northern Sweden, including a nitrogen fertilised and a reference plot. The resulting annual and daily GPPiso/SF estimates agreed well with PRELES, in the fertilised plot and the reference plot. We discuss the GPPiso/SF method as an alternative which can be widely applied without terrain restrictions, where the assumptions of EC are not met. This article is protected by copyright. All rights reserved.

3 Some studies using  13 C to estimate WUEi ( 2006). Because gm is finite, assuming that it is infinite leads to an overestimation of WUEi (Seibt et al.,111 2008; Wingate et al., 2007). Considering gm associated with δ 13 Cp measurements would considerably 112 improve GPP estimates, especially for conifers, which have relatively low gm (Rascher,Máguas,& 113 Werner, 2010). There is as yet no agreement about how to model gm, but it has often been estimated 114 from gS (Warren, 2008). 115 116 We should therefore give similar results. The GPPiso/SF method can also provide information on how GPPiso/SF 126 responds to fertilisation in terms of assimilation and gS. 127 128 A boreal forest is particularly suited for such a method comparison because of its simple species 129 composition (Hänninen, 2016;Högberg, 2007). Moreover, because this biome is strongly nitrogen (N)-130 limited (Du et al., 2020), N additions induce a strong response in terms of growth and C fluxes (see 131 reviews and references therein in Högberg, 2007 andTamm, 1991). These increases should be captured 132 by all methods. However, a positive N-fertilisation effect on GPP was not always observed. At our site, 133 previous studies showed no effect of N supply on GPP when measured from biometrics (Lim et al.,134 2015) or shoot-scale gas exchange (Tarvainen, Räntfors, Näsholm, & Wallin, 2016), but Tian et al. 135 (2020), who used eddy covariance data to parametrise a model, did find higher GPP in the fertilised plot 136 than in the reference plot. Thus the GPP results have been mixed, depending on which method was used. 137 The method we propose in this paper aims to provide an alternative stand-scale estimate of GPP that is 138 independent of eddy covariance. Our first objective here was to compare estimates of GPP based on 139 stable isotopes and sap flux against GPP based on PRELES, a process-based model parameterised with 140 eddy covariance data. The second objective was to determine how fertilisation treatment influenced the 141 canopy GPP with the sap flux/isotope method. Finally, the third objective explores alternative methods 142 for incorporating an empirical gm estimate and how these alternatives influence the GPP estimate. The study was carried out in a mature ~90 year-old Scots pine forest (Pinus sylvestris L.) at Rosinedal,148 near Vindeln in northern Sweden (64°10' N, 19°45' E) in 2012 and 2013. The site was an even-aged 149 and mono-specific stand, located on sandy soil. Two 15-ha plots were studied; a fertilised plot (F) and 150 a reference plot (R). In both plots, the sparse understory was dominated by Ericaceous shrubs, esp. was applied annually in mid-June to the fertilised plot (F) at a rate of 10 g N m -2 yr -1 , but reduced to 5 g 155 N m -2 yr -1 in 2012 and thereafter, using Skog-Can fertiliser (Yara, Sweden), containing NH4 (13.5%), 156 NO3 (13.5%), Ca (5%), Mg (2.4%), and B (0.2%) (Lim et al. 2015). 157 4 158

| Environmental data 159
Environmental data included half-hourly relative humidity (RH, %), photosynthetic photon flux density 160 (PPFD, µmol m -2 s -1 ), ambient temperature (Ta, °C) and soil water content (SWC, m 3 m -3 ), and 161 precipitation (mm) ( Figure S1). PPFD was measured at the R plot only and precipitation came from 162 Svartberget station, which is located about 8 km from the study site . During the period 1981-2010, mean  163 annual temperature and precipitation at Svartberget was 1.8 °C and 614 mm, respectively (Laudon et 164 al., 2013). Gaps in the meteorological data, due to instrument failure, were filled using measurements 165 from the Svartberget forest. All abbreviations, their units, and values of constants are summarised in 166 Table 1.  167  168 The temperature data were used to define the "thermal growing season" which estimates the period 169 theoretically suitable for vegetation growth for a given year (Cornes, van der Schrier, & Squintu, 2019; 170 Linderholm, 2006 2015), but we reduced the summer filter threshold to the minimum that would allow us to keep as many 283 data as possible. We filled the resulting GPP gaps using a predictive model (gC = a × Â +b) with a and 284 b determined for each combination of treatments. We replaced the GPPiso/SF outliers (beyond the 95% 285 confidence interval of the predicted values) and filtered values by the predicted functions only during 286 the thermal growing season. We did this because the common gapfill functions are based on EC data 287 and we wished to maintain our independence from EC data. The gaps were much larger outside the 288 thermal growing season than within it; because tree photosynthesis is reduced during that time we chose 289 not to fill these gaps. 290 291 2.5 | Carbon, discrimination, intrinsic water use efficiency and GPP calculations 292 Using the phloem samples collected between October 2011 and September 2012, we estimated isotopic 293 discrimination against 13 C (Δ, ‰), assuming it was mainly constituted from photosynthetic 294 carbohydrates. It was calculated as follows: 295 whole year at stand level, we argue that this short term variability would average out over the growing 303 season. The literature also describes variation in Δ between leaves and phloem contents and among 304 compounds in the phloem; we address this variation in the Discussion. 305 7 306 The intrinsic water use efficiency for the stand (WUEi) was then inferred from the following equation, 307 in each plot: 308 where Ca is the atmospheric CO2 concentration (µmol mol -1 ), r the ratio of diffusivities of water vapour 310 relative to CO2 in air (1.6), b the fractionation during carboxylation ( Gross primary production (g C m -2 ground area d -1 ) was then calculated from Eqn. 10 and Eqn. 12: 340 Eqn. 14 341 with MC the molar mass of C (12 g mol -1 ). The definitions of Wohlfahrt and Gu (2015) distinguish 342 between "canopy net photosynthesis," which includes carboxylation, respiration, and photorespiration, 343 "canopy apparent photosynthesis," which includes only carboxylation and photorespiration, and "true 344 photosynthesis," which includes carboxylation only. They point out that the flux-partitioning algorithms 345 used to calculate "GPP" with eddy-flux data are intended to estimate apparent photosynthesis (Wohlfahrt 346 & Gu, 2015). The sap flux/isotopic estimate also provides an estimate of canopy apparent 347 photosynthesis, at least in theory, because respiration is not allowed to influence the photosynthate pools 348 loaded into the phloem. However, Wohlfahrt and Gu (2015) go on to note that the flux-partitioning used 349 with eddy-flux data is inexact because it neglects the reduction in leaf respiration in the light. It is beyond 350 the scope of the current manuscript to solve that problem, but we hope that it can be addressed in the 351 near future. There was no replicate of the R and F treatments so it was impossible to perform analyses of variance 386 to infer any fertilisation effect. However, we could not ignore the effect of the fertilisation on the F plot 387 (Lim et al., 2015). We therefore presented the plot differences recognising that they may include a pre-388 existing plot effect as well as a fertiliser effect. 389 390 However, because 15 trees were sampled at each site for δ 13 Cp estimate, we did analyse a 'plot effect'. 391 We Errors in those inputs were assumed to follow normal distributions or truncated normal distributions 414 (see Table S1). The 95% confidence intervals were calculated to illustrate the predictive uncertainty in 415 our GPPiso/SF estimate ( Figure S3). The Sobol indices (Saltelli et al., 2008) were also calculated to 416 partition the variance into these uncertainty sources (Table S1) All analyses were conducted with R software, version 3.5.1 (R Core Team, 2016). 430 During this period the ratio Â/Â max was close to 1 ( Figure S4) meaning that photosynthetic capacity had 449 reached its seasonal maximum (Mäkelä et al., 2008). Conductance fell through September and October, 450 returning to zero in both plots on the 25 th of October in 2012 and the 4 th of December in 2013 ( Figure  451 1). 452

| Isotopic data 453
Isotopic data from the atmosphere and from the phloem were also used to infer WUEi. We observed 454 strong, but different, patterns of seasonal variation for atmospheric δ 13 Ca and for phloem contents 455 (δ 13 Cp). From January to the beginning of February, δ 13 Ca decreased to a minimum of -9.2‰ (Figure 2-456 10 A). Then δ 13 Ca increased rapidly, by about 1‰, during the initial weeks of high photosynthesis in late 457 June and early July. The main peaks of δ 13 Ca occurred during the thermal growing season, when canopy 458 conductance was also the highest. It then stabilised until late September, when it again began to fall 459 (Figure 2-A). In contrast, the phloem data ( Figure. 2-B) did not simply track the atmosphere. Instead 460 they showed a steep drop only at the beginning of the thermal growing season. The δ 13 Cp value depended 461 significantly on the date (p-value < 0.0001, df = 1, F = 53.09, Figure 2-B). It was significantly higher in 462 the F plot (-27.5‰) than in the R plot (-28.0‰) (p-value < 0.0001, df = 1, F = 76.96) as well. 463

| Comparison of GPP estimates 474
Our first objective was to compare GPPiso/SF to GPPPRELES for 2012 and 2013. To simplify the figure, we 475 chose to represent only the assumption that gm/g CÂ = 2.67 (Figure 4), which allows gm to vary during 476 the season. The seasonal GPP patterns were similar between PRELES and the sap flux/isotopic method 477 (Figure 4). Recall that GPPPRELES included understorey vegetation. Correlation coefficients among 478 methods and plots were all high, with minimum r = 0.91 ( Figure S5). However, the fit was nonlinear; in 479 2012 and 2013 GPPiso/SF approached an asymptote at high levels of GPPPRELES ( Figure S5). The highest 480 GPPPRELES values did not match with the highest GPPiso/SF values; the peak of GPPiso/SF occurred earlier 481 in the season than those of GPPPRELES. Interestingly, confidence intervals for GPPiso/SF and GPPPRELES 482 overlapped most of the time, even during the fall, when the offset was bigger than the rest of the year. 483 However, the VPD filters removed many values during the fall, which allowed us to draw a confidence 484 interval only during small periods at that time. As previously mentioned, the GPP values were gapfilled 485 to draw a complete seasonal pattern, at least during the thermal growing season. The resulting annual 486 sums were higher for GPPiso/SF than for PRELES on the control plot, but not on the fertilised plot ( Figure  487 5-A). 488 489 3.6 | Fertilisation effect 490 Our second objective was to assess the effect of fertilisation on GPP. Using the annual sums, neither 491 GPPiso/SF nor GPPPRELES was significantly different between the F and the R plots ( Figure 5-A). However, 492 there were consistent trends; GPPiso/SF was higher by 10% in the F plot than in the R plot and GPPPRELES 493 was higher by 16% ( Figure 5-A). Using the daily data corrected for autocorrelation, we found a 494 significant increase in the F plot; GPPiso/SF was higher by 8% and GPPPRELES was higher by 16% (Table  495 2 and see Figure S6). 496 497 3.7 | Mesophyll conductance assumptions 498 The third objective was to compare GPPiso/SF using different methods of estimating gm. Globally, there 499 was a significant effect of 'plot' (p-value = 0.007, df = 5, F = 19) and 'gm assumptions' (p-value = 500 0.0002, df = 5, F = 75). Focusing on one plot at a time, we found a significantly lower GPPiso/SF in the 501 control plot estimates when using gm/g CÂ = 2.67 as compared to the others. In the fertilised plot, we 502 found significantly lower GPPiso/SF of gm/g CÂ = 2.67 compared to gm∞. The F plot was not significantly 503 different from the R plot by any of these methods (Figure 5-B). 504 505 4 | DISCUSSION 506 Our study provided a new and simple method of independently estimating GPP and compared it to 507 estimates from PRELES, a model parameterised with EC data. The two methods yielded similar 508 estimates for both annual totals and seasonal patterns. We then used the two methods to compare a 509 11 fertilised to an unfertilised plot. Both methods detected higher GPP on the F plot, but only when using 510 the more abundant daily estimates ( based on earlier measurements of NPP at our site, was ~1000 g C m -2 y -1 (Lim et al., 2015). We compared 517 our GPPiso/SF estimate minus our standard deviation for the reference plot (1350 -43 = 1303 g C m -2 y -518 1 ) to the mean of these published values plus the standard deviation (1007 + 43 = 1050 g C m -2 y -1 ) and 519 found that the published values were consistently lower than our GPPiso/SF estimate. We next discuss the 520 strengths and weaknesses of each method. 521 522

| Strengths and weaknesses of PRELES 523
The key advantage of PRELES is that it is a compromise between predictive accuracy and model 524 complexity. It can be calibrated with a few variables derived from EC measurements. Once it is 525 calibrated, it can be run with an even smaller set of environmental variables (VPD, PPFD, precipitation 526 and air temperature). The required EC data are available from many sites around the world (Baldocchi Although the availability of EC data is an advantage for PRELES, EC data must be viewed with caution. 532 In particular, at our site, preliminary analyses of the data revealed significant problems in the data despite 533 the flat ground surface, uniform canopy, and low leaf area index. A careful study of the problem revealed 534 significant decoupling of the above-and below-canopy air masses, which often led to advection (Jocher  535  et  The key advantage of the sap flux/isotopic approach is that it is independent of eddy covariance. 549 Moreover, it leans on two methods, sap flux (Poyatos et al., 2007) and isotopic measurements (Bowling,550 Pataki, & Randerson, 2008 and references therein) that have been widely used at many sites by 551 ecosystem ecologists. The sap flux/isotopic approach combines them to estimate GPP at the tree scale, 552 which can then be scaled up to the stand. In simple stand structures, that scaling is relatively easy. We 553 used a model of sapflux based on measurements at our site scaled up in this way. It provided a simple 554 method to estimate tree GPP that, in combination with measurements of ground vegetation GPP, yields 555 an alternative estimate for comparison with GPP estimated by EC. 556 557 One critical advantage of the sap flux/isotopic method for estimating GPP is that its requirements for 558 the terrain and atmospheric conditions are less restrictive than for EC measurements. It thus provides an 559 empirical method that can be applied in hilly topography, complex canopy structure, and non-turbulent 560 atmospheres. 561

| Phloem contents and isotopic interpretation 562
The sap flux/isotopic method also has several important limitations. The literature describes several 563 post-photosynthetic modifications in the isotopic composition of the carbon that appears in the phloem 564 12 (Cernusak et al., 2009) related to transport and they are largely insoluble in water. In fact, their near absence would suggest 572 that phloem contents are less likely to show evidence of such modifications than bulk tissue. In this 573 sense, theory would suggest that phloem contents provide a better estimate of  of raw photosynthate 574 than does bulk leaf tissue. 575 576 Another set of post-photosynthetic modifications have been attributed to transport into and out of the 577 phloem during downward vertical transport. If these modifications reflected additions of photosynthate 578 from shaded branches, they might improve our estimates of whole-canopy photosynthesis. However, if 579 they were due to leakage and refilling with an isotopic fractionation, then they would degrade our 580 estimates (Gessler et al., 2009). In a detailed analysis of vertical changes in phloem composition in 581 Scots pine at our site, we were unable to detect a vertical δ 13 C gradient (data not shown). This argues 582 that the isotope signal is preserved during transport. 583 Post-photosynthetic modifications may also result from chloroplast starch hydrolysis and phloem 584 loading. Starch hydrolysis leads to diurnal changes in the isotopic composition of the sugars derived 585 from it (Gessler et al., 2009). In one study, the sugars leaving the leaf in the phloem had nearly the same 586 isotopic composition as the starch being hydrolyzed (Gessler et al., 2007). This result suggests that the 587 photosynthates were not substantially altered upon phloem loading. On the other hand, some authors 588 have found differences between δ 13 C of leaf soluble organic matter and the sugars in the phloem 589 (Brandes et al., 2006). This latter comparison assumes that the entire pool of leaf soluble organic matter 590 is available for export and that insoluble compounds, like starch, are not used as substrate for export. If 591 the assumption is true, it would suggest fractionation upon loading into the sieve tubes (Hobbie &  592 Werner, 2004). 593 594 Isotopic changes in phloem contents could also arise from compound-specific isotopic signatures in the 595 phloem. Such differences among compounds have been observed in phloem contents (Smith,Wild,596 Richter, Simonin, & Merchant, 2016) and they were especially noteworthy in the polyols in Douglas-fir 597 xylem sap (Bögelein, Lehmann, & Thomas, 2019), which represented 37% of the phloem solutes and 598 were approximately 2‰ more depleted than sucrose. It is not clear where the heavy carbon would go at 599 polyol synthesis, but one might expect that it is retained in the substrates. Similarly, phloem sap contains 600 N-compounds (e.g., amino acids and polyamines) as well. The δ 13 C analysis of phloem contents allowed 601 us to determine a C:N ratio, which was 119 ± 32 (SD) and 42 ± 20 in the R plot and the F plot, 602 respectively. On both plots, the values were high enough to consider that non sugar compounds would 603 have a small effect on the global isotopic signature. We acknowledge that a more detailed analysis would 604 improve our predictions. In the meantime, we have assumed that the bulk fractionation is negligible. 605 Phloem contents must be used carefully before photosynthesis begins in spring. During this period 606 before photosynthesis has begun, the phloem must contain, mobilized C reserves ( that the higher discrepancies between the GPPiso/SF and GPPPRELES in the fall and to a lesser extent in the 678 spring occurred because the constant ratio did not adequately account for seasonal dynamics in gm. The 679 need to refine our description of gm is confirmed by the uncertainty analysis (Table S1 and Figure S3) 680 The Sobol indices, which describe sources of uncertainty, showed that almost 75% of the GPPiso/SF 681 uncertainty came from the gm/g CÂ estimate. 682 683 4.4 | Difference between fertilisation treatments 684 We found a slightly higher GPP in the fertilised plot than in the reference plots with the sap flux/isotopic 685 method. Indeed, WUEi in the F plot was higher than in the R plot, although g CÂ was not different. The difference between the F and the R plots was only significant at the daily time scale, perhaps because 692 of the large number of repeated measurements (Table 2, Figure S5). However, this sap flux/isotopic 693 result, corrected for autocorrelation, was validated with the daily PRELES estimates (Table 2, Figure  694 S5). However, it should be recognized that these daily estimates are not independent and may exaggerate 695 our ability to detect a difference. In contrast, the annual sums did not detect a difference ( Figure 5), 696 perhaps because we were able to compare only two years, limiting the power of ANOVA. Thus, our 697 annual sums did not find a significantly higher GPP in the F plot compared to the R plot, agreeing with 698 previous studies focused on photosynthetic activity at shoot (Tarvainen et al., 2016) and stand scale 699 (Lim et al., 2015). The daily estimates did not agree. Based on these mixed results, we suggest that GPP 700 under the F treatment might be slightly higher, but that a replicated study would be necessary to settle 701 this question. 702 703 However, the magnitude of the GPP increase differed between PRELES and sap flux/isotopic methods. 704 The 8% increase in GPPiso/SF due to fertilisation was nearest to Lim et al (2015), who inferred a 3% 705 difference in GPP between the same F plot and the R plot based on biometric measurements. In contrast, 706 the GPPPRELES value in the F plot was 16% higher than in the R plot, almost twice the increase estimated 707 from GPPiso/SF and five times higher than in Lim et al. (2015). 708 709 4.5 | Role of understorey species 710 A key difference between the GPP methods is that GPPiso/SF quantified GPP of the trees only whereas 711 GPPPRELES quantified GPP of the whole ecosystem, which included understorey GPP. Understorey GPP 712 was 41 g C m -2 in a 120-year-old Scots pine boreal forest (Kulmala et al., 2011)  The GPPiso/SF method provides an alternative empirical method to estimate forest stand GPP that is 726 independent of eddy covariance (EC). We compared GPPiso/SF estimates from PRELES, a semi-empirical 727 model parameterised with EC data. When annual means were compared across two years, the GPP 728 15 estimates from the two methods were not significantly different. Moreover, the annual means showed 729 no effect of the fertiliser treatment. However, when compared using daily estimates, the fertilized plot 730 was 8% higher than the reference plot. The annual comparison agrees with previous estimates on this 731 site, the daily comparison does not. Future work will continue to explore this question (e.g., Tian et al.,732 in review). Adjusting GPPiso/SF for gm was necessary; we explored three alternatives for doing so. The 733 inclusion of mesophyll conductance provides an empirical/mechanistic means of connecting isotopic 734 measurements to gas-exchange measurements and GPPiso/SF provides a means of scaling from individual 735 trees to tree stands and canopies. Finally, a critical advantage of the sap flux/isotope based method for 736 estimating GPP is that its requirements for the terrain and atmospheric conditions are less restrictive 737 than for EC measurements. It can be applied in complex terrain, complex canopy structure, and non-738 turbulent atmospheres. 739  Table 2 Coefficients of the linear regression (corrected for autocorrelation) between daily GPP in the fertilised plot vs daily GPP in the reference plot across 1211 methods and years (slope (a) and intercept (b) ± SE).·,*, **, *** correspond to p < 0.1, 0.05, 0.01 and 0.001, respectively, after t-test when comparing the slope 1212 of the regressions with 1:1 regression.  , assuming a gm∞ assumption (green), a gm/g CÂ = 2.67 assumption (yellow) or a gm = 0.31 mol CO2 m -2 s -1 assumption (blue). Grey areas represent the thermal growing season Statistical results comparing WUEi between fertilised and reference plots: *** correspond to p < 0.001, respectively, after anova.