INTRODUCCIÓN
La descripción correcta del patrón de crecimiento en altura dominante a través de modelos matemáticos permite clasificar el potencial productivo (calidad de estación) de rodales y bosques, al determinar el índice de sitio (IS), lo que representa una herramienta fundamental para el manejo y conservación de los recursos forestales. El IS se define como la altura de los árboles dominantes o codominantes a una edad base o de referencia (Clutter et al., 1983; Skovsgaard y Vanclay, 2008).
La calidad de estación se refiere a la capacidad de un sitio forestal (rodal) para sustentar el crecimiento de cierta cantidad de árboles o de cualquier otro tipo de vegetación (Van Laar y Akça, 2007), donde la mejor calidad acumula mayor biomasa, razón por la que el potencial de un rodal está estrechamente relacionado con el volumen maderable medido en la corta final (Clutter et al., 1983).
La evaluación de la calidad del sitio requiere de un modelo que no sólo represente la relación altura-edad, sino el comportamiento de esta relación y que permita calificar la productividad de los sitios con base en su productividad mediante los patrones de crecimiento derivados de las curvas de crecimiento obtenidas a partir de esta relación funcional (Cieszewski y Bailey, 2000; Torres-Rojo, 2001).
El uso de la altura dominante ha demostrado ser uno de los elementos más idóneos para la determinación del potencial productivo de un rodal, debido a la escasa influencia que ejercen la densidad y los tratamientos silvícolas intermedios (Álvarez et al., 2004). Para cuantificar la productividad existen métodos directos e indirectos donde el IS, como método indirecto, utiliza un patrón de crecimiento de altura dominante a una edad de referencia o edad base establecida (Martín-Benito et al., 2008; Quiñonez-Barraza et al., 2015; Vanclay, 1994).
En la modelación forestal se han propuesto diversas metodologías para la construcción de funciones de altura dominante e IS. Bailey y Clutter (1974) presentaron una técnica para derivar ecuaciones dinámicas de crecimiento conocida como método de diferencia algebraica (Algebraic Difference Approach, ADA por sus siglas en inglés), que consiste en la sustitución de un parámetro del modelo base para expresarlo como una función del sitio y con esto construir curvas anamórficas y polimórficas. Esta metodología permite que los modelos sean invariantes de la edad base y del camino de simulación, con propiedades lógicas en las proyecciones del crecimiento (Bailey y Clutter, 1974).
Cieszewski y Bailey (2000) propusieron una generalización a las ecuaciones ADA, denominada método de diferencia algebraica generalizada (Generalized Algebraic Difference Approach, GADA por sus siglas en inglés), en el cual, la ecuación base puede ser expandida con más de un parámetro dependiente de la calidad de sitio y las curvas de crecimiento obtenidas son más flexibles. Las ecuaciones dinámicas tienen la estructura general Y = f (t, t0, Y0, β), donde Y0 y Y son los valores de la condición inicial y de proyección en función de los valores de la edad t (edad de proyección) y t0 (edad inicial) y β es el vector de parámetros de la ecuación (Cieszewski, 2001; 2003).
Los componentes de un sistema de predicción y proyección del área basal, diámetro, volumen y mortalidad se relacionan de manera directa con el crecimiento en altura dominante (Hall y Bailey, 2001). Se debe tener en cuenta que el patrón de crecimiento en altura dominante afectará de forma instantánea la proyección de las variables y la precisión de las estimaciones del sistema (Santiago-García et al., 2015); por tanto, una correcta estimación de la altura dominante, basada en ecuaciones de IS, no sesgadas y precisas, es fundamental para modelar el crecimiento y rendimiento de las masas forestales (Vargas-Larreta et al., 2013).
En el programa de manejo forestal de Ixtlán de Juárez autorizado en 2015 (STF, 2015) la clasificación de la productividad de los rodales de Pinus patula Schiede ex Schltdl. & Cham. se realizó con el modelo anamórfico de Chapman-Richards, el cual fue ajustado con información de sitios temporales de muestreo y el método de la curva guía; sin embargo, para lograr una estimación más eficiente del potencial productivo de los rodales se requieren datos remedidos en sitios permanentes o de análisis troncales que permitan utilizar enfoques para modelar de manera dinámica el crecimiento de la altura dominante y las tasas de cambio específicas de cada sitio forestal o por árbol individual.
Pérez-López et al. (2017) y Castillo-López et al. (2018) generaron en el área de estudio ecuaciones dinámicas de IS para P. patula con datos de sitios permanentes y de análisis troncales, respectivamente. La ecuación desarrollada por Pérez-López et al. (2017) fue generada con datos de dos inventarios consecutivos en 66 parcelas permanentes medidas en 2015 y 2016. Durante 2017 se realizó una tercera remedición de las parcelas, lo que permite ampliar el intervalo de aplicación del modelo y flexibilizar las curvas de crecimiento.
El objetivo de este estudio fue ajustar ecuaciones dinámicas de crecimiento en altura dominante e índice de sitio para rodales coetáneos de P. patula en Ixtlán de Juárez, Oaxaca y realizar una comparación de la mejor ecuación obtenida con las ecuaciones previas desarrolladas para la especie en el área de estudio. La hipótesis fue que un patrón de crecimiento con polimorfismo asintótico describe mejor el crecimiento en altura dominante de P. patula, dada la variabilidad natural de los rodales con manejo forestal.
MATERIALES Y MÉTODOS
Área de estudio
El estudio se realizó en rodales de P. patula localizados en bosques de la comunidad de Ixtlán de Juárez, Oaxaca, entre las coordenadas geográficas 17° 23’ 0.50’’ y 17° 23’ 0.58’’ latitud N y 96° 28’ 45’’ y 96° 28’ 53’’ longitud O, a una altitud promedio de 2780 m. Pinus patula es la especie que tiene mayor distribución en el predio, con aproximadamente 5000 ha como especie dominante, es de crecimiento rápido y posee un interés comercial alto para la región (Castellanos-Bolaños et al., 2008). El predio, gestionado de forma comunitaria, abarca 19,310.14 ha de superficie.
Datos dasométricos
La base de datos que se utilizó fue generada a partir de tres mediciones en 66 sitios permanentes de muestreo de 400 m2, las cuales se realizaron en 2015, 2016 y 2017; además, se utilizó una base de datos de 80 sitios de intervalo establecidos en 2007 y remedidos en 2008, con la finalidad de tener mayor representatividad del área de estudio y ampliar el intervalo de edades, que varió de cinco a 76 años.
Las variables registradas en campo fueron diámetro normal de todos los árboles dentro de la parcela con una forcípula (Mantrax Blue®, Haglöf, Sweden) y la altura total de al menos ocho individuos por sitio con clinómetro electrónico (HCC Clinometer Compass®, Haglöf, Sweden), de los cuales cuatro se identificaron como dominantes, de acuerdo con el concepto de altura dominante que emplea los 100 árboles más gruesos por hectárea (Alder, 1980; Assmann, 1970).
Los 146 sitios permitieron obtener un intervalo amplio de distribución de la muestra (5-76 años) con heterogeneidad alta, lo cual permitió mejor apalancamiento de los modelos. El conjunto de datos se caracterizó por tener: una edad promedio de 28.1 ± 16.5 años, con valores mínimo de cinco y máximo de 76 años, la altura dominante promedio fue de 21.5 ± 8.9 m, con un mínimo de 4.2 y de 39.7 m como altura máxima.
Ecuaciones dinámicas
En este estudio se utilizaron ecuaciones polimórficas con asíntotas comunes y variables (tipo ADA y GADA). Los modelos base que se ensayaron han sido documentados en investigaciones sobre biometría forestal; son flexibles, tienen un punto de inflexión y tendencia a alcanzar una asíntota horizontal (Hernández-Cuevas et al., 2018; Kiviste et al., 2002; Krumland y Eng, 2005; Rojo-Alboreca et al., 2017). En los Cuadros 1 y 2 se presentan las ecuaciones ajustadas a los datos.
Modelo base | Modelo de proyección (ADA)/Referencia |
---|---|
Korf (Lundqvist, 1957) |
Polimórfico (β1) (M1)/(Rojo-Alboreca et al., 2017) |
Hossfeld IV (Hossfeld, 1822) |
Polimórfico (β2) (M2)/(Hernández-Cuevas et al., 2018) |
Levakovic II (Kiviste et al., 2002) |
Polimórfico (β1) (M3)/(Pérez-López et al., 2017) |
Y0: altura dominante (m) a la edad t0; Y: altura dominante (m) a la edad t; t0: edad inicial (años); t: edad de proyección (años); βi: parámetros a estimar.
Modelo base | Parámetro relacionado al sitio | Solución para X0 con valores iniciales (Y0, t0)/Referencia |
---|---|---|
Bertalanffy (1949; 1957) - Richards (1959) |
|
Ecuación dinámica (M4)/(Krumland y Eng, 2005) Donde |
Hossfeld IV (Hossfeld, 1822) |
|
Ecuación dinámica (M5)/(Hernández-Cuevas et al., 2018) |
Korf (Lundqvist, 1957) |
|
Ecuación dinámica (M6)/(Rojo-Alboreca et al., 2017) |
Y0= altura dominante (m) a edad t0, Y: altura dominante (m) a la edad t, t0= edad inicial (años), t: edad de proyección (años), ai= parámetros en la ecuación base, bi= parámetros globales en la ecuación dinámica (formulación GADA).
Ajuste de los parámetros
La técnica de ajuste empleada en la estimación de parámetros de los modelos fue un método invariante con respecto a la edad de referencia (base-age invariant, BAI por sus siglas en inglés), denominado método iterativo anidado (Tait et al., 1988), el cual se utilizó como criterio para lograr la estabilización de los parámetros globales, de tal forma que el error medio cuadrático del modelo entre dos iteraciones fuera menor de 0.0001 (Vargas-Larreta et al., 2010). El método estima los efectos específicos del sitio y asume que siempre existen errores aleatorios y de medición en los datos, los cuales deben ser modelados (Cieszewski, 2003). El ajuste de las ecuaciones se realizó con el procedimiento MODEL de SAS/ETS® 9.3 (SAS Institute, 2011).
Para reducir los efectos por heterocedasticidad, la varianza del error fue asumida como una función de potencia de la altura predicha (Diéguez-Aranda et al., 2005; Huang et al., 2000). El factor de ponderación utilizado fue w i = Ĥk 1, donde k es una constante que tomó el valor de 0.5, para lograr la homogeneidad de los residuales y la consistencia de los estimadores de los parámetros (Quiñonez-Barraza et al., 2015).
Indicadores de ajuste
La bondad de ajuste de las ecuaciones se analizó mediante el análisis numérico y gráfico de los residuales. Los estadísticos utilizados fueron coeficiente de determinación ajustado por el número de parámetros (R 2 adj ), raíz del error cuadrático medio (RECM), sesgo promedio absoluto (Ē) y criterio de información de Akaike (AIC) (Basuki et al., 2009; Zhao et al., 2007); la autocorrelación de los residuales de los modelos fue medida con el estadístico Durbin-Watson (dw) (Durbin y Watson, 1951).
La significancia de los parámetros se evaluó mediante la prueba t-Student, bajo H 0 : βi = 0, bi = 0 y H 1 : βi ≠ 0, bi ≠ 0 con un valor de significancia de α = 0.05 (Infante y Zárate, 2012) y se realizó un análisis gráfico de las ecuaciones ajustadas con la finalidad de verificar su capacidad predictiva (Goelz y Burk, 1992; Sharma et al., 2011).
RESULTADOS Y DISCUSIÓN
Ajuste de modelos
Los estadísticos de bondad de ajuste para las seis ecuaciones dinámicas mostraron resultados satisfactorios, las ecuaciones explican 99.8 % de la varianza total observada en la altura dominante, con un error aproximado de 0.33 y 0.34 m para las ecuaciones ADA y GADA, respectivamente; asimismo, todos los estimadores de los parámetros resultaron diferentes de cero a un nivel de significancia del 5 % (P < 0.0001) (Cuadro 3).
Modelo | Valor | SD | Valor t | P > |t| | RECM | R2adj | AIC | Ē | DW | |
---|---|---|---|---|---|---|---|---|---|---|
ADA.M1 | β0 | 66.0688 | 0.5143 | 128.45 | < 0.0001 | 0.30 | 0.9988 | -850.36 | 0.00005 | 2.3 |
β2 | 0.7075 | 0.0051 | 138.59 | < 0.0001 | ||||||
ADA.M2 | β0 | 39.4207 | 0.0865 | 455.88 | <0.0001 | 0.38 | 0.9982 | -697.80 | 0.01194 | 2.5 |
β1 | 4.0882 | 0.0069 | 589.37 | <0.0001 | ||||||
ADA.M3 | β0 | 54.2913 | 0.1614 | 336.35 | <0.0001 | 0.30 | 0.9988 | -852.76 | 0.00032 | 2.3 |
β2 | 3.3155 | 0.0688 | 48.18 | <0.0001 | ||||||
GADA.M4 | b1 | 0.0899 | 0.0008 | 112.26 | < 0.0001 | 0.33 | 0.9986 | -781.77 | 0.00125 | 2.3 |
b2 | -8.3678 | 0.2118 | -39.5 | < 0.0001 | ||||||
b3 | 36.6950 | 0.7849 | 46.75 | < 0.0001 | ||||||
GADA.M5 | b1 | 16.2696 | 0.513 1 | 31.72 | < 0.0001 | 0.34 | 0.9986 | -778.34 | 0.00076 | 2.3 |
b2 | 2939.081 | 84.3011 | 34.86 | < 0.0001 | ||||||
b3 | 1.8165 | 0.0088 | 207.03 | < 0.0001 | ||||||
GADA.M6 | b1 | -11.2988 | 1.0641 | -10.62 | < 0.0001 | 0.36 | 0.9984 | -738.48 | 0.00279 | 2.4 |
b2 | 88.9546 | 4.4425 | 20.02 | < 0.0001 | 9 | |||||
b3 | 0.9191 | 0.0075 | 122.37 | < 0.0001 |
SD: error estándar, RECM: raíz del error cuadrático medio, R2 adj: coeficiente de determinación ajustado por el número de parámetros, Ē: sesgo promedio absoluto, AIC: criterio de información de Akaike, DW: estadístico de Durbin-Watson.
Los estadísticos de bondad de ajuste (RECM, R 2 adj , AIC) mostraron que el modelo de Korf (M1) polimórfico (β 1 ) es el más apropiado para modelar la altura dominante bajo el enfoque ADA; sin embargo, las curvas de crecimiento generadas (Figura 1a) no cubren en su totalidad los datos en las calidades de estación altas.
A pesar del ajuste preciso, los modelos Korf (M1) polimórfico (β 1 ) y Hossfeld IV (M2) polimórfico (β 2 ) presentaron algunas inconsistencias en la estimación de las alturas a edades jóvenes y, particularmente, en calidades de sitio altas (M2 = IS: 39 m). Adicionalmente, las curvas trazadas no se apegaron a la tendencia de los datos observados; por tal razón, se descartaron estas ecuaciones para modelar el crecimiento en altura dominante en los rodales de P. patula (Figuras 1a y 1b). Resultados similares fueron obtenidos por Corral-Rivas et al. (2004) y Diéguez-Aranda et al. (2005) al emplear estas ecuaciones en especies como P. cooperi C. E. Blanco, P. durangensis Martínez, P. engelmannii Carrière y P. sylvestris L.
Las curvas obtenidas con la ecuación de Levakovic II (M3) polimórfica (β 1 ) fueron mejores al tener una tendencia similar a la observada y características deseables en una curva de crecimiento (Figura 1c); por esta razón, en el enfoque ADA se eligió dicha curva como la más apropiada para representar el patrón de crecimiento en altura dominante de P. patula. Estos resultados son similares a los reportados por Santiago-García et al. (2017) y Elfving y Kiviste (1997), quienes destacaron que el modelo es flexible y predice de manera adecuada los datos observados en un sistema de crecimiento y rendimiento maderable. El estadístico de Durbin-Watson mostró que el modelo M3 (dw = 2.3) bajo el enfoque ADA no presenta problemas de autocorrelación; es decir, los estimadores de los parámetros son eficientes e insesgados; además, las curvas de crecimiento presentaron una tendencia biológica más realista al cubrir adecuadamente la variabilidad de la muestra.
Las curvas obtenidas con los modelos bajo el enfoque GADA mostraron un patrón de crecimiento más apegado a la trayectoria de los datos observados que las ecuaciones ADA. Por otro lado, la distribución de los residuales estandarizados presentó un comportamiento adecuado para ambos enfoques (Figura 2); sin embargo, los indicadores de bondad de ajuste indicaron que los modelos de Bertalanffy-Richards (M4) y de Hossfeld IV (M5) fueron los mejores al presentar características razonables y una tendencia biológica deseable en las curvas, además de cubrir completamente el intervalo de alturas de la muestra (Figuras 1d y 1e); no obstante, el modelo de BertalanffyRichards (M4) mostró mayor capacidad para modelar el patrón de crecimiento en altura dominante, y la corrección de la heterocedasticidad generó como resultado que los residuales estandarizados fueran homogéneos (Figura 2d); además, la prueba de Durbin-Watson reveló que los residuales no se encuentran correlacionados (dw = 2.3) (Cuadro 3). Estos resultados son similares a los reportados en los trabajos de Corral-Rivas et al. (2004), Vargas-Larreta et al. (2013) y Rojo-Alboreca et al. (2017), quienes obtuvieron curvas ADA y GADA con tendencia biológica deseable para estimar la productividad forestal en especies como P. cooperi C. E. Blanco, P. pseudostrobus Lindl., P. halepensis Mill., entre otras.
En los últimos años las metodologías ADA y GADA se han empleado para la modelación del crecimiento en altura dominante en el área de estudio; Hernández-Cuevas et al. (2018) ajustaron modelos de crecimiento ADA y GADA a partir de datos de análisis troncal de P. ayacahuite Ehrenb. ex Schltdl., y concluyeron que los modelos de Hossfeld IV bajo el enfoque GADA y Chapman-Richards polimórfico tipo ADA (β 2 ) fueron los más apropiados para clasificar la productividad de esta especie. Por otro lado, Castillo-López et al. (2018) generaron curvas dinámicas de índice de sitio tipo GADA con el modelo de Von Bertalanffy (1949; 1957) para P. oaxacana Mirov, P. douglasiana Martínez, P. patula Schiede ex Schltdl. & Cham. y P. pseudostrobus Lindl. a nivel de la Unidad de Manejo Forestal (UMAFOR) 2001 que abarca la región de la Sierra Norte del estado de Oaxaca, México. Estos estudios demuestran que el enfoque GADA es una de las metodologías más utilizadas en la actualidad para modelar el crecimiento en altura dominante; asimismo, el enfoque ADA ha brindado resultados satisfactorios porque las ecuaciones dinámicas construidas con esta metodología para P. patula han permitido construir tablas de rendimiento de densidad variable y simular infinidad de escenarios silvícolas (Santiago-García et al., 2017).
Las curvas generadas con los modelos M3 (ADA) y M4 (GADA) describen el crecimiento en altura dominante de manera apropiada porque cubren adecuadamente la trayectoria de los datos observados; además, las curvas generadas bajo el enfoque GADA son curvas polimórficas más flexibles, con tasa de crecimiento y potencialidades variables (Figura 3).
Con base en los indicadores de bondad de ajuste, las curvas de crecimiento de IS generadas, los gráficos de residuales, y por la estructura compacta de la ecuación, el modelo M4 es el más apropiado para modelar el patrón de crecimiento en altura dominante de los rodales de P. patula en los bosques de Ixtlán de Juárez, Oaxaca (Cuadro 3; Figuras 1d y 2d).
De esta manera, las curvas de índice de sitio generadas con el procedimiento GADA (Figura 1d) permitirán clasificar los bosques naturales de P. patula de acuerdo con su potencial productivo, y planificar las actividades silvícolas (McKenney y Pedlar, 2003) al hacer posible la predicción del crecimiento y rendimiento maderable a nivel de rodal, por categoría diamétrica y árbol individual (Santiago-García et al., 2014; Santiago-García et al., 2017).
Comparación de ecuaciones de altura dominante para Pinus patula
La obtención de datos recientes permitió desarrollar una ecuación dinámica de crecimiento de altura dominante para clasificar el potencial productivo de los rodales de P. patula en Ixtlán de Juárez, Oaxaca. La ecuación propuesta por Pérez-López et al. (2017) (Modelo previo), que expresa el índice de sitio, utiliza una edad base de 40 años, por lo que para efectos comparativos en este estudio se utilizó la misma edad de referencia; la comparación gráfica entre el modelo previo y el modelo M4 (actual) se muestra en la Figura 4. Al respecto, el modelo actual generó un patrón de crecimiento más coherente con los datos observados que el modelo previo (Figura 4).
Las curvas polimórficas y con polimorfismo asintótico generadas por Pérez-López et al. (2017) y Castillo-López et al. (2018), respectivamente, donde las primeras subestiman y las segundas sobreestiman la altura de los árboles antes y después de la edad base (Figura 4). Esto era de esperarse para el modelo previo, porque tales curvas se construyeron con menor información bajo la metodología ADA, en tanto que en el modelo de Castillo-López et al. (2018) se empleó el enfoque GADA con información de análisis troncales, pero un hecho importante a destacar es que se desarrolló a nivel de UMAFOR, por lo que la escala espacial es mucho mayor que la usada en este estudio; por tanto, se esperaría que una ecuación a nivel local permita modelar mejor las condiciones específicas del sitio.
Es importante definir la ecuación que describa mejor el comportamiento del crecimiento en altura dominante, puesto que, al interrelacionarse con un sistema de crecimiento y rendimiento permite predecir la edad óptima de cosecha; así mismo, cuando se combina con herramientas de manejo de la densidad permite definir posibilidades maderables en aclareos y optimizar el turno de los rodales (Santiago-García et al., 2013).
Con el modelo de Castillo-López et al. (2018) hay una sub- y sobreestimación a edades jóvenes; sin embargo, el modelo propuesto por Pérez-López et al. (2017) es ligeramente mejor en calidades bajas y altas (Figura 4). Este tipo de situaciones es común en la modelación forestal, debido a que los extremos climáticos anormales pueden influir más en la altura total de árboles jóvenes que en árboles adultos, y con frecuencia es determinado por factores distintos a la calidad de sitio (Huang, 1999); al respecto, Vargas-Larreta et al. (2013) señalaron que no es deseable seleccionar una edad de referencia demasiado joven para clasificar la calidad del sitio de los rodales.
La ecuación actual (M4) genera predicciones mejores en la amplitud de la base de datos utilizada, en comparación con las ecuaciones previas contrastadas (Castillo-López et al., 2018; Pérez-López et al., 2017) (Figura 4). El propósito mayor de modelar el crecimiento en altura dominante no es la predicción de la altura, sino la selección de un patrón de crecimiento de la altura dominante que el rodal puede seguir durante todo su desarrollo (Clutter et al., 1983).
CONCLUSIONES
La ecuación propuesta en este estudio bajo el enfoque GADA describe de manera adecuada el patrón de crecimiento en altura dominante para los datos de Pinus patula. Al clasificar la productividad de rodales de P. patula en los bosques de Ixtlán de Juárez, Oaxaca se recomienda el uso del modelo Bertalanffy-Richards tipo GADA (M4); esta herramienta ayudará a definir estrategias en el corto y largo plazo para planificar mejor los tratamientos silvícolas del manejo forestal y con ello lograr la conservación de los recursos forestales.