Introducción
En los ecosistemas terrestres la captura y almacenamiento del carbono atmosférico es realizado principalmente por las plantas, organismos que mediante la fotosíntesis transforman el dióxido de carbono (CO2) en compuestos orgánicos (Farquhar et al. 2001). El almacenamiento de carbono, también denominado producción primaria o secuestro de carbono (Farquhar et al. 2001, Brown 2002, Phillips et al. 2004) es un proceso relevante para sustentar el crecimiento y la formación de tejidos vegetales aéreos (tallos, ramas, gajos, hojas, flores y frutos) y subterráneos (raíces gruesas y finas, rizomas y tallos subterráneos) (Brown 2002), así como para disminuir la concentración atmosférica de Gases de Efecto Invernadero (GEI) (Farquhar et al. 2001, Fonseca-González 2017, IPCC 2019b) y mantener y regular las funciones ambientales de los ecosistemas (Millennium Assessment 2005).
Si bien la captura y almacenamiento de carbono en los ecosistemas terrestres es un proceso constante en el tiempo y espacio (Phillips et al. 2004), anualmente a nivel mundial se estima que 12 (± 2.9) Gt CO2e (valor estimado de Dióxido de Carbono equivalente entre el periodo 2007 a 2016) son liberados a la atmósfera (gradual o abruptamente) como consecuencia de la deforestación, degradación de la vegetación y el posterior uso de la tierra (IPCC 2019b); lo cual significa la pérdida de reservorios y sumideros de carbono (Farquhar et al. 2001, Brown 2002, Saatchi et al. 2011, Baccini et al. 2012, Baccini & Asner 2013), así como de otras funciones ambientales (Millennium Assessment 2005).
De acuerdo con la Global Forest Watch (2021), durante la gestión 2019 el 66.2% de la pérdida de la cobertura vegetal a nivel mundial se concentró en 10 países, entre los cuales se destaca Bolivia, ya que la disminución de la cobertura alcanzó aproximadamente a 861 mil ha de superficie, lo cual es equivalente a la pérdida de más de 101 millones de toneladas de biomasa vegetal aérea.
Tradicionalmente, en los ecosistemas terrestres el inventario y/o cuantificación de la biomasa y carbono almacenado en la vegetación es realizado mediante la implementación de parcelas de muestreo, el levantamiento de datos dasométrios y la posterior aplicación de modelos alométricos (Brown 2002, Phillips et al. 2004, Chave et al. 2005, Mostacedo et al. 2010, Feldpausch et al. 2012, Chave et al. 2014, Phillips et al. 2016). Este método, previamente aplicado en algunas regiones de Bolivia (Araujo-Murakami et al. 2006, 2014, 2016, Mostacedo et al. 2010, Morandi et al. 2018), tiene la ventaja de generar estimaciones locales más precisas (Brown 2002, Chave et al. 2005, 2014, Galindo et al. 2011). Sin embargo, cuando se requiere cubrir grandes extensiones de terreno, en términos de tiempo e inversión financiera su aplicación se torna poco viable (Brown 2002, Fonseca-González 2017) y no confiable, pues, de por sí, la simple extrapolación de los valores obtenidos a partir de los muestreos no incluye la heterogeneidad espacial del ambiente [p.e. Dauber et al. (2002)].
Para contrarrestar este problema, uno de los métodos más utilizados es la aplicación de algoritmos matemáticos que combinan y/o relacionan los inventarios de campo con variables ambientales generadas mediante la percepción remota, lo cual, además de disminuir el sesgo de la heterogeneidad ambiental del terreno, también permite realizar extrapolaciones a diferentes escalas (espacial o temporal); pudiendo destacarse los modelos pantropicales de densidad de biomasa y carbono generados por Baccini & Asner (2013) y Avitabile et al. (2016); y los modelos globales propuestos por Saatchi et al. (2011) y el Woods Hole Research Center (https://www.woodwellclimate.org/); este último utilizado como referencia por Andersen et al. (2016) y Cuéllar & Larrea-Alcázar (2016) para estimar las emisiones de CO2 producto de la deforestación y cambio de uso de suelo en Bolivia.
Si bien en la actualidad los modelos globales y pantropicales son reconocidos por el IPCC (2019a) como una fuente de información válida para realizar inventarios de emisiones de Gases de Efecto Invernadero (GEI), estos tienden a perder precisión e incrementan la incertidumbre de los resultados cuando son empleados en la estimación de depósitos y pérdida de biomasa y carbono a escalas geográficas menores (regional o local), principalmente por el efecto del tamaño del pixel y el enmascaramiento de la heterogeneidad ambiental y paisajística del terreno (Saatchi et al. 2011), tornándose necesario la generación de modelos ajustados a escalas regionales y/o locales, tales como los producidos por Galindo et al. (2011) para Colombia y por Durán et al. (2014) y Csillik et al. (2019) para Perú.
Considerando la constante disminución de la cobertura vegetal en Bolivia, así como la incerteza de los valores de densidad de biomasa y su distribución espacial, para el presente estudio nos planteamos los siguientes objetivos: i. modelar la distribución espacial de la densidad de la biomasa vegetal aérea a partir de la relación de datos de campo y variables ambientales (predictoras); ii. evaluar la eficiencia de predicción del modelo de biomasa vegetal aérea y compararlo con modelos previamente propuestos; y iii. estimar la cantidad de biomasa vegetal aérea almacenada en las áreas protegidas nacionales.
Métodos
Inventarios de campo y base de datos
A partir de la contribución de datos de inventarios de campo generados por los autores de esta investigación se construyó una base de datos conformada por 504 parcelas de muestreo (372 de tipo temporal y 132 de tipo permanente); las mismas que fueron establecidas y/o evaluadas entre el año 2000 y 2020, y están distribuidas sobre ocho de los nueve departamentos de Bolivia, así como en 11 de sus 13 ecorregiones [límites modificados a partir de Ibisch et al. (2003); Fig. 1].
Los elementos que se integraron en la base de datos son: i. localización geográfica de la parcela, expresada según el sistema geodésico de coordenadas geográficas WGS 84; ii. código de la parcela, siendo este de tipo alfanumérico y designado en función del lugar de procedencia o según la fuente de información; iii. superficie de muestreo de la parcela, habiéndose incluido para la estimación de la biomasa leñosa a parcelas de 0.1 a 1 ha de superficie, y de 1 m2 para la cuantificación de la biomasa gramíneo- herbácea; iv. tipo de fisionomía vegetal donde fue instalada la parcela, distinguiéndose según Villarroel et al. (2016) a fisionomías boscosas (dominada por una cobertura arbórea continua y/o ≥ 70%), sabánicas (con un estrato gramíneo- herbáceo continuo y con árboles y arbustos dispersos) y campestres (dominado por un estrato gramíneo-herbáceo, pero ocasionalmente con algunos arbustos y subarbustos, rara vez con árboles); v. identidad taxonómica de los individuos registrados en las parcelas, considerando las jerarquías de familia, género y especie, las mismas que fueron corregidas y estandarizadas en base a The Plan List v.1.1 mediante el algoritmo TPL del paquete “Taxonstand - Taxonomic Standardization of Plant Species v.2.2” (Cayuela et al. 2019); vi. diámetro del fuste de los individuos registrados en las parcelas, habiéndose incluido a todos aquellos individuos ≥ 10 cm de DAP (diámetro a la altura del pecho, 1.3 m del suelo), pero excluyendo a las lianas; vii. densidad de la madera de los individuos registrados en las parcelas, valores que fueron obtenidos de la Global Wood Density Database (Zanne et al. 2009); y viii. altura total de los individuos registrados en las parcelas, estimado según el modelo alométrico propuesto por Chave et al. (2014), cálculo que fue realizado mediante el algoritmo modelHD del paquete “BIOMASS - Estimating Aboveground Biomass and Its Uncertainty in Tropical Forests v.2.1.3” (Réjou-Méchain et al. 2017).
La ubicación geográfica de las 504 parcelas fue revisada y/o corregida (cuando fue necesario) en base a la información proporcionada por los investigadores. Esta revisión y/o corrección fue realizada mediante la plataforma Google Earth Pro v.7.3.3.7786; con lo cual también se verificó que dichos puntos de muestreo aún estuvieran en áreas con cobertura vegetal.
Estimación de la biomasa vegetal aérea (BA)
La BA de todos los individuos ≥ 10 cm de DAP (biomasa leñosa) fue estimada mediante las ecuaciones alométricas propuestas por Chave et al. (2005) para bosques secos y bosques húmedos (Tabla 1), las cuales incluyen como variables predictoras la densidad de madera (p), el DAP (D) y la altura total del individuo (H). La BA leñosa fue estimada principalmente para los puntos de muestreo localizados en fisionomías boscosas y sabánicas, así como para algunos puntos de muestreo implementados sobre fisionomías campestres (presencia de individuos ≥10 cm de DAP es rara o nula). En las fisionomías campestres la BA gramíneo-herbácea fue calculada en función del peso seco de toda la materia vegetal extraída de los cuadrantes de 1 m2 (Tabla 1).
Ecuación alométrica aplicada | Ecorregión | Fisionomía vegetal | Total Muestreo | ||
---|---|---|---|---|---|
Bosque | Sabana | Campo | |||
Bosque húmedo BA= 0.112 x (ρD2H)0.916 | Bosque del sudoeste de la Amazonia | 79 | 6 | 85 | |
Bosque tucumano - boliviano | 102 | 2 | 5 | 109 | |
Yungas | 26 | 3 | 29 | ||
Bosque seco chiquitano | 32 | 32 | |||
Bosques secos interandinos | 47 | 3 | 50 | ||
Chaco serrano | 61 | 6 | 4 | 71 | |
Bosque seco BA= exp(-2.977 + ln(ρD2H) | Gran Chaco | 46 | 4 | 50 | |
Abayoy | 3 | 11 | 14 | ||
Cerrado | 17 | 23 | 11 | 51 | |
Llanos de Moxos | 2 | 4 | 6 | ||
Pantanal | 4 | 3 | 7 | ||
Total muestreo | Total muestreo | 417 | 56 | 31 | 504 |
Las estimaciones de la BA leñosa y la BA gramíneo- herbácea fueron extrapoladas y/o uniformizadas a unidades equivalentes a toneladas por hectárea de superficie (t ha-1).
Variables predictoras
Con base a datos e información espacial producida a partir de la percepción remota, tres tipos de variables predictoras fueron obtenidas y/o generadas para Bolivia, las cuales son:
Variables bioclimáticas, seleccionando a cinco (bio 1 = temperatura promedio anual; bio 4= estacionalidad de la temperatura; bio 12= precipitación anual; bio 14= precipitación de trimestre más seco; bio 15= estacionalidad de la precipitación) de las 19 variables generadas por Fick & Hijmans (2017). Estas cinco variables bioclimáticas fueron descargadas de la plataforma Worldclim (www.worldclim.org, versión 2, 1970-2000, resolución de pixel: 1 km). La importancia de estas variables bioclimáticas para la estimación espacial de la BA fue destacada por Baccini et al. (2012), Chave et al. (2014) y Perea-Ardila (2018). Por otro lado, Saatchi et al. (2011) y Simard et al. (2011), también resaltaron la relación y el efecto que estas variables ejercen sobre la estructura y cobertura de la vegetación, características que resultan relevantes para la estimación de la BA (Saatchi et al. 2011, Baccini et al. 2012, Avitabile et al. 2016).
Variable geofísica, basada en el modelo de elevación digital MERIT [Multi-Error-Removed Improved-Terrain DEM; resolución espacial de tres segundos de arco (~ 90 m)] propuesto Yamazaki et al. (2017). Índices espectrales, habiéndose seleccionando un conjunto de métricas previamente utilizadas para el modelamiento espacial de la biomasa a nivel global y Pantropical (Saatchi et al. 2011, Baccini et al. 2012, Baccini & Asner 2013), así como en algunos países de América del Sur (Galindo et al. 2011, Baccini & Asner 2013, Durán et al. 2014). Las métricas seleccionadas fueron: i. Índice de Vegetación de Diferencia Normalizada (NDVI; Rouse et al. 1973); ii. Índice de Vegetación Mejorado (EVI; Haute et al. 2002); iii. Índice de Vegetación Ajustado al Suelo (SAVI; Huete 1998); y iv. Índice de Agua de Diferencia Normalizada (NDWI; Gao 1996). Todas estas métricas fueron calculadas y promediadas en base a series temporales de las imágenes Landsat 8 Surface Reflectance Tier 1 correspondientes a la estación seca (mayo a septiembre) de los años 2013-2018, adicionando, además, la métrica de NDVI promedio anual. La corrección atmosférica de las métricas obtenidas se realizó con el algoritmo LaSRC (Land Surface Reflectance Code); y el enmascaramiento de nubes, sombra, agua y nieve, mediante el algoritmo CFMASK (The C Function of Mask; https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC08_C01_T1_SR). La generación, corrección y descarga de las métricas fue realizada mediante la plataforma de Google Earth Engine (Gorelick et al. (2017); https://earthengine.google.com).
Una vez obtenidas las 11 variables ambientales, estas fueron reajustadas y uniformizadas a pixeles de 100 m, para así finalmente evaluar, identificar y excluir las variables predictoras con altos niveles de colinealidad (Factor de Inflación de Varianza; VIF≥10) y/o correlación (coeficiente de correlación de Pearson; r≥0.8), tal como sugiere Chen et al. (2018). Estos análisis fueron realizados mediante el algoritmo VIF del paquete “usdm - Uncertainty Analysis for Species Distribution Models v.1.1-18” (Naimi 2017).
Datos de calibración y evaluación
Generalmente los datos de calibración y evaluación para la construcción de modelos espaciales son seleccionados de forma aleatoria, y la cantidad de éstos es definida en función de porcentajes establecidos por los investigadores [p.e. Saatchi et al. (2011), Baccini et al. (2012), Baccini & Asner (2013), Avitabile et al. (2016)]. Sin embargo, considerando a alta variabilidad del contenido de BA según las fisionomías vegetales (formaciones campestres vs. formaciones sabánicas vs. formaciones boscosas), así como la heterogeneidad paisajística que caracteriza a las ecorregiones de llanura [p.e. Llanos de Moxos, el Pantanal y el Cerrado; Ibisch et al. (2003), Villarroel et al. (2016)] y montaña [p.e. Yungas, Bosque tucumano-boliviano; Ibisch et al. (2003)], la aplicación de este método podría generar o incrementar el sesgo de muestreo (selección de puntos con cantidades de BA similar y que no representan el amplio rango de valores de BA) y disminuir la representatividad espacial y ambiental (selección de puntos agregados con condiciones ambientales similares). Por tanto, siguiendo la propuesta de Saatchi et al. (2011), la base de datos fue dividida inicialmente en nueve clases o intervalos de densidad de biomasa (0.1-10 t ha-1, 10-25 t ha-1, 25-50 t ha-1, 50-75 t ha-1, 75-100 t ha-1, 100-150 t ha-1, 150-200 t ha-1, 200-300 t ha-1, >300 t ha-1). Posteriormente, a partir de cada clase y considerando una distancia mínima de 5 km entre los puntos de muestreo, se seleccionó aleatoriamente las muestras de calibración y evaluación (proceso de selección repetido 999 veces), que fue realizado mediante el algoritmo thin del paquete “spThin - Implements random spatial thinning algorithm v.0.2.0” (Aiello-Lammens et al. 2015).
Generación del modelo de BA
El modelamiento espacial de la densidad de la BA fue realizado mediante el método de aprendizaje automatizado Random Forest (Breiman 2001), método estadístico no paramétrico que cuantifica la relación de las variables predictoras (métricas ambientales) y la variable respuesta (BA; datos de calibración) mediante la generación aleatoria de múltiples árboles de desiciones, obteniendo como resultado un árbol que representa los mejores atributos de la clasificación y predicción.
Este análisis fue realizado mediante el algoritmo randomForest y predict del paquete “randomForest: Breiman and cutler's Random Forests for classification and regression v. 4.6-14” (Liaw & Wiener 2018), para lo cual se ajustó los siguientes argumentos: i. mtry, valor que representa el número adecuado de varibles predictoras a ser seleccionadas aleatoriamente para la generación de los múltiples árboles de decisiones, y con lo cual se minimiza el sesgo o error de predicción del modelo (OOB= Out-of-bag), el valor de mtry es obtenido a partir de la relación del OOB promedio y el incremento progresivo del número de variables incorporadas en los modelos; ii. nodesize, atributo que determina el número mínimo de muestras (variables respuestas) que deben ser incluidos en los nodos terminales de los árboles (nodos pequeños generan árboles con divisiones más profundas pero tienden a incrementar la varianza y sobre ajustar los modelos; nodos grandes generan árboles poco profundos y subjetivos, inducen el riesgo de sobre dimensionar las predicciones) y cuyo valor es obtenido en función de la relación del OOB promedio y el incremento progresivo del número de muestras; y iii. ntree, argumento que define el número mínimo de árboles que deben ser generados para obtener el mejor árbol de concenso, siendo obtenido en función de la relación del OOB promedio y la generación progresiva de árboles de predicción ajustados según los valores previamente calculados de mtry y ntree.
Evaluación del modelo
El rendimiento del modelo de predicción de la BA fue determinado mediante el valor del coeficiente de determinación y ajuste (R2), y el error cuadrático medio (RMSE), esperando obtener un alto valor de R2 y un bajo valor de RMSE.
Los valores de la densidad de la BA cuantificada en campo (muestras de evaluación) y la BA estimada en este estudio a partir de los datos de calibración (biomasa modelada) fueron comparadas estadísticamente mediante el test de Wilcoxon utilizando la función wilcox.test del paquete “Stats v.4.0.2” (R Core Team & contributors worldwide 2020). La evaluación del poder de predicción del modelo de BA fue analizada mediante el coeficiente de correlación de Pearson, contrastando los valores de la densidad de BA cuantificada en campo vs. la densidad de BA estimada.
La incertidumbre espacial del modelo de BA fue calculada en función de la diferencia de valores expresada en porcentaje de los valores de la BA modelada a partir de los datos de calibración vs. los valores de la BA modelada mediante los datos de evaluación.
Así también, para corroborar y comparar el desempeño del modelo de la BA con relación a otros modelos espaciales, se correlacionó (coeficiente de correlación de Pearson) los valores de la densidad de la BA cuantificada en campo vs. la densidad de la MBA estimada según los modelos propuestos por Baccini et al. (2012) y Avitabile et al. (2016).
El análisis de correlación de Pearson fue realizado mediante la función cor.test del paquete “psych: Procedures for Psychological, Psychometric, and Personality Research v.2.1.6” (Revelle 2021).
Todos los paquetes utilizados para la realización de los análisis previamente indicados fueron ejecutados sobre la plataforma R versión 4.0.2 (R Core Team 2020).
Resultados
Datos generales
De los 504 puntos de muestreo, 315 fueron empleados para la construcción del modelo de distribución espacial de la BA (datos de calibración; min= 1.9 t ha-1; max= 375.4 t ha-1) y los 189 puntos restantes para la validación del modelo (datos de evaluación; min= 2.2 t ha-1; max= 357.5 t ha-1). Por otro lado, considerando los umbrales de correlación (r≥0.8) y colinealidad (VIF≥10), sólo cinco de las 11 variables ambientales fueron conservadas para la construcción del modelo de BA, ya que éstas demostraron contener información con baja redundancia espacial (Tabla 2).
Ecuación alométrica aplicada | Ecorregión | Fisionomía vegetal | Total Muestreo | ||
---|---|---|---|---|---|
Bosque | Sabana | Campo | |||
Bosque húmedo BA= 0.112 x (ρD2H)0.916 | Bosque del sudoeste de la Amazonia | 79 | 6 | 85 | |
Bosque tucumano - boliviano | 102 | 2 | 5 | 109 | |
Yungas | 26 | 3 | 29 | ||
Bosque seco chiquitano | 32 | 32 | |||
Bosques secos interandinos | 47 | 3 | 50 | ||
Chaco serrano | 61 | 6 | 4 | 71 | |
Bosque seco BA= exp(-2.977 + ln(ρD2H) | Gran Chaco | 46 | 4 | 50 | |
Abayoy | 3 | 11 | 14 | ||
Cerrado | 17 | 23 | 11 | 51 | |
Llanos de Moxos | 2 | 4 | 6 | ||
Pantanal | 4 | 3 | 7 | ||
Total muestreo | Total muestreo | 417 | 56 | 31 | 504 |
Modelo de la biomasa vegetal aérea
El modelo de la distribución espacial de la BA alcanzó un alto nivel de predicción (R2= 0.789; RMSE= 39.4 t ha-1), siendo las métricas espectrales del EVI sec (r= 0.681±0.023) y el NDVI med (r= 0.306±0.011), y el modelo bioclimático bio 15 (r= 0.286±0.004) las variables ambientales más importantes para su construcción (Fig. 2).
La BA estimada mediante el modelo de predicción (
De acuerdo con el modelo, hasta el 2018 la cobertura vegetal remanente en Bolivia habría almacenado aproximadamente 9.77 Gt de BA (Fig. 4a), siendo los departamentos de Santa Cruz (3.51 Gt), Beni (2.58 Gt), Pando (1.41 Gt) y La Paz (1.34 Gt) donde se concentró el 90.1% del total de la BA estimada para el país.
La mayor cantidad de BA (4.33 Gt) estuvo almacenada sobre el 15.6% del territorio nacional (177.000 km2) en zonas que contienen densidad de 200-300 t h-1 (Figs. 4a y 5), las cuales se distribuyen principalmente en los departamentos del Beni (1.37 Gt), Pando (1.07 Gt), Santa Cruz (0.89 Gt) y La Paz (0.77 Gt). Por otro lado, la menor concentración de BA (0.14 Gt) estuvo en zonas cuya cobertura vegetal contienen densidades de 0 - 10 t h-1, las cuales se extienden sobre el 22.8% del territorio nacional (258 mil km2; Figs. 4a y 5), principalmente en los departamentos de Potosí (0.07 Gt), Oruro (0.03 Gt) y La Paz (0.02 Gt); así como en zonas con densidades de 10-25 t h-1 [0.11 Gt de BA distribuida sobre el 5.9% del territorio nacional (66 mil km2); Figs. 4a y 5] en los departamentos del Beni (0.04 Gt) y Santa Cruz (0.03 Gt).
Evaluación del modelo
Los valores de la BA estimada mediante el modelo presentaron una alta concordancia con relación a los valores de la BA observada en campo (r= 0.977; intervalos de confianza del valor de r entre 0.972 a 0.980; p-value <0.001), indicando así que el modelo generado posee un óptimo desempeño de predicción (Fig. 6a).
Espacialmente la incertidumbre de predicción del modelo cambió contrastantemente en función de los valores de densidad de la BA estimada (Fig. 4b; Tabla 3); siendo mayor en zonas que almacenan entre 0 - 10 t ha-1 [85.2% de su extensión total (258 mil km2) presentó una incertidumbre de predicción >50%]; y menor en zonas con densidades de BA superior a 50 t ha-1 (el 49.1% del territorio nacional presentó una incertidumbre de predicción <20%).
Superficie (%) por clases de biomasa aérea (t ha -1) | |||||||||
---|---|---|---|---|---|---|---|---|---|
Incertidumbre (%) | 0-10 5.8625 km2 | 10-25 6.6613 km2 | 25-50 83.800 km2 | 50-75 66.373 km2 | 75-100 55.114 km2 | 100-150 157.334 km2 | 150-200 95.567 km2 | 200-300 177.069 km2 | 300-360 10.274 km2 |
0-10 | 2.91 | 23.36 | 23.52 | 37.02 | 46.12 | 72.78 | 63.17 | 58.72 | 4.26 |
10-20 | 2.43 | 19.92 | 19.07 | 29.91 | 31.17 | 19.04 | 32.05 | 36.87 | 75.78 |
20-30 | 2.21 | 13.96 | 13.73 | 14.67 | 11.97 | 6.77 | 4.25 | 4.21 | 19.96 |
30-40 | 3.39 | 10.24 | 13.49 | 7.48 | 7.56 | 1.04 | 0.48 | 0.20 | 0.00 |
40-50 | 3.86 | 8.72 | 10.39 | 4.03 | 2.56 | 0.31 | 0.05 | 0.01 | 0.00 |
>50 | 85.20 | 23.79 | 19.81 | 6.88 | 0.63 | 0.06 | 0.01 | 0.00 | 0.00 |
La capacidad de predicción del modelo de BA propuesto en este estudio fue mayor (Fig. 6a) que el alcanzado por los modelos pantropicales propuestos por Baccini et al. (2012) (r= 0.555; intervalos de confianza del valor de r entre 0.477 a 0.624; p-value <0.001; Fig. 6b) y Avitabile et al. (2016) (r= 0.655; intervalos de confianza del valor de r entre 0.600-0.704; p-value <0.001; Fig. 6c); ya que ambos modelos expresan valores de densidad de BA poco concordantes con relación a los valores de la BA observada en campo.
Biomasa vegetal aérea en áreas protegidas
Del total de la BA cuantificada para Bolivia (9.77 Gt), sólo Gt (20.7%) se encuentran almacenadas dentro de las 22 áreas protegidas nacionales (185 mil km2). El 93.2% de las 2.02 Gt de BA está concentrado en 10 áreas protegidas (Fig. 7), siendo el PN y ANMI Madidi (0.34 Gt), el ANMI San Matías (0.33 Gt), el PN Noel Kempff Mercado (0.28 Gt) y el TI y PN Isiboro Sécure (0.23 Gt) los que almacenan la mayor proporción.
Discusión
El mapa de BA generado en este estudio se constituye en la primera iniciativa por cartografiar el almacenamiento de este atributo biológico y ecosistémico a nivel nacional, siendo así una fuente de información relevante y de alta precisión (nivel 2 y nivel 3) para la realización de inventarios (nacionales, regionales y/o locales) de emisiones de gases de efecto invernadero como consecuencia de la deforestación, degradación de la vegetación y el posterior uso de la tierra (IPCC 2019a); pues, por ejemplo, de acuerdo con la tercera comunicación nacional del Estado Plurinacional de Bolivia referente a “los avances y logros en la implementación de medidas que aportan al cumplimiento de los compromisos del país ante la Convención Marco de las Naciones Unidas para el Cambio Climático (CMNUCC) y el inventario de emisiones de Gases de Efecto Invernadero (GEI)” (Autoridad Plurinacional de la Madre Tierra 2020), indica que, en el 2008 el país genero aproximadamente 0.069 Gt de CO2, de los cuales, el 80.9% (0.056 Gt) fueron producidas por el sector Uso de la Tierra y Cambio del Uso de la Tierra y Silvicultura (LULUCF, por sus siglas en inglés). Sin embargo, metodológicamente, este inventario de GEI, al igual que otros previamente realizados [p.e. Apaza et al. (2009)] en Bolivia no especifican la fuente de información y/o detallan los procedimientos aplicados para la obtención de dichos resultados. Por lo que, para los futuros informes, el modelo de BA que generamos podrá ser utilizado como una fuente de información de nivel 2 (IPCC 2019a).
Si bien existen modelos pantropicales internacionalmente aceptados como una fuente de información referencial (nivel 1) para estimar (escala nacional o regional) la densidad de la BA y su relación con la pérdida de reservorios de carbono y emisiones de CO2e (Plan Vivo 2015, IPCC 2019a), éstos tienden a subestimar o sobrestimar los resultados cuando son empleados a escalas geográficas menores, ya que la resolución del tamaño de pixel generalmente no logra capturar la heterogeneidad de la cobertura vegetal ni ambiental (Saatchi et al. 2011) (Fig. 8). Por ejemplo, respecto al modelo de BA generado en este estudio (9.77 Gt;
Sin duda alguna, a nivel nacional e internacional la deforestación es la principal causante de la pérdida de los sumideros de carbono, así como es una de las principales fuentes generadoras de GEI (IPCC 2019b, Global Forest Watch 2021). De acuerdo con la ABT (2018), hasta el 2017, los cuatro departamentos identificados en este estudio como los que almacenan la mayor proporción de la BA en Bolivia (90.1%; Fig. 4a), son los que históricamente concentraron el 91% de la deforestación a nivel nacional (70.042 km2), siendo el departamento de Santa Cruz el más afectado (55.928 km2 de deforestación).
La deforestación en Bolivia se ha extendido gradualmente al interior de las áreas protegidas nacionales (resguardan 2.02 de las 9.77 Gt de la BA total), habiéndose cuantificado hasta 2010 la pérdida de aproximadamente 1201.5 km2 de bosques naturales (Servicio Nacional de Áreas Protegidas 2013) y de 395.9 km2 en el periodo 2012-2017 (ABT 2018), afectando principalmente a las áreas protegidas que almacenan una alta cantidad de BA (Fig. 7), tales como el PN y ANMI Kaa-Iya, ANMI San Matías y PN y ANMI Madidi hasta 2010, mientras que el PN y ANMI Amboró, ANMI San Matías y PN Carrasco hasta 2017.
Si bien las zonas que almacenan una densidad de BA entre 0-10 t ha-1 poseen una incertidumbre de predicción espacial mayor al 50%, éstas no llegan a afectar significativamente el valor de la BA total estimada para el país, ya que sólo representa el 1.5% de las 9.77 Gt de BA almacenadas en la cobertura vegetal remanente. El alto valor de incertidumbre de predicción en estas zonas concentradas principalmente en la ecorregión puneña de los departamentos de Potosí, Oruro y el suroeste de La Paz (Figs. 1 y 4b) es atribuido a la ausencia de puntos de muestreo de BA. Por lo que, la necesidad de desarrollar estudios de BA y carbono en estas regiones es altamente relevante para mejorar el modelo de BA y disminuir la incertidumbre de predicción.
Conclusiones
Si bien las 11 variables predictoras seleccionadas inicialmente en este estudio pueden llegar a representar la heterogeneidad ambiental del terreno en una amplia escala geográfica, no todas llegaron a ser relevantes y/o se tornaron redundantes para el modelamiento de la biomasa vegetal aérea, por lo que, la selección variable es un aspecto significativamente importante para reducir al mínimo necesario el número de variables predictoras que permitan optimizar el rendimiento e interpretación del modelo.
La combinación de inventarios de campo y variables ambientales espectrales, analizadas mediante el algoritmo de Random Forest permitió capturar las variaciones espaciales de la distribución de la densidad de la biomasa vegetal aérea, obteniendo un alto valor de rendimiento y predicción, pero con niveles de incertidumbre que cambiaron principalmente en función de la representatividad del número de puntos de muestreo.
En la actualidad, las áreas protegidas nacionales conservan una baja proporción de los reservorios y sumideros de carbono, la cuales, conjuntamente con el resto de la cobertura vegetal remanente están cada vez más amenazadas por la deforestación.
Este mapa, el cual representa el primero de su tipo para Bolivia permitirá identificar y estimar los impactos de la conservación y/o restauración/reforestación como parte de los esfuerzos nacionales para cumplir los compromisos asumidos ante la Convención Marco de las Naciones Unidas para el Cambio Climático (CMNUCC); así como también, establece las bases para cuantificar y comprender el impacto de la deforestación como el primer paso para el desarrollo económico del país vs. la pérdida de las funciones ambientales.