Biomass burning-agriculture coupling in the Orinoco savannas—Particulate matter emission scenarios

The Colombian Orinoco savannas (254 thousand km²), also known as Orinoquia or Llanos, have been steadily transformed into pastures for more than a century, and since the 1990s, into commodity crop intensified production. The cropland area expanded at 12% yr⁻¹ during the 2007–2018 period (65% larger than in 1996–2007). Yet, we estimate that cattle ranching occupied ten times more area (34%) than cropland (3.2%) in 2018. The rest of Orinoquia, including indigenous reservations and protected areas, was in a semi-natural state, although also exposed to seasonal fire. The three main crops, oil palm, corn, and rice (72% of the sown area in 2017), accounted for 68% of the expansion, with permanent crops expanding two times faster (18% yr⁻¹) than short-cycle crops. An extrapolation of trends indicates that the cultivated area will double by 2040 (reaching 20 thousand km²), with oil palm as the dominant crop. Satellite measurements show that 7% of Orinoquia burned every year during the 1997–2016 period, yet with large spatial and interannual variations (±26%), and significant decrease trends (up to −4% yr⁻¹). Up to 40% of the burned area (BA) interannual variability was linked to irregular rainfall and drought. The areas with the larger fractional BA were also those with the least fractional cropland cover. A model developed to describe this coupling, along with rainfall and other effects, successfully explained most of Orinoquia’s BA variability (r² = 0.93). The fitted model indicates that each sown hectare reduced the BA by 0.17 ha. This model predicts that the combination of cropland expansion and independent BA decline will lead to a fourfold reduction of Orinoquia’s BA by 2040 referred to 1997. Orinoquia’s crop production generated 3 Gg of PM10 (particulate matter <10 µm) in 2016, mostly from short-cycle crops, while biomass burning generated 57 Gg, i.e., 95% of the combined emissions. These are expected to halve during the 2017–2040 period, despite an 83% increase in crop production emissions, as total and seasonal emissions will remain controlled by biomass burning. Such a large pollution burden reduction should have tremendous positive impacts on public health in Orinoquia and the Andes.

1 Introduction

The increase in population, per capita income and consumption, and the market economy, along with the search for energy independence through the production of biofuels, has inevitably led to the incorporation of new lands into agriculture and its intensification (Lark et al., 2015). In the current globalized economy, these land use changes, particularly in developing countries, are frequently aimed at the supply of global demands, often driven by urbanization and economic expansion. Many countries, particularly in Southeast Asia, the Middle East and North Africa, some of them densely populated and rapidly industrializing, have practically run out of new arable land. FAO (2018) estimates that harvested areas will globally increase by 238 million hectares [Mha] between 2012 and 2050, mostly in Latin America and the Caribbean (20% increase) and Sub-Saharan Africa (62% increase). South America has been the region with the largest land use change for commodity production during the last decades. Pasture, cropland and tree plantation areas increased by 13%, 10% and 288%, respectively, from 1985 to 2018 (Zalles et al., 2021). Based on legal, institutional and economic reforms, infrastructure development, and better agricultural practices, including the use of agricultural inputs, Brazil, Argentina and other South American countries have dramatically increased their agricultural production since the 1990s (Schnepf et al., 2001; Graesser et al., 2018). Brazil, Argentina, Paraguay and Bolivia are currently the second, third, sixth and eighth largest soy producers in the world, respectively. These four countries jointly produced more than half of the 333 million tons of soybean consumed worldwide in 2019 (https://www.fao.org/faostat/en/#data/QCL). One very concerning factor is the flattening of the agricultural yield increase rate, which in 1985–2005 was only one third of that achieved in 1965–1985 (Foley et al., 2011). Moreover, global temperature increase and regional climate change are expected to reduce crop productivity and increase plant pathogen virulence and infection rates, further reducing crop yield (Anderson et al., 2020). This could hamper the reach of agriculture intensification.

Similar intensive agriculture development processes are taking place in the Colombian Orinoco River basin, also known as the Colombian Orinoquia or Llanos, whose agricultural transformation processes (García Durán, 2009; Etter et al., 2010) are similar to those carried out in the Brazilian Cerrado since the 1960s (Schnepf et al., 2001; Klink and Moreira, 2002), particularly the transformation of pastures for extensive cattle ranching into intensive agriculture (Quezada et al., 2022).

The clearing of savannas and forests for pastures and the replacement of pastures by crops (Quezada et al., 2022), along with other land use changes and the ensuing agricultural activities, cause multidimensional environmental impacts (Foley et al., 2011). These include biodiversity reduction and loss of habitats, unsustainable and conflicting water use, contamination of water and soils, including eutrophication (Tilman, 2001), changes in albedo, evapotranspiration, carbon inventories in soils, and emissions of greenhouse gases (GHG) (Silva et al., 2011) and air pollutants, including particulate matter (PM) (Aneja et al., 2008, Aneja et al., 2009; Zhang et al., 2008). Lumped together, agriculture (6.2 ± 1.4 Pg CO₂eq year⁻¹) and land use change and forestry (5.8 ± 2.6 Pg CO₂eq yr⁻¹) emitted a quarter of the anthropogenic long-lived GHGs between 2007 and 2016 (IPCC, 2019).

Here we estimate and analyze atmospheric emissions from agricultural activities and burning of savannas, and explore their coupling. We focus on emissions of PM10, i.e., of suspended particulate matter with aerodynamic diameter lower than 10 µm that can pass the larynx and reach the lungs. Agricultural emissions occur infrequently upon land use change and recurrently once crop production is established (Sharratt and Auvermann, 2014; Maffia et al., 2020). We focus on operational emissions from crop production, disregarding those from cattle production for reasons that are explained below. Routine emissions are proportional to either sown or harvested area, but depend heavily on several factors, including land preparation practices and harvest frequency.

Land clearing, preparation, tillage and harvesting emit dust (Sharratt and Auvermann, 2014; Sharratt and Van Pelt, 2020), i.e., crustal (soil) particles, that are mostly coarse (>2.5 µm). These and other pollutant emissions depend strongly on agricultural practices, soil humidity, edaphological factors, crop type and weather conditions. Emissions from dry soil traditional tillage (up to ∼35 kg PM10 ha⁻¹) can be three orders of magnitude higher than on wet soils (Wang et al., 2010). These aerosol flux levels are comparable to those in polluted urban (Pachón et al., 2018) or industrial areas. Finer particles, usually measured as PM2.5 (aerodynamic diameter <2.5 µm), are also emitted by operations on soils, diesel combustion for mechanization, and burning of biomass, including for land preparation and pre-harvest and crop residue management practices (Mateus-Fontecha et al., 2022). Diesel and biomass burning particles are mostly carbonaceous (soot and organic) and frequently toxic.

Agriculture intensification generally requires increased inputs per unit area, e.g., capital, agrochemicals, and energy (including mechanization), compared to traditional agriculture. This implies that also pesticides, nitrogen oxides (NOx), ammonia (NH₃), reduced sulfur compounds (e.g., H₂S) and odors (e.g., organic acids) are released to the atmosphere (Tilman, 2001). Some of these primary pollutants are also precursors of secondary pollution, e.g., secondary inorganic PM2.5 can form from NH₃. The estimation of these emissions and the resulting generation of secondary inorganic and organic aerosols (SIA and SOA, respectively) are beyond the scope of this paper.

Orinoquia, which a few decades ago was mainly an extensive cattle ranching region, is today experiencing rapid agricultural expansion and intensification. The cultivated area and agricultural production increased by more than 50% between 1996 and 2007 (Viloria de la Hoz, 2009). At the beginning of the 2010s, agricultural and industrial investors acquired over 2,000 km² (200 thousand hectares [kha]) and invested more than 1.4 billion dollars in intensive agriculture projects aimed at commodity crop production, including oil palm and sugarcane as biofuel feedstocks and corn and soybean for pig and poultry farming (Dinero, 2011). The characteristics of the region and the Colombian public policy that promotes the development of Orinoquia’s High Plains (DNP, 2014) and the production of biofuels (DNP, 2008), suggest that this rapid transformation pace will be sustained in the decades to come.

Romero-Ruiz et al. (2010) estimated that 18.9% of the area transformed to livestock and agriculture in the High Plains in 2007 was from forest ecosystems. Also, Etter et al. (2010) estimated that the rate of transformation of native savannas into pastures, tree plantations and crops would reach 2% per year in the 2020s (equivalent to >2,000 km² yr⁻¹). Unless riparian forests and savannas are effectively conserved and agriculture conservation practices are implemented (Quezada et al., 2022), these transformations will turn Orinoquia from a small GHG sink into a significant source (Etter et al., 2010; Morales-Rincon et al., 2021). Similarly, air pollutant emissions are expected to significantly increase but the information on emissions and air quality is very limited (Hernandez et al., 2019). Also, no projection has been made so far on changes in air pollutant emissions due to land use changes and intensive agriculture in Orinoquia. Although our focus is on atmospheric impacts, it must be stressed that land use and land use changes cause multidimensional impacts, including on Orinoquia’s mega-biodiversity (Lasso et al., 2010).

In a previous investigation, we found that biomass burning (BB) was by far the largest PM10 emission source in the Colombian Orinoquia (Hernandez et al., 2019). At the peak of the dry season, the transport of BB aerosols from Venezuela and within Colombia triggered multi-day pollution episodes with very high PM10 concentrations in two Orinoquia cities (Hernandez et al., 2019; Rodríguez-Gómez et al., 2021). Likewise, it has been demonstrated that these BB plumes were frequently transported to the Colombian Andes Eastern Cordillera (see Figure 1), also causing episodes in Bogota (Mendez-Espinosa et al., 2019; Ballesteros-González et al., 2020; Rincón-Riveros et al., 2020) and other Andean cities.

FIGURE 1

FIGURE 1. A map of the Orinoco River savannas, also known as Orinoquia or Llanos, depicting landscapes, provinces and cities (in Colombia only). Roughly, the Meta River splits the Colombian Orinoquia into the High Plains to the south and a seasonally floodable region to the north.

Etter et al. (2010) hypothesized that agriculture expansion in Orinoquia would cause “reduction in fire frequency”. Based on satellite-derived products and after correcting for precipitation variation effects, Andela et al. (2017) reported that the number, mean area of fires and the aggregated burned area (BA) were decreasing worldwide, mainly due to agriculture expansion and increasing population densities. The estimated reduction trend for South America ranged from −1.4% yr⁻¹ to −2.7% yr⁻¹, depending on the satellite product (from −0.5% yr⁻¹ to −1.8% yr⁻¹ after adjusting for precipitation variability effects). They proposed a conceptual model that accounted for the agriculture-fire coupling complexity, in which BB increases or decreases depending on the level of sophistication and investment in agriculture. Thus, cropland expansion could have contradictory effects on PM emissions: while crop production PM10 emissions are relatively small compared to those of BB, the rapid expansion of intensified agriculture might suppress BB. Subsequent studies have confirmed the BB trend, its beneficial impact on carbon uptake (Arora and Melton, 2018), and shown progress towards BB area prediction using mechanistic (Pechony and Shindell, 2009; Knorr et al., 2016; Lasslop and Kloster, 2017; Rabin et al., 2017; Venevsky et al., 2019) and neural networks (Joshi and Sukumar, 2021) models. Here we use the available data to estimate the coupling and combined emissions scenarios.

Using time series of sown and harvested areas, precipitation and BA, we aimed at answering the following research questions: What are the sown area trends and their variability in the Colombian Orinoquia, disaggregated by crop and region? How time and spatially variable and dependent on precipitation are the burned areas? How to tell apart the effects of cropland expansion on BB from those of other drivers in the Colombian Orinoquia? What are the combined effects of crop production and BB on PM emissions when their coupling is considered? For this purpose, we estimated and extrapolated the cropland expansion, implemented a detailed emission model, and estimated the resulting cropland PM emissions. Then, we developed a model that enlightens the agriculture-BB coupling, looked for and analyzed precipitation-BA and then sown-burned areas relationships, and calculated BA and emissions scenarios for the Colombian Orinoquia based on variations of the coupling coefficient. It should be noted that, in the tropics, the same physical area can be sown over 2-3 cycles per year. The available information is not detailed enough for precise discrimination between physical and sown area, nor does it allow determining what percentage of cropland expansion is due to displacement of the agricultural frontier, i.e., of the interface between agricultural areas and areas with no direct economic purpose, typically native savanna or riparian forests. Also, available data on livestock production land use are currently not accurate enough so they were not explicitly included in our analyses.

2 Theoretical foundations

2.1 Agriculture expansion modeling

The total area available and suitable in a region for sowing a given crop is limited. This includes areas already used for agriculture and “new” land suitable for use change, e.g., from degraded or improved pastures to commercial crops (Castiblanco et al., 2013; Quezada et al., 2022). The limitation is also determined by meteorological variables, particularly rainfall and photosynthetic active radiation (Morales-Rincon et al., 2021), soil and terrain characteristics, location and infrastructure (Etter et al., 2010; Castiblanco et al., 2013), technology and security. If the maximum area that can be sown (i.e., the area asymptotic or saturation limit, also known as “carrying capacity”) is large, as it is in Orinoquia, the initial expansion rate is expected to be slow and linear, then faster around the expansion mid-point (exponential phase), and afterwards again slow as the sown area approaches its carrying capacity. This sown area ($A_S$) sigmoidal curve behavior is appropriately represented by the logistic model (Carrillo and González, 2002),

$A_S=\frac{A_{S,\infty }}{1+e^{-\left(t-t_{1/2}\right)/\tau }}(1)$

where $A_{S,\infty }$ is the carrying capacity, $t_{1/2}$ is the expansion midpoint time, i.e., time at which $A_{S,\infty }/2$ is reached, and $\tau$ is the saturation time constant, which is proportional to the exponential phase duration and inversely proportional to its normalized rate.

Fitting sown area to the logistic model at the beginning of the expansion process, when only a small fraction of the carrying capacity has been developed, is likely to produce misleading, highly biased parameter estimations and forecasts, particularly if sown areas are highly variable. In the initial expansion phase, a more robust representation and forecast can be obtained by fitting a simple normalized linear model,

$A_S/A_S^{\ast }=a_0+a_1\left(t-t_{ref}\right)(2)$

where $A_S^{\ast }$ is a suitable-for-normalization sowing area (for instance the sown area at the beginning of the time series), $t_{ref}$ is a suitable reference time, $a_0$ is the fitted normalized area at the reference time, and $a_1$ is the normalized linear growing rate referred to $A_S^{\ast }$. Eq. 2 was individually fitted to the sown area of each crop and crop group (as explained below), the sown area in the four Orinoquia provinces, and also used to estimate the observed BA decreasing rate using the mean BA as a normalization parameter. The regressions of crops and BAs were all referred to the year 2000 ($t_{ref}=2000$).

2.2 Estimation of crop production and biomass burning emissions

We included only direct crop production emissions in our inventory, i.e., only those directly attributable to crop production in the sown area. This excluded, e.g., emissions due to transportation of crop inputs and harvest products to and from farms. The full crop lifecycle was considered in the calculations to account for both recurring and infrequent emissions. Three activity groups were included: A) Land preparation and sowing (infrequent for permanent crops; rather recurrent for short-cycle (SC) crops); B) Harvesting (recurrent—only after long maturation periods for some permanent crops); and C) Machinery (diesel and aerial spraying aviation gasoline combustion emissions—recurrent). Land preparation emissions were thus “diluted” over their lifecycle, which, e.g., for oil palm, implies diluting month-long emissions over ∼25 years. Each activity group was further disaggregated into their constituting processes, e.g., land preparation and sowing included land planning, subsoiling, discing, harrowing, weeding, and mechanized sowing operations, the emissions of which depended on the number of passes applied. We first calculated process-based emission factors (EFs) for each of the 69 crops, and then, based on their similarity, we grouped similar crops and calculated weighted crop-group equivalent EFs. For this, we prepared a binary matrix that described the activities required by each crop. The information on crop-dependent practices and their frequency was obtained from crop production professionals and practitioners with experience in Orinoquia. The activities in the groups A and B mainly generate dust, i.e., mostly coarse mode crustal particles (>2.5 µm). The machinery group mainly generates fine mode (<2.5 µm) carbonaceous particles. Further details are provided in the Supplementary Material. The resulting crop group EFs are presented in Table 1.

TABLE 1

TABLE 1. Crop production emission factors disaggregated by activity and crop group.

Biomass burning emissions were bottom-up calculated by multiplying the mass-based particulate matter (PM) emission factor by the biomass burned. This last one was estimated as the product of burned area, biomass density, and combustion efficiency. Mass-based emission factors, biomass density and combustion efficiency, all depend on the type of vegetation burned. For Orinoquia savannas, the biomass density (“fuel load”) ranges from 6 to 7 metric tons (ton) per hectare (ha) and the combustion efficiency from 0.70 to 0.85 (Etter et al., 2010). In their compilation of emission factors, Akagi et al. (2011) report a savanna burning PM emission factor of 7.17 ± 3.14 g (kg dry matter)⁻¹. This paper also indicates that most of the BB emissions are in the fine mode, i.e., BB PM2.5 and PM10 emissions are nearly identical. The resulting area-based emission factor (EF) for PM10 is 36.12 kg (burned ha)⁻¹.

2.3 An agriculture-biomass burning coupling model

Vegetation fires require fuel and ignition, and once set, they can spread or be suppressed. Most fire predictive models (Rabin et al., 2017) estimate fractional BA (e.g., Knorr et al., 2016) (or number of ignitions per unit area—Venevsky et al., 2019) by non-linear factoring available vegetation, flammability, and human population density. Some models use adjoint biogeochemical models to estimate available biomass. As modeled, anthropogenic fire frequency peaks at mid population densities as ignitions are infrequent with sparse populations and fire is actively suppressed near urban and suburban settlements. Despite their complexity, only a few models explicitly incorporate cropland fraction as a strong fire suppression factor, most notably JSBACH-SPITFIRE (Lasslop and Kloster, 2017). The conceptual model below was developed to better understand the linkages between land use, land-use heterogeneity, and BB, rather than for vegetation fire mechanistic prediction.

Orinoquia’s traditional cattle ranchers have viewed burning as a native pasture “renovation” practice (López-Hernández et al., 2005; García Durán, 2009; Etter et al., 2010). As this is a common practice, and cattle occupies a large fraction of Orinoquia, it is reasonable to assume that a significant portion of the cattle grazing native pastures are burned on a regular basis. Unfortunately, there are presently no statistics available on this practice. On the other hand, due to the perceived high value of commercial crops and the soils that sustain them, particularly if conservation practices are in place, areas sown with improved pastures or commercial crops are expected to be less susceptible to burning than other areas, e.g., those dedicated to extensive cattle ranching or native savanna. We thus hypothesize that the net effect of cropland expansion in Orinoquia has been a reduction of the burned area (variable $A_B$), as observed in other regions with intensive and mechanized agriculture (Rosan et al., 2019).

As per the considerations above, only three lumped land uses were included in the model: 1) cropland (subscript S, which stands for sown); 2) livestock production (subscript L), mainly extensive cattle ranching with low stocking rates; and 3) areas mainly in a near natural state, protected or not, including indigenous people reservations (subscript 0). The total area, $A_{tot}$, is thus $A_{tot}=A_0+A_L+A_S$, where $A_S$, $A_L$, and $A_0$ are the cropland, livestock and “natural” areas, respectively. The annual fractional BA of each of these k-th land use types, $f_k$ [yr⁻¹], later referred to as normalized burning rates, can be calculated as the product of the number of annual fire events per unit area, $r_k$ [km⁻² yr⁻¹], a random variable, and the burned area per event, $s_k$ [km²], a highly variable parameter, i.e., $f_k=r_k{s}_k$. Thus, the annual burned area expected value, $X\left(A_B\right)=\overline{A_B}$ (hereafter denoted as $A_B$), is simply the sum of the areas on which burning can occur, $A_k$, “weighted” by their annual fractional BA, $A_B={\displaystyle \sum }_k{A}_k{f}_k$.

The high value of crops and soils also implies that farmers actively suppress fire not just on their farms but also around them. This perimeter fire suppression notion is based on various conversations, including with a representative of the agricultural community fire brigade at Puerto López (Meta, Colombia). The brigade, which started with just four people, is currently about 60 strong and supported by local and national authorities, including the Civil Defense and Emergency Committees. Besides setting 4–6 m wide buffer land strips (firebreaks) around their croplands, brigade farmers communicate and suppress any proximity fires. Moreover, some farmers now deploy ISO water tanks near their farm fields and own water tankers. Overall, this suggests that cropland fire suppression spreads to neighboring areas, i.e., livestock and “natural” areas near cropland areas are less likely to burn than their inland counterparts. Therefore, $A_0$ and $A_L$ are further divided into their inland and perimeter components (subscripts I and p, respectively), so that $A_0=A_{0,I}+A_{0,P}$, and $A_L=A_{L,I}+A_{L,P}$. The cropland-perimeter surrounding area, $A_{k,P}$, can be easily estimated as the area defined by a fire suppression width, $w$, around the cropland perimeter, $P$ ($w$ is likely far wider than the above-mentioned physical firewall width). See Supplementary Figure S1 for a schematic. The inland area is simply the difference, $A_{k,I}=A_k-A_{k,P}$. The functional relation between the cropland area, $A_S$, and its perimeter, $P$, is described by Mandelbrot’s equation (Feder, 1988), $P=\gamma {\left(\sqrt{A_S}\right)}^D$, where $D$ is the fractal dimension and $\gamma$ is a constant, both of which increase with landscape fragmentation (Blanco Jorge and Garcia, 1997). Intensive agriculture fields are often far larger in Orinoquia than elsewhere in Colombia. Furthermore, mechanization usually necessitates fields with basic shapes. Accordingly, and for simplicity, we assume that $D\approx 1$. Therefore, the total cropland-perimeter surrounding area can be calculated as

$A_{0,P}+A_{L,P}=wP\approx \gamma w\sqrt{A_S}(3)$

It seems reasonable to assume that the normalized burning rates of livestock and “natural” areas near the cropland perimeter are similar, i.e., $f_{0,P}\approx f_{L,P}=f_P$, where $f_P$ is the mean near-perimeter normalized burning rate. Hence, the burned area expected value is

$A_B\approx \left(f_{0,I}{A}_{0,I}+f_{L,I}{A}_{L,I}\right)+f_S{A}_S+\gamma w{f}_P\sqrt{A_S}(4)$

where $f_{0,I}$, $f_{L,I}$, and $f_S$ are the normalized burning rates of inland “natural”, inland livestock and cropland areas, respectively. As the total area, $A_{tot}$, is constant, the cropland area expands at the expense of livestock or “natural” areas. Therefore, the term in brackets is not constant under cropland expansion, but it likely diminishes over time provided the normalized burning rates remain constant. To further simplify the model, based on information for oil palm (Castiblanco, 2014; Furumo and Aide, 2017; Quezada et al., 2022), we assume that most of the new croplands in Orinoquia are subtracted from previously cattle ranching areas, i.e., $\Delta A_S=-∆A_L$. This is clearly not the case in Orinoquia’s bordering and transition areas with Amazonia, where forest fires are linked to illegal land grabbing (Armenteras et al., 2019). The resulting expression connects the inland livestock area with the cropland area,

$A_{L,I}\approx \left(A_{tot}-A_0\right)-A_S-w{\gamma }_L\sqrt{A_S}(5)$

where ${\gamma }_L$ is the livestock-cropland perimeter specific parameter. Finally, assuming that Orinoquia’s “natural” areas have not been significantly, or at least not rapidly, transformed into agriculture (crops or livestock) in recent times, we derive an equation for the biomass burning—cropland coupling coefficient, $K_S$,

$K_S=\partial A_B/\partial A_S\approx \left(f_S-f_{L,I}\right)+\frac{w}{2}\left(\gamma f_P-{\gamma }_L{f}_{L,I}\right)A_S^{-1/2}(6)$

As per its assumptions, this equation describes the effect on BA of replacing cattle grazing areas by croplands. The first term to the right describes the replacement direct effect owing to the change in normalized burning rate upon land use change. The rightmost term describes the perimeter fire suppression second order effect due to change in cropland perimeter.

We know that currently in the Colombian Orinoquia $A_0>A_L\gg A_S$, and that burning of permanent crops (e.g., plantations like oil palm and rubber tree) is actively inhibited and suppressed, thus very infrequent. On the other hand, information provided by commercial crop producers indicates that local communities rather frequently set harvest residues on fire, particularly those of corn and during the dry season. Thus $f_S$ is expectedly a small non-zero rate. The same holds true for $f_P$. On the other hand, $f_{L,I}$ is a much larger rate due to the degraded pasture renovation practice discussed above. Therefore, the direct and second order terms of the coupling coefficient, $K_S$, are expected to be negative. Romero-Ruiz et al. (2010) reported a satellite-derived annual fractional burned area for “ranches” of $f_L=0.19\pm 0.03$ a decade ago. If the cropland normalized burning rate is very small, i.e., $f_S\approx 0$, and the effect of its perimeter expansion is negligeable, i.e., the rightmost term approaches zero, then one would expect $K_S\approx -0.2$.

3 Region, datasets and methodology

3.1 The Colombian Orinoquia

The Orinoco River basin (ORIB) is the 8^th largest watershed in area and the 3rd highest in water discharge to the ocean in the world. Its total area of 981,446 km² is shared by Venezuela (65%) and Colombia (35%, i.e., 343,506 km²). This investigation is about the Colombian Orinoquia, or Colombian Llanos, i.e., the Colombian share of the mostly-savanna fraction of ORIB. There are several definitions of Orinoquia, and hence its boundaries and areas vary (Hernández Villamizar, 2019). Here we use the definition of Trujillo González et al. (2010), which equates the Colombian Orinoquia to the provinces of Arauca, Casanare, Meta, and Vichada. These 4 provinces add up to 254,271 km² (25.4 Mha) of mostly savannas (Supplementary Table S1). The remaining 26% of the Colombian ORIB is made up of transition zones with neighboring basins and regions (Andes, Amazonia, and Guayana Shield).

The Colombian Orinoquia is a complex set of 5 geomorphological landscapes (Piedmont, Alluvial Plains, Eolian Plains, High Plains, and Dissected High Plains), covering 22.3% of the Colombian continental territory, with altitudes between 80 and 500 m above the sea level, and ecological gradients between the Colombian Andes (Eastern Cordillera) and Amazonia. Their soils are mostly of low fertility, acidic, with low organic carbon content and toxic levels of iron and aluminum (Viloria de la Hoz, 2009). Among other factors, this is the result of seasonal exposure to floods, droughts, burning, and high winds. Because they are not subject to seasonal flooding, and due to their topographic characteristics, the flat High Plains (Altillanura in Spanish) are being rapidly transformed into mechanized agriculture. DNP (2014) estimated that 28 thousand km² of the 42 thousand km² of the High Plains had potential for crop and cattle production, and forestry.

The Colombian Orinoquia is distributed among the provinces (Departmentos) of Arauca (23,818 km²), Casanare (44,640 km²), Meta (85,635 km²), and Vichada (105,947 km²). Officially, it also includes smaller areas of the provinces of Guaviare, Guainía, Boyacá and Cundinamarca (Viloria de la Hoz, 2009). The four main provinces were home to 1.3 million inhabitants (Mhab) in 2018 (57% in Meta and 23% in Casanare). Although Orinoquia is one of the least populated regions of Colombia, its population grew 30% faster than Colombia’s during the 1985–2015 period, mainly due to internal migration. During the same period, the fraction of rural population in the Colombian Orinoquia decreased from 43% to 30% (https://microdatos.dane.gov.co/index.php/home).

Publicly available information on land use in the Colombian Orinoquia is currently deficient. There are accurate data on 1) indigenous peoples’ reservations (40,759 km²—Díaz Lemus et al., 2017), and 2) national parks and other public and private protected areas (24,897 km²—https://www.parquesnacionales.gov.co/portal/es/sirap-orinoquia/areas-protegidas-del-sirap-orinoquia/). Based on detailed permanent and SC crop sown area (see below), and assuming that SC crops are harvested 2 times a year and undergo fallow every 4 years, and that infrastructure (e.g., internal roads and storage facilities) occupies 10% of the cropland area, we reckoned the 3) cropland physical area at ∼8,192 km² in 2018, i.e., 95% of the reported total sown area. Based on cattle statistics (https://www.fedegan.org.co/estadisticas/inventario-ganadero) and stocking rate estimations per province (∼0.36 animal units—adults—per ha, weighted average), and information on calf share (41.4%), we gauge the 4) livestock area at ∼86,058 km² (see Supplementary Table S1 for details and references). The remainder ∼94,365 km² out of the 254,271 km² of the Colombian Orinoquia are likely non-official “natural” areas, as defined by the coupling model above. Thus, we estimate that “natural”, livestock, and cropland areas used 63%, 34% and 3.2%, respectively, of the Colombian Orinoquia in 2018, which is consistent with Etter et al. (2010) estimation of savanna clearing (37%) by 2020. The number of cattle in the Colombian Orinoquia increased almost linearly by ∼60% during the 2001–2021 period (https://www.fedegan.org.co/estadisticas/inventario-ganadero). If this trend continues and stocking rates are not improved, ∼1,886 km² yr⁻¹ will most certainly be incorporated for livestock in the years to come, likely by subtraction from “natural” areas.

3.2 Sown and burned areas

We retrieved the Colombian Orinoquia sown and harvested areas [ha yr⁻¹], and their production [ton yr⁻¹] and yields [ton ha⁻¹] from the Colombian Ministry of Agriculture (MinAgricultura) AGRONET database (https://www.agronet.gov.co). The available, most recent dataset covers the 2007–2018 period at annual resolution, with 69 crops disaggregated by 59 municipalities in the 4 provinces. We further reaggregated the crops into 11 groups and by province for simplicity. After considering the characteristics of the crops, the primary criterion for crop grouping was the similarity of their effective emission factors. Oil palm was not clustered with any other crop considering its relevance and emissions characteristics. The crop groups were the following (in descending order according to sown area in 2017): oil palm, rice (upland and paddy), corn, soybean, musaceae (banana plantain included), fruit trees (cocoa included), sugarcane, rubber tree, tubercles (and root crops, cassava included), others, and coffee. “Others” is a rather heterogeneous group, but with a very small aggregated sown area. Unfortunately, there is no information on the area sown with improved (“exotic”) pastures as these are not officially accounted as crops.

We analyzed the cropland expansion by crop and by province using the logistic and linear models above. Many permanent crops take years to reach maturity so they can be harvested. Due to the implicit value of these long payback periods, one would expect them to be more actively protected from fire than SC crops. Moreover, SC crops generate more PM10 emissions than permanent crops (Table 1). Therefore, we also analyzed the trends of the permanent (oil palm, musaceae, fruit trees, sugarcane, rubber tree, coffee) and short-cycle (rice, corn, soybean, tubercles) crop categories.

We obtained the BA time series for the polygon-defined areas of the 4 Colombian Orinoquia provinces for the longest, most recently available time period (1997–2016), from the Global Fire Emissions Database (GFED4, Giglio et al., 2013) site (http://globalfiredata.org/pages/data/). Although GFED data has been extensively validated, it is well known that GFED and other satellite products systematically miss small fires (Ramo et al., 2021). As Orinoquia savanna fires are predominantly small (Hernandez et al., 2019; Hernández Villamizar, 2019), GFED’s burned areas are likely underestimated yet suitable for the purposes of this paper. BA linear trends were estimated using Eq. 2.

3.3 Burned area dependance on precipitation

Vegetation fires require biomass to accumulate over time and then to dry enough to burn. The accumulation process is complex and non-linearly controlled mainly by water availability and solar radiation (Morales-Rincon et al., 2021). Also, nominal drought conditions must persist for several days for biomass to be set on fire. These two processes are thus mainly ruled by precipitation. Therefore, the BA is expected to be significantly, although probably also non-linearly, correlated with precipitation, particularly under severe drought conditions, e.g., due to El Niño (ENSO) events (drought in Orinoquia). Our research showed that the monthly and daily precipitation levels played a more important controlling role in the Orinoco savannas’ growth than the accumulated annual precipitation (Morales-Rincon et al., 2021).

For correlation purposes and to ensure full spatial coverage, we used reanalysis ERA5 (Hersbach et al., 2019) monthly precipitation data, spatially averaged for the 4 Colombian Orinoquia provinces. These data were retrieved from the Copernicus Climate Data Store (CDS) (https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels-monthly-means?tab=form). We validated the reanalysis data by comparison with IDEAM’s (Colombian Institute of Hydrology, Meteorology and Environmental Studies) raingauge network measurements (http://dhime.ideam.gov.co/), which are unevenly distributed and scant in a large portion of Orinoquia. See Supplementary Figure S2. We also analyzed the 40-year monthly precipitation reanalysis time series, particularly for trend detection using simple linear regression and a LOESS (locally estimated scatterplot smoothing)-based seasonal trend decomposition procedure (STL; Cleveland et al., 1990).

We correlated annual province burned area and its anomaly (burned area residual after linear trend removal) with statistical descriptors of province annual precipitation level (monthly mean, maximum precipitation, annual sum of monthly positive anomalies), annual variability (annual range, monthly standard deviation), and annual drought (monthly minimum precipitation, annual sum of monthly negative anomalies), and linear, exponential and logarithmic combinations of non-covariant descriptors. We selected a correlation function based on goodness of fitting, simplicity and plausibility.

3.4 Cropland-burned area coupling, emissions and scenarios

Satellite measurements show that BB is decreasing worldwide for reasons that are not fully understood (Knorr et al., 2016; Andela et al., 2017; Joshi and Sukumar, 2021). In many regions of the world, including Orinoquia, the declining trend is blurred by high time and spatially variability. BB in Orinoquia is a very complex phenomenon, mostly anthropogenic (Hernández Villamizar, 2019), but strongly modulated by meteorological conditions. Our fitting and scenario calculation model assumes that the BA can be reasonably explained by the linear combination of the following four terms: A) A constant term equal to the BA at zero relative time corrected for rainfall variability and cropland coupling effects; B) A rainfall-dependent term expected to explain most of the BA interannual variability as sown area changes slowly over time; C) A cropland area—BA coupling term; and D) A time-decreasing trend term intended to explain BA decline unrelated to cropland expansion, e.g., due to increasing landscape fragmentation and woody fraction (at least partially attributed to CO₂ fertilization—Rosan et al., 2019). Orinoquia’s sown area is currently much smaller than its BA. Thus, one would expect the term “C” to be comparatively smaller than the terms “B” and “D”. For this reason, we obtained sequential, least-squares approximations to the model and fitting coefficients per province and Orinoquia’s as follows: 1) BA linear trend estimates (terms “A” and “D”); 2) Rain descriptors and trends as explanatory variables (terms “A”, “B”, and “D”); 3) Full model.

The agricultural PM10 emissions per crop group and province were calculated using the sown area based effective emission factors presented in Table 1. These emissions were also aggregated by province and crop type (permanent or SC). BB emissions were estimated using the effective emission factor as discussed in Section 2.2.

PM10 emissions due to crop production in Orinoquia are much lower, less variable and have a much clearer trend than BB emissions. Thus, for building emission scenarios to year 2040, we disregarded the uncertainty in the sown area trends, and instead focused on the BA uncertainty. Various studies indicate that climate change could shift the current anthropogenic-driven fire regime to a temperature-driven regime from the mid-21st century on (Pechony and Shindell, 2010; Knorr et al., 2016). Incorporating climatic and meteorological trends or shifts into the Orinoquia BA scenarios is beyond the scope of this paper. Therefore, we based the province-level and Orinoquia-aggregated BA scenarios on the coupling coefficient uncertainty. Besides the mean scenario, i.e., projected province and Orinoquia-aggregated burned areas calculated with a least-squares coupling coefficient ($K_S=\widehat{K_S}$), we calculated “no-coupling” ($K_S=0$), and “highly-positive” ($K_S=\widehat{K_S}+3{{\sigma }_K}_S$) coupling scenarios, where ${\sigma }_{K_S}$ is the standard deviation of $\widehat{K_S}$. These scenarios were compared to the BA linear trend.

4 Results and discussion

4.1 Agriculture trends and spatial distribution

The time-dependent logistic model in Eq. 1 severely underestimated the cropland expansion in Orinoquia. Moreover, in the case of Meta, the logistic model was unreliable, as evidenced by the large discrepancy between forecasts that included or excluded the last data point (Supplementary Figure S3). These undesirable features were associated with the sown area high variability and reduced number of data points, and the fact that the cropland area is currently quite far from its “carrying capacity” level. Cropland currently occupies only ∼10 thousand km² out of the 254 thousand km² of the Colombian Orinoquia. Using improved pastures and integrated crop-cattle production systems, cattle raising could significantly decrease its footprint. A 30% increase in stocking rates could “free” ∼20 thousand km², which is twice as much as the “new” cropland area required by 2040. It must be mentioned that the resource-regulated logistic models applied to Orinoquia by other investigators (Etter et al., 2010; Castiblanco et al., 2013) apparently do not suffer from the problems reported here.

The Colombian Orinoquia sown area increased by 494 ± 47 km² (49 ± 5 kha; plus/minus is 1-sigma) on average every year during the 2007–2018 period. See Figure 2. This is ∼65% faster in relative terms than during the 1996–2007 period (Viloria de la Hoz, 2009). Almost all (86%) of this expansion took place in Meta (275 km² yr⁻²) and Casanare (148 km² yr⁻²), the provinces with the largest cropland areas (6,172 km² yr⁻¹ and 2,707 km² yr⁻¹, respectively, out of the 10,028 km² yr⁻¹ sown in 2017). Arauca and Vichada had much smaller cropland areas but comparable or much larger expansion rates in relative terms (10 ± 3% yr⁻¹ and 152 ± 43% yr⁻¹, respectively). Vichada, a 100 thousand km² province, i.e., slightly larger than Hungary, had less than 500 km² of crops during the 2010s, yet this was a ∼20-fold increase compared to the sown area in 2007 (22 km²). Agriculture has been rapidly but also irregularly expanding in Vichada, the province with the most fragile soils in Orinoquia. Also, it is likely that a significant fraction of the cropland expansion in Vichada has occurred at the agricultural frontier, i.e., involving clearing of native savannas.

FIGURE 2

FIGURE 2. Sown area in Colombia’s Orinoquia, broken down by province. The dashed lines are the extrapolated linear fits (Eq. 2) for the provinces and their aggregate for Orinoquia.

A large fraction of the sown area growth came from the sustained expansion of oil palm (207 km² yr⁻²) and corn (66 km² yr⁻²), and the much variable growth of rice (62 km² yr⁻²). See Figure 3. Moreover, crops that were minor in 2007, particularly sugarcane and rubber tree, grew at staggering rates of 46 ± 8% yr⁻¹ and 89 ± 3% yr⁻¹, respectively (Supplementary Figure S4). In relative terms, permanent crops have expanded more than twice faster and more steadily (18 ± 1% yr⁻¹) than short-cycle crops (7 ± 2% yr⁻¹). The steady growth of permanent crops is to some extent tied to their more stable nature and likely to their higher market stability and lower price volatility. On the opposite side, the impact of higher volatility is evident for rice, which still is the second most important crop in Orinoquia but with a comparatively lower and highly variable growth rate (4 ± 2% yr⁻¹).

FIGURE 3

FIGURE 3. Sown areas and their growth rates in Colombia’s Orinoquia, broken down by crop group.

4.2 Burned area and precipitation trends and correlation

GFED estimated that 16,806 ± 4,328 km² (1.68 ± 0.43 Mha) burned annually in Orinoquia during the 1997–2016 period. This means that 6.6% of Orinoquia burned annually, with large differences among provinces, from 3% in Casanare to 9% in Vichada. This is a large burned area, both in absolute and relative terms (a factor ∼2 larger than the cultivated area in Orinoquia and just a factor ∼2 smaller than the area burned annually by wildfires in the United States, a country that is ∼40 times larger than Orinoquia). See Supplementary Figure S5. Half of the burned area occurred in Vichada and another ∼30% in Meta. The provinces with the larger fractional burned areas (9% in Vichada and 6% in Arauca) were also the provinces with the least fractional cropland cover. Moreover, an analysis of the normalized BA revealed that despite the large absolute differences in BA among the provinces, the detrended normalized BA anomalies (BA minus trend) are significantly and positively correlated, particularly between Vichada and Casanare (r² = 0.62), and to a lesser extent, between Meta and Arauca (r² = 0.28). See Figure 4. This suggests that the burned area is strongly modulated by weather effects. For instance, the 2014–2016 regional BB peaking periods (positive anomalies) might be linked to a strong El Niño (ENSO) episode (Navarro-Monterroza et al., 2019). Also, significant BA decreasing trends were observed in Meta (−3.7 ± 1.1% yr⁻¹), Vichada (−1.3 ± 1.1% yr⁻¹), and Arauca (−1.0 ± 1.0% yr⁻¹). These are comparable to the South America trends reported by Andela et al. (2017). Casanare’s BA trend (−0.2 ± 0.8% yr⁻¹) was indistinguishable from its variability. Meta, the province with the largest permanent (3,451 km² yr⁻¹ in 2017) and SC (2,578 km² yr⁻¹ in 2017) crop areas, showed the largest decreasing trend. On the other hand, with less than ∼500 km² of crops, the decreasing trend in Vichada cannot be explained by cropland “inhibition” alone.

FIGURE 4

FIGURE 4. Mean-normalized burned areas in Colombia’s Orinoquia, disaggregated by province. Dashed lines are linear fits (Eq. 2).

The raingauge-validated, 40-year long ERA5 precipitation dataset (Supplementary Figure S2) showed that the annual rainfall in the High Plains provinces of Vichada (2760 ± 195 mm yr⁻¹) and Meta (2630 ± 225 mm yr⁻¹) was climatologically ∼10% higher than in the seasonally flooded provinces of Arauca (2557 ± 296 mm yr⁻¹) and Casanare (2357 ± 264 mm yr⁻¹). Although all of these represent comparatively high precipitation levels, it must be stressed that Orinoquia’s annual precipitation cycle encompasses a ∼3-month long dry season, which can be severe enough to provide conditions for a large number of fires. The dry season severity was higher in the seasonally flooded provinces.

Both linear regression and STL revealed quite concerning long-term precipitation decrease trends in the four provinces, particularly from the beginning of the 1990s on (Supplementary Figure S6). On average, the annual rainfall in Arauca, Casanare, Meta, and Vichada decreased by 17, 16, 10, and 7 mm per year, respectively, during the 1979–2019 period. The much larger decreasing trend in the seasonally flooded provinces has accumulated over time, so annual rainfall was 18% lower in Arauca and Casanare compared to the High Plains provinces during the 2016–2019 period. Although biomass accumulates less if rainfall diminishes, a standing reduction in precipitation likely implies that the fire frequency and possibly also the burned area per fire event will increase over time. Instead, the burned area has diminished in most of Orinoquia, which suggest that effects other than climatological have played a controlling role in BB (Figure 5). These might include landscape fragmentation, woody vegetation encroachment, and cropland expansion. Explicit modeling of the first two effects is beyond the scope of this paper, so we simply modeled them as a weather-uncorrelated time trend.

FIGURE 5

FIGURE 5. Normalized linear trends for rainfall, burned area, and burned area corrected for rainfall variability during the 1997–2016 period.

The function we found that described with the highest statistical significance the burned area-rain relation was the logarithm of the ratio of rainfall annual standard deviation to rainfall minimum (“rainfall function” hereafter). This was also a biogeochemically sound function as BA is expected to increase during severe dry seasons or drought periods and when rainfall becomes very irregular. The BA normalized trend and the rainfall function were simultaneously fitted for each province using the following model,

$A_B/\overline{A_B}=b_0+b_{1}t+\lambda \ln \left[SD\left(L\right)/\min \left(L\right)\right](7)$

where $\overline{A_B}$ is the time series mean annual burned area, $L$ is the monthly precipitation [mm month⁻¹] (or rainfall, as there is no snow and almost no hail in Orinoquia), SD is the annual standard deviation, min denotes annual minimum, $\lambda$ is an adimensional, expectedly positive rainfall function constant, $b_1$ is a weather-uncorrelated time trend constant, and $b_0$ an adimensional offset. The rainfall function reasonably reproduced the observed BA interannual variability in the provinces of Meta and Arauca (Supplementary Figure S7). The fitted function indicates that ∼40% of the BA in these two provinces was associated with drought and irregular rainfall conditions. On the other hand, the rainfall function was unable to reproduce the unusually positive 2004 and 2007 BA anomalies in Vichada and Casanare, which might imply that factors other than meteorological were at play during those 2 years. The Pearson coefficients ranged from r² = 0.56 for Meta, to r² = 0.35 for Arauca, down to r² = 0.21 for Vichada and r² = 0.06 for Casanare. The first two might be considered acceptable given the model simplicity. The rainfall function constant was higher, very similar and of low uncertainty for Meta and Arauca (λ = 0.2), compared to lower and highly uncertain values for Vichada (λ = 0.1) and Casanare (λ = 0.05). The removal of the rainfall-related interannual variability revealed that the BA has decreased faster than reckoned with a plain BA trend (Figure 5). It is also noticeable that the BA in Meta decreased ∼3 times faster (at −4.5% yr⁻¹) than elsewhere in Orinoquia, which could be associated to the large cropland expansion in Meta compared to the other three provinces.

4.3 Cropland—burned area coupling coefficient

As explained in Section 3.4, we adopted a sequential fitting approach (plain trends → rainfall effects → coupling coefficient) given the fact that the cropland area in Orinoquia is still relatively small. The province-level BA full model thus incorporates rainfall, cropland fire suppression and unidentified trend effects as follows,

$A_{B,j}=c_{0,j}+c_{1,j}t+{\mathrm{\Lambda }}_{\mathrm{j}}\ln \left[SD\left(L_j\right)/\min \left(L_j\right)\right]+K_S{A}_{S,j}(8)$

All parameters but the coupling coefficient, $K_S$, are province specific (subscript j). $A_{B,j}$ is the annual BA [km² yr⁻¹], $c_{0,j}$ is a multilinear intercept value [km² yr⁻¹], $c_{1,j}$ a weather and cropland uncorrelated time trend constant [km² yr⁻²], and ${\Lambda }_{\mathrm{j}}$ is the absolute constant [km² yr⁻¹] to the above-described rainfall function. Thus, the model that described BA in Orinoquia had 13 parameters (3 per province plus the coupling coefficient) to be estimated by least-squares applied to the 4-province dataset at once for statistical robustness. After testing several fitting variants, we chose to keep the $c_{0,j}=b_{0,j}\overline{A_{B,j}}$ and ${\Lambda }_j={\lambda }_j\overline{A_{B,j}}$ parameters from the individual province fitting (Eq. 7), so as to allow the uncorrelated and coupling coefficient terms to “compete” for the remaining variance. We also tested multilinear correlations with permanent and SC crop areas as independent variables per province and aggregated. We found that permanent and SC crop areas are colinear enough to degenerate the multilinear correlations, so we kept the sown areas per province as fitting variables. The simultaneous fitting of Eq. 8 to the 4-province dataset explained a large fraction of the absolute BA variability and decline (r² = 0.93), much more than the normalized rainfall model. See Figure 6 and Supplementary Figure S8. The least-squares estimated coupling coefficient was $K_S=-0.17\pm 0.09$ (1-sigma), which implies that each sown hectare reduced the burned area by 0.17 ha (±55%). This value is very close to what was theoretically anticipated (Section 2.3). The least squares brought Arauca’s and Casanare’s $c_1$ to zero and significantly lowered Meta’s, which indicates that all the BA decline in the first two provinces and a significant fraction of Meta’s could be attributed to cropland expansion. On the contrary, almost all the BA reduction in Vichada was uncorrelated to cropland expansion (${c_1\approx b}_1\overline{A_B}$).

FIGURE 6

FIGURE 6. Observed and projected aggregated sown and burned areas in Orinoquia from 1997 to 2040. The projection included the following three scenarios: coupling as least-squares fitted (K_S = −0.17), no coupling (K_S = 0), and highly-positive coupling (K_S = 0.11).

4.4 Cropland expansion, and burned areas and emissions scenarios

The cropland area was projected by aggregating the estimations based on both the province and the crop-group trends (Eq. 2), reaching perfect agreement. Both estimate that the sown area will double by 2040, reaching 20,425 km² yr⁻¹ (2.04 Mha yr⁻¹). See Figures 2, 6. The linear trends predict that oil palm will become the dominant crop in Orinoquia with 39% of the sown area by 2040, compared to 32% in 2017. Most of oil palm’s gained share is estimated to come from rice, the participation of which is estimated to drop from 28% in 2017 to 18% by 2040. Estimated changes in participation for other crops are minor. Also, Casanare’s and Vichada’s participation will increase, at the expense of Meta’s, from 27% to 2% in 2017 to 29% and 5% in 2040, respectively. It is worth emphasizing that these simple linear models disregard crop prices, market trends, meteorological factors, and public policies, which likely account for most of the sown area change and variability.

The full model (Eq. 8) with the least-squares fitted coupling coefficient ($K_S=-0.17$) predicts that, on average, the BA will keep on declining in the four provinces (Figure 6 and Supplementary Figure S8). The model results indicate that a large fraction of the excursions around these mean trends can be attributed to weather interannual variability, which cannot be predicted. Thus, the BA scenarios were all calculated with the province-specific climatological-mean (1979–2019) rainfall functions, $\overline{\ln \left[SD\left(L\right)/\min \left(L\right)\right]}$, which ranged from 1.0 for Meta to 2.3 for Casanare. According to the fitted model ($K_S=-0.17$), the mean aggregated BA in Orinoquia will drop ∼4 fold by 2040, from ∼20 thousand km² (∼2 Mha) in 1997. The no-coupling scenario ($K_S=0$), which is indistinguishable from the plain linear trend, predicts a less rapid BA decline down to ∼6 thousand km² by 2040. Even the positive coupling scenario ($K_S=0.11$) predicts a largely reduced BA by 2040 (∼8 thousand km²). Overall, this indicates that cropland expansion accelerates BA reduction but also that a large fraction of the BA drop comes from Vichada’s and Meta’s unexplained “residual” decline trends, i.e., uncorrelated to sown area or weather. These are possibly related to unaccounted drivers such as landscape fragmentation, woody vegetation encroachment, cultural practices changes, and the replacement of native by exotic pastures. The uncorrelated trends, along with climate change effects, are by far the largest sources of uncertainty in the model.

Our emission factor (EF) calculations show that SC crops emitted a factor ∼3 more PM10 per sown hectare than permanent crops (Table 1). Short-cycle crops require more land preparation (∼7 times higher EFs), for which they use more diesel machinery (∼3 times higher EFs) than permanent crops. This is further amplified by the fact that SC crops may be harvested up to three times per year (accounted for as sown area). In turn, the sown area weighted crop production EF, 3.3 kg PM10 (sown ha)⁻¹ in 2016 (3 kg PM10 (sown ha)⁻¹ in 2018), was an order of magnitude smaller than the biomass burning EF.

We estimate that crop production in Orinoquia generated 3 Gg (kton) of PM10 in 2016, 71% of which came from SC crops (Figure 7). In turn, biomass burning generated 57.5 Gg of PM10 that same year, i.e., 95% of the combined emissions. The BB emissions were not just ∼20 times larger but also significantly more toxic. We predict that the sown area will increase by a factor of 2 by 2040 referred to 2017. In turn, crop emissions will increase only by 83%, due to the faster expansion of oil palm compared to SC crops. In contrast, the BB emissions (with $K_S=-0.17$) will drop by a factor of 2.5 and the combined emissions by a factor of 2 by 2040 referred to 2017 (down to 17.1 and 23.3 Gg yr⁻¹, respectively), despite the increase in crop emissions. With a coupling coefficient of $K_S=-0.17$ burned ha (sown ha)⁻¹, the cropland expansion fire suppression effects (−0.17 burned ha (sown ha)⁻¹ ⋅ 36 kg PM10 (burned ha)⁻¹) would outweigh the crop emissions increase (3 kg PM10 (sown ha)⁻¹). It must be emphasized that the coupling coefficient is still a very uncertain value. Despite this, the other two scenarios also predict substantial PM10 emission reductions. The no-coupling scenario ($K_S=0$), which is basically identical to the plain trend, and the highly positive coupling scenario (with $K_S=0.11$) predict 23% and 53% higher combined emissions by 2040, respectively, than the fitted $K_S$ scenario. This persistent reduction is to a large extent explained by the uncorrelated BA decline trends of Vichada and Meta, the drivers of which are not understood.

FIGURE 7

FIGURE 7. Observed and projected crop production, biomass burning, and combined PM10 emissions until 2040. The projection included the same three coupling scenarios as in the previous figure.

The combined BB and crop production annual PM10 emissions in the Colombian Orinoquia in 2016 (60.5 Gg yr⁻¹) were very similar to the total PM10 emissions in Bogota (68.2 Gg yr⁻¹—Pachón et al., 2018), a 9.3 Mhab metropolitan area. The main difference between two is their timing: 95% of Orinoquia’s emissions are due to BB, most of which occurs during the 3-month dry season (Hernandez et al., 2019). Thus, PM10 fluxes from the Colombian Orinoquia are seasonally ∼4 times larger than Bogota’s urban fluxes. These emissions are transported to Bogota and other Andean cities causing major air pollution episodes (Mendez-Espinosa et al., 2019; Ballesteros-González et al., 2020; Rincón-Riveros et al., 2020). If the uncorrelated component of Orinoquia’s BA decline trend holds, a sustained BA decrease in Orinoquia will have major positive impacts on air quality and public health, both in Orinoquia (1.3 Mhab) and the Colombian Andes (>10 Mhab).

5 Conclusion

We estimate that 34% of Orinoquia’s 254 thousand km² (25.4 Mha) was used for cattle ranching (mainly extensive) in 2018, which, based on cattle statistics and stocking rate data, expands at 2% yr⁻¹ (1,886 km² yr⁻¹). Indigenous reservations, national parks (and other protected areas), and crops occupied 16%, 9.8%, and 3.2% of Orinoquia’s, respectively. The remaining 37% of the region is likely still in a semi-natural state, including fire-susceptible native savanna at risk of incorporation into agriculture, particularly at the agricultural frontier.

Colombia’s Orinoquia cropland area (10 thousand km² in 2017) has been rapidly and steadily expanding since the 1990s. The current (2007–2018) area-averaged linear expansion rate (12% yr⁻¹) is 65% larger than that observed during the 1996–2007 period, which indicates that cropland expansion in Orinoquia has been long-term accelerating. Most (86%) of the cropland expansion in the 2007–2018 period took place in the Meta and Casanare provinces, which share 89% of sown area. Although oil palm, corn, and rice accounted for 68% of the sown area growth rate (494 km² y⁻²), permanent crops expanded more than twice as fast (18% yr⁻¹) and more steadily than short-cycle (SC) crops.

According to GFED, on average 7% of Orinoquia (17 thousand km²; 1.7 Mha) burned every year during the 1997–2016 period, although with large differences among provinces and significant interannual variability (±26%). This is a large burned area (BA), both in absolute and relative terms. Up to 40% of the BA interannual variability was linked to irregular rainfall conditions and seasonal drought, as modeled by the “rainfall function”, $\lambda \ln \left(SD\left(L\right)/\min (L\right))$ (L is monthly rainfall and SD is standard deviation). The provinces with the larger fractional BA were also those with the least fractional cropland cover. Also, significant BA decrease trends were observed in Meta (−4% yr⁻¹), Vichada, and Arauca.

A BA model incorporating rainfall effects and uncorrelated decline, along with a linear BA ($A_B$)—sown area ($A_S$) coupling term, $K_S{A}_S$, successfully explained most of Orinoquia’s BA variability and decline (r² = 0.93). The fitted coupling coefficient, $K_S=\partial A_B/\partial A_S=-0.17$, which means that each sown hectare reduced the BA by 0.17 ha, was very close to the theoretically anticipated value. It also confirms and quantifies Etter et al. (2010) and Andela et al. (2017) cropland-biomass burning coupling hypothesis.

We estimate that the sown area in Orinoquia will double during the 2017–2040 period, with a much more rapid expansion of oil palm compared to rice. The BA full model predicts that this expansion coupling effect, combined with the uncorrelated BA decline, will cause a factor 2.5 reduction in Orinoquia’s mean burned area by 2040 compared to 2017. The largest sources of uncertainty in this prediction are the uncorrelated BA decline trends, i.e., those not explained by weather or cropland expansion, and expected climate change effects.

We calculate that crop production in Orinoquia generated 3 Gg of PM10 in 2016 (71% from short-cycle crops), while biomass burning generated 57 Gg, i.e., 95% of the combined emissions that year. Due to the increasing share and lower emissions of oil palm compared to short-cycle crops, the aggregated crop production emissions are expected to increase by a factor of 1.8 during the 2017–2040 period. This increase is small considering the factor 2.5 decrease in biomass burning emissions during the same period, resulting in a factor two drop in the combined emissions (23 Gg of PM10 by 2040). A pollution burden reduction of this magnitude would have tremendous positive effects on public health, not just in Orinoquia (1.3 Mhab) but also in its direct outflow region, the Colombian Andes (>10 Mhab).

Data availability statement

The data used in this publication will be made available on request.

Author contributions

RJ and EL-V conceived and managed the project. AH, AA, AV-B, and RJ reviewed the literature, obtained and analyzed the data, and implemented the emission models. RJ developed the rainfall and coupling models and 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 Orinoquia,” and by Universidad EAN.

Acknowledgments

We gratefully acknowledge the information on agricultural and fire suppression practices provided by colleagues and crop production professionals and practitioners. We are also thankful to Glenn Hyman and two reviewers for their detailed and insightful revision of the manuscript.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fenvs.2022.689844/full#supplementary-material

References

Akagi, S. K., Yokelson, R. J., Wiedinmyer, C., Alvarado, M. J., Reid, J. S., Karl, T., et al. (2011). Emission factors for open and domestic biomass burning for use in atmospheric models. Atmos. Chem. Phys. 11, 4039–4072. doi:10.5194/acp-11-4039-2011

CrossRef Full Text | Google Scholar

Andela, N., Morton, D. C., Giglio, L., Chen, Y., Van Der Werf, G. R., Kasibhatla, P. S., et al. (2017). A human-driven decline in global burned area. Science 356, 1356–1362. doi:10.1126/science.aal4108Available at: http://science.sciencemag.org/content/356/6345/1356.abstract.

PubMed Abstract | CrossRef Full Text | Google Scholar

Anderson, R., Bayer, P. E., and Edwards, D. (2020). Climate change and the need for agricultural adaptation. Curr. Opin. Plant Biol. 56, 197–202. doi:10.1016/j.pbi.2019.12.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Aneja, V. P., Schlesinger, W. H., and Erisman, J. W. (2009). Effects of agriculture upon the air quality and climate: Research, policy, and regulations. Environ. Sci. Technol. 43, 4234–4240. doi:10.1021/es8024403

PubMed Abstract | CrossRef Full Text | Google Scholar

Aneja, V. P., Schlesinger, W. H., and Erisman, J. W. (2008). Farming pollution. Nat. Geosci. 1, 409–411. doi:10.1038/ngeo236

CrossRef Full Text | Google Scholar

Armenteras, D., Schneider, L., and Dávalos, L. M. (2019). Fires in protected areas reveal unforeseen costs of Colombian peace. Nat. Ecol. Evol. 3, 20–23. doi:10.1038/s41559-018-0727-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Arora, V. K., and Melton, J. R. (2018). Reduction in global area burned and wildfire emissions since 1930s enhances carbon uptake by land. Nat. Commun. 9, 1326. doi:10.1038/s41467-018-03838-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Ballesteros-González, K., Sullivan, A. P., and Morales-Betancourt, R. (2020). Estimating the air quality and health impacts of biomass burning in northern South America using a chemical transport model. Sci. Total Environ. 739, 139755. doi:10.1016/j.scitotenv.2020.139755

PubMed Abstract | CrossRef Full Text | Google Scholar

Blanco Jorge, L. A., and Garcia, G. J. (1997). A study of habitat fragmentation in Southeastern Brazil using remote sensing and geographic information systems (GIS). For. Ecol. Manage. 98, 35–47. doi:10.1016/S0378-1127(97)00072-8

CrossRef Full Text | Google Scholar

Carrillo, M., and González, J. M. (2002). A new approach to modelling sigmoidal curves. Technol. Forecast. Soc. Change 69, 233–241. doi:10.1016/S0040-1625(01)00150-0

CrossRef Full Text | Google Scholar

Castiblanco, C., Etter, A., and Aide, T. M. (2013). Oil palm plantations in Colombia: A model of future expansion. Environ. Sci. Policy 27, 172–183. doi:10.1016/j.envsci.2013.01.003

CrossRef Full Text | Google Scholar

Castiblanco, C. (2014). Scenarios of the future expansion of oil palm in Colombia: Impacts generated by the biofuels sector. Bogota (Colombia): Pontificia Universidad Javeriana. Doctoral dissertation. Available at:https://repository.javeriana.edu.co/handle/10554/15737 (Accessed November 19, 2022).

Google Scholar

Cleveland, R. B., Cleveland, W. S., McRae, J. E., and Terpenning, I. (1990). Stl: A seasonal-trend decomposition procedure based on loess. J. Off. Stat. 6, 3–73.

Google Scholar

Díaz Lemus, E., Moreno Gaviria, O., Urango Petro, Á., and Castro Sabbagh, C. I. (2017). Los recursos de la Asignación Especial del SGP para Resguardos y el reto para los pueblos indígenas en el ejercicio de la función pública. Bogotá: Contraloría General de la República. Available at: https://www.contraloria.gov.co/.

Google Scholar

Dinero (2011). Los nuevos llaneros. Rev. Dinero (February 21) (367). Available at: https://www.semana.com/caratula/edicion-impresa/articulo/los-nuevos-llaneros/113492/.

Google Scholar

DNP (2008). Lineamientos de política para promover la producción sostenible de biocombustibles en Colombia. Bogotá: Contraloría General de la República. Available at: https://colaboracion.dnp.gov.co/CDT/Conpes/Económicos/3510.pdf.

Google Scholar

DNP (2014). Política para el desarrollo integral de la Orinoquia: Altillanura – fase I. Bogotá: Contraloría General de la República. Available at: https://colaboracion.dnp.gov.co/CDT/Conpes/Económicos/3797.pdf.

Google Scholar

García Durán, O. (2009). Carimagua: La investigación y el desarrollo en ecosistemas de baja fertilidad. Rev. Colomb. Ciencias Pecu. 22, 74–78.

Google Scholar

Etter, A., Sarmiento, A., and Romero, M. H. (2010). “Land use changes (1970–2020) and carbon emissions in the Colombian Llanos,” in Ecosystem function in savannas. Editors M. J. Hill, and N. P. Hanan (Boca Raton: CRC Press). doi:10.1201/b10275

CrossRef Full Text | Google Scholar

FAO (2018). The future of food and agriculture – alternative pathways to 2050. Available at: http://www.fao.org/3/I8429EN/i8429en.pdf (Accessed November 19, 2022).

Google Scholar

Feder, J. (1988). The perimeter-area relation. Fractals, 200–211. doi:10.1007/978-1-4899-2124-6_12

CrossRef Full Text | Google Scholar

Foley, J. A., Ramankutty, N., Brauman, K. A., Cassidy, E. S., Gerber, J. S., Johnston, M., et al. (2011). Solutions for a cultivated planet. Nature 478, 337–342. doi:10.1038/nature10452

PubMed Abstract | CrossRef Full Text | Google Scholar

Furumo, P. R., and Aide, T. M. (2017). Characterizing commercial oil palm expansion in Latin America: Land use change and trade. Environ. Res. Lett. 12, 024008. doi:10.1088/1748-9326/aa5892

CrossRef Full Text | Google Scholar

Giglio, L., Randerson, J. T., and van der Werf, G. R. (2013). Analysis of daily, monthly, and annual burned area using the fourth-generation global fire emissions database (GFED4). J. Geophys. Res. Biogeosci. 118, 317–328. doi:10.1002/jgrg.20042

CrossRef Full Text | Google Scholar

Graesser, J., Ramankutty, N., and Coomes, O. T. (2018). Increasing expansion of large-scale crop production onto deforested land in sub-Andean South America. Environ. Res. Lett. 13, 084021. doi:10.1088/1748-9326/aad5bf

CrossRef Full Text | Google Scholar

Hernandez, A. J., Morales-Rincon, L. A., Wu, D., Mallia, D., Lin, J. C., and Jimenez, R. (2019). Transboundary transport of biomass burning aerosols and photochemical pollution in the Orinoco River Basin. Atmos. Environ. X. 205, 1–8. doi:10.1016/j.atmosenv.2019.01.051

CrossRef Full Text | Google Scholar

Hernández Villamizar, A. J. (2019). Assessment of the impact of biomass burning on air quality in the Colombian Orinoco River Basin. Bogota (Colombia): Universidad Nacional de Colombia. Doctoral dissertation. Available at: https://repositorio.unal.edu.co/handle/unal/76979.

Google Scholar

Hersbach, H., Bell, B., Berrisford, P., Biavati, G., Horányi, A., Muñoz Sabater, J., et al. (2019). ERA5 monthly averaged data on single levels from 1979 to present. Copernic. Clim. Chang. Serv. Clim. Data Store. Data repository. doi:10.24381/cds.f17050d7

CrossRef Full Text | Google Scholar

Joshi, J., and Sukumar, R. (2021). Improving prediction and assessment of global fires using multilayer neural networks. Sci. Rep. 11, 3295–3314. doi:10.1038/s41598-021-81233-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Klink, C. A., and Moreira, A. G. (2002). “Past and current human occupation, and land use,” in The Cerrados of Brazil ecology and natural history of a neotropical savanna. Editors P. S. Oliveira, and R. J. Marquis (New York, United States: Columbia University Press), 69–88. doi:10.7312/oliv12042.8

CrossRef Full Text | Google Scholar

Knorr, W., Arneth, A., and Jiang, L. (2016). Demographic controls of future global fire risk. Nat. Clim. Chang. 6, 781–785. doi:10.1038/nclimate2999

CrossRef Full Text | Google Scholar

Lark, T. J., Meghan Salmon, J., and Gibbs, H. K. (2015). Cropland expansion outpaces agricultural and biofuel policies in the United States. Environ. Res. Lett. 10, 044003. doi:10.1088/1748-9326/10/4/044003

CrossRef Full Text | Google Scholar

Lasslop, G., and Kloster, S. (2017). Human impact on wildfires varies between regions and with vegetation productivity. Environ. Res. Lett. 12, 115011. doi:10.1088/1748-9326/aa8c82

CrossRef Full Text | Google Scholar

Lasso, C. A., Usma, J. S., Trujillo, F., and Rial, B., A. (2010). Biodiversidad de la cuenca del Orinoco: Bases científicas para la identificación de áreas prioritarias para la conservación y uso sostenible de la biodiversidad. Bogotá (Colombia):Instituto de Investigación de Recursos Biológicos Alexander von Humboldt, WWF Colombia, Fundación Omacha, Fundación La Salle e Instituto de Estudios de la Orinoquia Universidad Nacional de Colombia. Available at http://repository.humboldt.org.co/handle/20.500.11761/34982 (Accessed November 19, 2022).

Google Scholar

López-Hernández, D., Hernández-Hernández, R. M., and Brossard, M. (2005). Historia del uso reciente de tierras de las sabanas de américa del sur. Estudios de casos en sabanas del orinoco. Interciencia 30, 623–630. Available at: http://ve.scielo.org/scielo.php?script=sci_arttext&pid=S0378-18442005001000008.

Google Scholar

Maffia, J., Dinuccio, E., Amon, B., and Balsari, P. (2020). PM emissions from open field crop management: Emission factors, assessment methods and mitigation measures – a review. Atmos. Environ. X. 226, 117381. doi:10.1016/j.atmosenv.2020.117381

CrossRef Full Text | Google Scholar

Mateus-Fontecha, L., Vargas-Burbano, A., Jimenez, R., Rojas, N. Y., Rueda-Saa, G., et al. (2022). Understanding aerosol composition in a tropical inter-Andean valley impacted by agro-industrial and urban emissions. Atmos. Chem. Phys. 22, 8473–8495. doi:10.5194/acp-22-8473-2022

CrossRef Full Text | Google Scholar

Mendez-Espinosa, J. F., Belalcazar, L. C., and Morales Betancourt, R. (2019). Regional air quality impact of northern South America biomass burning emissions. Atmos. Environ. X. 203, 131–140. doi:10.1016/j.atmosenv.2019.01.042

CrossRef Full Text | Google Scholar

Morales-Rincon, L. A., Hernandez, A. J., Rodriguez-Hernandez, N. S., and Jimenez, R. (2021). Carbon exchange and accumulation in an orinoco high plains native savanna ecosystem as measured by eddy covariance. Front. Environ. Sci. 9, 1–14. doi:10.3389/fenvs.2021.673932

CrossRef Full Text | Google Scholar

Navarro-Monterroza, E., Arias, P. A., and Vieira, S. C. (2019). El Niño-Oscilación del Sur, fase Modoki, y sus efectos en la variabilidad espacio-temporal de la precipitación en Colombia. Rev. Acad. Colomb. Cienc. Exactas Fis. Nat. 43, 120. doi:10.18257/raccefyn.704

CrossRef Full Text | Google Scholar

Pachón, J., Galvis, B., Lombana, O., Carmona, L., Fajardo, S., Rincón, A., et al. (2018). Development and evaluation of a comprehensive atmospheric emission inventory for air quality modeling in the megacity of bogotá. Atmos. (Basel). 9, 49. doi:10.3390/atmos9020049

CrossRef Full Text | Google Scholar

Pechony, O., and Shindell, D. T. (2010). Driving forces of global wildfires over the past millennium and the forthcoming century. Proc. Natl. Acad. Sci. U. S. A. 107, 19167–19170. doi:10.1073/pnas.1003669107

PubMed Abstract | CrossRef Full Text | Google Scholar

Pechony, O., and Shindell, D. T. (2009). Fire parameterization on a global scale. J. Geophys. Res. 114, D16115. doi:10.1029/2009JD011927

CrossRef Full Text | Google Scholar

Quezada, J. C., Guillaume, T., Poeplau, C., Ghazoul, J., and Buttler, A. (2022). Deforestation-free land-use change and organic matter-centered management improve the C footprint of oil palm expansion. Glob. Chang. Biol. 28, 2476–2490. doi:10.1111/gcb.16069

PubMed Abstract | CrossRef Full Text | Google Scholar

Rabin, S. S., Melton, J. R., Lasslop, G., Bachelet, D., Forrest, M., Hantson, S., et al. (2017). The Fire Modeling Intercomparison Project (FireMIP), phase 1: Experimental and analytical protocols with detailed model descriptions. Geosci. Model Dev. 10, 1175–1197. doi:10.5194/gmd-10-1175-2017

CrossRef Full Text | Google Scholar

Ramo, R., Roteta, E., Bistinas, I., van Wees, D., Bastarrika, A., Chuvieco, E., et al. (2021). African burned area and fire carbon emissions are strongly impacted by small fires undetected by coarse resolution satellite data. Proc. Natl. Acad. Sci. U. S. A. 118 (9), 1–7 (e2011160118). doi:10.1073/pnas.2011160118

PubMed Abstract | CrossRef Full Text | Google Scholar

Rincón-Riveros, J. M., Rincón-Caro, M. A., Sullivan, A. P., Mendez-Espinosa, J. F., Belalcazar, L. C., Quirama Aguilar, M., et al. (2020). Long-term Brown carbon and smoke tracer observations in bogotá, Colombia: Association with medium-range transport of biomass burning plumes. Atmos. Chem. Phys. 20, 7459–7472. doi:10.5194/acp-20-7459-2020

CrossRef Full Text | Google Scholar

Rodríguez-Gómez, C., Echeverry, G., Jaramillo, A., and Ladino, L. A. (2021). The negative impact of biomass burning and the orinoco low-level jet on the air quality of the Orinoco River basin. Atmósfera 35 (3), 497–520. doi:10.20937/atm.52979

CrossRef Full Text | Google Scholar

Romero-Ruiz, M., Etter, A., Sarmiento, A., and Tansey, K. (2010). Spatial and temporal variability of fires in relation to ecosystems, land tenure and rainfall in savannas of northern South America. Glob. Chang. Biol. 16, 2013–2023. doi:10.1111/j.1365-2486.2009.02081.x

CrossRef Full Text | Google Scholar

Rosan, T. M., Aragão, L. E. O. C., Oliveras, I., Phillips, O. L., Malhi, Y., Gloor, E., et al. (2019). Extensive 21st-century woody encroachment in south America’s savanna. Geophys. Res. Lett. 46, 6594–6603. doi:10.1029/2019GL082327

CrossRef Full Text | Google Scholar

Schnepf, R. D., Dohlman, E., and Bolling, C. (2001). Agriculture in Brazil and Argentina: Developments and prospects for major field crops. Washington, DC: U.S. Department of Agriculture, Market and Trade Economics Division, Economic Research Service, Agriculture and Trade Report WRS-01-3. Available at: https://www.ers.usda.gov/webdocs/outlooks/40339/15081_wrs013_1_.pdf?v=7292 (Accessed November 19, 2022)

Google Scholar

Sharratt, B., and Auvermann, B. (2014). “Dust pollution from agriculture,” in Encyclopedia of agriculture and food systems (Amsterdam, Netherlands: Elsevier), 487–504. doi:10.1016/B978-0-444-52512-3.00089-9

CrossRef Full Text | Google Scholar

Sharratt, B. S., and Van Pelt, R. S. (2020). “Erosion by wind: Source, measurement, prediction, and control” in Managing Soils and Terrestrial Systems. Editor B. D. Fath, and S. E. Jorgensen (Boca Raton (FL): CRC Press). Chapter 37, 321–341. doi:10.1201/9780429346255-41

CrossRef Full Text | Google Scholar

IPCC (2019). Climate change and land: An IPCC special report on climate change, desertification, land degradation, sustainable land management, food security, and greenhouse gas fluxes in terrestrial ecosystems. Editors J. M. P. R. Shukla, J. Skea, E. Calvo Buendia, V. Masson-Delmotte, H.-O. Pörtner, D. C. Robertset al. (Geneva, Switzerland: Intergovernmental Panel on Climate Change). Available at: https://www.ipcc.ch/srccl/.

Google Scholar

Silva, J. M. N., Carreiras, J. M. B., Rosa, I., and Pereira, J. M. C. (2011). Greenhouse gas emissions from shifting cultivation in the tropics, including uncertainty and sensitivity analysis. J. Geophys. Res. 116, D20304–D20321. doi:10.1029/2011JD016056

CrossRef Full Text | Google Scholar

Tilman, D., Fargione, J., Wolff, B., D'Antonio, C., Dobson, A., Howarth, R., et al. (2001). Forecasting agriculturally driven global environmental change. Sci. (80-.) 292, 281–284. doi:10.1126/science.1057544

PubMed Abstract | CrossRef Full Text | Google Scholar

Trujillo González, F., Beltrán Gutiérrez, M., Diaz-Pulido, A., Ferrer Pérez, A., and Payan Garrido, E. (2010). “Mamíferos,” in Biodiversidad de la cuenca del Orinoco: Bases científicas para la identificación de áreas prioritarias para la conservación y uso sostenible de la biodiversidad. Editors C. A. Lasso, J. S. Usma, F. Trujillo, and A. Rial (Bogota (Colombia): Instituto de Investigación de Recursos Biológicos Alexander von Humboldt, WWF Colombia, Fundación Omacha, Fundación La Salle, Instituto de Estudios de la Orinoquia Universidad Nacional de Colombia), 311–336. Available at: http://repository.humboldt.org.co/handle/20.500.11761/34982 (Accessed November 19, 2022).

Google Scholar

Venevsky, S., Le Page, Y., Pereira, J. M. C., and Wu, C. (2019). Analysis fire patterns and drivers with a global SEVER-FIRE v1.0 model incorporated into dynamic global vegetation model and satellite and on-ground observations. Geosci. Model Dev. 12, 89–110. doi:10.5194/gmd-12-89-2019

CrossRef Full Text | Google Scholar

Viloria de la Hoz, J. (2009). Geografía económica de la Orinoquia. Colombia: Cartagena.

Google Scholar

Wang, J., Miller, D. R., Sammis, T. W., Hiscox, A. L., Yang, W., and Holmén, B. A. (2010). Local dust emission factors for agricultural tilling operations. Soil Sci. 175 (4), 194–200. doi:10.1097/SS.0b013e3181dae283

CrossRef Full Text | Google Scholar

Zalles, V., Hansen, M. C., Potapov, P. V., Parker, D., Stehman, S. V., Pickens, A. H., et al. (2021). Rapid expansion of human impact on natural land in South America since 1985. Sci. Adv. 7, eabg1620–12. doi:10.1126/sciadv.abg1620

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Wu, S. Y., Krishnan, S., Wang, K., Queen, A., Aneja, V. P., et al. (2008). Modeling agricultural air quality: Current status, major challenges, and outlook. Atmos. Environ. X. 42, 3218–3237. doi:10.1016/j.atmosenv.2007.01.063

CrossRef Full Text | Google Scholar