An empirical model for estimating daily atmospheric column-averaged CO2 concentration above São Paulo state, Brazil.

BACKGROUND: The recent studies of the variations in the atmospheric column-averaged CO2 concentration ([Formula: see text]) above croplands and forests show a negative correlation between [Formula: see text]and Sun Induced Chlorophyll Fluorescence (SIF) and confirmed that photosynthesis is the main regulator of the terrestrial uptake for atmospheric CO2. The remote sensing techniques in this context are very important to observe this relation, however, there is still a time gap in orbital data, since the observation is not daily. Here we analyzed the effects of several variables related to the photosynthetic capacity of vegetation on [Formula: see text] above São Paulo state during the period from 2015 to 2019 and propose a daily model to estimate the natural changes in atmospheric CO2.
RESULTS: The data retrieved from the Orbiting Carbon Observatory-2 (OCO-2), NASA-POWER and Application for Extracting and Exploring Analysis Ready Samples (AppEEARS) show that Global Radiation (Qg), Sun Induced Chlorophyll Fluorescence (SIF) and, Relative Humidity (RH) are the most significant factors for predicting the annual [Formula: see text] cycle. The daily model of [Formula: see text] estimated from Qg and RH predicts daily [Formula: see text] with root mean squared error of 0.47 ppm (the coefficient of determination is equal to 0.44, p < 0.01).
CONCLUSION: The obtained results imply that a significant part of daily [Formula: see text] variations could be explained by meteorological factors and that further research should be done to quantify the effects of the atmospheric transport and anthropogenic emissions.

Background

Understanding the variability of atmospheric carbon dioxide (CO2) concentration in time and space is a crucial task so that we can adopt mitigation strategies. In this sense, several studies analyze the average concentration of this greenhouse gas not only on a global scale [1, 2] but also to estimate anthropogenic emissions in urban centers [3, 4]. Other studies focus on understanding the column-averaged of carbon in the atmosphere () above tropical forests [5], or above agriculture crops in different seasons of the year [6, 7].

In a recent regional study, da Costa et al. [7] analyze the spatio-temporal variability of in a sugarcane-producing area in the southeast region of Brazil. They observed an important inverse relationship between the average carbon concentration in the atmosphere with climatic and vegetative variables. Concluding that the dependence of the natural carbon cycle is related to the predominant agriculture crop in the region and how Global Radiation (Qg), relative humidity (RH), and the Sun Induced Chlorophyll Fluorescence (SIF) was related to this behavior. Similarly, Morais Filho et al. [6] conducted a study that analyzed three different crops and the temporal variability of and SIF in these environments, they also found a significant negative correlation between these variables.

However, there is still a temporal gap in the data collected by remote sensing, since the measurements are not daily [8, 9]. This type of measurement is important to several factors, such as, estimate the potential capability of atmospheric CO2 assimilation by vegetation, establishing public strategies at local levels for climate adaptation and mitigation, and even in economy incorporating daily trends in the carbon market and ecosystems services payments [10-16].

Daily CO2 measurements can be made using the Eddy Covariance technique [17-19], although this has the disadvantage of being a point (local) study. In this sense using orbital data, such as Orbiting Carbon Observatory-2 (OCO-2), has become more common [1, 2]. Remote sensing data also can be used to estimate the daily variations of different aspects (e.g., climate, meteorological, land-use changes, ecosystems services) for a larger area [16, 20–23].

Several studies confirm that photosynthesis is the main regulator of atmospheric carbon sinks [10, 24–26]. However, photosynthesis is a process sensitive to climatic variations such as relative humidity [27], precipitation [28], evapotranspiration [29], and incident solar irradiance [30].

Therefore, the natural cycle of CO2 is dependent on several aspects, such as vegetation and climate, being necessary data from several different bases for understanding this dynamic [7], turning pre-processing techniques and analysis of autocorrelations necessary, since the multicollinearity introduces an uncertainty due to the model overfit [31, 32]. In this sense, we aim to model the atmospheric CO2 cycle above the state of São Paulo to estimate the time changes on a daily scale, based on vegetative and climatic variables retrieved from different orbital platforms, applying a technique to remove the collinearity and after employing a stepwise forward selection, improving in this way the regional understanding of CO2.

One of our assumptions, is given that we detrend the to maintain only the variability related to the natural interactions [6, 33], the transport by the wind in the atmosphere is not significant, other studies such Hakkarianen et al. [34], that proposed an anomaly model of , also disregard the atmospheric and wind transport in their study, however, this introduces a limitation of our approach that not account this aspect [35]. In the same way, the trend and increase due to anthropogenic sources are also simplified by this detrend.

Results

Variance Inflation Factor (VIF) analysis (Table 1) shows it was possible to reduce the number of variables related to (according to the adopted criterion, VIF < 10) as shown comparing Fig. 1a with b, before and after the selection, respectively, and therefore reducing the overfit source of uncertainty. Despite wind speed (Ws) hade a VIF < 10, the Pearson’s correlation was not significant (p > 0.05). Variables most related to were the Global Radiation (Qg), Sun-Induced chlorophyll Fluorescence at 757 nm (SIF 757), and Relative Humidity (RH).

Table 1

Variance Inflation Factor (VIF) of the studied variables

Variable	VIF
Qg	9.35
RH	5.10
SIF 757	1.81
Prec	10.13
Temp	21.43
Ws	4.06
LST	19.54
NDVI	22.15
LAI	87.21
Fpar	65.40
ET	33.47

Qg Global radiation, RH Relative humidity, SIF 757 Solar-Induced Chlorophyl Fluorescence at 757 nm, Prec Precipitation, Temp Temperature at 2 m, Ws Wind Speed, LST Land Surface Temperature (MODIS), NDVI Normalized Difference Vegetation Index, LAI Leaf Area Index, ET Evapotranspiration

Fig. 1

Heatmap of the Pearson’s correlation matrix, where: a before the Variance Inflation Factor (VIF) selection and b after the selection by Variance Inflation Factor (VIF)

Regarding the temporal variability of , the maximum mean for the analyzed period was 393.09 ± 0.17 ppm and occurred in October 2019, while the minimum average was in November 2018, being 390.11 ± 0.15 ppm (Fig. 2a). Meanwhile, the Qg (Fig. 2b) ranged between 24.3 ± 0.09 and 13.07 ± 0.04 (MJ m−2 day−1), with the maximum average occurring in December 2018 and the minimum in June of the same year.

Fig. 2

Monthly variability of (a), Qg (b), RH (c), and SIF 757 (d) over the period from January 2015 to December 2019. Where column average of carbon dioxide in the atmosphere (ppm), Qg global radiation (MJ m−2 day−1), RH relative humidity (%), SIF 757 sun-induced chlorophyll fluorescence at 757 nm (Wm−2 sr−1 μm−1)

SIF 757 (Fig. 2d) had the highest average recorded in the period in November 2015 [1.1 ± 0.05 (Wm−2 sr−1 μm−1)] and the lowest in September 2017 [0.3 ± 0.06 (Wm−2 sr−1 μm−1)], while the Relative Humidity (Fig. 2c) ranged from 84.86 ± 0.07 to 70.44 ± 0.19%, where the highest mean was observed in March 2016 and the lowest in October 2019.

Regarding SIF 757 the minimum averages occurred in June of 2015 and 2016, September 2017, November 2018, and July 2019, ranging from 0.3 to 0.46 Wm−2 sr−1 μm−1 (Fig. 2d). The minimum Qg averages vary between May and June for the entire series approximately between 13.07 and 14.71 MJ m−2 day−1 (Fig. 2b). Maximum Qg averages are concentrated between December and January of each year, reaching 24 MJ m−2 day−1 in those months (Fig. 2b). The maximum average of SIF 757 occurs between November and February of each year, ranging from 0.8 to 1.1 Wm−2 sr−1 μm−1 (Fig. 2d).

The stepwise forward selection method, with multiple cross-validation, had the best result with two variables, with a root mean squared error (RMSE) of ~ 0.60 ppm in the training sample (Fig. 3), the selected variables being Qg and RH, respectively (Eq. 1).

Fig. 3

RMSE score for the training sample

The model built in the training (Eq. 1) was applied in the test sample of the variables cited (Qg and RH), and from the cross-validation of the estimated data with the observed data, we observe an R2 of 0.44, the values of the metrics MSE, RMSE, and MAE were 0.22, 0.47, and 0.37 (ppm) respectively, and for MAPE we found a value of 1.54% (p < 0.01) (Fig. 4a), with this we were able to reduce the time scale of the OCO-2 satellite from every 15 days to a daily scale (Fig. 4b).

Fig. 4

a Cross-validation between estimated by stepwise and observed by OCO-2 and b Daily downscale of natural using Eq. 1 and daily measurements of NASA/POWER from January 2015 to December 2019

Discussion

The natural annual cycle of is affected by factors related to climate and vegetation aspects [6, 36, 37]. Due to the VIF analysis, we were able to summarize three main factors for São Paulo state: Global Radiation (Qg), Relative Humidity (RH) and Sun-Induced chlorophyll Fluorescence at 757 nm (SIF 757), reducing the uncertainties in the model formulation since we removed the overfit caused by multicollinearity [31, 32]. Several studies have already been conducted using this method to identify which variables select for ecological studies [38], computational studies [39], and remote sensing studies [40].

Except for wind speed (Ws), all variables studied correlated negatively with (Fig. 1), hence, related to the sink of atmospheric CO2. The non-significant correlation between and Ws could be related to the detrending of the atmospheric CO2 concentration (see Methods section), which removes the transport effect and simplify the variability only for the biochemical cycle [6, 7, 33]. In general, the highest concentrations of are observed in the months corresponding to the Brazilian autumn and winter (April to August) and lowest in the summer, from December to February. Studies such as by Siabi et al. [41] and Falahatkar et al. [42] reported how the different seasons affect the average CO2 concentration in the atmosphere.

Recently, researches were conducted at regional scales in Brazil such as by Morais Filho et al. [6] and da Costa et al. [7], indicating negative correlations between and SIF over agricultural areas, approximately − 0.5 and − 0.8, respectively. SIF is a variable directly related to the photosynthesis of plants, laboratory-scale experiments have demonstrated this relation [43], and remote sensing studies at the canopy and global level reported positive relations between SIF and Gross Primary Production, and also a negative correlation between SIF and the [5, 44–46].

As a result of photosynthesis, it is expected that SIF increases during summer [7, 41], as in this season, higher precipitation events and higher temperatures are observed [47]. Our results show higher SIF average values in the months when summer occurs in the São Paulo state, and an inverse relationship between SIF and . The lowest average values of usually occur during the summer period in the study region. This is due to plant CO2 assimilation [48], printing quasi-periodical , and SIF time changes as well as observed in other studies [5, 6, 41, 49].

Most of São Paulo’s state has a wet summer and dry winter [47] resulting in a positive correlation between precipitation and SIF (Pearson’s correlation = 0.61 and p < 0.05), while negative with (r = − 0.49, p < 0.05) (Fig. 1a). Precipitation is a photosynthetic control factor, so the greater availability of water that exists in the summer in São Paulo’s state induces plants to perform more photosynthesis through primary productivity, which leads to a reduction of atmospheric CO2. The opposite is observed in the dry winter because water availability is lower resulting in less photosynthesis, or less CO2 assimilation by plants, either in natural or agricultural areas [7, 28, 50].

Another effect observed during summer in the region is the increase of relative humidity (RH), which reduces the water transfer between soil or plant to the atmosphere [51], inducing plants to keep their stomata open, where CO2 assimilation occurs [52]. Studies have already shown the relationship of stomata opening in periods with good water availability is related to plant growth [53, 54]. Thus, establishing the negative relationship between RH and , also previously observed by Golkar et al. [27].

In the same way, another requirement for photosynthesis occurs is sunlight, which is the source of energy to carry out the biochemical processes of this phenomenon. Therefore, as the amount of radiation (Qg) is absorbed by the plant, photosynthesis tends to increase, and consequently higher CO2 assimilation, decreasing in this way the concentration of this greenhouse gas in the atmosphere [7, 30]. We can observe these relationships in our results (Fig. 3b), Qg correlates positively with SIF, and these variables relate negatively with .

Since we are dealing with the natural annual cycle of CO2 the main factor of the higher concentrations of this gas in the atmosphere is due to the lowest photosynthetic absorption by plants. The autumn and winter have low available water and sunlight for plants, leading to a decrease in photosynthesis, also another important factor is that the annual calendar for agriculture in the state of São Paulo has harvest periods between these seasons [55], and as consequence decreasing the cover area by vegetation. Shekhar et al. [56] show how the crop’s grown in summer decrease the values of over the Nile Delta and when the harvest starts the values of are higher, also, they found that SIF values are higher in the grown season.

Our model was based on Qg and RH, which are two variables related to the CO2 assimilation process, or CO2 sink. The model has lower RMSE values than have been reported in previous studies, such as by Guo et al. [57] where the values of this metric ranged from 0.7 to 1.1 ppm. In a more recent study by Taylor et al. [58] when evaluating initial OCO-3 data results from the globe and model-related errors, they found an RMSE between 1 and 2 ppm. Another important measure is the MAPE, which shows in percentage how much we are getting wrong, studies with remote sensing have already demonstrated errors below 10% as being considered extremely low for predicting plant and climate aspects [59, 60]. With this, we can evaluate that the performance of the model proposed in this work presents a very low error.

The coefficient of determination (R2) was 0.44, an increment of almost 20% from the simple linear fit with Qg alone, which has a higher importance in the model. Although the R2 is moderate, studies using other orbital sensors such as MODIS to model the average CO2 concentration in the atmosphere have reported similar results [23]. In addition, we should consider that although OCO-2 and NASA-POWER are two high quality and validated databases [8, 9, 61], the difference between grids and spatial resolution (see Table 2 in Methods and Fig. 5b) cannot be disregarded, as it is an aspect that can influence these results, leading us to consider the coefficient of determination observed in this study as being high.

Table 2

Studied variables, data base, temporal and spatial resolution

Variable	Data base	Temporal resolution	Spatial resolution
GHG
X_CO_2 (ppm)	OCO-2 “OCO-2 Data product user’s guide, 2016” V9	16 days	1.29 km × 2.25 km
Climate
Surface solar shortwave irradiance (Global radiation, Qg) (MJ m−2 day−1)	FLASH Flux Version 3 (A, B,C) NASA/POWER	Daily	111.3 km × 111.3 km
Average air temperature at 2 m (Temp) (ºC)	GEO-5 FP-IT (NASA/POWER)	Daily	111.3 km × 111.3 km
Land surface temperature (LST) (ºC)	MOD11A1.006 V6 MODIS-TERRA	Daily	1200 km × 1200 km
Wind speed at 10 m (WS) (m s−1) (ºC)	GEO-5 FP-IT (NASA/POWER)	Daily	111.3 km × 111.3 km
Relative humidity (RH) (%)	GEO-5 FP-IT (NASA/POWER)	Daily	111.3 km × 111.3 km
Precipitation (Prec) (mm day−1)	GEO-5 FP-IT (NASA/POWER)	Daily	111.3 km × 111.3 km
Vegetation
SIF 757	OCO-2 “OCO-2 Data product user’s guide, 2016” V9	16 days	1.29 km × 2.25 km
LAI (m2 m−2)	MCD15A2H.006 V6 MODIS- CFPAR	8 days	500 m × 500 m
Fraction of Photosynthetically Active Radiation (Fpar) (%)
Evapotranspiration (ET) (kg m−2 day−1)	MOD16A2.006 V6 MODIS-TERRA	8 days	500 m × 500 m
NDVI	MOD13A1.006 V6 MODIS-TERRA	16 days	500 m × 500 m

Fig. 5

Flowchart of data acquisition, processing, analysis (a), and sounding map of the satellite observations in the study region (b). Where red dots represent the OCO-2 soundings, the black dots represent the NASA-POWER and the dark green represent the MODIS

These differences between the databases can be suppressed by the greater temporal coverage of NASA-POWER, allowing us to estimate the daily temporal variability of the natural CO2 cycle in the atmosphere for the state of São Paulo, besides reducing in the future the spatial scale of obtained from OCO-2 and gaining greater spatial resolution cover. Other vegetation index-based models aimed at reducing the spatial sampling of OCO-2 data, but focused on SIF, as is the case of Zhang et al. [62] and Yu et al. [63].

Despite the errors associated with the model and the uncertainty measures due to the difference in satellite resolution, an advantage of using models similar to the one proposed here is being able to have a daily measure of the variability of atmospheric CO2 and how the climate parameters affect this dynamic, also serving as an indirect indicator of how is the daily assimilation capacity of this gas in a region.

Conclusions

In summary, the cycle of in the state of São Paulo has higher average values during April to October, periods of lower intensity of rainfall, and is considered as the winter in the state, in the other hand the lowest averages of were usually observed between December to March, this period corresponds to the summer, and the inverse behavior was observed for SIF 757, global radiation (Qg) and relative humidity (RH). This pattern is due to the relationship between photosynthesis and Carbon assimilation, given that photosynthesis is a process sensitive to climate variation and a process that depends on water and light, in summer this process tends to be greater, leading to a decrease in CO2.

Concerning the daily model presented, it performed well when we looked at the set of metrics presented. Given this, we were able to estimate the daily behavior of natural in general for the state of São Paulo, a semi-periodical wave with a maximum peak between March and July, and a minimum peak between December to February. There are still challenges in this aspect, such as the transport process in the atmosphere, which was simplified due to the detrend in the dataset, that also remove the anthropogenic sources in the CO2 cycle, however, this study was capable in advance in the temporal gap, and properly address how to estimate the natural behavior of this gas in a synthetic way using daily meteorological open access data, establishing a low-cost approach, and we believe that this study will serve as a basis for further implementations.

We suggest that for future work, the relationship between soil respiration and factor controlling organic matter decay in soil with the would be needed to better understand CO2 dynamics, as well the addition of variables related to activities, such as in transports or the data of fossil fuel consumption, in big cities to improve predictions, as well the atmospheric transport.

Methods

Study region

The state of São Paulo (SP) (Fig. 5b), southern Brazil, has approximately 249 × 103 km2 and 645 municipalities, with a demographic density of 179.84 habitants/km2 [64] being one of the main agricultural hubs of Brazil, regarding the production of sugarcane and citrus [65]. According to Rolim et al. [47] the climate of the state, in general, has its areas characterized by a humid subtropical climate with dry winter, followed by humid tropical dry winter and sub-humid tropical dry winter, according to the climate classification proposed by Camargo [66].

Products of remote sensing: acquisition and processing

Greenhouse gas, climate, and vegetation data were collected from different satellites (Table 2) for a time series from 2015 to 2019 and were aggregated on a monthly scale. The primary product of the Orbiting Carbon Observatory-2 (OCO-2) consists of georeferenced estimates of the mean atmospheric CO2 concentration (), in addition, the Sun Induced Chlorophyll Fluorescence (SIF), retrieved due to the overlap that occurs in the SIF wavelengths with the O2 absorption wavelength (680–850 nm) [8, 9, 43]. Data from this satellite have already been validated by Crisp et al. [8] and, according to O'Dell et al. [9], this satellite provides about 65,000 quality observations per day worldwide.

Here we used the version 9 of the OCO-2 with a bias-correction and considered only the measurements with the best quality flag (quality flag = 0, meaning that has no cloud cover) [67, 68], also, we do not consider the data with more than 12 alert level at nadir viewing [33, 69]. Concerning the SIF, we take into account only the SIF at 757 nm, this was due to previous studies that exploited the relationship in the São Paulo’s State [6, 7] and, also because this wavelength is closer to the far-red peak (~ 740 nm) in the whole SIF signal [43].

MODIS sensor data were extracted from the “Application for Extracting and Exploring Analysis Ready Samples” (AppEEARS). This application allows users to obtain subsets of large databases using spatial and temporal parameters. Two types of sample requests are available: point samples by entering geographic coordinates and area samples using vector polygons. Sample requests submitted to AppEEARS provide users with not only data values but also associated quality data values. Interactive visualizations with summary statistics are provided for each sample within the application, which allows users to view and interact with their samples before downloading the data [70].

Nasa Power data (https://power.larc.nasa.gov ) consists of precipitation (mm), surface solar shortwave irradiance (MJ m−2 day−1), average air temperature (ºC), and relative humidity at 2 m (%). This platform consists of a NASA project entitled: Worldwide Energy Resource Forecast (POWER) and was initiated to enhance the current renewable energy dataset and create new datasets from new satellite systems [71].

To minimize the differences between the spatial and temporal resolutions of the different orbital sensors used in this study, the process described in Fig. 5a was employed, which establishes a standard for the acquisition of data from the coordinates obtained in the OCO-2 platform (Fig. 5b). We emphasize that several studies have been conducted using different time and spatial scales [6, 7, 27].

Pre-process of the data

Using the regression method proposed by Gujarati and Potter [72], we removed the trend from data, in order to understand the natural and regional variability of and its relationships with other factors [6, 7, 33]. The other variables were standardized using the function scale from the R language [73].

Variance Inflation Factor (VIF)

Variance Inflation Factor (VIF) analysis was performed. This analysis is a method of detecting multicollinearity within a database since the relation between the predictors for a multi-regression model can affect the estimative and the standard errors associated with the regression model [31]. The VIF is based on the R2 value (Eq. 2), and should not be greater than 10, however, this can vary according to the study [31, 32]. where R2 is the coefficient of determination.

Temporal variability, Pearson’s correlation, and dependency analysis

The data was processed using month averages for the analysis period, except precipitation, which consists of monthly sums for the entire state of SP (ST.1). The means were subjected to analysis of variance (F-test) to obtain the mean standard errors. Simultaneously, the basic assumptions of analysis of variance and, normality of errors, and homogeneity of variances were tested for the selected variables by VIF analysis. To understand the variation of with the other variables, Pearson correlation analyses were performed. More about the descriptive statics of selected variables in VIF, such as the number of observations (soundings) for each month, can be found in Additional file 1: Table S2

Stepwise: forward selection

The stepwise method used in this study was the forward selection method being performed in R language [73], as can be seen in the flow chart (Fig. 6), the variables selected in the VIF analysis were separated into a training and test samples (70% and 30% of the dataset respectively). The training sample was submitted to the train () function of the caret package, using repeated cross-validation (cv) method. This technique consists in randomly splitting the training dataset into k-subsets, one of them is reserved and the model is trained with the others, and after is validated with the reserved subset, this process is repeated until each subset serves as a test sample, finally, the average error is how the performance is given [74]. The model is based on the lowest Root Mean Squared Error (RMSE) and, from variables selected in training, the generated model is applied to the test sample defined at begging for estimating the with these independent data. Finally, cross-validation between the estimated data and observed data in the test sample was performed and from this, we derive the metrics Mean squared error (MSE), Root Mean Squared Error (RMSE), Mean Absolute Error (MAE), R2, and Mean absolute percentage error (MAPE).

Fig. 6

Flowchart of the stepwise construction

Additional file 1: Table S1. Variables used in the VIF analyses. These variables were standarlized by scale() in R language. Table S2. Descreptive statistics.