Área de estudio

El sitio se localiza en el Eje Neovolcánico Transversal entre las latitudes N 19° 21’ - 19° 34’ y longitudes O 102° 08’ - 102° 17’ (Figura 1), con altitud promedio de 2 550 m, clima templado húmedo, temperatura promedio anual de 18 °C y precipitación de 1 600 mm. La CINSJP tiene 18 318 ha de superficie total, de las cuales 9 914 ha son de producción maderable, esta última comprende 8 927.5 ha en bosque natural y 986.5 ha de plantaciones comerciales. El manejo se hace bajo el Método de Desarrollo Silvícola (MDS) con regeneración de árboles semilleros y el Método Mexicano de Ordenación de Bosques Irregulares (MMOBI) con regeneración por selección. Se tienen 10 anualidades de producción, nueve de ellas se manejan con el MDS y una con MMOBI (^{Dirección Técnica Forestal, 2017}).

Figura 1 Ubicación geográfica del área de estudio en la Comunidad Indígena de Nuevo San Juan Parangaricutiro (CINSJP), Michoacán, México. 

Datos

Las observaciones de densidad-tamaño provienen de 9 559 parcelas circulares temporales de 1 000 m² establecidas en el 2016 en bosque natural (Cuadro 1), en las que se estimó el área basal (G), densidad en número de árboles por hectárea (N), volumen y diámetro cuadrático medio DCM=40,000πGN ).

En el inventario forestal se identificaron 26 especies que fueron organizadas en cuatro grupos, con base en el género y su tolerancia a la sombra: Pinus, Quercus, otras coníferas y latifoliadas. Nueve especies de Pinus se localizaron en 98.1 % de las parcelas: Pinus pseudostrobus Lindl., Pinus montezumae Lamb., Pinus leiophylla Schltdl. & Cham., Pinus devoniana Lindley, Pinus douglasiana Martínez, Pinus teocote Schiede ex Schltdl., Pinus lawsonii Roezl., Pinus ayacahuite Ehren. y Pinus patula Schltdl. & Cham.; siete de Quercus (72.8 % de las parcelas): Quercus rugosa Neé, Quercus laurina Bonpl., Quercus candicans Neé, Quercus castanea Neé, Quercus dysophylla Benth., Quercus obtusata Bonpl. y Quercus magnoliifolia Neé; dos agrupadas en ‘otras coníferas’ (17.4 %): Abies religiosa (Kunth) Schltdl. & Cham. y Cupressus lusitanica Klotzsch ex Endl.; y ocho latifoliadas (34.6 %): Alnus jorullensis Kunth, Alnus acuminata H. B. K., Ternstroemia pringlei (Rose) Standl., Clethra mexicana DC., Arbutus xalapensis Kunth, Carpinus caroliniana Walter, Ilex pringlei Standl. y Tilia mexicana Shltdl.

En un análisis preliminar se hizo una selección de parcelas con máxima densidad. El método consiste en obtener el IDR de Reineke para cada sitio y después definir un porcentaje para seleccionar las parcelas con mayor ocupación a través del IDR (i. e. 10 %) (^{Solomon & Zhang, 2002}). En este estudio se eligieron 956 sitios con la técnica anterior, de los cuales 91 % está bajo manejo con MDS y 9 % con MMOBI. Las parcelas cubrieron una variedad de combinaciones de mezcla de especies. De acuerdo con la proporción de área basal en cada parcela, los grupos de Pinus, latifoliadas, Quercus y otras coníferas dominaron 84.73 %, 9.21 %, 3.24 % y 2.82 % del total, respectivamente. En promedio hubo de tres a cuatro especies en las parcelas seleccionadas, con un intervalo desde una hasta ocho especies.

Estimación del límite superior de la relación densidad-tamaño máxima

Se usaron dos enfoques, el primero consideró la composición del rodal con la estimación de parámetros globales, y el segundo incorporó la composición como proporción por grupo de especies. Con el primer enfoque, la RDTM se estimó con la ecuación linealizada de Reineke o modelo potencial:

ln⁡N=β0+β1∙ln⁡DCM+ϵ

donde, N es el número de árboles, DCM es el diámetro cuadrático medio, ln es el logaritmo natural, β₀ = ln (β₀) y ε es el error.

Con el mismo enfoque se usó también el modelo exponencial en su forma linealizada propuesto por ^{Quiñonez-Barraza y Ramírez-Maldonado (2019)}:

ln⁡N=β0+β1∙DCM+ϵ

El segundo enfoque se basó en el modelo potencial con proporción de especies (MP_SP) de ^{Torres-Rojo y Velázquez-Martínez (2000)} con dos supuestos:

El intercepto varía de acuerdo con la proporción por grupo de especies y presenta una misma pendiente (MP_SP1):

ln⁡N=α1∙PS1+α2∙PS2+α3∙PS3+α4∙PS4+β∙ln⁡DCMw+ϵ

donde, α _i es el intercepto ponderado por PS _i , i es el i-ésimo grupo (1 = Pinus, 2 = Quercus, 3 = otras coníferas y 4 = latifoliadas), PS _i es la proporción de área basal por grupo con respecto al área basal total de la parcela, y DCM _w es el diámetro cuadrático medio ponderado por la proporción PS _i de cada grupo de especie DCMw=DCM∙PS1+DCM∙PS2+DCM∙PS3+DCM∙PS4=DCM.

Intercepto y pendiente varían de acuerdo con la proporción o agrupación (MP_SP2); los estimadores se obtienen con la ecuación siguiente:

ln⁡N=α1∙PS1+α2∙PS2+α3∙PS3+α4∙PS4+(β1∙PS1+β2∙PS2+β3∙PS3+β4∙PS4)∙ln⁡DCMw+ϵ

Técnicas de ajuste

Los modelos base y las variantes del MP_SP se ajustaron con dos técnicas. Una de estas es la regresión de frontera estocástica (RFE) (^{Aigner et al., 1977}; ^{Meeusen & van den Broeck, 1977}) que incorpora un error compuesto (ε) con dos términos aleatorios v∼N(0,σv2) y u para el cual se usaron las formas seminormal RFE-SN ui∼N+(0,σu2) y normal-truncado RFE-NT ui∼N+0,σu2 (^{Meeusen & van den Broeck, 1977}). Los valores de σ2=σv2+σu2 y γ=σv2σ2 se obtuvieron con la parametrización de ^{Battese y Corra (1977)}. Con una prueba de razón de verosimilitud con X ² y bajo H ₀ : µ = 0 se evaluó si la distribución de SN era más apropiada que NT para µ (^{Bi, 2004}).

La otra técnica es la regresión cuantílica (RC), la cual usa cuantiles τ∈(0,1) para describir posiciones no-centrales de una distribución (^{Koenker & Bassett Jr., 1978}). La estimación de los parámetros resuelve un problema de optimización minimizando una función asimétrica con la ecuación β^τ=minβ∈R⁡∑i=1nρτ∙N-N^, donde ρτ∙N-N^=N-N^∙τ-1∙N-N^<0 es una función de verificación, de manera que:

ρτ∙N-N^=N-N^∙(τ-1),  si N-N^<0N-N^∙(τ),  si N-N^≥0;  τ=∈(0.95, 0.975, 0.99)

La pendiente estimada con los tres cuantiles usados en la RC se comparó con una segunda prueba de razón de verosimilitud.

Los ajustes se realizaron con el paquete estadístico R^® (^{R Core Team, 2020}), evaluando la significancia de los parámetros con α = 0.05. Para RFE se usó el paquete ‘frontier’ con la técnica de máxima verosimilitud, y para RC ‘quantreg’ con el algoritmo ‘simplex’. Los errores estándar e intervalos de confianza para RC se obtuvieron con el método de Hall-Sheather (^{Koenker, 2019}). El ajuste con RFE se evaluó con el criterio de información de Akaike (AIC) y el logaritmo de la verosimilitud (LogLik) (^{Quiñonez-Barraza et al., 2018}), y RC se evaluó con la Pseudo R², el AIC y el LogLik (^{Condés et al., 2017}). Además, se evaluó visualmente la dispersión de N vs. DCM con relación a la línea de ajuste respectiva.

El mejor modelo del primer enfoque se utilizó para construir dos guías de densidad que relacionaron la densidad (N) y el tamaño promedio (DCM) en escala logarítmica. A partir del valor del IDR de máxima densidad (IDR_max) del modelo potencial (IDR_R) y exponencial (IDR _Q ) se trazaron tres isolíneas que identifican tres etapas del desarrollo del rodal: autoaclareo o mortalidad inminente, crecimiento constante y crecimiento libre. Las isolíneas se proyectaron al 70 %, 40 % y 20 % del IDR_max con las ecuaciones siguientes para los modelos potencial y exponencial, respectivamente (^{Quiñonez-Barraza & Ramírez-Maldonado, 2019}):

IDRR=NDCMRDCMβ^1 IDRQ=Ne-β^1DCM-DCMR

donde, IDR es el índice de densidad del rodal, DCM_R es el diámetro cuadrático medio de referencia (i. e. DCM_R = 25 cm), e es el exponencial y β^1 es el parámetro de la pendiente del modelo potencial o exponencial según sea el caso.

Límite superior de la relación densidad-tamaño máxima sin considerar la composición en el rodal

El Cuadro 2 indica que los parámetros estimados para los modelos potencial y exponencial con las dos técnicas de regresión fueron significativos (P < 0.0001), si bien la pendiente del modelo potencial no fue diferente de -1.605 de acuerdo con la prueba t-student (α = 0.05). El intercepto y pendiente, estimados con RFE, fueron similares para las formas seminormal (SN) y normal-truncado (NT). Los valores de σv2 fueron 0.015 (potencial) y 0.018 (exponencial), con σu2 cercana a cero en ambos modelos; RFE-SN fue mejor que RFE-NT para µ de acuerdo con la prueba de razón de verosimilitud (P > 0.05). Las estimaciones de µ con RFE-NT no fueron significativas (P > 0.05) y sí cercanas a cero, por lo que su ajuste no es apropiado (^{Weiskittel et al., 2009}).

Tanto en el modelo potencial como exponencial, los interceptos y pendientes fueron mayores en RC que en RFE; al aumentar τ en el ajuste con RC, los parámetros también incrementaron en ambos modelos (Cuadro 2). La prueba de razón de verosimilitud no detectó diferencias significativas en la pendiente acorde con los cuantiles utilizados (P > 0.05).

De acuerdo con el Cuadro 3, el modelo potencial fue el mejor para analizar la RDTM con valores mayores de Loglik y Pseudo R² en los ajustes con RFE y RC, respectivamente. RFE-SN fue mejor que RFE-NT con base en el Loglik y AIC, mientras que, en la RC, los estadísticos con τ = 0.95 fueron mejores que los otros dos cuantiles tanto en el modelo potencial como en el exponencial.

La tendencia de la RDTM fue similar en RFE-SN y RFE-NT; sin embargo, las líneas de ajuste de RC se situaron por arriba de las anteriores y, conforme el valor de τ fue mayor, estas se acercaron al límite superior de los datos (Figura 2). El ajuste de RC con τ = 0.99 fue adecuado tanto para el modelo potencial como el exponencial, por su mayor acercamiento al límite de la RDTM de los sitios seleccionados con máxima densidad (puntos oscuros) y del total de parcelas (puntos grises); nueve y cuatro parcelas sobrepasaron el límite de la RDTM estimado con los modelos potencial y exponencial, respectivamente.

Figura 2 Trayectoria de la relación densidad-tamaño máxima obtenida con los modelos potencial y exponencial ajustados con regresión de cuantiles (RC) y regresión con frontera estocástica seminormal (RFE-SN) y normal-truncado (RFE-NT). DCM: diámetro cuadrático medio. 

Límite superior de la relación densidad-tamaño máxima considerando la proporción de especies

Los estimadores para MP_SP1 y MP_SP2 fueron diferentes de cero (P < 0.0001) con excepción de β^2 para MP_SP2 usando RC con τ = 0.99 (Cuadro 4). Los estimadores obtenidos con las dos variantes de RFE fueron similares, aunque diferentes a los de RC. Las pendientes β^'s ) del MP_SP1 fueron más pronunciadas en RC que en RFE. La ocupación inicial α^i en términos de N, cuando DCM se aproxima a cero, fue menor en los grupos de Pinus y Quercus que en latifoliadas y otras coníferas, para RFE y RC con τ = 0.95 y τ = 0.975; sin embargo, en RC con τ = 0.99 sucedió lo contrario. Los estimadores del intercepto del MP_SP2 fueron semejantes a los de MP_SP1, aunque con tasas de mortalidad variable por grupos de especies. Las pendientes de MP_SP2 fueron más pronunciadas para latifoliadas y otras coníferas con respecto a Pinus y Quercus.

La estimación del MP_SP2 fue mejor que MP_SP1 con un AIC más bajo con la RC, Loglik superior en RC y RFE y Pseudo R² mayor en los tres τ de RC (Cuadro 5). También, RFE-SN fue mejor que RFE-NT, y RC con τ = 0.95 mostró los mejores estadísticos de bondad de ajuste con respecto a los otros dos cuantiles para MP_SP1 y MP_SP2. Conforme a los resultados anteriores, la RC con τ = 0.95 del MP_SP2 es superior para definir la RDTM con proporción de especies.

Límite superior de la relación densidad-tamaño máxima sin considerar la composición en el rodal

En este análisis, la RC estimó adecuadamente la RDTM para los modelos potencial y exponencial. Las bondades de la RC con respecto a RFE se centran en los supuestos de la distribución del error, ya que permite hacer inferencias de forma correcta con otros cuantiles (^{Salas‐Eljatib & Weiskittel, 2018}). No obstante, la RC puede ser subjetiva porque es necesario seleccionar previamente un τ para su ajuste (^{Condés et al., 2017}; ^{Tian et al., 2021}). A pesar de lo anterior, en este estudio se buscó un límite superior que cubriera la mayor parte de las parcelas con máxima densidad, por lo que la RC se considera apropiada para definir la RDTM. La RC se ha utilizado con el modelo potencial en otros trabajos y se han obtenido resultados satisfactorios en la estimación de la RDTM, cuando se usan valores de τ cercanos a uno (^{Aguirre et al., 2018}; ^{Andrews et al., 2018}; ^{Condés et al., 2017}).

El uso de los modelos potencial y exponencial en bosques con mezcla de especies se ha analizado y comparado en pocos trabajos (^{Quiñonez-Barraza & Ramírez-Maldonado, 2019}). Debido a las diferencias en las estructuras matemáticas de ambos modelos, el valor de la pendiente varía considerablemente y resulta infructuoso hacer una comparación estadística. No obstante, a nivel biológico, los modelos muestran formas distintas que definen la dinámica alométrica entre el tamaño promedio de los árboles y el espacio de crecimiento (^{Condés et al., 2017}; ^{Pretzsch & Biber, 2005}), así como la representación de la competencia interespecífica (^{Quiñonez-Barraza & Ramírez-Maldonado, 2019}; ^{Weiskittel et al., 2009}).

Los resultados del modelo potencial mostraron tasas de mortalidad constantes durante el desarrollo del rodal con mezcla de especies. La RDTM del modelo potencial es fundamental en la identificación del tamaño promedio máximo que los rodales alcanzan antes del autoaclareo (^{Zeide, 1985}), independientemente de la densidad con que se originan (^{VanderSchaaf, 2010}). Por ello, la línea de máxima densidad del modelo potencial solo es representativa para un intervalo específico de la dinámica del rodal (^{VanderSchaaf, 2010}; ^{Zeide, 2010}), principalmente cuando ocurre el cierre completo de las copas (^{Quiñonez-Barraza et al., 2018}; ^{Zeide, 2005}). En el presente estudio, la RDTM estimada con el modelo potencial cubrió un intervalo de 25 cm a 50 cm del tamaño promedio de los árboles (Figura 2), lo cual muestra que los rodales mezclados con densidades altas alcanzan más rápido la línea de densidad máxima (^{Ningre et al., 2016}).

La estimación promedio de la RDTM con el modelo exponencial mostró un patrón distinto a la del modelo potencial, y su trayectoria cóncava cubrió la mayoría de las parcelas bajo estudio. El valor del intercepto estimado con el modelo exponencial presentó un valor más realista que el modelo potencial y describió las etapas iniciales del rodal de mejor forma para el área de estudio. La línea de máxima densidad del modelo exponencial puede considerarse como una representación promedio que caracteriza los diversos eventos de competencia (^{Quiñonez-Barraza & Ramírez-Maldonado, 2019}), conforme a la combinación de las observaciones de densidad-tamaño evaluadas.

El comportamiento curvilíneo del modelo exponencial mostró pendientes cercanas a cero cuando el rodal tiene densidades altas y un tamaño promedio pequeño, mientras que la pendiente se vuelve más inclinada conforme el tamaño de los árboles aumenta (^{Zeide, 2010}). Este comportamiento sugiere que los rodales más viejos o con árboles maduros no son tan densos como los rodales jóvenes debido a la creciente proporción de claros (^{Zeide, 1995}). La misma tendencia se encontró cuando el modelo de Reineke se modificó para describir la RDTM en forma curvilínea para rodales de Pinus elliotti var. elliotti Engelm. (^{Cao et al., 2000}), Pinus tadea (L.) (^{VanderSchaaf & Burkhart, 2008}) y Picea abies (L.) Karst. y Fagus sylvatica L. (^{Schütz & Zingg, 2010}). No obstante, existen aún pocos estudios que apliquen estas nuevas formulaciones para rodales con mezcla de especies.

Con la pendiente estimada por ambos modelos se obtuvieron los valores del IDR_max con un DCM_R igual a 25 cm (^{Quiñonez-Barraza & Ramírez-Maldonado, 2019}; ^{Reineke, 1933}). En el modelo potencial, el IDR_max fue 883 y en el modelo exponencial fue 800, valores que junto al del modelo respectivo sirvieron para delimitar las zonas de crecimiento: mortalidad inminente (IDR = 618 y 560), crecimiento constante (IDR = 309 y 280) y crecimiento libre (IDR = 177 y 160) para los modelos potencial y exponencial, respectivamente (Figura 3). Las diferencias entre las guías de densidad radican en las proyecciones lineal y curvilínea de las familias de curvas de los respectivos IDR. En el modelo potencial, las proyecciones de N decrecen linealmente conforme aumenta el tamaño promedio de los árboles, mientras que tal tendencia es curva en el modelo exponencial (^{Cao et al., 2000}; ^{Schütz & Zingg, 2010})

Figura 3 Guías de densidad basadas en el modelo potencial (a) y exponencial (b) para bosques mezclados de Nuevo San Juan Parangaricutiro, Michoacán, México. DCM: diámetro cuadrático medio 

Límite superior de la relación densidad-tamaño máxima considerando la proporción de especies

La determinación de la RDTM en rodales mezclados es compleja (^{Quiñonez-Barraza & Ramírez-Maldonado, 2019}; ^{Rivoire & Moguedec, 2012}), debido a que tienen que considerarse todas las combinaciones posibles de especies que se encuentren en un lugar en específico (^{Torres-Rojo & Velázquez-Martínez, 2000}). En este trabajo se seleccionaron parcelas que cumplieran con dicha condición; sin embargo, la carencia de observaciones que mostraran tamaños de DCM aproximados a 0 cm hasta 15 cm puede ser un factor para no tener una mejor representación de la RDTM (^{VanderSchaaf & Burkhart, 2012}).

Los parámetros estimados con el mejor modelo con proporción de especies describieron las diferencias entre las tasas de mortalidad, tolerancia y la competencia intra e intererespecífica (^{Aguirre et al., 2018}; ^{de Prado et al., 2020}; ^{Torres-Rojo & Velázquez-Martínez, 2000}), así como la relación alométrica para alcanzar el potencial máximo de una combinación de especies (^{Weiskittel et al., 2009}).

A nivel general, los grupos de otras coníferas y latifoliadas tienden a tener densidades altas en la fase inicial del rodal, pero muestran una tasa de mortalidad alta conforme se desarrollan, debido a que demandan más espacio de crecimiento en la fase madura (^{Pretzsch & Biber, 2005}), así como un desarrollo robusto del dosel (^{Zeide, 1985}). Los grupos Pinus y Quercus mostraron los valores más bajos en el intercepto, lo cual indica que en la etapa madura del rodal dominarán con una menor demanda de espacio de crecimiento (^{Pretzsch & Biber, 2005}). Estos grupos muestran capacidad alta de crecimiento y sobrevivencia en presencia de competencia interespecífica (^{de Prado et al., 2020}).

Las pendientes β _i variaron considerablemente para cada grupo de especies. Los grupos otras coníferas y latifoliadas fueron más tolerantes a la sombra, debido a que mostraron valores más negativos con respecto a Pinus y Quercus. Esto sugiere que cada una de las especies cuenta con su propia relación alométrica densidad-tamaño, demanda de recursos y espacio de crecimiento (^{Pretzsch & Biber, 2016}). ^{Weiskittel et al. (2009)} analizaron la RDTM en plantaciones y bosques naturales de tres especies y encontraron que Tsuga heterophylla (Raef.) Sarg. tuvo una pendiente más negativa con mayor tolerancia a la sombra respecto a Alnus rubra (Bong.) Carr. y Pseudotsuga mensiezii var. menziesii (Mirb.) Franco. En otros trabajos también llegaron a la conclusión que las especies con pendientes β _i pronunciadas, en rodales con mezcla de especies, son más tolerantes a la sombra (^{Andrews et al., 2018}; ^{de Prado et al., 2020}; ^{Torres-Rojo & Velázquez-Martínez, 2000}).

Las especies tolerantes a la sombra no tienen garantizada su sobrevivencia en etapas maduras del rodal, ya que esta puede ser influenciada por la autotolerancia de la especie, principalmente por la capacidad del crecimiento y la supervivencia de los árboles que compiten en el rodal (^{Zeide, 1985}). Los resultados indican que las especies menos tolerantes a la sombra pueden ser más autotolerantes que las tolerantes, característica que puede representar una ganancia por la permanencia en el rodal con menores tasas de mortalidad (^{Zeide, 2005}). No obstante, es importante analizar otras características como la copa, la productividad del sitio y la gravedad específica de la madera para entender las diferencias específicas de la relación del autoaclareo en rodales mezclados, así como su parámetro alométrico (^{Woodall et al., 2005}).