An emergent constraint on the thermal sensitivity of photosynthesis and greenness in the high latitude northern forests

Despite the general consensus that the warming over the high latitudes northern forests (HLNF) has led to enhanced photosynthetic activity and contributed to the greening trend, isolating the impact of temperature increase on photosynthesis and greenness has been difficult due to the concurring influence of the CO2 fertilization effect. Here, using an ensemble of simulations from biogeochemical models that have contributed to the Trends in Net Land Atmosphere Carbon Exchange project (TRENDY), we identify an emergent relationship between the simulation of the climate-driven temporal changes in both gross primary productivity (GPP) and greenness (Leaf Area Index, LAI) and the model’s spatial sensitivity of these quantities to growing-season (GS) temperature. Combined with spatially-resolved observations of LAI and GPP, we estimate that GS-LAI and GS-GPP increase by 17.0 ± 2.4% and 24.0 ± 3.0% per degree of warming, respectively. The observationally-derived sensitivities of LAI and GPP to temperature are about 40% and 71% higher, respectively, than the mean of the ensemble of simulations from TRENDY, primarily due to the model underestimation of the sensitivity of light use efficiency to temperature. We estimate that the regional mean GS-GPP increased 28.2 ± 5.1% between 1983–1986 and 2013–2016, much larger than the 5.8 ± 1.4% increase from the CO2 fertilization effect implied by Wenzel et al. This suggests that warming, not CO2 fertilization, is primarily responsible for the observed dramatic changes in the HLNF biosphere over the last century.

Figure S3 The spatially-derived relationship between GS-GPP and GS-temperature.The points include grid cells located north of 50 degrees latitude with at least 40% tree cover, with colorcoded by latitudes.A. OCO-2 SIF-constrained GPP; B. FLUXCOM GPP based on ANN algorthm; C. FLUXCOM GPP based on random forest (RF); D. FLUXCOM GPP based on MARS algorithm.X-axis: temperature (unit, °C), and y-axis is GPP with unit of gC/m 2 /day.The black curve shows exponential fitting between temperature and GPP while the gray line shows linear fitting.S1.
Figure S7 The spatial covariation between growing-season GPP (GS-GPP) and growing-season temperature during one of the 20-year groups (1971)(1972)(1973)(1974)(1975)(1976)(1977)(1978)(1979)(1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990).The points include grid cells located north of 50 degrees latitude with at least 40% tree cover listed in Table S1.The color bar represents latitudes of those points.On top of each panel is the name of TRENDY models.(1971)(1972)(1973)(1974)(1975)(1976)(1977)(1978)(1979)(1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990) for the models that were not selected in this study due to its low spatial coherence with temperature spatial gradient.The points include grid cells located north of 50 degrees latitude with at least 40% tree cover listed in Table S1.The color bar represents latitudes of those points.On top of each panel is the name of TRENDY models.Figure S12 The same as Figure S7, except that all the models using the same condensed MODIS tree cover fraction.We have excluded grid points with near zero GS-GPP values and grid cells that do not have clearly-defined growing season.Figure S19 The relationships between temporal percentage change of GPP and temperature increase due to CO2 and climate effect (blue) and climate effect only (black) for each of the selected TRENDY models.n is from 2 to 10. DGPP(CO2+T), DGPP(T), and are defined in section 3 in Method.

Figure
Figure S1 Workflow to derive the observation-constrained photosynthesis/greenness -climate feedback factors  !"" #$% and  &'( #$% , which are then used to predict changes of GPP and LAI from historical increase in temperature over the high latitude northern forests.S2 and S1 represent TRENDY S2 and S2 runs respectively.S2 runs have varying climate and CO2, and S1 runs have varying CO2.

Figure
Figure S4 The spatially-derived relationship between GS-LAI to GS-temperature.The points include grid cells located north of 50 degrees latitude with at least 40% tree cover.The dots are color-coded by latitudes.A: GIMMS MODIS LAI; B. GIMSS AVHRR LAI; C. NOAA AVHRR LAI; D. MODIS LAI; X-axis: temperature (unit, °C), and y-axis is LAI with the unit (m 2 / m 2 ).
Figure S5 The spatially-derived relationship between GS-fPAR (fraction of Photosynthetic Active Radiation) and GS-temperature.The points include grid cells located north of 50 degrees latitude with at least 40% tree cover.The color bar represents latitudes of those points.A: GIMMS MODIS; B. GIMSS AVHRR; C. NOAA AVHRR; D. MODIS MCD15 product; X-axis: temperature (unit, °C), and y-axis is fPAR with the unit (%).
Figure S6 The S1, S2, and S3 runs have similar magnitude of spatially-derived sensitivity of growing season LAI (A) and GPP (B) to temperature.The error bars are the standard deviations of the corresponding spatially-derived sensitivity among the 10 groups in each model.The model names corresponding to each model ID are listed in TableS1.

Figure
Figure S8The same as FigureS7, except this is for LAI with unit of m 2 /m 2 .

Figure
Figure S10 Growing season definition based on SIF-constrained GPP.

Figure
Figure S11 The relationship between growing season mean photosynthetic active radiation (PAR) (W/m 2 ) and growing season mean temperature (°C) colored by latitude (degrees).

Figure
Figure S13The same as FigureS12, except this is for LAI.
Figure S14 Same as Figure 1 in the main text, except that all models use the condensed MODIS tree cover fraction.

Figure
Figure S15 A. The relationship between the spatially-derived and the temporally-derived sensitivities of LAI to temperature in the simulations ( &'( ! ).Different from Figure 2A in the main text, we added model D with tundra included (labeled as D+).Model D simulates time-invariant LAI over tundra.

Figure
Figure S18 Temperature driven change of NH CO2 seasonal cycle amplitude.Solid black line: the observed CO2 seasonal cycle amplitude (SCA) changes between IGY (1958-1963) and HIPPO (2009-2011) aircraft campaigns at 500 hPa.Dashed black line: the approximate CO2 SCA change between IGY and HIPPO at 500 hPa attributed to the HLNF, which is the multiplication between observed CO2 SCA and the ratio between model simulated CO2 SCA forced by forest NEE over 50-75°N and the observed CO2 SCA (Liu et al., 2020).Red: the calculated CO2 SCA changes due to the temperature increase of the HLNF.Blue: The simulated CO2 SCA changes due to climate using the TRENDY S1 and S2 runs.The grey and red lines are the same as Figure 4 in Liu et al., 2020.

FigureFigure
Figure S20The same as the FigureS19, except this is for LAI.

Table S1 .
Summary of TRENDY v6 models.Only models that have monthly GPP and LAI from both S1 and S2 runs, and either R 2 (LAI, T) or R 2 (GPP, T) larger than 0.2 were selected.For models without land cover fraction information, the compressed IGBP data from MODIS was used to identify grid with at least 40% tree cover.