Background & Summary

The eddy-covariance flux stations provide reliable ecosystem-scale measurements of the carbon and energy fluxes at a high temporal resolution1. They have become crucial tools to generate observation datasets to verify and benchmark the Earth surface models2,3. In particular, it is possible to construct a carbon-water flux simulation model from the station-scale to the regional- or global-scale by means of a large-scale eddy covariance4 measurement network (e.g. Fluxnet, AmeriFlux and ChinaFlux). However, the existing flux stations are sparsely and unevenly distributed and yield rather discontinuous observation data1. This restricts studies on the carbon-water fluxes at a large-scale3, for example in Eurasia, where a strong spatial heterogeneity is exhibited on complex terrains. The meteorological stations, in contrast, are densely spread around the world with long-term continuous observation data5, which could have great potential to mine the more extensive carbon-water flux information, particularly combined with machine learning (ML) and remote sensing (RS). This could greatly offset the limitations of the flux station-based observations.

Machine learning is increasingly used to extract the patterns and insights from big geospatial data6. Many studies have focused on the comparative evaluation of different ML algorithms and have found the accuracy performance of the same algorithm varies in different research contexts7,8,9,10. The data-driven ML algorithms are similar to the encapsulated complex empirical algorithms, which demonstrate a high simulation accuracy3,11. But the ML algorithms are still influenced by the quality, processing methods, and spatio-temporal representativeness of the data12,13,14. Compared with the process-based land surface or ecosystem models, the ML has a higher carbon-water flux simulation accuracy1,6. However, when transferred to other site or regional or spatial (grid) scales, the applicability of both the ML models and process models need to be evaluated due to the distinct spatio-temporal heterogeneity. That is to say, there is no guarantee that these models are applicable to all sites, grids or regions. If this evaluation of the model applicability is not considered, the simulation results will generate new uncertainties. This issue has become a major problem affecting the simulation accuracy of the carbon-water fluxes at different scales.

In this study, ML with flux observations was used to build carbon-water flux simulation models (random forest model, RFM) to simulate the carbon-water fluxes of the meteorological stations in Eurasia. We proposed a framework to evaluate the applicability of the model transfer and to build a bridge from the flux stations to the meteorological stations. We used this framework to generate four carbon-water flux datasets for the Eurasian meteorological stations. Due to the high precision, these datasets could be regarded as quasi-observational at the site level, which might be used to assess the simulation accuracy of the regional- or global-scale ecosystem carbon-water fluxes based on the ecosystem or land surface or remote sensing or atmospheric inversion models. Our study can, therefore, benefit terrestrial water management and carbon dynamic assessments.

Methods

The RFM was constructed based on the Eurasian flux stations. We built a total of 3,600 RFMs at site scale in accordance with the classification of the flux stations. The simulation accuracy of these RFMs at each flux station in the test set was validated by the spatial cross-validation, thus generating thousands of determination coefficients (R2) at test stations. According to the third law of geography15, the factors (variables) similarity between the test station and the training set of the RFM determines the similarity between their fluxes, that is, the R2 of the RFM at the test station are determined. The similarity between the datasets composed of the same factors could be characterized by the Euclidean distance. Based on the R2 and Euclidean distance, the R2 simulation model (RSM) was built by using multiple linear regression (MLR) to evaluate the applicability of RFM on meteorological stations. So that the RFMs can be reasonably transferred to meteorological stations to simulate the carbon-water fluxes. Figure 1 shows the detailed flowchart of the data processing, RFM construction and RFM transfer to the meteorological stations.

Fig. 1
figure 1

Research framework. R2, determination coefficient; MLR, multiple linear regression. N/A (not applicable) indicates that the RFM could not be transferred to the specific meteorological stations. NEE-RS: net ecosystem carbon dioxide exchange (NEE) based on the RFM with remote sensing (RS); WF-RS: water flux (WF) based on the RFM with RS; these explain the fact that the RS variables were used in the RFM construction. NEE-WRS: NEE based on the RFM without RS; WF-WRS: WF based on the RFM without RS; these demonstrate that the RS variables were not applied in the RFM construction. RS variables include the fraction of the photosynthetically active radiation, enhanced vegetation index, land surface water index and surface reflectance for the Moderate Resolution Imaging Spectroradiometer bands 1–7.

Data processing

We selected 156 flux stations in Eurasia from five different landscape types (Fig. 2a), as detailed in the flux station information file16. For the flux stations from the National Tibetan Plateau Data Center (NTPDC)17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133 and European Fluxes Database Cluster (EFDC)134,135,136,137 (http://www.europe-fluxdata.eu/home), the flux data from one hour before (and after) rainfall were excluded. The data collected at 10-min or 30-min intervals were interpolated using the marginal distribution sampling (MDS) method in REddyProc138. All final data were converted into daily data. For the flux stations from FLUXNET139,140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159,160,161,162, the data were extracted with quality control values ≥ 0.8 for the net ecosystem carbon dioxide exchange (NEE) and latent heat fluxes (LE). The water fluxes (WF) were converted from LE (W/m2) with a conversion factor of 0.408 × 10−6 163,164,165. For the selected 6856 meteorological stations from the Global Surface Summary of the Day dataset in the National Centers for Environmental Information (https://www.ncei.noaa.gov/metadata/geoportal/rest/metadata/item/gov.noaa.ncdc%3AC00516/html#), the vapour pressure deficit variable was calculated using the air temperature and dew point temperature. The downward shortwave radiation (DSR) of meteorological stations for 2002–2020 and 1983–2018 were extracted from the GLASS dataset166,167 and the dataset of high-resolution global surface solar radiation168,169 from the NTPDC, respectively. For the remote sensing (RS) variables (including the fraction of the photosynthetically active radiation extracted from the MCD15A3H data170, enhanced vegetation index, land surface water index and surface reflectance for the Moderate Resolution Imaging Spectroradiometer bands 1–7 extracted from the MOD09GA data171), a linear interpolation was carried out for the missing data with continuous missing days < 8165,172. Terrain and soil variables were extracted from the MERIT DEM data173 and the HWSD data174, respectively. High quality RS variables, terrain variables and soil variables averaged over a 500-meter spatial extent centered on the station were integrated into the flux stations and meteorological stations (Table S1).

Fig. 2
figure 2

Study area and the accuracy of carbon-water flux simulation models (random forest model, RFM). (a) Distribution of the 156 Eurasian flux stations covering five main landscape types. (b) The accuracy assessments of the RFM based on the Eurasian flux stations in the framework of the 10-time 10-fold cross-validation. The figure shows the percentage of the RFMs with R2 ≥ 0.5 tested in the test sets for nine categories.

Due to the significant spatial heterogeneity of the earth’s surface, the flux stations and meteorological stations were divided into nine categories according to the following four strategies. The first is based on the International Geosphere-Biosphere Programme classification from the MCD12Q1 data175, including Wetland (i.e. permanent wetlands), Cropland (i.e. croplands and cropland/natural vegetation mosaics), Grassland (i.e. grasslands, savannas, woody savannas) and Forest (i.e. evergreen needleleaf forests, evergreen broadleaf forests, deciduous needleleaf forests, deciduous broadleaf forests, and mixed forests). The second is based on the continents, e.g. Asia and Europe. The third is the arid and non-arid regions classified by the dryland dataset176. This dataset identified areas with an aridity index (AI) less than 0.65 as drylands, which were described as arid in this study, and the remaining areas (i.e. AI ≥ 0.65) were classified as non-arid regions. The fourth comprises entire Eurasia, that is, overall. We used the datasets of nine categories from flux stations as input to RFM, and the detail is presented in the division of flux stations file177.

RFM construction

In this work, the random forest method6,178 was used to construct the RFM using the scikit-learn library in Python 3.7.6 (https://pypi.org/project/scikit-learn/1.0.1/). The RFM is a combination model based on independent regression trees179, of which the predictions were made by averaging the results across all regression trees. The random-search optimizer180 was applied to identify the optimal hyperparameter settings for the RFM (Table S2). In addition, the most important step in the RFM construction is the k-fold cross-validation (CV)181. Suppose the flux station dataset is composed of data from n flux stations and it might generate k sub datasets (D1, D2, …, Dk) if it is equally divided into k parts. Here, each sub dataset is a test set, which is composed of data from the m = n/k flux stations, while the remaining sub datasets constitute training sets, which are composed of data from the (n-m) flux stations. Thus, each training set could be used to establish the RFM for the flux simulation; consequently, a total of k models (RFM1, RFM2, …, RFMk) were built. Each RFM is tested and verified at each flux station in the test set. Furthermore, the R2 is calculated for each case. Hence, the R2 amount that could be generated after each RFM training and test is given by k × m = k × n/k = n. The flux station dataset (n) requires multiple (p) repetitions of the k-fold division to avoid the contingency of the station division. In this way, k × p RFMs could be constructed and the number of R2 that could be generated is n × p.

According to the above-mentioned principles, we used the 10-time (split) 10-fold CV (i.e. p = 10 and k = 10) to set up 100 RFMs for each of the nine categories under two scenarios: with and without RS variables in the model building (i.e. RFM-RS and RFM-WRS), respectively. That is, we constructed 900 RFM-RS models to simulate NEE and WF, respectively. And we also constructed 900 RFM-WRS models to simulate NEE and WF, respectively. Each RFM was validated at each flux station in the corresponding test set, and the R2 was generated to represent the validation accuracy of the RFM. The R2 also represents the applicability of these models in the test flux stations. The higher the R2 of the model on the test flux station, the more suitable the model is for the specific flux station, that is, the more similar the data characteristics of the training set for building the model are to the data characteristics of the flux station according to the third law of geography.

RFM transfer to the meteorological stations

We screened available RFMs for meteorological stations by using the RSM, which was used to evaluate the RFM applicability on the meteorological stations. The framework for the evaluation was designed (Fig. 1) as follows:

① Euclidean distances of the influencing factors between test flux stations and training sets of RFMs.

The R2 of each RFM is determined by the similarity of the influencing factors between the training set and the test set181. This could be characterized by the Euclidean distance182.

For a specific factor affecting the RFM, the Euclidean distance (ds) between a flux station in the test set (test station) and a flux station in the training set (training station) is expressed as:

$$ds(x,y)=\sqrt{{({x}_{1}-{y}_{1})}^{2}+{({x}_{2}-{y}_{2})}^{2}+\cdots +{({x}_{t}-{y}_{t})}^{2}}=\sqrt{{\sum }_{(i=1)}^{t}{({x}_{i}-{y}_{i})}^{2}},t\le 365$$
(1)

where x represents one of the factors influencing the carbon-water fluxes in the training station (Table S1), y denotes the corresponding influencing factor in the test station, and t is the sample size of the factor. The data of the training station and test station were averaged on the same day (day of the year, DOY) for multiple years, respectively. Then, the two stations could be matched day by day based on DOY to ensure that these have the same daily data sample size.

For the factor j influencing the RFM, the final Euclidean distance (d) is the average of all n-m ds between the test station and each of the n-m training stations, which is the Euclidean distance (dj) of the factor j between this test station and the training set of the RFM (Eq. 2). In the same way, the Euclidean distances of all influencing factors are produced and denoted as d1, …, dw-1, dw. In this way, the R2 of the RFM tested in the test station from the test set and the Euclidean distances d1, …, dw-1, dw between this test station and the training set for building the RFM constitute a complete data sample (Fig. S3). Similarly, all test stations could generate such samples, which constitute a dataset with a quantity equal to 10 × n × 9. Samples of the same category are integrated into one dataset, that is, nine datasets produced under nine categories, i.e. Dataset 1 in Fig. 1.

$$d(x,y)=\frac{d{s}_{1}(x,y)+d{s}_{2}(x,y)+\cdots +d{s}_{(n-m)}(x,y)}{(n-m)}=\frac{{\sum }_{(i=1)}^{(n-m)}d{s}_{i}(x,y)}{(n-m)}$$
(2)

where n represents the number of the flux stations of each category and m illustrates the number of the flux stations in each test set.

② Construction of the RSM.

Based on Dataset 1, the RSM is constructed using MLR183 for nine categories under NEE and WF scenarios, of which each one is expressed as:

$${R}^{2}={a}_{0}+{a}_{1}{d}_{1}+{a}_{2}{d}_{2}+\cdots +{a}_{(w-1)}{d}_{(w-1)}+{a}_{w}{d}_{w}$$
(3)

where a0, a1, a2, …, aw-1, aw are regression coefficients and d1, d2, …, dw-1, dw indicate the Euclidean distances of the factors influencing the carbon-water fluxes between the test station and the training set.

③ Euclidean distances of the influencing factors between meteorological stations and training sets of RFMs.

The same processes of ① are applied to the meteorological stations so as to calculate the Euclidean distance for the influencing factors between each RFM training set and meteorological station for each category under two scenarios, which yields a large dataset, i.e. Dataset 2 in Fig. 1.

④ Prediction of the R2 of the RFM transfer to the meteorological stations.

Before a RFM is transferred to a specific meteorological station, the RSM could predict the R2 value on the station using Dataset 2 in Fig. 1. Only if predicted R2 ≥ 0.5, its corresponding RFM might be transferred to the corresponding meteorological stations. Otherwise, the RFM was assumed to be not applicable to the meteorological station. If there was more than one RFM applicable to a meteorological station, the RFM corresponding to the maximum predicted R2 was screened as the model that could be linked to the meteorological station. Not all meteorological stations did have a RFM which is applicable to the meteorological station.

⑤ Carbon-water flux simulation of the meteorological stations.

For the meteorological stations in Eurasia that could be linked with an applicable RFM, the corresponding RFM can be used to simulate the daily carbon-water fluxes and to build high-precision carbon-water flux datasets of the Eurasian meteorological stations to analyze the carbon-water dynamics. These datasets184 consist of two essential building blocks: (i) datasets related to remote sensing, including the net ecosystem carbon dioxide exchange (NEE-RS) and water fluxes (WF-RS) simulated by the RFM-RS; (ii) the net ecosystem carbon dioxide exchange (NEE-WRS) and water fluxes (WF-WRS) simulated by the RFM-WRS.

Data Records

Our carbon-water flux datasets184 are available at figshare (https://doi.org/10.6084/m9.figshare.21347721.v3). The data record contains two daily carbon dioxide flux datasets (NEE-RS and NEE-WRS) and two daily water flux datasets (WF-RS and WF-WRS) of the Eurasian meteorological stations. The coverage period of the NEE-RS and WF-RS has been recorded from 2002 to 2020 and the one of NEE-WRS and WF-WRS from 1983 to 2018. The data of each meteorological station was deposited separately in the CSV file format under the dataset folders. The file name indicates the identification number of the meteorological station in the meteorological station information file185 (https://doi.org/10.6084/m9.figshare.23695920.v2). The list of flux stations used in this study and the details of flux station division used for the RFM construction are shown in the flux station information file16 (https://doi.org/10.6084/m9.figshare.23899701.v1) and the division of flux stations file177 (https://doi.org/10.6084/m9.figshare.23695980.v2) stored at figshare, respectively. In addition, the details of the RSM construction are presented in the RSMs information file186 (https://doi.org/10.6084/m9.figshare.23899785.v1) deposited at figshare. The file specific fields are as follows:

Carbon-water flux datasets file (.csv)

  1. (1)

    id: Identification of the meteorological station.

  2. (2)

    lon: Longitude of the meteorological station.

  3. (3)

    lat: Latitude of the meteorological station.

  4. (4)

    year: Year of the data record.

  5. (5)

    month: Month of the data record.

  6. (6)

    day: Day of the data record.

  7. (7)

    doy: Day of the year.

  8. (8)

    NEE: Net ecosystem carbon dioxide exchange (g C m−2 d−1).

  9. (9)

    WF: Water fluxes (mm d−1).

Meteorological station information file (.xlsx)

  1. (1)

    Identification of meteorological station.

  2. (2)

    Station name: Name of the meteorological station.

  3. (3)

    Longitude: Longitude of the meteorological station.

  4. (4)

    Latitude: Latitude of the meteorological station.

  5. (5)

    Elevation: Elevation of the meteorological station (m).

  6. (6)

    Continent: Continent of the meteorological station.

  7. (7)

    Drought situation: Drought situation of the meteorological station.

  8. (8)

    Landscape: Landscape of the meteorological station.

  9. (9)

    Data source: Data source of the meteorological station.

  10. (10)

    Classification of simulated NEE-RS: Accuracy classification of NEE-RS for meteorological stations (1: low quality, R2 < 0.5; 2: moderate quality, 0.5 ≤ R2 < 0.7; 3: high quality, R2 ≥ 0.7).

  11. (11)

    Classification of simulated NEE-WRS: Accuracy classification of NEE-WRS for meteorological stations (1: low quality, R2 < 0.5; 2: moderate quality, 0.5 ≤ R2 < 0.7; 3: high quality, R2 ≥ 0.7).

  12. (12)

    Classification of simulated WF-RS: Accuracy classification of WF-RS for meteorological stations (1: low quality, R2 < 0.5; 2: moderate quality, 0.5 ≤ R2 < 0.7; 3: high quality, R2 ≥ 0.7).

  13. (13)

    Classification of simulated WF-WRS: Accuracy classification of WF-WRS for meteorological stations (1: low quality, R2 < 0.5; 2: moderate quality, 0.5 ≤ R2 < 0.7; 3: high quality, R2 ≥ 0.7).

Flux station information file (.xlsx)

  1. (1)

    Identification of flux stations.

  2. (2)

    Flux station: Name of the flux station.

  3. (3)

    Longitude: Longitude of the flux station.

  4. (4)

    Latitude: Latitude of the flux station.

  5. (5)

    Elevation: Elevation of the flux station (m).

  6. (6)

    Continent: Continent of the flux station.

  7. (7)

    Drought situation: Drought situation of the flux station.

  8. (8)

    Landscape: Landscape of the flux station.

  9. (9)

    Study period: Study period of the flux station used in this study.

  10. (10)

    Data source: Data source of the flux station.

Division of flux stations file (.xlsx)

  1. (1)

    Categories: Category of the flux station.

  2. (2)

    Number of flux station: Number of flux stations under each category.

  3. (3)

    Split: Identification of divisions for 10-fold cross-validation on flux stations.

  4. (4)

    Fold: Identification of folds for cross-validation on flux stations.

  5. (5)

    Identification of flux station: List of identifications for flux stations under each fold.

RSMs information file (.xlsx)

  1. (1)

    Models: Name of the RSM.

  2. (2)

    Categories: Category of the RSM.

  3. (3)

    N: Number of samples used by the RSM.

  4. (4)

    R2rsm: Determination coefficient of the RSM.

  5. (5)

    Adj. R2rsm: Adjusted determination coefficient of the RSM.

  6. (6)

    F-statistic: F-statistic of the RSM.

  7. (7)

    P value: Significance probability of the RSM.

  8. (8)

    RSMs: Equation of the RSM.

Technical Validation

Model validation

The R2 and RMSE (root mean square error) were used to evaluate the performance of the RFM to simulate the NEE-RS (NEE simulated by RFM-RS), WF-RS (WF simulated by RFM-RS), NEE-WRS (NEE simulated by RFM-WRS) and WF-WRS (WF simulated by RFM-WRS)184. The model’s simulation accuracy for the WF was much higher than for NEE under each category and the performance of the RFM-RS was also better than the RFM-WRS (Fig. 2b, Table S3). For the WF simulation of the RFM-RS and RFM-WRS under each category, the percentage of the models with R2 ≥ 0.5 in the test sets was larger than 70%, while for the NEE simulation, it was lower than 50% (Fig. 2b). For the category ‘overall’, the RFM generally indicated a high simulation accuracy (Fig. 2b, Table S3). The simulation accuracy of the RFM was generally higher in Asia and the arid regions than in Europe and the non-arid regions. For the cropland and forest, the RFMs demonstrated the highest simulation accuracy under the scenarios of NEE-RS and NEE-WRS; for the grassland and wetland, the RFMs demonstrated the highest simulation accuracy under the WF-RS and WF-WRS scenarios.

The box plots (Fig. 3, Fig. S1) present the simulation performance of the RFM for NEE and WF in 10-time 10-fold CVs, with each box representing the R2 distribution for the test flux stations in each split (time). The simulation accuracy of the same RFM for different test flux stations varied widely, indicating that the RFM cannot be applied to all flux stations and that not all stations could match at least one available RFM model. The maximum R2 distribution for each flux station was observed in a 10-time 10-fold CV (Table S4). The proportion of the flux stations with R2 ≥ 0.5 of the RFM test measured 60.9%, 46.2%, 89.7% and 88.5% under the NEE-RS, NEE-WRS, WF-RS and WF-WRS scenarios, respectively.

Fig. 3
figure 3

The accuracy performance of the carbon-water flux simulation models (random forest model, RFM) at test flux stations. The NEE (net ecosystem carbon dioxide exchange) and WF (water flux) R2-based accuracy performance of the RFM of each split of the 10-time 10-fold cross-validation for (a) Overall with 156 stations, (b) Wetland with 16 stations, (c) Cropland with 23 stations, (d) Grassland with 47 stations and (e) Forest with 64 stations. The box plots show the R2 distribution of each flux station of the test set for different categories, in which the whiskers indicate the 1.5 times’ interquartile range.

Accuracy pre-assessment of the carbon-water flux simulation at the meteorological stations

Using the MLR model in which the Euclidean distances were independent variables and R2 concerned a dependent variable, the RSMs were constructed for different categories under different scenarios, as shown in the RFMs information file186. By using the RSM to simulate the R2 of the RFM at the test flux stations, the overall accuracies (of the RSM) for a correct classification of R2 under nine categories amounted to 80.1%, 84.0%, 91.0% and 89.1% for the NEE-RS, NEE-WRS, WF-RS and WF-WRS, respectively (Fig. S2). This might prove that the RSM is reliable and could be utilized to predict the accuracy of the RFMs applied to the meteorological stations.

Finally, the RFMs were transferred to all meteorological stations in Eurasia and the R2 was predicted for every meteorological station for each category under two scenarios (Fig. 4a–d, Table S5). In this study, the criteria for screening the RFM imply that the RFM corresponding to the highest predicted R2 of a given meteorological station and its R2 ≥ 0.5 was screened as the simulation model of the carbon-water fluxes for the meteorological station. The percentages of the meteorological stations in Eurasia were 84.5%, 68.2%, 99.1% and 98.7% for the NEE-RS, NEE-WRS, WF-RS and WF-WRS, respectively, in which the RFMs met the above-mentioned criteria (Fig. 4e). The RFMs have much higher applicable percentages and seem more accurate for the WF simulation than the NEE simulation at the meteorological stations. The RFM models of the forest and grassland categories were highly applicable and more accurate regarding the NEE and WF simulation of the meteorological stations than for cropland or wetland (Table S5).

Fig. 4
figure 4

The distribution of the R2 predicted by the R2 simulation model (RSM) at the meteorological stations. Spatial distribution of the R2 at (a) 4466 meteorological stations under the scenario of NEE-RS, (b) 6849 meteorological stations under the scenario of NEE-WRS, (c) 4466 meteorological stations under the scenario of WF-RS and (d) 6849 meteorological stations under the scenario of WF-WRS, respectively. (e) The percentage distribution of R2 < 0.5, 0.5 ≤ R2 < 0.7 and R2 ≥ 0.7 of the meteorological stations in the different scenarios.

The input data for the RFM were primarily derived from the meteorological stations’ observations and remote sensing data. Moreover, machine learning models (such as RFM) have the advantage and the predictive ability in the non-linear relation fitting and have been proven in the application research of relevant geoscience6, which is generally superior to linear regression, ecosystem process models, remote sensing inversion models, etc. Therefore, the carbon-water flux datasets of the meteorological stations generated in this study demonstrate a relatively high accuracy and constitute the attribute of quasi-observation, which might be considered as a quasi-observational dataset. They could be applied as benchmark data to verify the simulation results produced by the process-based models or remote sensing inversion models related to the carbon-water fluxes, which overcomes the challenge of insufficient observational data on the carbon-water fluxes165,187.

Spatio-temporal patterns of the Eurasian NEE and WF

We have further investigated the spatio-temporal distribution of the mean daily values of NEE and WF simulated by the RFM-RS during the period March-November from 2003 to 2020 and simulated by the RFM-WRS from 1984 to 2018. The meteorological stations with at least 30 data volumes of each spring, summer and autumn were used for statistical analysis. The mean daily values of the NEE-RS, NEE-WRS, WF-RS and WF-WRS at the meteorological stations are −3.9~0.7 g C m−2 d−1, −2.6~0.4 g C m−2 d−1, 0.8~3.8 mm d−1 and 0.5~4.3 mm d−1, respectively (Fig. 5a,b,d,e). The spatial distribution of these mean daily NEE fluxes reveals that the ecosystem carbon dioxide loss had increased in Eurasia during 2003–2020, with 457 carbon dioxide loss stations during this period, which means 178 more than from 1984 to 2018 (Fig. 5a,b). The daily average NEE (generally presented as net carbon dioxide uptake) has shown an increasing trend from 1984 to 2002, while a slightly decreasing tendency from 2003 to 2020, with slight fluctuations during these two periods (Fig. 5c). The temporal variation of the WF has demonstrated a rising trend with a distinct fluctuation from 1984 to 2020 (Fig. 5f). The differences between RS and WRS products might be caused by the differences in the input DSR dataset, the RFMs, and the number of meteorological stations (Fig. 5c,f).

Fig. 5
figure 5

Spatio-temporal variations of the carbon-water fluxes at the Eurasian meteorological stations. Spatial distribution of the mean daily values during the period March-November of (a) NEE-RS from 2003 to 2020 for 3436 meteorological stations and (b) the NEE-WRS from 1984 to 2018 for 4352 stations. (c) The annual temporal variation of the mean daily NEE (net ecosystem carbon dioxide exchange) values for the meteorological stations and the corresponding 95% confidence interval shown as a shaded band. Spatial distribution of the mean daily values during the period March-November of (d) WF-RS from 2003 to 2020 for 3990 stations and (e) WF-WRS from 1984 to 2018 for 6302 stations. (f) The annual temporal variation of the mean daily WF (water flux) values for the meteorological stations and the corresponding 95% confidence interval shown as a shaded band.

Comparison of the NEE and WF with other carbon-water flux products

We further compared the NEE and WF datasets with those from FLUXCOM1,3, GOSAT L4A (https://data2.gosat.nies.go.jp/) and MODIS (MOD16A2 Version 6)188 (Fig. 6, Table S6). The WF from the other products were converted from the LE. Our NEE and WF datasets, as well as the fluxes from the other products, were converted into a monthly scale. The months with the same number of stations for each product were selected for a comparison. All data show a similar seasonal variation, with high carbon-water fluxes in summer and low during winter (Fig. 6c,d). The distributions of the carbon-water fluxes from FLUXCOM and NEE from GOSAT are relatively discrete (Fig. 6a,b). The multi-year monthly averages of the NEE (NEE-RS = −0.31 g C m−2 d−1 and NEE-WRS = −0.34 g C m−2 d−1) and WF (WF-RS= +1.57 mm d−1 and WF-WRS= +1.48 mm d−1) simulated herein are less than those from FLUXCOM (NEE = −0.61 g C m−2 d−1, and WF= +1.79 mm d−1), whereas the averages are larger than those from GOSAT (NEE = −0.20 g C m−2 d−1) and MODIS (WF= +1.51 mm d−1). The WF from MODIS were almost consistent with our results (Fig. 6d). Because of the fact that the carbon-water flux datasets of the meteorological stations (generated by the RFM in this study) could be considered as “quasi-observational data”, the Eurasian carbon-water fluxes from FLUXCOM may be overestimated, while the NEE from GOSAT could rather be underestimated (Fig. 6c,d).

Fig. 6
figure 6

Comparison of the monthly NEE (net ecosystem carbon dioxide exchange) and WF (water flux) in this study with those from GOSAT, MODIS and FLUXCOM during the period 2010–2013. The box plots of the monthly values (black dots) for (a) NEE and (b) WF, respectively, in which the whiskers indicate the 1.5 times’ interquartile range. The monthly changes in (c) NEE and (d) WF and the corresponding 95% confidence interval shown as a coloured line and shaded band, respectively. RFM-RS: NEE or WF based on the RFM with remote sensing (RS), representing the fact that the RS variables were used in the RFM construction. RFM-WRS: NEE or WF based on the RFM without RS, illustrating that the RS variables were not used in the RFM construction. The RS variables include a fraction of the photosynthetically active radiation, enhanced vegetation index, land surface water index and surface reflectance for the Moderate Resolution Imaging Spectroradiometer bands 1–7. GOSAT, the GOSAT L4A data; MODIS, the MOD16A2 Version 6 data; FLUXCOM, an initiative to upscale the biosphere-atmosphere fluxes from the FLUXNET sites to the continental and global scales.