Carbon Exchange and Accumulation in an Orinoco High Plains Native Savanna Ecosystem as Measured by Eddy Covariance

Savanna ecosystems cover ∼20% of the total land surface and account for ∼30% of the terrestrial global net primary production. They are also highly sensitive to climate change, since their carbon (C) sink capacity may decline under rising temperatures and irregular rainfall. These responses, which will define the future climate role of the savanna ecosystems, are currently not well understood. The Colombian Orinoco River basin (“Llanos”) natural savannas are being rapidly converted to agriculture. The impact of this transformation on C fluxes and accumulation is not clear. It is thus urgent to understand the Llanos natural savanna ecosystem services, including their C cycle and underlying mechanisms. Here we report and analyze 2 years of measurements of carbon dioxide fluxes from a naturally-restored (secondary) Llanos High Plains savanna ecosystem, using eddy covariance. Meteorological conditions, particularly rainfall, were quite variable during the measurement period. During the first year of measurements, the savanna was a weak carbon source (35 gC m−2), while during the second year, the system was a comparatively strong carbon sink (−273 gC m−2), despite receiving less rainfall than during the first year. As expected, the savanna net ecosystem exchange (NEE) was highly dependent on global solar radiation, soil water content, and ecosystem respiration. We found that after ∼10 days of nominal drought, i.e., with less than ∼5 mm/day of precipitation, the NEE became highly dependent on drought duration. The ecosystem reached a critical condition of low photosynthetic activity after ∼60 days of nominal drought. Based on these findings, we developed and applied a simple standard meteorology-based model that properly reproduced the observations. Our results indicate that a shift to a climate with similar total precipitation but split into extreme dry and wet seasons might eventually suppress the savanna C uptake capacity.


INTRODUCTION
Terrestrial ecosystems sink about one third of the total anthropogenic greenhouse gas (GHG) emissions from industrial activity and land-use changes (Keenan and Williams, 2018). The urgency of climate change mitigation requires these ecosystem services to be understood, quantified, monitored, modeled, and enhanced when feasible. This is particularly important in the tropics, where soil and aboveground carbon accumulation, and other ecosystem services are not sufficiently well understood. Remaining natural ecosystems are severely threatened by pressures of land-use changes, including those related to food and industrial feedstock production. Understanding the functioning and services of natural and restored tropical ecosystems is of paramount importance to conservation efforts and for rational mitigation decision-making.
Savanna ecosystems cover ∼20% of the total land surface and account for ∼30% of the global net primary production (NPP) (Lipsett-Moore et al., 2018). Therefore, they play a very important role in the C cycle; but they are also highly sensitive to climate change, since their C sink capacity is expected to decline under rising temperatures (Bond-Lamberty et al., 2018) and irregular rainfall (Fei et al., 2017). The vegetation of tropical savannas consists of a continuous herbaceous cover of mostly C4 grasses and sedges along with significant but discontinuous cover of C3 shrubs and trees (Walker, 1987). These ecosystems are also characterized by strong rainfall seasonality, which implies periods of significant water stress (Frost et al., 1985). Conversion of savannas into agricultural lands may result in a reduction of its C sink capacity to the extent of changing the ecosystem into a C source. Currently, the global rate of savanna loss is not well-established, but it may exceed 1% per year, approximately twice as fast as that of rainforests. This likely imply carbon dioxide (CO 2 ) emissions at least as large as those arising from rainforest deforestation (Grace et al., 2006).
Savannas cover an area of ∼270 million ha (Mha) in South America, mainly distributed in Brazil (76%), Venezuela (9%) and Colombia (8%). The large areas of savannas are seen as the main alternative to avoid agricultural expansion in tropical areas of greater ecological fragility, such as the Amazon rainforests (López-Hernández et al., 2005). Although most of their soils suffer from physical and chemical limitations to agricultural production (López-Hernández et al., 2005), large areas of South American savannas have been turned into monocultures or pastures (Borghetti et al., 2019). The Colombian Orinoco River basin savannas, also known as "Llanos", particularly, the flat High Plains region, has also been rapidly transformed into mechanized agriculture. The agriculture expansion of this region is expected to continue (Jimenez et al., this issue) since market demand and government policy for the development of the Llanos High Plains create conditions that encourage investment to take advantage of the agricultural and agrobusiness potential of the region. The Colombian government estimates that 2.8 Mha of the High Plains have potential for agricultural, livestock and forestry use (DNP, 2014).
This agricultural expansion in the Llanos High Plains creates a concern about its impact on the carbon cycle. There are very few reports (Etter et al., 2010) on the C role of the Colombian Llanos savanna ecosystems and a critical knowledge gap on their controlling factors. San José's et al. carbon flux measurements in the seasonallyflooded Venezuelan savannas found that their C content increases was mainly due to increases in their woody fraction (San José, 1991;San José et al., 1998a;San José, 2001;San José and Montes, 2007;San José et al., 2014). Regarding the Colombian High Plains savannas ( Figure 1A), we found only two studies on GHG fluxes, specifically on soil respiration (Lavelle et al., 2014;Agrosavia, 2017). We found no reports on Net Ecosystem Exchange (NEE) over this essential ecosystem, nor about its response and main environmental controlling factors. Our research addresses the knowledge gap on the Colombian savanna C fluxes by answering the following question: What environmental factors control CO 2 exchange in the High Plains savannas?
In this study, we used the eddy covariance technique to measure the exchange of CO 2 between the atmosphere and a naturally-restored tropical, Northern South America savanna ecosystem located in Llanos High Plains. We conducted measurements during two years of very distinct meteorological conditions. The accumulated rainfall during 2016 was significantly higher than during 2017 but the rainfall distribution during 2016 was very anomalous, with an extremely low-precipitation dry season and higher-thanaverage precipitation during its wet season. Based on our findings, we developed and applied a simple, standard meteorology-based model that properly reproduces the observed CO 2 fluxes.

Site Description
The Llanos comprises four subregions (piedmont, high plains-flat and dissected, alluvial overflow plains, and aeolian plains), each having different landforms, soil, savanna formations and vegetation patterns. The Llanos westerly wind patterns are synoptically controlled by a Walker-type circulation (Mesa-Sánchez and Rojo-Hernández, 2020), and its climate is highly influenced by the position of the Inter-Tropical Convergence Zone. The precipitation regime is primarily unimodal in most of the region, except for the piedmont, where the regime is bimodal due to the Andes influence (Pulwarty et al., 1998;Pacheco and León-Aristizábal, 2001). There is a decreasing precipitation gradient through the Llanos with a west-to-east axis, with annual precipitations from ∼3,000 mm toward the southwestern Andes foothills in Colombia to 900-1,300 mm in the Orinoco River delta in Venezuela (San José et al., 1998b;Rippstein et al., 2001). Consistently, the dry season length varies from 1 or 2 months in the southwest to 5-6 months in the northeast (Rondón et al., 2006).
The Colombian Llanos cover around 23 Mha (20% of the Colombian continental territory) with altitudes between 80 and 500 m a.s.l., and ecological gradients between the Colombian Andes and the Amazon (Figure 1). Its soils are mostly of low fertility, acidic, with low organic carbon content and toxic levels of iron and aluminum (De la Hoz, 2009). Among other factors, this is the result of their seasonal exposure to floods, droughts, burning and strong winds. The Llanos High Plains, located south of the Meta River in Colombia (∼42% of the Colombian Llanos) are slightly more elevated than the alluvial and aeolian plains, making them much less subject to seasonal flooding. In spite of their low fertility, the Colombian Llanos High Plains, particularly the flat High Plains, are being rapidly transformed into mechanized agriculture.
Measurements were conducted over a unique, large patch of seminatural (native but secondary) savanna located in the Colombian High Plains (Figure 1), specifically at Agrosavia's Taluma Experimental Station (EC-TS, 4°22′24.10″N, 72°13′19.23″W, 173 m a.s.l.). The patch, naturally transformed back into a High Plains native savanna, is flat, homogeneous, and with enough area to guarantee the quality of the measured fluxes (more than 4 ha). The site has been actively protected from fires, agricultural, and grazing activities for at least the last 20 years. The savanna site is under a monomodal precipitation regime, in which the dry season typically runs from December to mid-March with late March and November as transition time periods (Figure 2). The mean annual rainfall in the region during the last 15 years (2005-2020) was 2,294 mm/yr (1836 mm/yr to 2,607 mm/yr), with almost all the precipitation (∼90%) occurring between april and november (at Hacienda Las Margaritas station; 4°20′26.00″N, 72°9′22.90″W). The average daily temperature at Taluma was 26.4°C (22.6-31.2°C) and the According to studies conducted at Taluma by Agrosavia, the predominant soils are Oxisols (Tropeptic Haplorthox), with a loamy-sandy texture and an average pH of 4.3. These soils are naturally compressed, with an apparent density of 1.53 g cm −3 . They are susceptible to compaction and erosion, have low infiltration, low water storage capacity, and are highly vulnerable to mechanization. Macropores account for 3-5%, while mesopores add up to 9-15% of the porous space. Measured calcium (Ca) content of the soils was within 0.05-0.10 cmol kg −1 , magnesium (Mg) 0.02 to 0.01 cmol kg −1 at the same depth, and there was a low content of available phosphorus, 1.3 mg kg −1 . Minor elements like boron (Bo), copper (Cu) and zinc (Zn) have been found in concentrations lower than 0.5 mg kg −1 (Agrosavia, 2017).
The herbaceous layer in the Colombian High Plains savannas is typically 0.3-1.0 m high (Rippstein et al., 2001). Shrubs are abundant at our site (see Supplementary Figure S1). When these are included, the average canopy height at our measurement footprint was 1.6 m on average. The families Poaceae (grasses), Fabaceae (legumes) and Cyperaceae are the most representative plant families and are endemic to this type of savannas. They are followed by the families Asteraceae (Compositae), Rubiaceae, Labiatae (Lamiaceae), and Melastomataceae. The Poaceae family is the richest in species and its dominant genera are Paspalum, Panicum, Axonopus, and Andropogon. In the Fabaceae family, the genera are better distributed, and in the family Cyperaceae, the dominant genus (in species) is Rhynchospora (Rippstein et al., 2001).

Eddy Covariance and Ancillary Environmental Measurements
Fluxes of sensible heat (H), latent heat (λE), and CO 2 were measured with an eddy covariance (EC) system composed of a CSAT3 sonic anemometer (Campbell Sci., Logan, UT, United States) and an open-path Li-7500A gas analyzer (Li-Cor, Lincoln, NE, United States), hereafter known as the EC-TS. The gas analyzer was set-up with a separation of −10 cm northward, −5 cm eastward and −10 cm vertical. The EC system was mounted on a tripod at 3.1 m height. The gas analyzer was calibrated using two standard gases to reset the CO 2 zero (Linde Zero-Air gas) and another for CO 2 span (518 ppmv CO 2 source with an accuracy of ±1 ppm). A Biomet 101 system (Li-Cor, Lincoln, NE, United States) was used to measure ancillary biometeorological variables at EC-TS. This equipment was available for a limited time (from December 2015 to May 2016) at this site. The measured variables included soil heat fluxes, soil water content, soil temperature, air temperature, humidity, and global, net, and photosynthetic active radiation (PAR). These data were used to validate the simpler Taluma Weather Station (TWS-iMetos 3.3 Pressl instrument, Werksweg, Austria) data. The period when Biomet and TWS measured simultaneously was used for developing correction factors to TWS (see Supplementary Table S1). When Biomet data was not available, the main environmental parameters were measured with TWS.
The variables needed for EC CO 2 flux measurements, i.e., the vertical and horizontal wind speed components and the CO 2 mixing ratio, along with the virtual sonic temperature, were measured at 10 Hz. The companion, low-frequency meteorological measurements were averaged and recorded at a time resolution of 1 min. Sensors were set up following Munger et al. (2012) guidelines. Measurements were conducted from 2015-12-20 to 2017-12-14. The time series had two large data gaps due to logistical issues: 1) from 2016-10-08 to 2016-12-08 and 2) from 2017-03-08 to 2017-05-05.
The Mauder and Foken (2004) micrometeorological quality assessment evaluates whether the theoretical requirements of EC are fulfilled, using tests for stationarity and integral turbulent characteristics. The results from these two tests were combined to derive a data quality classification into three groups: flag 0, flag 1, and flag 2, which correspond to high-quality, intermediatequality, and low-quality data, respectively. We only considered the high and intermediate quality data for evaluating fluxes. The data collected during precipitation events were also discarded.
A storage term was added to the CO 2 fluxes to account for changes in concentrations within the measurement volume, i.e., under the savanna canopy during the 30 min averaging period (Aubinet et al., 2012). The CO 2 storage was estimated from CO 2 concentrations and based on a 1-point profile. This technique is based on the assumption that the 30 min changes in concentration at the instrument level are representative of the whole air layer underneath the EC system. EC measurements can be affected by high-frequency spikes due to fast changes in turbulence conditions, changes in footprint, water droplets on sonic anemometers, and other factors. Although most of the spikes are removed in the high frequency data processing, some anomalous data get into the half-an-hour time series from which the half hourly fluxes are estimated, producing false fluxes. Outlier detection techniques are applied to prevent this from happening (Papale et al., 2006). Two additional data assessment criteria were applied to identify these outliers: 1) coherence or plausibility, and 2) consistency or step test (Vejen et al., 2002;Zahumenský, 2004). For both tests, the evaluation was performed using a customized moving average filter that used a time window size and a signal "amplification" factor as parameters. This filter is described in Section 3 of the Supplementary Material. For coherence, the thresholds for the variables were defined using a long moving time window (7 days) and a vertical amplification factor of 3 times the standard deviation (SD) of data within the window. Data above and below the 3 SDs within the sliding window were flagged as outliers. All fluxes flagged as outliers and their associated variables were discarded (e.g. for H, H storage, flux error, and variance of sonic temperature were discarded). This procedure was performed twice. For consistency, the thresholds for the variables were defined by applying the moving average filter using a short horizontal window (2 h before and 2 h later of the evaluated point), and a vertical amplification factor of 3 times the standard deviation (SD) of data in the horizontal window. Data out of these thresholds were also discarded.

Nighttime Flux Estimation
During nighttime, eddy covariance provides direct measurements of nighttime ecosystem respiration. However, carbon flux measurements can be severely hindered by the low turbulence conditions prevailing at night. The first step to amend this was to identify and discard biased flux data. The filtering methodology proposed by Aubinet et al. (2012) was applied to nighttime data (19:00-05:30 LT), but using the standard deviation of vertical velocity (σ w ) as turbulence parameter (Acevedo et al., 2009;Foken, 2017). The σ w threshold was set to 0.07 m s −1 as per out nighttime turbulence decay analysis and following Pirvulescu (2013) (Supplementary Figure S3).
Once Aubinet's filter was applied, concerns remained on EC flux measurements of nighttime ecosystem respiration, R eco . Under low wind conditions, CO 2 might be stored beneath the measurement plane within the canopy. Some of the stored CO 2 might be horizontally transported by advection. Typically, most of the remaining accumulated CO 2 is vertically released by turbulence bursts or upon turbulence onset after sunrise. On the other hand, early evening measurements are unlikely to be affected by advection (Van Gorsel et al., 2007), and can be considered representative of nighttime fluxes. San José et al. (2014) chamber measurements over a similar Llanos savanna ecosystem showed that nighttime soil respiration fluxes did not change much along the night, and that the daytime ones were comparable in magnitude although more variable along the day. We did not find any significant correlation between nighttime respiration fluxes and soil temperature. In absence of better information, we assumed that respiration did not vary within a day, and that this flux was identical to the nighttime one. Following San José et al. (2014) and as per our discussion above, R eco was assumed to be constant and estimated as the median of the early night (19:00-20:30 LT) fluxes corrected by storage. At our Equatorial site, sunset and sunrise times, and thus the length of the day, change by less than 30 min along the year. Therefore, we defined nighttime as the time period from 19:00 to 05:30. Solar radiation during this time period is well below 10 W m −2 throughout the year (see Supplementary Figure S4).
As valid nighttime data was sometimes unavailable, and to ensure the representativeness of the nighttime respiration median, a wing-variable time window (from 3 to 13 days) was used to guarantee having at least 30% of the possible data of the window. Also, in the processed CO 2 flux time series, the nighttime period (19:00-5:30 LT) data were substituted by the calculated, constant R eco . The uncertainty of the respiration and nighttime fluxes was estimated as the standard deviation (SD) of the data used to calculate the flux.

Seasonal and Daily Cycle of Carbon Fluxes
We used the Mitscherlich light response curve (LRC, Eq. 1; Falge et al., 2001) to investigate the savanna ecosystem CO 2 flux variability and controlling factors. The LRC fitting for each day included the daytime (06:00-17:00 LT) NEE data of that day along with that of the immediately previous and next days. For R eco , we used nighttime respiration of that day derived from its early-evening fluxes as explained above. Occasionally, R eco was directly estimated from the LRC model fitting when the early evening method failed to produce a R eco estimate.
Since PAR measurements were available at the site for a limited period of time only, global radiation (R g ) was used instead of PAR. These two variables were highly correlated (see Supplementary Figure S6), which allowed for the replacement of the original variables α (µmol CO 2 µmol (photons) −1 ) and PAR (µmol (photons) m −2 s −1 ) for β (µmol CO 2 J −1 ) and R g (W m −2 ), respectively. CO 2 fluxes during the 2016-03-15-2016-03-19 period were excluded from the LRC analysis as they did not follow the LRC expected behavior (low or positive NEE when negative was expected).
F c−sat is the asymptotic NEE (µmol CO 2 m −2 s −1 ) at light saturation and R eco is the dark ecosystem respiration (autotrophic and heterotrophic) that we estimated from the early evening fluxes. LRC was further refined to account for the NEE limitation due to stomatal closure, modeled using Eq. 2, where the vapor pressure deficit (VPD) was calculated from RH and T a , according to Campbell and Norman, (1998). We initially used a VPD 0 of 10 hPa following Lasslop et al. (2010). Our measurements indicated that a VPD 0 of 17.25 hPa is better suited for the Llanos savannas (see Supplementary Figure  S7). F c−sat and β were constrained to positive values.

Gap Filling
Besides statistical analysis, the main reason for filling gaps was the proper estimation of carbon accumulation or loss by the ecosystem over extended periods of time. Daytime and nighttime gaps were tackled separately. For nighttime, the gaps in the R eco time series (as described in )Section Nighttime Flux Estimation were filled as described in Section Seasonal and Daily Cycle of Carbon Fluxes. A five-method, gap-length dependent recipe was implemented for daytime gaps: (1) Linear interpolation for gaps shorter than 1 h (Falge et al., 2001).
(3) The method proposed by Reichstein et al. (2005) was used for similar meteorological conditions (SMC) (only the filling quality A was used). (4) A modification of the SMC method (SMCM) that included time and vertical moving windows and used the median as a central tendency estimator, was also used. Time windows of 7 or 15 days were used. The criterium for not increasing the time window length was to have at least 5 data points to calculate statistics. Vertical windows were 20, 40, 60, 80 and 100% of the following defined thresholds: R g 50 W m −2 , air temperature (T a ) 2.5°C, VPD 5 hPa. (Reichstein et al., 2005). Frontiers in Environmental Science | www.frontiersin.org June 2021 | Volume 9 | Article 673932 (5) A modification of the mean diurnal variation method (MDV, Falge et al., 2001) that used a moving time window to guarantee at least 30% of the possible data within the window, was also used. These procedures were run twice to fill long gaps.
Data gap filling methods were selected as per the following sequence: method 1); method 2); the least uncertain resulting from methods 3) and 4); and finally, method 5). The flux uncertainty was calculated for each method as follows: linear interpolation of the flux random error for method 1); residual standard error of the LRC model for method 2); and standard deviation of the data within the time window for methods 3), 4) and 5).

Meteorological Conditions
The meteorological conditions, particularly rainfall, were quite variable during the measurement period ( Figure 2). Our measurements began in the 2016 dry season, during which precipitation was 93% lower than the long-term average (34 mm/month for the 2005-2020 period), as estimated from a meteorological station located at ∼7 km east from our site (Hacienda Las Margaritas). Moreover, the precipitation during the 2016 dry season was 73% lower than the average for the five driest years of the 2005-2020 time period (9.3 mm/month). No rain (0 mm) was observed during January 2016. Only 5 mm were observed during February 2016. Conversely, during the 2016 wet season the precipitation was 8% higher than the long termaverage (259.9 mm/month). This likely implies that the ecosystem experienced high stress during both seasons, due to severe drought during the dry season, and higher-than-normal precipitation during the wet season.
Precipitation during the 2017 dry season (30.5 mm/month) was very close to the long-term average, but during the wet season was 22% lower than the long-term average (259.9 mm/month) ( Figure 2B). Figure 3 shows the half-hourly resolved, 3-consecutive-day CO 2 flux time series during three contrasting periods: 1) the beginning of the dry season ( Figure 3A), 2) wet season ( Figure 3B) and 3) extreme dry season ( Figure 3C). During the first period, the ecosystem served as a net carbon sink (daily NEE ranged from −0.77 to −1.40 gC m −2 d −1 ). Measured nighttime ecosystem fluxes were highly variable, likely due to intermittent turbulence CO 2 flushing. Nighttime fluxes estimated from periods of nearly constant flux early in the night, showed ecosystem respiration fluxes of ∼4.6 µmol CO 2 m −2 s −1 . During the second period, the ecosystem also served as a slightly less effective net carbon sink with daily NEE ranging from −0.74 to −1.03 gC m −2 d −1 . Ecosystem respiration during this period was slightly higher than respiration during the first period. During the third period, the ecosystem was a carbon source with daily NEE ranging from 0.49 to 0.82 gC m −2 d −1 . The ecosystem respiration (∼1.8 µmol CO 2 m −2 s −1 ) was almost a factor ∼3 smaller during this period than during the two other periods. Figure 3 also shows the fitted LRC model for each day with Eqs 1, 2, and without VPD adjustments (i.e. with Eq. 1 only). Although both models represent quite well the NEE daily cycle, the model that considers stomatal closure, based on VPD, represents better the NEE peak, especially during the beginning of the dry season. In contrast, the VPD-adjusted model underperformed the simpler LRC model, i.e., without stomatal closure, for most part of the late 2016 dry season. Figure 4 shows the weekly-averaged CO 2 fluxes, ecosystem respiration (R eco ), daily carbon accumulation (i.e., cumulative NEE), and precipitation. There were significant differences between the monthly NEE of 2016 and 2017. Overall, the largest continued carbon uptake occurred during the beginning of the dry seasons, mainly during December and January. NEE interannual differences among late dry seasons and the wet seasons appear to be highly influenced by differences in precipitation distribution between the two years.

Seasonal and Daily Cycle of Carbon Fluxes
During 2016, as the dry season progressed with almost no rain between January and middle March, a systematic decrease in NEE was observed. During February 2016, the ecosystem became a carbon source and maintained that condition until the end of april, almost one month after the beginning of the wet season. After the first precipitation events in middle March, there was an increase in ecosystem respiration that resulted in NEE levels within 1.5-2.5 gC m −2 d −1 during the two weeks that followed the precipitation events. Between May and September 2016, the ecosystem acted as a carbon sink during most of the weeks, but had much less carbon uptake than during the same period of 2017. During late September, the ecosystem decreased its carbon uptake, turned into a carbon source, and maintained that condition until middle November. No data was obtained during the period October 8-December 8 due to logistic problems. Data presented for that period corresponded to estimations based on gap filling procedures taking into account the data trends near that gap. During 2017, the ecosystem was a carbon sink during most of the weeks. Variability in ecosystem respiration significantly contributed to NEE variability. The savanna ecosystem acted as a small carbon source during 2016 (NEE 35 gC m −2 ), and as a significant carbon sink during 2017 (NEE −273 gC m −2 ; from 2017-01-01 to 2017-12-14, i.e. 17 days short). Figure 5 shows the 2016 and 2017 monthly averaged daily CO 2 flux and VPD profiles. The monthly precipitation data for each year (P 16 and P 17 ) and the accumulated precipitation since January of the year (P C16 and P C17 ) are also shown. Large differences are observed between the 2016 and 2017 CO 2 flux profiles during the dry season with substantially lower CO 2 uptake during the drier months of 2016 (January-February) compared to 2017. The anomalous CO 2 fluxes during the 2016 dry season are also consistent with larger vapor pressure deficits that year compared to 2017. As the dry seasons progressed, the carbon uptake decreased in spite of decreasing VPD in march. Conversely, precipitation during april was 32% greater in 2016 compared to 2017. This resulted in significantly higher respiration in april 2016 compared to april 2017. At the beginning of the wet season, although the diurnal carbon uptake increased, there was also an increase in ecosystem respiration that resulted in positive carbon fluxes during april. Precipitation during the wet season was higher in 2016 than in 2017. With the exception of May and September, all the wet season months were wetter in 2016 compared to 2017. The 2016 annual accumulated precipitation ended up being 428 mm higher in 2016 compared to 2017. As the wet season progressed the differences between 2016 and 2017 diurnal carbon fluxes started to increase, due to a systematic decrease in the 2016 ecosystem carbon uptake.

Environmental Controls on CO 2 fluxes
The results above clearly show that the savanna CO 2 fluxes and cumulative NEE are highly dependent not just on the total annual and monthly precipitation, but particularly on its temporal distribution, and that this dependence is a nontrivial one. While precipitation in 2017 was 19% lower than in 2016, the net C uptake in 2017 was far higher than in 2016. To understand the impact of drought severity on the functioning of the savanna ecosystem, we first explored the CO 2 flux dependency on accumulated precipitation over variable accumulation period lengths. No clear correlation Frontiers in Environmental Science | www.frontiersin.org June 2021 | Volume 9 | Article 673932 8 was found. Instead, we found a very compact dependence of NEE on an indicator of nominal drought. This indicator was defined as the consecutive number of days without effective precipitation. We found that a minimum precipitation of 5 mm/day minimizes the spread of the NEE-nominal drought indicator relationship. The drought indicator is thus the consecutive number of days with less than 5 mm/day of precipitation. Figure 6 shows the daily CO 2 fluxes as a function of this indicator. It clearly shows that the carbon uptake decreases monotonically but non-linearly from 10 days of nominal drought on. The ecosystem becomes a carbon source after ∼40 days of nominal drought.
The effects of drought on the ecosystem were also evident on the light-response curve parameters. Figure 7 shows the NEE at light saturation (F c−sat ) and the photosynthetic efficiency (β) as a function of the "drought duration," i.e., the number of consecutive days with less than 5 mm/day of precipitation, starting at 10 days. Both, β and F c−sat , decreased exponentially with the drought duration. β decreased by a factor of ∼3 from about 0.075 µmol CO 2 J −1 at 10 days of drought to 0.025 µmol CO 2 J −1 at more than 60 days of drought. The effect on F c−sat , the NEE at light saturation was even more pronounced. It decreased from 20 µmol CO 2 m −2 s −1 at 10 days of nominal drought to about 2 µmol CO 2 m −2 s −1 at more than 60 days of nominal drought. The modeled parameters were restricted to minimum values to properly represent the basic functioning of the vegetation.
There was also high variability in ecosystem respiration. During the dry season, R eco displayed a statistically significant dependence on VPD, with lower respiration rates at larger VPD values (Figure 8). This dependency was represented by a logistic curve. The R eco dependence on drought could be better described by soil water content (SWC), but this variable was unavailable for much of the measurement campaign. These interrelations are further confirmed by the good correlation between VPD and SWC during the dry season (Supplementary Figure S8).

Effect of Meteorological Factors on Carbon Fluxes
We observed high variability in monthly carbon fluxes and ecosystem respiration during our two years of measurements. During the two wet seasons and the early dry seasons, i.e. when the ecosystem had enough water, the maximum hourly NEE levels were rather similar, and the carbon fluxes were controlled by solar radiation. As the dry seasons progressed, long periods of low precipitation were experienced leading to decreasing water availability. These reduced the ecosystem gross primary productivity (GPP) but also its respiration. The reduction of respiration was insufficient to counterbalance the GPP reduction, thus the drought resulted in a systematic decrease of NEE, which seasonally turned the ecosystem into a carbon source. Light response curves properly described the savanna carbon fluxes at the daily scale. This representation of diurnal fluxes by the LRC model was improved when F c−sat was made parametrically constrained by VPD. We also found that, besides the VPD effect on carbon fluxes at the daily scale, the ecosystem response was also highly dependent on the "recent" (i.e., monthly) history of precipitation, which determines the amount of water actually available for its functioning. The dependence of carbon fluxes, NEE at light saturation (F c−sat ), and photosynthetic efficiency (β) on the number of days of nominal drought, reveals that 1) "drought" is better defined as precipitation below an effective precipitation value (∼5 mm/day) rather than the total absence of precipitation; and 2) the water stress effects on the ecosystem functioning show up after ∼10 days of nominal drought. The ecosystem seems to function normally at shorter drought periods; 3) the ecosystem fluxes become almost fully determined by the drought duration during an extreme dry season.
Ecosystem respiration was highly variable and dependent on VPD, particularly during the dry season. At the end of the dry season, heterotrophic respiration was inhibited by the very low soil water contents. Then, we observed a significant increase in ecosystem respiration during the dry to wet season transition. This is likely explained by the occurrence of CO 2 respiration FIGURE 7 | Dependence of (A) NEE at light saturation (F c−sat ), (B) photosynthetic efficiency (β) on drought duration. Coefficient uncertainties are ±1 standard deviation (68% CI). bursts upon moistening of dry soil, known as the Birch effect (Birch and Friend, 1956;Griffiths and Birch, 1961). This is a complex phenomenon that results from rapid stimulation of microbial activity, reproduction and growth from sudden increase in soil water (Ma et al., 2012). The Birch effect has been previously reported in African (Williams et al., 2009) and other Orinoco savannas (San José et al., 2014). Based on these findings, we developed a simple NEE lightresponse curve (LRC) model, all the variables of which can be calculated from basic meteorological measurements (precipitation, global radiation, temperature, and relative humidity). Functions of these variables properly represent the effects of stomatal closure and drought duration on the ecosystem NEE and respiration. After testing various modeling variants, we found that the incorporation of the drought effects into the model was essential to properly represent the ecosystem C accumulation. We also found that, if ecosystem functioning was not limited by water availability, its C exchange could be reasonably well-represented by a single set of LRC parameters, i.e., for the non-drought conditions experienced during the wet and early dry seasons. An ad hoc model is then used to incorporate the water limitation effects on NEE, using the drought duration models explained above.
The proposed LRC-drought model used two sets of parameters. The first set of parameters (F c−sat 11.42 μmol CO 2 m −2 s −1 , β 0.0679 µmol CO 2 J −1 , and R eco 4.24 µmol CO 2 m −2 s −1 ) were derived from and applied to non-drought conditions (November 2016 to December 2017), which represent the usual, non-drought conditions of the ecosystem. The second parameter set allows simulating the ecosystem under nominal drought conditions (defined as more than 10 consecutive days with precipitations under 5 mm/day). These parameters were obtained by fitting of 2016 severe dry season data to the droughtduration exponential decay curves shown in Figure 7. These parameters were applied only if they led to lower uptake that those used for non-drought conditions. R eco was calculated as a function of VPD (Figure 8).
It is very important for a C flux model, either simple, as the one described above, or complex as DNDC (DeNitrification- DeComposition model) to be as bias-free as possible. Otherwise, even a small bias at the hourly or daily time scales may lead to a significantly biased estimation of the C accumulation over extended periods of time (years to decades). Therefore, although the droughtno-drought model presented above reproduces the hourly C fluxes quite well, it is necessary to assess its ability to reproduce the longer-term C accumulation. Figure 9 shows the measured carbon accumulation (black trace), along with several estimated variants: 1) "Daily LRC":LRCs fitted on a daily basis using VPD when was possible (red trace); 2) "Generalized LRC": single set of parameter for LRC ignores drought and stomatal closure effects (blue trace); 3) "Generalized LRC with VPD": single parameter set for LRC incorporates stomatal closure effects (VPD function) but ignores drought effects (orange trace); 4) "Generalized LRC with drought effects": drought-no-drought LRC-based model explained above incorporates drought instead of stomatal closure effects (green trace). None of the generalized LRC models (i.e. models 2-4) was able to properly reproduce NEE during the 2016 dry to wet season transition nor during the 2016 extreme precipitation period (see Figure 9A). This suggests that soil processes play a major yet complex role in the ecosystem C exchange during rapid transitions to high soil humidity conditions. Due to their complexity, these periods were excluded from the modeled C accumulation comparison. Incorporating stomatal closure effects into the light-response curve ("Generalize LRC with VPD") improves the ability of the model to simulate the accumulated carbon dynamics, particularly during the dry season; but the model has a significant bias and poor performance. The "Generalized LRC with drought effects" model constitutes a significant improvement with both a better dynamic representation and a much smaller bias. This model was able to reproduce carbon accumulation during the extreme 2016 dry season, the 2017 dry season, and the 2017 wet season.

Carbon Storage and Climate Change Mitigation
There were large differences in net carbon fluxes between 2016 and 2017. During 2016 the ecosystem acted as a net source (annual NEE of 35 gC m −2 ) during a year characterized by a severe dry season followed by extreme precipitation during the wet season. During 2017, the savanna served as a carbon sink with an annual NEE of −273 g C m −2 . Although annual precipitation during 2017 was lower than during 2016, the average precipitation during the 2017 wet season was high enough to provide ecosystem water requirements, which resulted in high ecosystem productivity. This change of carbon sink to carbon source was also observed in a grazed savanna ecosystem in South Africa (Räsänen et al., 2017), where the annual NEE changed from -85-139 g C m −2 and the magnitude and direction of NEE was reported to be related to water stress and precipitation distribution rather than annual precipitation. Specifically, NEE was found to be highly dependent on the number of days when soil moisture was higher than 7%. The annual C uptake measured during 2017 was similar to the annual NEE of Campo Sujo savannas (−288 g C m −2 ), an open shrub savanna with scattered shrubs similar to the Llanos High Plain savannas (Santos et al., 2003). Th annual C uptake reported for the Venezuelan Llanos were lower than our measured 2017 NEE (−6, −116 and −139 gC m −2 yr −1 for an herbaceous savanna, a tree savanna and a woodland savanna, respectively; San José et al., 2008). The savanna vegetation in our study could be considered similar to the herbaceous and tree savanna measured by San José et al. (2008). However, the dry season is longer in the Venezuelan savannas than in Colombian High Plains and Campo Sujo, which may explain the lower carbon uptake in the Venezuelan savannas.
The Colombian Llanos High Plains (9.7 million hectares-Mha) have undergone significant transformation. Etter et al. (2010) estimate that 21% of the High Plains savannas had been cleared for agriculture by 2010. The same study also estimates that an additional ∼10% of the High Plains savannas would be cleared by 2020. The clearing of the Dissected High Plains savannas was estimated to be lower (∼10% until 2010). Overall, more than 5 Mha of savannas are currently reported as being used for cattle grazing in Meta and Vichada, the two Colombian Provinces (Departmentos) that host the High Plains region. Colombia's Nationally Appropriate Mitigation Actions (NAMA) for the livestock sector have set an emissions reduction goal of 11.15 Tg CO 2 -eq/year (11.15 Mton CO 2 -eq/yr) by 2030. One of the mitigation measures is land restoration. During our 2 years of observations, including a very anomalous one, the native savanna accumulated ∼2.5 ton C/ha. Considering this carbon sequestration capacity, we feature a simple scenario in which 10% of the area currently used for cattle grazing is left alone, so it naturally restores. We estimate that this action alone might sink ∼1 Tg CO 2 -eq/yr, i.e., 10% of the NAMA's emission reduction goal by 2030. This estimation is based on the assumption that a 10-year period would be composed of two cycles as follow: 1-year cycles of benign meteorological conditions, meteorologically and C accumulationwize similar to 2017, one year of extreme weather, similar to 2016, and three 1-year periods of mild meteorological conditions with 70% of the carbon accumulation measured in 2017. We also considered that during the restoration process, the carbon uptake capacity of the ecosystem would increase from 10% of the abovementioned capacity to 55% during the 10th year (a 5% per year increase rate). Based on these assumptions, the restored ecosystems could uptake 1.8 Tg CO 2 -eq/year by 2050.

CONCLUSION
This research reported carbon dioxide flux measurements in a Northern South America tropical naturally restored savanna ecosystem located in the Llanos High Plains, using the eddy covariance technique. Besides quantifying the CO 2 flux in this ecosystem (not previously reported in the literature), we identified the main relationships between fluxes and meteorological conditions.
The savanna ecosystem annual carbon uptake was highly dependent on climate conditions, particularly the length and intensity of the dry and wet seasons. Despite receiving more precipitation during 2016, the ecosystem went from a minor C source in 2016 (annual NEE 35 g C m −2 ) to a significant C sink in 2017 (annual NEE −273 g C m −2 ). The savanna poor performance in 2016 was mainly due to 2 factors: 1) increasing C release during its long and severe dry season; and 2) enhanced respiration at the beginning of the wet season (Birch effect). This indicates that the C performance of Orinoco High Plains savannas is highly sensitive to uneven distributions of rainfall. This savanna ecosystem is thus especially vulnerable to the intensification of climate variability, which might turn savannas into net C sources.
Our results suggest that this savanna ecosystem became waterstressed when precipitation went below ∼5 mm/day. After ∼10 days under this nominal drought condition, the savanna ecosystem C fluxes became controlled by the nominal drought duration rather than by solar radiation. The nighttime ecosystem respiration was also dependent on drought conditions, however its dependency was better described by vapor pressure deficit rather than by drought duration.
We developed a simple NEE light-response curve (LRC) model, all the variables of which can be calculated from basic meteorological measurements (precipitation, global radiation, temperature, and relative humidity). Functions of these variables properly represent the effects of drought duration on ecosystem NEE and respiration. The model was able to adequately reproduce carbon accumulation during the extreme 2016 dry season, the 2017 dry season, and the 2017 wet season. The model did not represent well the dry to wet season transition, when the ecosystem was under extreme wet conditions, thus likely flooded. Long-term CO 2 flux measurements, along with ancillary variable measurements, are required to better understand and represent the ecosystem respiration, particularly during the dry to wet season transition, over which heterotrophic respiration may become dominant.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
RJ conceived and managed the project. LM-R, AH, and RJ planned the observations, developed the models and analyzed the data. LM-R, AH and NR-H conducted the field measurements. LM-R and RJ prepared the manuscript with significant contributions from all the authors.

FUNDING
This research was funded by Colciencias (Colombian Administrative Department of Science, Technology and Innovation), through grant number 1101-569-35161 for the project "Atmospheric emissions and impact of land use change and intensive agriculture on air quality in the Colombian Llanos."