Observational Constraints on the Response of High‐Latitude Northern Forests to Warming

Since the 1960s, carbon cycling in the high‐latitude northern forest (HLNF) has experienced dramatic changes: Most of the forest is greening and net carbon uptake from the atmosphere has increased. During the same time period, the CO2 seasonal cycle amplitude (SCA) has increased by ~50% or more. Disentangling complex processes that drive these changes has been challenging. In this study, we substitute spatial sensitivity to temperature for time to quantify the impact of temperature increase on gross primary production (GPP), total ecosystem respiration (TER), the fraction of Photosynthetic Active Radiation (fPAR), and the resulted contribution of these changes in amplifying the CO2 SCA over the HLNF since 1960s. We use the spatial heterogeneity of GPP inferred from solar‐induced chlorophyll Fluorescence in combination with net ecosystem exchange (NEE) inferred from column CO2 observations made between 2015 and 2017 from NASA's Orbiting Carbon Observatory‐2. We find that three quarters of the spatial variations in GPP can be explained by the spatial variation in the growing season mean temperature (GSMT). The long term hindcast captures both the magnitude and spatial variability of the trends in observed fPAR. We estimate that between 1960 and 2010, the increase in GSMT enhanced both GPP and the SCA of NEE by ~20%. The calculated enhancement of NEE due to increase in GSMT contributes 56–72% of the trend in the CO2 SCA at high latitudes, much larger than simulations by most biogeochemical models.


Introduction
The High-Latitude Northern Forest (HLNF, 50°N to 75°N) is composed of needle leaf trees (both evergreen and deciduous), deciduous broadleaf trees, and shrubland ( Figure 1). It is responsible for~20% of global gross primary productivity (GPP) (Jung et al., 2017(Jung et al., , 2020Tramontana et al., 2016). Over the region northward of 50°, precipitation substantially exceeds evaporation due to moisture convergence by the largescale eddy circulation and low net radiation, which results in a generally well-watered, thermally limited ecosystem (Nemani et al., 2003). The northern extent of the forest is limited by the Arctic ocean and/or where summertime temperatures are too low. Because of its seasonality (days with nearly 24 hr of sunlight followed by days with 24 hr of darkness), uptake and release of carbon from the atmosphere is responsible for a large fraction of the seasonal cycle amplitude (SCA) in atmospheric CO 2 . Without the HLNF, the SCA of CO 2 in the Northern Hemisphere (NH) would be only about half as large (Graven et al., 2013).
Since the 1960s, carbon cycling in the HLNF has experienced dramatic changes: Most of the forest is greening (Keenan & Riley, 2018;Zhu et al., 2016), and net carbon uptake from the atmosphere has ©2020 Jet Propulsion Laboratory. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
increased . Although the underlying processes driving these trends are uncertain (Bastos et al., 2019;Forkel et al., 2016;Gray et al., 2014;Piao et al., 2018;Randerson et al., 1997;Zeng et al., 2014), it is generally agreed that both the increase in temperature and CO 2 are responsible (Bastos et al., 2019;Forkel et al., 2016;Piao et al., 2018;Randerson et al., 1997). Global mean CO 2 concentration has increased from 320 ppm in 1958 to~400 ppm in 2010s (https://www.esrl.noaa.gov/gmd/ccgg/trends/). Biogeochemical processes ranging from photosynthesis, net primary production (NPP), to the net carbon balance in general respond positively to the CO 2 increase, which is called CO 2 fertilization (Kimball, 1983;Long et al., 2004). The optimal temperature for photosynthesis is currently higher than the mean annual temperature over high-latitude biomes (Huang et al., 2019); thus, we expect that increase in temperature accelerates GPP, unlike tropics where the increase in temperature can suppress growth (Sullivan et al., 2020). With the increase of GPP and temperature, total ecosystem respiration (TER) is also expected to increase (Baldocchi, 2008;. The temperature increase over the NH high latitudes is nearly double that of the mean temperature increase over land (IPCC AR5). Quantifying the individual role of temperature and CO 2 increases on photosynthesis, respiration, and net carbon uptake is crucial for projecting future climate change (Friedlingstein et al., 2014). However, because the temporal changes in temperature and CO 2 have been so correlated, disentangling their individual contribution with observations is difficult, especially over HLNF where conventional observations are sparse (Schimel et al., 2015). Previous studies primarily make use of factorial simulations of biogeochemical models to isolate the impact of change in climate from other factors (e.g., Bastos et al., 2019). Studies using observation constraint are scarce.
In this study, we provide observational constraints on the impact of temperature change on GPP, TER, and the resulting net ecosystem exchange (NEE) of the HLNF. We use the spatial heterogeneity of GPP and TER inferred from solar-induced chlorophyll fluorescence (SIF) (Frankenberg et al., 2011) and a combination of NEE and GPP, respectively. The NEE is inferred from column CO 2 observations retrieved from spectra of NASA's Orbiting Carbon Observatory-2 (OCO-2) Eldering et al., 2017) between 2015 and 2017.
Substituting the observed spatial sensitivity of GPP and TER to temperature for time, we estimate the change in carbon cycling in the HLNF between 1960s and 2010s due to the warming this region has experienced. We diagnose what contribution these changes have made to the observed increase in the CO 2 SCA between 1958-1963 aircraft campaign) (Keeling et al., 1968) and 2009-2011 aircraft campaign) (Wofsy et al., 2011) and over NH high latitude CO 2 flask sites with long record (>25 years) (Alert (Canada), Cold Bay (Alaska, USA), Barrow (Alaska, USA), and Shemya Island (Alaska, USA). GPP is a function of fraction of photosynthetic active radiation (PAR) absorbed by plants (fPAR) (Farquhar et al., 1980), which reflects forest structure (leaf area index, LAI) (Zhu et al., 2016). To evaluate whether space and time can be substituted, we compare the trend and interannual variability in fPAR and GPP calculated with this assumption to observations of the changes over space and time.

Materials and Methods
We calculate monthly mean GPP (2.1) and NEE (2.2) for 2015-2017 with OCO-2 SIF and column CO 2 retrievals, respectively. TER is calculated as the difference between NEE and GPP (2.2) ( Figure S1 in the supporting information, SI). We substitute seasonal (i.e., spring, summer, and fall (2.3)) spatial GPP-T and TER-T relationship for time to predict the temporal changes of GPP and TER and the resulting changes in NEE from increase in temperature for the time period of 1980-2011 and between 1958-1963 and 2009-2011 (2.4) ( Figure S1). To evaluate the space-for-time substitution, we compare predictions of the change in fPAR between 1982 and 2012 based on the modern spatial fPAR-T relationship with the changes in GIMMS3g fPAR (Zhu et al., 2013) (2.4). We further compare predictions of seasonal cycle and interannual variability in GPP with observations of SIF-constrained GPP and data from FLUXNET2015 data set (2.4). Finally, we calculate the CO 2 SCA changes due to changes in NEE between 1958NEE between -1963NEE between and 2009NEE between -2011 and for the time period of 1980-2011 (2.5) ( Figure S1).

SIF-Constrained Monthly Mean GPP Calculation
We use SIF retrieved from OCO-2 to calculate monthly mean GPP. SIF is fluorescence from plants that occurs during photosynthesis. It is sensitive to both plant structure and photosynthetic function and is linearly correlated with GPP under high-light conditions at seasonal time scale (Frankenberg et al., 2011;Sun et al., 2018). OCO-2 SIF retrievals have footprint size~1.2 × 2 km 2 . Each SIF retrieval is scaled based on solar zenith angle, date, latitude, and local passing time,~1:30 pm, to produce a daily mean (Frankenberg et al., 2011). The monthly mean value at each 4°× 5°grid is a simple averaging of the daily mean SIF values. The 4°× 5°spatial resolution is chosen for the consistency with the optimized NEE from top-down atmospheric CO 2 flux inversion. We use OCO-2 level 2 SIF observations from land nadir mode only between January 2015 and December 2017. The time frame is constrained by the availability of the FLUXCOM GPP product, an upscaling GPP product based on GPP from flux towers, ancillary satellite, and meteorology (Jung et al., 2017(Jung et al., , 2020Tramontana et al., 2016).
The growing season mean SIF value increases exponentially with temperature over the grid points with tree cover larger than 40% ( Figure S2). The SIF-temperature sensitivity is within the sensitivity range of the three FLUXCOM GPP products ( Figure S2). There is no obvious discontinuity in the sensitivities of both SIF and FLUXCOM GPP products to temperature, so we calculate monthly mean GPP over 50-75°N at 4°× 5°spatial resolution as the product of OCO-2 SIF and a scaling factor derived using the growing season (GS) mean GPP from three FLUXCOM products and mean SIF. The scaling factor, 22.9 ± 2.6 (gram carbon micron sr per day per watts, gC μm sr day −1 W −1 ), is consistent with (Sun et al., 2018). In using the same conversion factor across the region, we have assumed that different tree types have the same relationship between SIF and GPP (see section 4).

Top-Down NEE Estimation With CMS-Flux and TER Calculation During 2015 and 2017
We optimize monthly mean NEE with CMS-Flux system Liu et al., 2014Liu et al., , 2017Liu et al., , 2018 at each 4°× 5°grid with OCO-2 B9 land nadir observations (Eldering et al., 2017) during September 2014 to April 2018. In spite of increasing CO 2 observational coverage from OCO-2, the spatial resolution of NEE from top-down inversion is still limited by observations. For readers interested in more details on atmospheric CO 2 flux inversion, please refer to Liu et al. (2014Liu et al. ( , 2015Liu et al. ( , 2017. The prior fluxes include climatology terrestrial biosphere carbon fluxes from CASA-GFED 3 (Randerson et al., 1997), biomass burning from  (Table S1). GPP is scaled SIF.

AGU Advances
GFED4, fossil fuel emissions from ODIAC (Oda et al., 2018), and ocean carbon fluxes from ECCO-Darwin (Menemenlis et al., 2008). Except the climatology terrestrial biosphere fluxes from CASA-GFED3, the other types of carbon fluxes are yearly dependent.
We calculate TER as the difference between the optimized NEE and GPP between 2015 and 2017. The right three panels in Figure S3 show the dependency of TER on temperature during the spring, summer, and fall. The estimated TER shows strong dependence on temperature, which is expected, since TER's two components, autotropic and heterotrophic respiration, depend on the amount of carbon fixed through photosynthesis (Baldocchi, 2008). The temperature dependence of TER is not related to the amount of above ground biomass, which has a weak relationship with TER on seasonal timescale (not shown), with R 2 less than 0.07. In Figure S3, we remove points that have nonphysical negative TER values.
We validate the top-down flux inversion results by comparing the CO 2 SCA sampled at Mauna Loa against observed CO 2 SCA and comparing posterior CO 2 concentrations between April and September against independent aircraft CO 2 observations north of 30°N. The SCA of posterior CO 2 concentration sampled at Mauna Loa is within the uncertainty range of the observed CO 2 SCA ( Figure S4). The random error and bias of regional mean monthly posterior CO 2 concentration is 0.4 and 0.1 ppm, respectively ( Figure S4 in SI).

Definition of Growing Season
We define spring, summer, and fall with SIF-optimized GPP product (Parazoo et al., 2014), which optimally combines OCO-2 SIF observations with climatology mean state from TRENDY models (Sitch et al., 2015). This product is only used in defining seasons, since the product maintains essentially the same seasonality as original SIF observations during growing season, while filling the gaps of SIF observation during the late growing season and winter. At each grid point, the growing season is defined as the time period when the GPP values are larger than 20% of the maximum GPP at that grid point. Spring and fall are defined as the time period when GPP is between 20% and 75% of maximum GPP, and summer is the time period between spring and fall. We choose 20% as the threshold to ensure that the plants have started growing and the start of the spring is well defined. We find that 10% of grid points would have no spring if the threshold is 30% instead of 20%, since GPP increases sharply during spring, especially in Eurasia.

Prediction of Growing Season Monthly Mean GPP, TER, and fPAR
We substitute the spatial sensitivity to temperature for time to predict the temporal changes in GPP, TER, and fPAR due to changes in temperature. Note that the spatial CO 2 gradient is typically less than 15 ppm, and we find that there is no covarying relationship in space between GPP and CO 2 on seasonal time scale over the HLNF (not shown). Our approach of substitution of space for time to predict changes in carbon states due to climate has been used by others. For example, Sullivan et al. (2020) predicted changes in tropical biomass from changes in temperature using the spatial biomass-temperature sensitivity. Figure S5 shows that the fPAR-T spatial sensitivity is relatively constant, with variability less than 5%, implying that the spatial sensitivity to temperature represents a semiequilibrium state. We use the spatial-temperature sensitivity shown in Figures S3 and S6 to predict temporal changes of GPP, TER, and fPAR.
In predicting historical GPP, we use GPP/PAR-T spatial sensitivity ( Figure S3) to calculate the temporal changes of GPP/PAR with changes in temperature, and then multiply PAR to get final GPP prediction.
Since PAR is only a function of latitude and time of year, we use the same PAR averaged over 2015-2017 to calculate GPP.
With the spatial-temperature sensitivity derived from the observation-constrained GPP and TER between 2015 and 2017 and CRU monthly mean air temperature during 2009-2011 and 1958-1963, we predict the monthly GPP and TER during growing season of these two time periods ( Figure S7). In predicting the GPP for 2009-2011 and 1958-1963 time periods, we assume that the temporal range of spring, summer, and fall in each grid cell is the same as in 2015-2017, which is a reasonable assumption since the temporal resolution of our study is month and the growing season length only increases by~0.25 month for every~1°i ncrease of temperature ( Figure S8). We generate 10-member ensemble GPP and TER predictions with the 95% confidence interval shown in Figure S3. We calculate the uncertainty from the ensemble GPP and TER hindcasts shown in Figure S7.

AGU Advances
We calculate NEE during growing season as the difference between ensemble TER and GPP and assume that nongrowing season NEE is proportional to the growing season NEE (Kauppi & Posch, 1985). The proportional assumption may be a conservative estimation of nongrowing season NEE changes (see section 4).
To calculate NEE between 1980 to 2011, we predict growing season GPP and TER with 4-year frequency, where we use monthly mean temperature averaged over every 4 years. We find that the results are not sensitive to the prediction frequency. In total, we have eight groups of NEE hindcast between 1980 and 2011. We choose 2011 as the ending year to be consistent with the time frame of HIPPO aircraft campaign.
With the same seasonally dependent GPP/PAR-T spatial sensitivity and monthly mean temperature, we predict monthly mean GPP between 1960 and 2014. Following the same method described above, we predict monthly mean fPAR with the seasonally dependent relationship calculated for 2013-2016 ( Figure S6). The sensitivity shown in Figure S6 is based on fPAR from GIMMS3g. Since fPAR from GIMMS3g has more than 30-year record, the trend calculation is more robust than that is based on MODIS. In the main text, we use fPAR from GIMSS3g in trend calculation but use fPAR from MODIS in Figure 2, since MODIS has higher spatial resolution than AVHRR used in the earlier record of GIMMS3g (https://nsidc.org/support/ 21524371-What-are-the-MODIS-sensor-s-strongest-attributes-and-how-do-they-compare-to-other-sensors).

CO 2 Simulations With GEOS-Chem Transport Model
We use GEOS-Chem transport model to simulate 3-hourly 10-member ensemble of CO 2 concentrations during 1958-1963 (IGY) and 2009-2011 (HIPPO) and from 1980 to 2011. The only surface boundary conditions are the predicted ensemble monthly NEE (i.e., TER-GPP). Therefore, the differences in simulated CO 2 are only due to the differences in the predicted monthly NEE between 50°N and 75°N. We assume the same diurnal biosphere fluxes during these two time periods, which are from CASA-GFED3.
We run GEOS-Chem transport model at 2°× 2.5°spatial resolution driven by MERRA2 reanalysis after regridding 4°× 5°NEE to 2°× 2.5°spatial resolution. The transport model has smaller transport errors over high latitudes at 2°× 2.5°resolution than at 4°× 5°resolution. For both IGY and HIPPO CO 2 simulations, we run ensemble transport model simulations continually from 2006 to 2011, and we then analyze CO 2 simulation results over the last 3-year simulation. The difference between IGY and HIPPO CO 2 simulations is only from the boundary forcing, not from the meteorology forcing. For each of the eight-group NEE from 1980 to 2011, we run ensemble transport model simulations continually for 6 years with meteorology from 2006 to 2011 and then analyze CO 2 simulation results for each group over the last 3-year simulation. As a result, the changes of CO 2 SCA are only due to changes of NEE, not due to changes of transport. We follow the same method described in Graven et al. (2013) (SI) to calculate CO 2 SCA for both surface sites and aircraft observations and the corresponding model simulations.

The Spatial GPP-Temperature Relationship
Over the HLNF, both the spatial distribution and seasonal evolution of GPP closely follow those of temperature ( Figure 1). From April to July, increases in GPP propagate from the west to the east of each continent following the spatial increase of monthly mean temperature ( Figure S9). From July to October, GPP decreases as cooler temperature propagates from the east to the west of each continent ( Figure S9). The spatial changes illustrate how the onset and length of GS are closely coupled to the surface temperature ( Figure S8). The growing season length increases~1 week with each degree increase in surface temperature, consistent with Park et al. (2016).
Quantitatively, the spatial variation in GSMT explains 70% of GPP/PAR (Photosynthetic Active Radiation) variability (Figures 1b, 1c, and 2a) (and 75% of the variability in GPP; Figure S10). The ratio between seasonally averaged GPP and the seasonally averaged PAR removes the influence of the somewhat larger solar flux at lower latitude. While PAR is a function of latitude over this domain, it contributes less than 10% to the GPP-temperature spatial relationship ( Figure S10). The correlation between GPP/PAR and GS mean temperature is not simply due to variations in latitude: Trees growing at the same latitude have large variability in productivity due to the longitudinal variation in GSMT. For example, at 62°± 2°N latitude (highlighted by 10.1029/2020AV000228

AGU Advances
closed circles in Figure 2a), GPP/PAR varies from 0.10 to 0.45 gram carbon per mega Joules (gC MJ −1 ), highly correlated with the variation in GSMT.
The distribution of vegetation cover is also highly correlated with GSMT ( Figure 2a). The mean GS temperatures are 9.2°C, 11.5°C, and 13.6°C for shrubland, needleleaf forest, and deciduous boreal forest, respectively. Between the transitional zones, mixed forest types coexist. The correlation illustrates the complex role that temperature likely plays on many aspects of the ecology of the HLNF. From this correlation, we might anticipate that warming would lead to succession of the dominant trees in a fashion consistent with the patterns observed spatially (e.g., evergreen forests will transition to deciduous forests as the GSMT reaches 13°C). Limited field measurements are consistent with this hypothesis (Bjorkman et al., 2018;Kharuk et al., 2007).
Why is GPP/PAR so sensitive to temperature in the HLNF? GPP/PAR is a product of fPAR and the light use efficiency (the efficiency of conversion of absorbed sunlight to assimilated carbon) (Farquhar et al., 1980). The fPAR is observed to increase with temperature in the HLNF (Figure 2b). The fraction of absorbed sunlight resulting in fixation of CO 2 (light use efficiency) depends on, among other things, plant functional types. Field experiments (Ahl et al., 2004) show that the light use efficiency of deciduous forest is~30% higher than needle leaf trees, consistent with the behavior seen in Figure 2c. In the absence of water or nutrient limitations, light use efficiency is also expected to increase with temperature in cold forests ( Figure 2c) (Horn & Schulz, 2011). Due to the high rates of precipitation, snow melt, and low net radiation, water limitations-especially earlier in the growing season-are not expected. As shown in Figure S11, GPP increases with the increase in vapor pressure deficit, opposite to the expected relationship for water-limited ecosystems. Because the rate of nitrogen (N) fixation also increases rapidly with temperature (Schimel et al., 1997;Welter et al., 2015), N limitation may not change as plants increase their demand as temperature increases. Thus, the sensitivity of GPP/PAR to GSMT results from the sensitivity of both forest structure (fPAR) and light use efficiency, with similar contributions (0.090 vs. 0.100) from each (Figures 2b and 2c).

Hindcast of the fPAR and GPP/PAR From the Changes in GSMT Since 1960
With more than 30 years (1982-2016) of fPAR record from GIMMS3g (Zhu et al., 2013), we can test how well the temporal change in fPAR is explained by the observed spatial temperature sensitivity. We find that the modern correlation of fPAR with GSMT accurately captures the variability and trend in fPAR over the past several decades. Remarkably, even the interannual variations in fPAR (Figure 3) are well predicted by the spatial fPAR-T correlations, with R 2 between observation and predicted fPAR equal to 0.62 (Figure 3d). The calculated temporal trend of fPAR during the period of fPAR observation is 0.20 ± 0.07%/year, same as the observed fPAR trend (0.20 ± 0.06%/year). Spring has the largest increasing trend, and R 2 between observation and predicted fPAR is the largest (R 2 = 0.84) ( Figure S12). The comparisons between fPAR observations and prediction indicate that the current spatial correlation between fPAR and GSMT can explain the historical changes in HLNF mean fPAR without invoking a significant role for CO 2 fertilization or other factors.
The hindcast fPAR and GPP show larger increasing trend during the later time period (after 1983). The fPAR trend is 0.05 ± 0.08%/year from 1960 to 1983, while it is 0.20 ± 0.07%/year during 1983-2014. The GPP increasing trend is 0.06 ± 0.19%/year and 0.45 ± 0.17%/year, respectively, during these two time periods. The larger increasing trend originates from the changing trend in the temperature. The trend in Figure 2. The spatial relationship between temperature and gross primary production (GPP) and its contributing factors between 50°N and 75°N. The x axis is the growing season mean temperature (°C), GSMT. (a) GSMT versus the ratio between GPP and photosynthetic active radiation (PAR) (unit: gram carbon per mega joules (gC MJ −1 ); (b) GSMT versus fraction of photosynthetic active radiation absorbed by the forest (fPAR) (%); (c) temperature versus light use efficiency (unit: gC MJ −1 ). Each point is an average over growing season at individual grid cell. Green: needleleaf forest; brown: deciduous boreal forest; blue: shrubland; the size of circles is proportional to the percentage of tree coverage. Filled circles are the points with latitude equal to 62°N. The minimum of tree coverage for each type is 40%. The uncertainty is 3σ of the fitting. The temperature 2-m air temperature from ERA-interim reanalysis (Table S1).

10.1029/2020AV000228
AGU Advances temperature is 0.01 ± 0.02°C/year during 1960-1983, while it is 0.04 ± 0.02°C/year during 1983-2014. The acceleration of temperature increasing rate after 1980s agrees with Sánchez-Lugo et al. (2018), which shows more than two times greater of temperature increasing rate over the high latitudes after 1980s.
The hindcast of fPAR also captures the spatial inhomogeneity in the observed fPAR changes. Between 2007Between -2011Between and 1983Between -1987, the northeast Eurasia and western Europe have experienced larger increases in GSMT than North America (NA) (Figure 3a). The changes in both the observed and hindcast GS fPAR show similar spatial pattern as the temperature differences between these two time periods (Figures 3b and 3c). This, again, is consistent with temperature playing the dominant role in controlling fPAR. For example, if CO 2 fertilization was responsible for enhanced fPAR, we would not expect the fractional change to be so well correlated with the spatial change in temperature given that the change in CO 2 was identical across both space and time. The overestimate in the calculated fPAR changes over NA may be caused by the land cover changes from disturbance  that the hindcast does not consider.
The trend in the hindcast GPP during GS is 0.36 ± 0.07%/year, about two times the hindcast fPAR trend during 1960-2014, consistent with Figure 2 showing that only half of the GPP/PAR-T spatial sensitivity is attributed to change in fPAR. The trend in GPP during spring is larger than summer and fall (0.43% ± 0.09% vs. 0.30% ± 0.08% and 0.29% ± 0.10%) (Figure S12), as a result of a larger increase in temperature in spring ( Figure S12) and larger PAR value in spring than late summer and fall. Consequently, the increasing trend in spring temperature has a larger impact on GPP and fPAR (Zhang et al., 2020).  (Table S1).

AGU Advances
Our analysis hinges upon the validity of space-for-time assumption. The comparison between fPAR observations and the corresponding hindcast provides compelling support for this assumption. In the SI (Peer Review History), we further test the space-for-time assumption with GPP constrained with SIF and GPP from 11 sites of FLUXNET2015 data set (https://fluxnet.org/data/fluxnet2015-dataset/). We compare observations and predictions in terms of seasonality and interannual variability (Figures S13-S15). Both tests support the space-for-time assumption.
We use fPAR-T spatial relationship and its temporal prediction to test whether the conclusion depends on spatial resolution. The fPAR at 1°× 1°resolution has almost the same spatial sensitivity to temperature as at 4°× 5°resolution at seasonal time scale (0.050 vs. 0.054, 0.074 vs. 0.079, and 0.056 vs. 0.055) ( Figure S16). As a result, the predicted fPAR trend is essentially identical to that at 4°× 5°resolution (0.021 ± 0.07%/year vs. 0.020 ± 0.07%/year) ( Figure S17). This test indicates our method and conclusion is robust at the spatial scale of hundred kilometers. An additional test with site-level data from FLUXNET data set ( Figures S14 and S15), which has spatial representativeness of 1-10 km 2 , further supports the robustness of the conclusion with respect to spatial resolution. Whether the method and conclusion are still valid at even smaller spatial scale needs further test with high-resolution data in the future.
We have assumed that the temporal range of spring, summer, and fall in each grid cell is the same between 1960 and 2014 and in 2015-2017. There is a reasonable possibility that a given month may change from "spring" in the 1960s to "summer" in later time period, since our temporal resolution is only monthly. To test the sensitivity of the prediction to the separation of seasons, we move the ending month of spring at each grid cell 1 month later between 1960 and 1979. Figure S18 shows that the trend of predicted fPAR between 1960 and 2014 is somewhat higher (0.16 ± 0.03%/year) but within the uncertainty of the control prediction (0.14 ± 0.03%/year) where the same season definition is used throughout the time period. With the available of high spatiotemporal SIF observations from TROPOspheric Monitoring Instrument (TROPOMI, launched in October 2017) (Köehler et al., 2018), future studies can increase the temporal resolution and take into account the changes of growing season length in the prediction.

The Impact of Increasing Temperature on the NH CO 2 Seasonal Cycle Amplitude
The increase of GPP will increase both net carbon uptake during growing season and net carbon release during nongrowing season (due to increased carbon pool size and increases in the rate of ecosystem respiration with temperature) (Zhao et al., 2014). As a result, it is expected that the increase in temperature will enhance the SCA of NEE. In turn, the increase of NEE SCA may have contributed to the documented large increase of NH CO 2 SCA (Graven et al., 2013). Between 1960 and 2010, the GPP and TER are calculated to have increased 20 ± 12% and 18 ± 12%, respectively, between April and September ( Figure S7). The NEE between April and Sep increased by 23 ± 12%. The NEE amplitude, which is defined as the difference between the net source and net sink, increases by 21 ± 8.0% over 50-75°N ( Figure S7). The net carbon sink increased by 0.25 ± 0.17 GtC between 1960 and 2010. The estimate of net carbon sink has, however, large uncertainties that require further study to evaluate.
We sample the CO 2 concentrations along the IGY and HIPPO aircraft locations, respectively, in the midtroposphere (~500 hPa) and those four high latitude sites. Figure S19 summarizes the CO 2 SCA changes from observations and predictions. The increase in GSMT is estimated to contribute 19-42% of the total observed CO 2 SCA increase north of 45°N between two aircraft campaigns and 34-76% of the CO 2 SCA increase due to changes in carbon cycling in the HLNF (Figure 4). The observed CO 2 SCA changes attributed to HLNF (dashed line in Figure 4) is calculated as the multiplication between observed CO 2 SCA and the ratio between model simulated CO 2 SCA forced by NEE over 50-75°N and the observed CO 2 SCA ( Figure S19).
The surface sites may more directly reflect the changes of NEE over the HLNF, and frequent measurements make it more precise compared to aircraft observations. At the four surface sites, increase in the GSMT is calculated to contribute 56-72% of the observed trends in the CO 2 SCA (Figure 4). The NEE over 50-75°N accounts for~60% of the CO 2 SCA at these surface sites. Assuming that a similar or even higher proportion (60-70%) of the trend in the observed CO 2 SCA can be attributed to HLNF and the rest is due to changes over other regions, then the increase in GSMT can explain almost the total observed CO 2 SCA trend attributed to HLNF, which agrees with the dominant temperature impact on the fPAR trend discussed earlier. Since our 10.1029/2020AV000228

AGU Advances
focus is the temperature impact on the carbon exchange over the HLNF, the exact proportion of the CO 2 SCA trend attributed to HLNF is out of scope of this study.

Discussion
Though the predicted changes in carbon cycling in the HLNF is supported by the fPAR and GPP observations, further study is needed to test some assumptions used in the study. We have assumed a single scaling factor between the observed SIF and GPP. While this assumption is supported by (Li et al., 2018;Sun et al., 2018), it requires additional testing within these biomes. In particular, about 50% of the observed correlation between SIF and temperature is not explained by fPAR but by the "light use efficiency." The estimated change in light use efficiency is correlated with biome (evergreen vs. shrub vs. deciduous trees; Figure 2c), and further work is needed to test whether all the differences in SIF of these biomes is indeed related to the differing photosynthetic efficiency of these tree types or may (partially) reflect differences in N content, or the optics of SIF within canopies differing in structure. Though the hypothesis that N fixation increases with warming has some support in the literature (e.g., Schimel et al., 1997), it is clear that more observations are needed to test such hypothesis.
The spatial covariation between TER and temperature may include some confounding factors, such as the available soil organic carbon and sensitivity of microbial activity to temperature (e.g., Nottingham et al., 2019). Here TER-temperature spatial sensitivity is a statistical relationship of the net impact of temperature on TER. A detailed process understanding on factor controls of TER still has large uncertainties and requires more research.
Though studies show that the magnitude of nongrowing season NEE is positively correlated with GS NEE (Yuan et al., 2011), the assumption that the change in nongrowing season NEE is proportional to the changes in GS NEE may not always be valid. In using this assumption, we have ignored the impact of changes of snow coverage and freezing/thaw on ecosystem respiration (Natali et al., 2019) and nonhomogeneous changes in temperature throughout the year. There is evidence that the rate of respiration outside the GS may have increased faster than during GS because of larger temperature increase during fall and winter (Schuur et al., 2015).  (2009( -2011( ) and IGY (1958( -1963 aircraft campaigns at 500 hPa. Dashed black line: the approximate CO 2 SCA change between HIPPO and IGY at 500 hPa attributed to the HLNF, which is the multiplication between observed CO 2 SCA and the ratio between model simulated CO 2 SCA forced by forest NEE over 50-75°N and the observed CO 2 SCA ( Figure S19) Remarkably, the hindcast based on the spatial temperature sensitivity captures both the long-term trend and some fraction of short-term interannual variability (Figures 3 and S12). We attribute this skill in part to the use of a single curve of the spatial data across multiple plant functional types and climate zones. The fit thus can capture possible species succession, which takes several decades, since the single curve assumes the same sensitivity to temperature among different forest types (Figure 2). At the same time, it also predicts with some skill the faster changes in GPP and other quantities associated with interannual temperature differences, because the single curve spans large range in temperature ( Figures S3 and S6) and the temperature anomaly of any single year likely sits within the range. This is especially true of springtime when the initiation of photosynthesis is be highly sensitive to small changes in temperature.
The rapid increase in GPP with GSMT is certainly limited: the data shown in Figure 2 extends only to a GSMT of 16°C. It is expected that at some point, photosynthetic, water, and other nutrient limitations will limit further increase in GPP with temperature, especially during summer. Even today there is evidence that by mid/late summer such limitations come into play (e.g., Buermann et al., 2018), and browning trend has been observed over limited regions (Beck & Goetz, 2011) Radiation may also limit the response of GPP to temperature increase in the future, especially in the fall (Zhang et al., 2020). In addition, the increase in abiotic (e.g., fire) and biotic (e.g., herbivory) disturbance with climate will also impact productivity (Sulla-Menashe et al., 2018).
Despite the limitations noted above, the spatial variations in fPAR and GPP in the HLNF are remarkably well correlated with temperature (Figures 2 and S2), consistent with a thermally limited ecosystem. The strong temperature control on changes in fPAR and GPP and its large contributions to changes in CO 2 SCA are in contrast to results obtained from most (Bastos et al., 2019;Piao et al., 2018) but not all biogeochemical models (Forkel et al., 2016), which generally suggest a dominant CO 2 fertilization effect on the changes of NEE over HLNF (and thereby the CO 2 SCA). While Zhu et al. (2016) show dominant climate impact on the greening trend (i.e., leaf area index) over high latitudes, later studies (Bastos et al., 2019;Piao et al., 2018) suggest a dominant CO 2 fertilization effect on CO 2 SCA changes with a similar set of TRENDY models. The high degree of spatial correlation between fPAR, GPP, and GSMT provides a key observational measure that can be used to evaluate and improve upon corresponding model simulated properties. A factorial simulation with changes only in temperature drivers for biogeochemical models will be able to isolate temperature effect from other factors and thereby test the model behavior against observations to identify possible model deficiencies.
That warming alone accounts for nearly all the observed fPAR trend is in agreement with prior observational studies that suggest CO 2 fertilization is not a major contributor to carbon exchange in mature forests (Jiang et al., 2020) and, more broadly, in the boreal (Girardin et al., 2011). Over mature forest, even under warmer climate (when CO 2 fertilization is expected to contribute much more than in cooler climates), GPP increases 12% with 150 ppm CO 2 enrichment (Jiang et al., 2020). Extrapolating the results in Jiang et al. (2020), thẽ 80 ppm increase in CO 2 between 1960s and 2010s implies a~6% increase of GPP over boreal forest, where a large proportion of forests is undisturbed, especially over Eastern and Northern Canada and Eurasia (https://daac.ornl.gov/NACP/guides/NA_Tree_Age.html). Since fPAR accounts for half of GPP increases, fPAR may only increases by~3% from CO 2 fertilization between 1960s and 2010, well within the uncertainty estimates shown in Figure 3. The more important role of temperature implied by the current spatial dependence of fPAR on GSMT is consistent with warming explaining most of the temporal trend in fPAR.
The ensemble of CMIP5 models suggest that by 2100, absent of efforts to mitigate climate forcing, these forests will experience an additional 2-5°C increase relative to 1986-2005 (RCP6.0, IPCC AR5) in GSMT. Figure 2 implies that such a change may dramatically alter the distribution of tree species over the region: evergreen trees will invade the shrublands all the way to the Arctic ocean, and much of the current evergreen forests will transition to deciduous trees. GPP over the region has the potential to double. Long-term consistent observation records are essential to monitor such possible changes over HLNF.

Conclusion
Three quarters of the spatial variations in gross primary production and fPAR over high-latitude northern forest can be explained by the spatial variation in growing season mean temperature. Substituting 10.1029/2020AV000228 AGU Advances observed spatial sensitivity to temperature for time, we estimate the change in carbon cycling in the HLNF between 1960s and 2010s and diagnose what contribution these changes have made to the observed increase in the CO 2 SCA. The hindcast captures both the magnitude and spatial variability of the trends in observed fPAR. We calculate that the observed high-latitude warming results in a 20 ± 12% and 18 ± 12% increase in GPP and TER, respectively, during growing season. The calculated changes in GPP and TER result in an increased seasonality of NEE by 21 ± 8% between the 1960s and 2010s. In turn, the increase of the NEE SCA contributes 56-72% of the observed increase in the CO 2 SCA at high latitude continental surface sites and 17-43% of the observed change in the SCA in the midtroposphere.