2.1 Input to the VODCA2GPP model

2.1.1 VODCA

The Vegetation Optical Depth Climate Archive (VODCA v1; Moesinger et al., 2020) consists of three single-frequency VOD products (Ku-, X-, and C-band), covering the period from 1987–2017 (Ku-band), 1997–2019 (X-band), and 2002–2019 (C-band). For VODCA2GPP, we used an updated VODCA version (VODCA v2 CXKu; Zotta et al., 2022) that merges all bands into a single dataset to obtain increased spatial and temporal coverage and reduce random errors compared to VODCA v1. VODCA v2 CXKu utilizes observations from the same sensors and frequencies as VODCA v1 (Table 1) to generate a single long-running, multi-frequency VOD time series. VODCA v2 CXKu is obtained by first scaling VODCA v2 observations from the C- and Ku-band to the X-band to remove systematic biases and then by computing a weighted average in order to fuse overlapping observations. The reference frequency for the scaling of the different frequencies is the X-band. VODCA v2 CXKu provides a single long-term vegetation metric covering over 30 years of observations (1988–2020) and thus exceeds the temporal length of the single-frequency multi-sensor products (VODCA v2 C, X, and Ku). VODCA v2 CXKu merges 15 passive nighttime VOD datasets retrieved from seven different sensors via the Land Parameter Retrieval Model (LPRM; Van der Schalie et al., 2017). The LPRM is based on radiative transfer theory introduced by Mo et al. (1982), uses forward modelling to simulate the top-of-atmosphere brightness temperatures under a wide range of conditions, and minimizes its difference with the actual satellite observation. Although primarily developed for soil moisture, it simultaneously solves for the VOD using an analytical solution by Meesters et al. (2005), utilizing the ratio between H- and V-polarized observations (Van der Schalie et al., 2017). The LPRM assumes that the soil and vegetation temperatures are equal, which may not be the case during the day due to uneven heating from solar radiation. VODCA v2 therefore uses only nighttime observations, which are assumed to be in thermal equilibrium (Owe et al., 2008). Scaling of the single-sensor VOD observations is achieved by means of cumulative distribution function (CDF) matching (Moesinger et al., 2020).

Table 1Input data for the merged-band VODCA v2 with the main sensor specifications: time periods used, local ascending equatorial crossing times (AECTs), and frequencies used. Table information is taken from Moesinger et al. (2020) and adapted for VODCA v2 CXKu.

Download Print Version | Download XLSX

The preprocessing of LPRM level-2 VOD data used in VODCA v2 follows the steps described in detail in Moesinger et al. (2020). These include projecting the data onto a $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ grid, using nearest-neighbour resampling, and selecting the closest nighttime value in a window of ±12 h for every 00:00 UTC (Zotta et al., 2022). Data are masked for radio-frequency interference (De Nijs et al., 2015), negative VOD retrievals, and temperatures. Differently from Moesinger et al. (2020), masking for low land surface temperature (LST<0 ^∘C), when the dielectric properties of water change drastically, is not based on Ka-band retrievals because these have high uncertainties over frozen land (Holmes et al., 2009). Instead, VODCA v2 uses the ERA5-Land (Muñoz-Sabater et al., 2021) soil temperature level-1 (stl1) data. To ensure that all observations taken under frozen conditions are masked, all observations with an associated surface soil temperature (stl1) below 3 ^∘C are masked (Zotta et al., 2022).

2.2 Reference datasets

3.1 VOD2GPP model

VODCA2GPP is based on the VOD-driven GPP estimation approach (the VOD2GPP model) introduced by Teubner et al. (2019, 2021). The VOD2GPP model describes the theoretical relationship between GPP and VOD. The biogeochemical basis of this model is the relationship between GPP and ecosystem net uptake of carbon (NPP) and autotrophic respiration (R_a) (Bonan, 2008):

where R_a can be again split into two terms: maintenance respiration and growth respiration (Bonan, 2008). The VOD2GPP model makes use of several VOD variables to represent the sum of NPP and R_a: the original VOD time series (VOD), which relates to maintenance respiration; temporal changes in VOD (Δ(VOD)), which relate to both growth respiration and NPP; and the temporal median of VOD (mdn(VOD)) derived from the complete time series, which serves as a proxy for vegetation density. Specifically, mdn(VOD) is incorporated to subtract larger structural vegetation components which makes the resulting model more closely related to biomass changes of smaller structural vegetation components such as leaves. It was shown that this increases model performance (Teubner et al., 2019). NPP is mostly represented by Δ(VOD), while R_a is represented by the original VOD signal and Δ(VOD). Thus, the VOD-based-only VOD2GPP model can be formulated as follows (Teubner et al., 2019):

where s() denotes the mapping function that maps the input variables to GPP. Δ(VOD) is derived for each pixel (x_i) by computing the difference between two consecutive VOD observations of the smoothed and 8 d aggregated VOD signal (Teubner et al., 2019):

The smoothing was performed in order to increase the robustness of the derivation and implemented using a Savitzky–Golay filter with a window size of 11 data points as suggested by Teubner et al. (2021).

Equation (2) represents a simplified model formulation connecting VOD to GPP but does not explicitly take into account the strong temperature dependence of autotrophic respiration (Wythers et al., 2013; Atkin et al., 2005; Tjoelker et al., 2001), which is mainly attributed to its maintenance part (Bonan, 2008; Ryan et al., 1997). Therefore, an improved formulation of the model was developed by considering the temperature dependence of maintenance respiration through a term representing the interaction between temperature (T2m) and VOD (Teubner et al., 2021):

The mapping and interaction functions were implemented using generalized additive models (GAMs). The usage of short time intervals (on the order of several days) for the computation of Δ(VOD) is crucial since it reduces the influence of larger vegetation components (e.g. stems) and makes the model more sensitive to changes in leaf biomass.

3.2 Generalized additive models

Generalized additive models (GAMs; Hastie and Tibshirani, 1990) are semi-parametric generalizations of linear models and combine properties of generalized linear models (GLMs) and additive models (Guisan et al., 2002). Link functions f() are trained and summed up for each predictor in order to relate the expected value of a response variable E(Y) to the explanatory variables x_i (Hastie and Tibshirani, 1990). The model can be written as

where β denotes a constant offset and n corresponds to the number of input predictor variables x_i. The link functions s_i() are implemented as smooth spline functions and allow the representation of non-monotonic and non-linear relationships, which give them a high degree of flexibility (Hastie and Tibshirani, 1990). Hence, the relationship between target and predictor variables does not require explicit a priori knowledge but can be estimated from the data themselves, which makes GAMs appropriate for the VOD2GPP model for which the exact relationship between VOD, air temperature, and GPP is difficult to determine (Teubner et al., 2019).

3.3 Preprocessing

The model input data (response variable, FLUXNET GPP; predictor variables, VODCA v2 CXKu, ERA5-Land T2m) were resampled from a daily to 8 d resolution using the 8 d means over the respective time period in order to reduce noise and computation times. This also means that the final VODCA2GPP represents the mean of daily GPP for an 8 d period with an estimate every 8 d. Since VODCA v2 CXKu already incorporates extensive quality flagging (e.g. for temperature), no additional data processing was necessary.

3.4 Model training and output

For each valid FLUXNET2015 in situ observation, the corresponding overlapping pixel values of VOD, Δ(VOD), mdn(VOD), and T2m were used to set up the GAM. Data from days with at least one invalid or missing observation were not considered for model training. While cross-validation was performed to evaluate the model (Sect. 3.5), all data were used for training of the final VODCA2GPP model. The GAM-based implementation of the VOD2GPP model is consistent with Teubner et al. (2021) and utilizes algorithms from the pyGAM Python package (Servén and Brummitt, 2018).

The trained VODCA2GPP model was applied to each pixel where all input variables (VOD, Δ(VOD), mdn(VOD), and T2m) were available. The result of this upscaling process is VODCA2GPP, which covers the period between January 1988 and July 2020. It has a spatial resolution of 0.25^∘, and its temporal resolution is 8 d. In rare cases (∼2.5 % of all data points) the machine learning model produced slightly negative results for GPP due to the extrapolation capacities of the trained GAM. As negative GPP is not possible, such estimates were set to zero.

3.5 Site-level evaluation and uncertainty assessment

The robustness of the model was evaluated based on a site-based cross-validation analysis during which the influence of the selection of available in situ stations on the GPP model was investigated. For the cross-validation 10 VODCA2GPP models were trained. Each of the 10 models was trained with 90 % of the available FLUXNET stations, while the remaining 10 % of the stations were retained for validation (Teubner et al., 2019). Every station was excluded exactly once, which is why this approach is classified as pseudo-random. Model performance was assessed at all sites that were omitted in the respective model run by computing root mean square errors (RMSEs), bias, and Pearson's r for different timescales (8 d, monthly, and yearly). In addition, the uncertainty in the VODCA2GPP model was then expressed through the minimum–maximum range as well as the standard deviation of the resulting 10 mean annual accumulated GPP estimates for each pixel. The standard deviation is also incorporated as an uncertainty map in the available dataset (layer name “Uncertainties”) to support users with an indicator for known uncertainties in VODCA2GPP.

3.6 Product evaluation and assessment

Mean annual, monthly, and 8 d GPP from VODCA2GPP, MODIS GPP, and FLUXCOM GPP was evaluated against GPP from FLUXNET. The error metrics used were RMSE, bias, and Pearson's r. Global spatial GPP patterns were compared against products by computing the mean annual GPP per pixel and the differences in mean annual GPP per pixel over the common observation period. Temporal agreements were tested by means of a Pearson correlation analysis for 8 d GPP. A correlation analysis of GPP anomalies was conducted for monthly GPP and also includes TRENDY-v7 GPP, which only provided monthly GPP data. Anomalies were derived by subtracting the long-term mean of the overlapping observation periods from monthly GPP estimates for each product.

Additionally, a trend analysis was conducted for all available GPP products in order to compare long-term changes in GPP. Trends in yearly median GPP were quantified using the Theil–Sen estimator (Theil, 1950; Sen, 1968), which calculates the slopes for each line between two points. The median of all computed slopes is then used for line fitting, making it insensitive to outliers and more robust than simple linear regression (Wilcox, 2010). Slopes were considered significant when the signs of the lower and upper 90 % confidence intervals were equal. For the trend analysis, yearly median GPP was used.

4.1 Spatiotemporal patterns in global GPP

The average annual GPP of VODCA2GPP exhibits spatial patterns similar to the remote-sensing-based MODIS and FLUXCOM products (Fig. 1). The agreement in annual GPP is high in northern latitudes (e.g. Europe, Russia, Canada), while there are relatively large differences in the Southern Hemisphere, especially in tropical and sub-tropical regions (Fig. 1d). The largest positive differences are found in the sub-tropics. Very arid regions (e.g. Australian deserts, Kalahari Desert) have low mean yearly productivity in all three datasets (Fig. 1a) but tend to be higher in VODCA2GPP compared to MODIS and FLUXCOM (Fig. 1b and c). The mean global total GPP as derived from VODCA2GPP amounts to $\mathrm{200}\pm \mathrm{2.2}\mathrm{Pg}\mathrm{C}{\mathrm{yr}}^{-\mathrm{1}}$. Comparison of the latitudinal distribution of FLUXNET stations shows that the closest agreement in yearly GPP is generally found in regions with a high density of FLUXNET in situ stations, while the largest discrepancies are found in regions with few or no FLUXNET stations.

Figure 1(a) Mean yearly aggregated GPP for the common observation period of the three products VODCA2GPP, MODIS, and FLUXCOM (2002–2016); (b, c) difference in mean annual GPP between VODCA2GPP and FLUXCOM and MODIS, respectively; (d) latitude plot of zonal means of mean annual GPP. The means were computed based on 8 d, 0.25^∘ sampling. The Min/Max area denotes the minimum–maximum latitudinal mean for the 10 model runs that were obtained with the site-based cross-validation. The dots represent the latitudinal location of the FLUXNET sites and their corresponding mean annual GPP. The brightness of the dots indicates the data availability for the respective FLUXNET station. Only pixels and observation dates that are available in all three datasets were used for these plots.

Similarly, uncertainties tend to be smaller in latitudes with a high density of FLUXNET stations (Fig. 1d). The lowest spread in the 10 models (i.e. the lowest uncertainty) is found north of 20^∘ N where the majority of FLUXNET GPP stations are also located. The Southern Hemisphere, where only a few in situ stations are located, generally exhibits a larger spread (higher uncertainty), indicating a considerable sensitivity of the model to the choice of stations. This emphasizes the need for a well-distributed network of in situ flux towers across all biomes. The uncertainty map (Fig. 2b) shows that arid regions (e.g. Sahara, Australian deserts, Arabian Peninsula) as well as various mountainous regions (e.g. Carpathians, Alps, Rocky Mountains, Andes) have the highest model uncertainties. Moderate to high model uncertainties are also exhibited for the tropics. Furthermore, significant uncertainties in VODCA2GPP are found in parts of eastern and western Siberia's boreal forests as well as in parts of southern China.

Figure 2(a) Mean annual GPP as derived from VODCA2GPP for the period 1988–2019. (b) Standard deviation of mean annual GPP (1988–2019) as obtained by the uncertainty analysis based on site-level cross-validation.

VODCA2GPP shows good temporal agreement with MODIS and FLUXCOM for most regions outside the tropics (Fig. 3). Pearson's r is the highest in regions with distinct interannual variability such as sub-arctic, temperate, and semi-arid regions and the lowest for dense tropical forests where even negative correlations occur. The median Pearson's r values between VODCA2GPP and the reference datasets MODIS and FLUXCOM are 0.77 and 0.75, respectively.

Figure 3Pearson's r between VODCA2GPP and MODIS GPP (a) and VODCA2GPP and FLUXCOM GPP (b). The correlations are based on the common observation period between 2002 and 2016 with a 0.25^∘ spatial and 8 d temporal resolution.

4.2 Site-level evaluation

VODCA2GPP's tendency towards a positive bias with respect to MODIS and FLUXCOM products is not mirrored in the comparison against FLUXNET GPP (Figs. 4 and A1). The bias with respect to FLUXNET site data is substantially smaller for VODCA2GPP than for MODIS and FLUXCOM. The RMSE and Pearson's r values of VODCA2GPP are slightly higher and lower, respectively, than for MODIS and FLUXCOM and are of the same magnitude for 8 d values and mean annual GPP. All three datasets underestimate productivity at high GPP values. A land-cover-based analysis (Fig. A3) shows that discrepancies in annual VODCA2GPP mostly occur in (semi-)arid environments (e.g. savannas, open shrublands, grasslands). VODCA2GPP performs best in temperate environments (e.g. wetlands, evergreen broadleaf forest, croplands). Wetlands and evergreen broadleaf forests exhibit the best performance for all products, while all three datasets underperform in open shrublands and deciduous broadleaf forest.

Figure 4GPP from FLUXNET plotted against GPP from VODCA2GPP, MODIS, and FLUXCOM for the period 2002–2016 with 8 d sampling.

Download

The site-based spatial cross-validation also exhibits only a small (negative) bias of VODCA2GPP (Fig. 5) for monthly and 8 d GPP values, while the bias for annual variations is positive and slightly higher. A high median Pearson's r for 8 d and monthly values indicates good model performance for seasonal variations. It is to be noted that there are only eight significant Pearson's r values for yearly sampling, which decreases the expressiveness of this value. This is explicable with the short observation period of most FLUXNET sites which might not exhibit interannual variability. The RMSE decreases with increasing observation length scales.

Figure 5Site-based cross-validation for 8 d, monthly, and yearly sampling of GPP from VODCA2GPP and FLUXNET. RMSE, bias, and Pearson's r were computed at each of the 10 % of FLUXNET sites that were omitted during the respective training run. Non-significant Pearson correlations (p value<0.1) were ignored. The boxplots for Pearson's r are based on the 71 (8 d), 66 (monthly), and 8 (yearly sampling) significant values for Pearson's r values. The whiskers of the boxplots extend to the 5th and 95th percentiles.

Download

4.3 Anomaly patterns in space and time

In terms of anomalies from the long-term climatology, VODCA2GPP shows good correlation with MODIS and TRENDY-v7 and weaker correlation with FLUXCOM (Table 2, Fig. 6). TRENDY-v7 correlates similarly well with VODCA2GPP, and MODIS and shows worse correspondence with FLUXCOM. The highest correlation is found between the two products based on optical remote sensing: MODIS and FLUXCOM.

Table 2Pearson's r correlation matrix for mean global monthly GPP anomalies between 2002 and 2016.​​​​​​​

Download Print Version | Download XLSX

The temporal evolution and spatial distribution of the anomalies exhibit similar patterns (Fig. A4). Several extreme events are captured in VODCA2GPP and in at least one of the other GPP datasets. An example of such GPP extremes is the strong positive anomalies between 2010 and 2011 at around 25^∘ S which were mainly caused by record-breaking rainfalls in Australia (Wardle et al., 2013). These positive anomalies are clearly visible in all examined GPP products apart from FLUXCOM (Fig. A4). Other GPP extremes that are noticeable in all products apart from FLUXCOM are the extremely low GPP in 2002, 2003, and early 2005 around 20^∘ S (Fig. A4). Both anomalies can be explained by extreme drought events that occurred in these years (Bureau of Meteorology, 2002, 2003, 2005; Horridge et al., 2005), which are associated with El Niño events (Taschetto and England, 2009). Also, the distinct drop in GPP in 2015 and 2016 in similar latitudes is likely linked to El Niño-related drought events (Malhi et al., 2018; Zhai et al., 2016). Generally, extreme events in VODCA2GPP are more pronounced than in the other datasets.

Figure 6Time series of mean global monthly GPP anomalies.

Download

4.4 Global GPP trends

Trends in global annual median GPP for the overlapping period between 2003 and 2015 are similar for VODCA2GPP, TRENDY-v7, and MODIS, and all show significant positive trends (Table 3 and Fig. 7). FLUXCOM does not exhibit a significant trend. The spatial distribution of GPP trends for the period 2003–2015 (Fig. 8) exhibits many similarities between all analysed products. Large patterns of strong positive trends are, for example, found in eastern parts of Siberia and China as well as in India and North America. Patterns of negative trends are found north of the Caspian Sea in all datasets. The remote-sensing-based products exhibit distinct patterns of declining GPP in central Siberia and significantly increasing GPP in western Russia. Generally, the trends of VODCA2GPP match better with MODIS and TRENDY-v7 than with FLUXCOM. While there are many similar patterns in the Northern Hemisphere, trends in the Southern Hemisphere do not match well or are even contradictory. Especially in the tropics, hardly any similarities are apparent. Note that the analysed time period is short and may be impacted by individual extreme events.

Table 3Theil–Sen trends in global yearly median GPP. Same signs of the upper and lower 90 % confidence interval indicate significant trends. The analysed periods are 2003–2015, which corresponds to the fully overlapping periods for all datasets (for the period 2002–2016 there were some data gaps in the MODIS and FLUXCOM data used at the very beginning/very end; since these data gaps could potentially impact the slope estimation, the slightly shorter period 2003–2015 was used); 1988–2016, which corresponds to the fully overlapping period of VODCA2GPP and TRENDY-v7; and 1988–2019, which corresponds to all available complete years of VODCA2GPP data.

Download Print Version | Download XLSX

Figure 7Time series of yearly median GPP with the regression lines as obtained by the Theil–Sen estimator. Areas around the regression lines indicate the 90 % confidence intervals.

Download

Figure 8Global map of trends in yearly median GPP for the period 2003–2015 for all analysed datasets. White indicates non-significant trends.

For the full time period (1988–2019), VODCA2GPP increases slightly on a global scale (Table 3), but this cannot be classified as significant due to contradictory upper and lower confidence intervals. The same is true for the slightly shorter period between 1988 and 2016 during which TRENDY-v7 detects a small significant positive trend on a global scale. The spatial distribution of long-term (1988–2019; Fig. 9) trends in VODCA2GPP is similar to that of the shorter period (2003–2015), but in general, long-term VODCA2GPP trends are less pronounced. The comparison of the fully overlapping period between VODCA2GPP and TRENDY-v7 (1988–2016; Fig. A5) shows that TRENDY-v7 exhibits weak but consistent positive trends for practically all biomes, while VODCA2GPP trends are spatially differing and for some regions even opposite in sign to the trends in TRENDY-v7.

Figure 9Global map of yearly median GPP trends for the period 1988–2019 for VODCA2GPP. White indicates non-significant trends.

A comprehensive comparison with in situ GPP trends is not possible because most FLUXNET time series are too short to derive reliable trends. However, trends that could be derived for a few stations with a long time series (Fig. A6) also suggest increasing GPP. The in situ analysis indicates that there is a comparatively good correspondence between VODCA2GPP and FLUXNET GPP trends. Together with the strong similarities between VODCA2GPP and MODIS and TRENDY GPP, this suggests that VODCA2GPP can provide a valuable contribution to the analysis of global GPP trends.

5.1 Uncertainties in the VODCA2GPP model

The results from the uncertainty analysis and the comparison with in situ GPP show that VODCA2GPP estimates can be viewed as reliable across most biomes. However, significant uncertainties were exhibited in some areas with extreme climatic or topographic conditions (e.g. deserts and mountain ranges). Also, parts of eastern and western Siberia and parts of southern China show a relatively large spread in predictions. The observed uncertainty patterns in Siberia might be associated with topography, land cover, generally lower data availability (due to frozen masking in VODCA), and a lack of FLUXNET stations. Complex topography is presumably also the main driver for uncertainties in southern China. The uncertainty analysis suggests that VODCA2GPP estimates tend to be too high in these regions and thus should be interpreted with caution. Furthermore, moderate uncertainties were also found for the tropics, which is likely due to the extremely low in situ data availability and higher absolute GPP than in mid-latitudes.

The comparison with in situ GPP shows clear differences in performance of the VODCA2GPP model across different biomes. High performance is achieved in densely vegetated biomes, while performance decreases in arid and less vegetated regions. A reason for the weaker performance in areas with less water availability might be adapted water regulation strategies of plants. Plants in drought-prone regions often reduce transpiration by limiting stomatal conductance to maintain a constant water potential even in times of extreme water scarcity (Sade et al., 2012). Since VOD is largely driven by the vegetation's water content, this isohydric behaviour of vegetation could at least partly explain relatively high VOD and consequently also overestimated GPP in those regions (Teubner et al., 2021).

Also, the observation bias which is introduced by unevenly distributed FLUXNET sites decreases the model's robustness. GPP is measured in situ only at a few locations, and these stations are mostly located in temperate regions (e.g. Europe and North America) while semi-arid and tropical forest regions are underrepresented in the training data.

A comparison of uncertainties between VODCA2GPP and GPP based on optical remote sensing (Xie et al., 2021) shows that in both cases topographic complexity decreases the reliability. Furthermore, the reliability of GPP estimates based on optical remote sensing is highly dependent on weather and illumination conditions. Clouds often contaminate or prevent the observations, which is presumably the main reason why the largest uncertainties for FLUXCOM and MODIS are found in the wet tropics where GPP is known to be underestimated (de Almeida et al., 2018; Jung et al., 2020). In contrast, VODCA2GPP shows good skill for densely vegetated areas, including broadleaf evergreen forests. On the other hand, the relatively high uncertainties in VODCA2GPP in water-limited regions have not been reported for FLUXCOM or MODIS GPP, indicating that these are VOD-specific and presumably caused by the above-mentioned isohydric behaviour of plants in arid regions.

5.2 Limitations in VODCA and their impact on VODCA2GPP

Certain limitations in the VODCA v1 product exist, as outlined in Moesinger et al. (2020), which are partly also evident in VODCA v2 (Zotta et al., 2022) and thus propagate to VODCA2GPP. A known issue of VODCA v2 is caused by an observation gap between October 2011 and July 2012 for AMSR-E and AMSR2 (Table 1), which prevents a direct bias removal between the sensors. However, scaling between the sensors is achieved by using TMI observations north/south of 35^∘ N/35^∘ S for the X- and Ku-band. Beyond these latitudes for the X- and Ku-band and globally for the C-band, AMSR-E data were matched directly to AMSR2 using the last 2 years of AMSR-E and first 2 years of AMSR2 as the reference period, under the assumption that trends between 2010–2014 are negligible (Moesinger et al., 2020). The result is that AMSR2 observations exhibit a slight positive bias in parts of North America, which is also evident in a spatial break in VODCA v1 X- and Ku-band trends (Moesinger et al., 2020). Although the impact of this procedure on VODCA2GPP trends is small and spatially limited, users are advised to keep the potential bias in mind when analysing VODCA2GPP data after 2012 for latitudes north/south of 35^∘ N/35^∘ S. Other limitations in VODCA concern the mixing of observations that were retrieved at different geometries (e.g. incidence angles) or observation times (Moesinger et al., 2020) and the data loss in certain regions, mostly in the Himalayas, which is caused by failure of the CDF-matching method due to insufficient input data (Moesinger et al., 2020). These issues, however, only have a small or spatially very limited influence on the final VODCA2GPP product. Furthermore, VOD retrievals exhibit a tendency towards saturation in regions with very dense vegetation, making it less likely to distinguish variability. A slight tendency towards saturation was also observed for VODCA2GPP, but the land-cover-based analysis exhibited a very high agreement between VODCA2GPP and in situ GPP, indicating high reliability of VODCA2GPP over densely vegetated regions.

Another limitation of VOD products in general, and thus also of VODCA2GPP, is the limited spatial resolution (0.25^∘). The lower spatial granularity achieved from passive microwave remote sensing is presumably another reason for the slightly weaker performance of VODCA2GPP in comparison with GPP derived from optical remote sensing since the VODCA grid cell might not always be well represented by the in situ measurements. Furthermore, the lower spatial resolution of VODCA2GPP is disadvantageous for the analysis of local GPP as small-scale variations in GPP might be hidden in VODCA2GPP. However, due to its long-term availability and generally high reliability, VODCA2GPP can still serve as a valuable source of data for various other applications (Sect. 5.4), especially concerning long-term climate-related studies and climate model evaluation.

5.3 Observed bias between VODCA2GPP and other remote-sensing-based GPP datasets

There is only a very small bias when comparing VODCA2GPP with eddy-covariance measurements from FLUXNET, but relatively large discrepancies in absolute GPP exist between VODCA2GPP and other remote-sensing-based products and to a lesser extent with process-based TRENDY-v7 models. In tropical regions, the positive bias between VODCA2GPP and MODIS and FLUXCOM GPP can be partly explained by a reported and observed tendency of FLUXCOM and MODIS to underestimate GPP in these regions (Turner et al., 2006; Wang et al., 2017; Figs. 1 and A3). Outside the tropics, discrepancies in absolute GPP among products might be caused by the assumed overestimation of VODCA2GPP in winter months (i.e. in times with very little or no primary productivity). This overestimation is explained by water content in vegetation that is also present in these dormant periods. The sensitivity of microwaves to this water content results in non-zero VOD and, consequently, non-zero GPP (Teubner et al., 2021). This effect is similar to the isohydricity effect discussed in Sect. 5.1, which is an explanation for overestimation of VODCA2GPP in arid regions.

Another potential explanation for the positive bias of VODCA2GPP compared to MODIS and FLUXCOM is the presence of surface water and its impact on VOD retrievals. The presence of surface water is known to decrease the brightness temperature of the earth's surface and thus significantly decreases VOD retrievals (Bousquet et al., 2021). The impact of surface water contamination is evident in VODCA pixels that partly contain waterbodies (e.g. lakes, rivers). These pixels exhibit systematically lower values than neighbouring pixels without waterbodies. On the one hand, this leads to underestimation in the VODCA2GPP model in pixels containing surface water. On the other hand, it also has an effect on model training. This effect is caused by FLUXNET stations located close to waterbodies, which hardly impacts in situ GPP retrievals but does cause erroneous VOD retrievals at the 0.25^∘ pixel scale. As a result, underestimated VOD is trained against unaffected in situ GPP, which causes a slight but systematic global overestimation. A potential solution would be the masking of water-contaminated VOD. However, due to the constraints with temperature in the interaction term (Eq. 3), this would strongly reduce the data available for training, which would potentially decrease the robustness of the VODCA2GPP model if the number of stations is not increased.

A general issue in the upscaling of GPP is the low availability of in situ GPP, which not only is problematic in model training but also hampers a fair evaluation and validation at the global scale. The remote-sensing-based reference products, FLUXCOM and MODIS, however, are also trained and calibrated using in situ GPP observations (Jung et al., 2020; Running et al., 1999) and can therefore not be viewed as fully independent from VODCA2GPP (Teubner et al., 2021). In contrast to observation-based GPP products, estimates from the TRENDY ensemble can be considered largely independent from VODCA2GPP.

5.4 Potential applications of VODCA2GPP

The validation results show that VODCA2GPP generally exhibits a high consistency with in situ GPP observations and global state-of-the-art GPP products, indicating that VODCA2GPP can be used complementarily to current global GPP products. For the analysis of global as well as regional GPP anomalies, VODCA2GPP can provide valuable insights that might be hidden in other observational products due to the fundamentally differing observation methods and the associated limitations related to saturation effects, cloud cover, and other atmospheric effects such as water vapour content or aerosols (Xiao et al., 2019).

Also, for the monitoring of global GPP trends, VODCA2GPP has the potential to serve as an independent and reliable source of data. The long-term trend analysis suggests that the majority of biomes have increased their primary productivity since 1988. There are several potential drivers for long-term increases in GPP, the most important ones being global warming, land-use changes, and elevated CO₂ concentration in the atmosphere (Piao et al., 2019). The observed long-term trends in GPP across the different products support the theory of elevated atmospheric CO₂ leading to an increased uptake of CO₂ (Haverd et al., 2020; Walker et al., 2020; Campbell et al., 2017; Schimel et al., 2015). The absence of trends in FLUXCOM does not contradict but rather supports this theory as FLUXCOM does not account for CO₂ fertilization effects (Jung et al., 2020). Due to the shortness of most in situ GPP time series, it is, however, difficult to draw final conclusions on the existence of, magnitude of, and reasons for long-term variations in GPP. Therefore, the global influence of atmospheric CO₂ on vegetation productivity remains uncertain, but VODCA2GPP allows us to gain new perspectives on long-term GPP trends and might help to identify and quantify driving factors for increasing long-term primary productivity.

Furthermore, VODCA2GPP can be used as a largely independent source of data for the intercomparison and validation of other existing or newly developed global GPP datasets and models. Currently, a multitude of global GPP products exists showing large inconsistencies among products (Zhang and Aizhong, 2021). Similarly to other global GPP datasets, VODCA2GPP cannot be seen as a true reference, but using it as an additional reference might help to acquire a more comprehensive picture of the performance of other datasets, especially in the context of long- and short-term variability in GPP.