Modelización del crecimiento y la producción de plantaciones de Eucalyptus globulus Labill. en el noroeste de España

TESIS DOCTORAL MODELIZACIÓN DEL CRECIMIENTO Y LA PRODUCCIÓN DE PLANTACIONES DE Eucalyptus globulus LABILL. EN EL NOROESTE DE ESPAÑA Juan Daniel García Villabrille R.D. 99/2011 Departamento de Enxeñaría Agroforestal Escola Politécnica Superior Lugo, 2015 TESIS DOCTORAL MODELIZACIÓN DEL CRECIMIENTO Y LA PRODUCCIÓN DE PLANTACIONES DE Eucalyptus globulus LABILL. EN EL NOROESTE DE ESPAÑA Fdo: Juan Daniel García Villabrille R.D. 99/2011 Departamento de Enxeñaría Agroforestal Escola Politécnica Superior Lugo, 2015 TESIS DOCTORAL MODELIZACIÓN DEL CRECIMIENTO Y LA PRODUCCIÓN DE PLANTACIONES DE Eucalyptus globulus LABILL. EN EL NOROESTE DE ESPAÑA Juan Daniel García Villabrille Memoria para optar al grado de Doctor realizada bajo la dirección de los Doctores del Departamento de Ingeniería Agroforestal de la Universidad de Santiago de Compostela, Dr. ALBERTO ROJO ALBORECA, Dr. ULISES DIÉGUEZ ARANDA y Dr. FELIPE CRECENTE CAMPO VºBº VºBº VºBº El Director de Tesis El Director de Tesis El Director de Tesis ALBERTO ROJO ALBORECA ULISES DIÉGUEZ ARANDA FELIPE CRECENTE CAMPO Lugo, 2015 Agradecimientos Quisiera dar las gracias a todas las personas que me han ayudado en la elaboración de esta tesis, ya que sin ellos no habría sido posible. A mis tutores, Alberto, Ulises y Felipe, por saber aclarar mis dudas y tener paciencia para hacerme entender las cosas. Son muchas y muy variadas sus aportaciones, y sin duda han sido los directores que cualquiera querría. Gracias por confiar en mí y por vuestra ayuda. A mi compañero de fatigas, Horacio, por hacer que el trabajo de campo sea un placer entre tanto tojal y silveiral, y aportar el punto de calidad que precisan los siempre necesarios datos. A Manuel (Rodil), por la inestimable ayuda con los códigos de R y aguantarme como compañero de laboratorio, te deseo lo mejor, aunque no necesites suerte porque eres un fenómeno. A Germán, Elías, Jorge y Campos, quienes ayudaron en los trabajos de campo y de laboratorio, facilitando esta ardua labor. A Juan Gabriel y Fernando Castedo, que siempre me echaron un cable en las dudas que les planteé. No quisiera tampoco olvidarme de Roque, que fue la primera persona con quien hice un trabajo de investigación, y que sin duda ha influido en mi carrera. A los amigos que habiendo terminado este camino antes que yo, siempre me animaron y me ayudaron, como César, Esteban, Fernando, Felipe o Edu. A Almudena, que sin ella sería imposible hacer todo y nada, y nos anima cada día. A todos los que pusieron su granito de arena cuando hizo falta, como Covelo, Fina, Klaus von Gadow, Ana Cabaneiro, Irene, Christine, Guillermo Riesco, Ana Daría, Javier Gorgoso, Sara, Beatriz Omil o Cristina Eimil. A Fernando, que me ayudó cuando me fui a Portugal y que ha hecho un gran trabajo con la aplicación EucaTool®. También quisiera agradecer a João Azevedo mi acogida en el Instituto Politécnico de Bragança, que puso todo los medios a mi disposición para facilitar mi estancia allí. Este trabajo ha sido financiado mediante el proyecto “Modelización del crecimiento y la producción en plantaciones de brinzales de Eucalyptus globulus Labill. del NO de España (AGL2010-22308-C02-01)”, concedido por el Ministerio de Ciencia e Innovación, financiado por la Unión Europea mediante fondos FEDER. También ENCE (Energía & Celulosa) financió los trabajos durante 6 meses, permitiendo que todas las mediciones se completasen. También se agradece a las empresas ENCE, Maderas Villapol y ASEFOR, a la “Asociación de productores de madera de Cerdido” (PROMACER), a la “Federación de Asociaciones de Productores de Madera de Galicia” (PROMAGAL) y a los propietarios privados de Lourenzá, Mondoñedo, Nois, Trabada, Valdoviño, y Xinzo de Pontareas su inestimable colaboración por aportar lugares para la instalación de las parcelas de inventario. Por supuesto quiero darle las gracias a mi familia, que siempre me ha apoyado desde la lejanía de mi tierrina, a mi novia y a mis amigos de Oviedo y Lugo. Muchas gracias a tod@s, y si me he olvidado de alguien, lo siento. A mi familia Resumen Resumen RESUMO En España hai preto de 640 000 ha do xénero Eucalyptus, sendo Eucalyptus globulus Labill. a especie máis importante, cunha extensión superior as 400 000 ha, principalmente no Norte e Noreste, e cun aproveitamento anual preto dos 4 millóns de metros cúbicos. Neste traballo preséntase un modelo de crecemento e produción para rodais da especie destinados principalmente á produción de celulosa. Os datos utilizados proveñen dunha rede de parcelas, instalada e diseñada ex professo, que trata de cubrir a variedade de materiais de reprodución, idades, densidades e calidades existentes para a especie na rexión. O modelo de crecemento inclúe funcións de transición para predecir o crecemento en altura dominante e área basimétrica, así como a evolución do número de pés. Tamén se inclúen tarifas de cubicación con e sen cortiza de árbore individual e de rodal, así como ecuacións de biomasa de árbore individual e de rodal. O último compoñente do modelo é unha función de desagregación para a predición de diámetros e o número de pés en cada clase diamétrica. Durante o desenvolvemento dos módulos avaliouse a necesidade de elaborar modelos específicos segundo os materiais de reproducción (brinzais, clons F0 e clons F1), e/ou solos (esquistos e granitos) e a inclusión de variables climáticas, resultando modelos únicos en tódolos casos. Os modelos desenvolvidos aportan información de precisión, así como novas ferramentas non dispoñibles para a especie ata a data, como o modelo dinámico de crecemento, a tarifa de cubicación sen cortiza (a única materia prima para a celulosa), e as ecuacións de biomasa de rodal. As proxeccións do modelo indican que as producións de calidades medias sitúanse preto dos 150 m3/ha cun turno de 13 anos. Tódalas ferramentas desenvolvidas foron implementadas na aplicación EucaTool®, de acceso gratuito (http://app.eucatool.com) que permite o empleo do modelo de forma sinxela. PALABRAS CHAVE Eucalyptus globulus, clons de eucalipto, modelo dinámico, volumen, biomasa RESUMEN En España hay cerca de 640 000 ha del género Eucalyptus, siendo Eucalyptus globulus Labill. la especie más importante, con una extensión superior a las 400 000 ha, principalmente en el Norte y Noroeste, y con un aprovechamiento anual cercano a los 4 millones de metros cúbicos. En este trabajo se presenta un modelo de crecimiento y producción para rodales de la especie destinados principalmente a la producción de celulosa. Los datos utilizados proceden de una red de parcelas, instalada y diseñada ex professo, que trata de cubrir la variedad de materiales de reproducción, edades, densidades y calidades existentes para la especie en la región. El modelo de crecimiento incluye funciones de transición para predecir el crecimiento en altura dominante y área basimétrica, así como la evolución del número de pies. También se incluyen tarifas de cubicación con y sin corteza de árbol individual y de rodal, así como ecuaciones de biomasa de árbol individual y de rodal. El último componente del modelo es una función de desagregación para la predicción de diámetros y número de pies en cada clase diamétrica. Durante el desarrollo de los módulos se evaluó la necesidad de elaborar modelos específicos según los materiales de reproducción (brinzales, clones F0 y clones F1), y/o suelos (esquistos y granitos) y la inclusión de variables climáticas, resultando modelos únicos en todos los casos. Los modelos desarrollados aportan información de precisión, así como nuevas herramientas no disponibles para la especie hasta la fecha, como el modelo dinámico de crecimiento, la tarifa de cubicación sin corteza (la única materia prima para la celulosa), y las ecuaciones de biomasa de rodal. Las proyecciones del modelo indican que las producciones de calidades medias se sitúan próximas a los 150 m3/ha con un turno de 13 años. Todas las herramientas desarrolladas fueron implementadas en la aplicación web EucaTool®, de acceso gratuito (http://app.eucatool.com) y que permite el empleo del modelo de forma sencilla. PALABRAS CLAVE Eucalyptus globulus, clones de eucalipto, modelo dinámico, volumen, biomasa ABSTRACT In Spain, trees belonging to the genus Eucalyptus cover about 640 000 ha of land. Eucalyptus globulus Labill. is the most important species, occupying more than 400 000 ha of land, mainly in the north and northwest, and yielding annual harvest volumes of almost 4 million cubic metres. This work presents a growth and yield model for stands of the species destined for pulp production. The data used were derived from a network of plots specifically designed for the purpose and installed with the aim of representing the existing range of reproductive material, ages, densities and site indexes for the species in the region. The model includes transition functions for predicting dominant height and basal area growth, as well as tree number reduction. Equations for estimating volume over and under bark of individual trees and stands are also included, along with equations for estimating biomass of individual trees and stands. The final component of the model is a disaggregation function for predicting mean diameter and number of trees in each diameter class. The need to develop specific models for different types of reproductive material (seedlings, clones clones F0 and F1) and/or soils (schists and granites) was Justificación y objetivos 3 9. Transferencia tecnológica a través de la implementación de las herramientas desarrolladas en una aplicación de acceso gratuito en la nube, y elaboración del correspondiente manual de usuario. Estos objetivos se han planteado siguiendo la publicación “Herramientas selvícolas para la gestión forestal sostenible en Galicia” (Diéguez-Aranda et al., 2009, actualizado en 2012), ya que con esta tesis también se pretende actualizar las herramientas existentes sobre E. globulus y dar cabida a los aspectos que aún no figuraban en dicho libro. Introducción 5 2. Introducción En este capítulo se realiza, en primer lugar, una descripción somera de la especie y de su importancia económica, centrando la atención en el caso de Galicia. A continuación, se exponen las prácticas selvícolas más comunes entre los productores de eucalipto. Finalmente, se presenta la clasificación más frecuente de los modelos forestales de crecimiento y se describen brevemente los modelos más destacables en España y en el mundo referentes a Eucalyptus globulus. 2.1. Distribución e importancia económica de Eucalyptus globulus Labill. La clasificación taxonómica de Eucalyptus globulus Labill. es la siguiente: Reino: Plantae División: Magnoliophyta o Angiospermae Clase: Magnoliopsida Orden: Myrtales Familia: Myrtaceae Género: Eucalyptus Especie: Eucalyptus globulus Labill. Nombres comunes: eucalipto blanco, eucalipto azul, blue gum (en inglés). La descripción que se presenta es una simplificación de la que se puede obtener en Paiva (1997), incluyendo algunos detalles que se consideran útiles, observados durante las diversas labores de campo realizadas en este trabajo. Eucalyptus globulus Labill. en plantaciones es un árbol esbelto que puede alcanzar grandes dimensiones. Su corteza es blanquecina, gris-azulada o verdosa y es muy característica su facilidad para desprenderse del tronco. Las hojas juveniles son ovadaslanceoladas, sésiles, de tonos verdes y glaucos; las adultas son lanceolado-falciformes, con largos peciolos y de color verde oscuro. Sus frutos son cápsulas de cuatro caras, a modo de pirámide invertida, rematadas en la parte superior por un botón. Generalmente se presentan aislados y sin pedúnculo, y son más grandes que los de otras especies del género, como E. nitens o E. camaldulensis, que presentan frutos más pequeños, agrupados y con pedúnculos. Esta especie, descubierta en 1792 y descrita en 1800 por el botánico Jacques Julien Houtton de Labillardière, se distribuye de forma natural en pequeñas manchas en las Juan Daniel García Villabrille 6 regiones de Victoria (sureste de Australia) y en el litoral del sudeste y sur de Tasmania, formando masas mixtas con otros eucaliptos, sobre relieves ondulados, con altitudes comprendidas entre los 0 y 550 m.s.n.m, y precipitaciones anuales entre 600 y 1500 mm, con pocos días de heladas y temperaturas estivales moderadas (Fernández-López, 1982). Dentro del área geográfica más próxima a su lugar de origen, las plantaciones de la especie tuvieron una fuerte expansión a partir de 1995 y hasta 2005 en el oeste y sur de Australia, la región de Victoria y Nueva Zelanda, llegando a ocupar unas 450 000 ha destinadas principalmente a la producción de celulosa (Wang y Baker, 2007). En el resto del mundo, y por lo tanto como especie exótica, se comenzó a plantar a finales del siglo XIX como pantalla cortavientos y combustible a lo largo de las vías de comunicación, pero las plantaciones productivas comenzaron de manera significativa a partir de 1955 (Turnbull, 1999). Actualmente se encuentra distribuida por América del Sur (Argentina, Brasil, Bolivia, Chile y Ecuador), Europa (España, Italia y Portugal), África (Argelia, Etiopía, Sudáfrica, Sudán y Túnez) y el sudeste asiático (China). En el caso de España, el eucalipto blanco fue introducido en el siglo XIX, pero su implantación generalizada comenzó en la década de los años 40 del siglo XX, lo que facilitó el impulso de las celulosas SNIACE y ENCE, fundadas en 1944 y 1956 respectivamente. Desde entonces, la superficie plantada con eucaliptos, principalmente con E. globulus, aumentó hasta llegar a 225 234 ha en la fecha de finalización de los trabajos del 1.er Inventario Forestal Nacional (IFN1) (Ministerio de Agricultura, 1980). En Galicia, donde se encuentra la mayor superficie de eucalipto de España, se estima que existen en la actualidad más de 400 000 ha de masas del género Eucalyptus. Según el IFN4 (MMAMRM, 2011), Eucalyptus globulus es la especie más importante en lo que se refiere a superficie y producción. Así, se extracta que hay 320 775 ha de masas puras de Eucalyptus globulus, 100 246 ha de masas mixtas con Pinus pinaster Ait. y 12 895 ha de masas mixtas con Quercus robur L., repartidas principalmente a lo largo de una franja costera de unos 45-50 km en las provincias de A Coruña, Lugo y Pontevedra, y sin apenas representación de la especie en la provincia de Ourense, donde las condiciones climáticas no son adecuadas para su desarrollo. Respecto a la producción, el volumen con corteza medio aprovechado anualmente de E. globulus en España en la década de los 80 fue de unos 500 000 m³, según se indica en el IFN2 (DGCN, 1998), cifra que aumentó a más de 2 M m³ en el IFN3 (DGCN, 2002). Respecto a Galicia, el volumen con corteza aprovechado fue de 3,5-4 M m³ durante los trabajos del IFN4 (MMAMRM, 2011). La Tabla 2.1 muestra el volumen total aprovechado en Galicia según los informes de la Confederación Española de la Madera entre 2008 y 2014 (CONFEMADERA, 2008- 2014). Estas cifras quedan resumidas en el Anuario de Estadística Forestal (MAGRAMA, 2012), donde se refleja que E. globulus proporcionó en Galicia aproximadamente el 50% de la madera total aprovechada en la región para todas las especies, casi el 75% del volumen de E. globulus en España y cerca Introducción 7 del 38% del volumen de todas las especies extraído en el Estado, lo que supone un aumento desde el 29,9% que representaba en el año 2010 (MAGRAMA, 2010). Tabla 2.1. Volumen con corteza (m3) de Eucalyptus globulus (diámetro en punta delgada, dpd = 7 cm) aprovechado en Galicia para el período 2008 – 2014. Fuente: CONFEMADERA (2008-2014). Año 2008 2009 2010 2011 2012 2013 2014 Volumen 3,7 3,1 3,6 4,1 3,9 4,1 3,7 Por último, la importancia de esta especie se acentúa por el hecho de que es la principal fuente de fibra de pasta para celulosa (y la de mayor calidad) en toda la Unión Europea, siendo España y Portugal los principales productores de BEKP (bleached eucalypt kraft pulp o pasta kraft blanqueada de eucalipto). 2.2. Selvicultura y aspectos productivos Aunque se puede encontrar una guía completa para los cuidados culturales de la especie en González-Río et al. (1997) o Ruiz et al. (2008), en este epígrafe se pretenden describir las prácticas que se han llevado a cabo por los productores, y que han sido adquiridas desde distintas experiencias. No obstante, también se incluyen aquellas prácticas realizadas que coinciden con las recomendaciones expuestas en los documentos mencionados anteriormente. Dado el carácter exótico de la especie, todas las masas presentes en Galicia provienen de repoblaciones (plantaciones) o de la regeneración natural de la especie por rebrote o por colonización mediante dispersión de semillas de nuevos terrenos, considerándose asilvestrada en Galicia (Fernández-López, 1982; Díaz-Balteiro et al., 2008). El método de repoblación más utilizado en el pasado, a mediados de siglo XX, consistía en la apertura de hoyos mediante azada, en los que se depositaban varias semillas. Esto implicaba que las densidades fuesen superiores a los 2000 e incluso a los 3000 pies por hectárea, como se puede comprobar en el trabajo de Pita (1966), cuyos datos procedían de plantaciones de la empresa SNIACE. El conocimiento de la especie y los distintos ensayos llevados a cabo a partir de ese momento sirvieron para establecer nuevos y mejores métodos de repoblación y densidades de plantación, de tal manera que las densidades de plantación publicadas en Fernández-López (1985) no difieren en gran medida de las empleadas actualmente en repoblación, entre 900 y 1500 pies/ha, si bien los marcos encontrados en las parcelas de esta tesis son muy variados. Por los datos recogidos para esta tesis, se observa que las menores densidades y mejores calidades corresponden a propietarios particulares, mientras que las densidades mayores y calidades inferiores a los montes de gestión empresarial. En general responden Juan Daniel García Villabrille 8 a densidades destinadas a la producción de celulosa entre 1100 y 1500 pies/ha. Por otra parte, la madera de sierra que se obtiene de la especie es el resultado de masas que han sobrepasado el turno clásico para estas plantaciones (Nutto y Touza, 2004), establecido en torno a los 15 años bajo un criterio de máxima renta en especie en volumen, con un intervalo entre los 13 y los 24 años (Álvarez-González et al., 2005a). Algunos propietarios privados, sin embargo, prefieren alargar los turnos dentro del rango descrito, unos porque aun teniendo extensas masas opinan que la densidad de la madera se incrementa a partir de los 15 años (más o menos), otros porque intentan obtener ganancia de la producción de superficies que en la mayoría de las ocasiones no supera los 3000 m2. Como dato, la densidad de la madera es mayor cuanto más externa está en el fuste (Gerard, 1994) y cuanto más alta se encuentra (Rego, 2008). Respecto a los tratamientos selvícolas, y teniendo en cuenta que el destino fundamental de la madera de E. globulus es la celulosa (trituración), las densidades de plantación son definitivas y no se aplican claras, mientras que las podas raramente se realizan ya que la especie tiene una poda natural intensa. En el caso de masas destinadas a la obtención de madera, se intuye la necesidad de ahondar en aspectos que han sido estudiados por Nutto y Touza (2004, 2006) y que cobran importancia ante la creciente demanda de madera estructural. Para ello se precisan una serie de tratamientos selvícolas que reduzcan la densidad de plantación y que no son frecuentes en Galicia, por lo que no han sido incluidas en esta tesis. En la práctica, por tanto, los principales cuidados de la especie se centran sobre todo en la fertilización y en la eliminación de la competencia. La fertilización generalmente se realiza en el momento de la plantación, de la reposición de marras (en propietarios privados) y de la propia eliminación de la competencia por matorral (géneros Ulex y Erica, mayoritariamente), gramíneas y helechos, a la cual es muy sensible la especie, siendo recomendables desbroces periódicos, en sus distintas modalidades, hasta que la masa domine a la competencia (González-Río et al., 1997). El hecho de eliminar la competencia, además, sirve de defensa en el caso de la prevención de incendios en una especie de alta inflamabilidad (durante todo el año) y alto poder calorífico (> 5000 kcal/kg) (Hernando et al., 2004). Por otra parte, las plantaciones destinadas a la obtención de biomasa, donde prima sobre todo la cantidad en peso, y que están suscitando de nuevo un gran interés, son escasas. Por ejemplo, estudios recientes llevados a cabo en Chile por la Universidad de Concepción, estableciendo ensayos de plantaciones en turnos muy cortos (4-5 años), tanto en el primer turno como en el rebrote posterior, indican que la fertilización, la densidad de plantación y la edad de la corta en cada turno son los factores más determinantes. Con respecto a plantaciones de celulosa en primer turno, transformadas a plantaciones energéticas en segundo y sucesivos turnos, aún están por conocerse los cuidados más apropiados. En este caso, los aspectos que se consideran importantes son la tipología de plantación previa, el momento de selección y el número de brotes a Introducción 9 mantener. No obstante, el control de las enfermedades que pueden afectar a la planta en edades jóvenes, p. ej. Mycosphaerella spp., presente en casi todas las masas en Galicia y que reduce el crecimiento inicial de la especie (Otero et al., 2007), también se considera un factor a tener en cuenta. Además de los problemas por hongos en edades jóvenes, existe otro problema que afecta a las plantaciones de Eucalyptus globulus en Galicia, como es la defoliación causada por el curculiónido Gonipterus scutellatus Gyll. Esta especie, que apareció en Europa en 1976 (Italia) y que se detectó en Pontevedra en el año 1991 (Mansilla, 1992), fue capaz de extenderse a pesar del freno que supusieron los programas de control biológico iniciados en 1994, llegando a afectar a la mayoría de las masas en el año 2001 (Pérez et al., 2003). Este fenómeno fue identificado en las muestras recogidas durante los trabajos de campo de esta tesis, siendo caracterizado para su posterior consideración en los análisis de algunas de las ecuaciones que se han ajustado. En lo referente al índice de Hart-Becking (ver Anexo I) en las masas de eucalipto, la inexistencia de claras hace que sus valores desciendan paulatinamente a lo largo del turno. Normalmente sus valores son muy altos en el primer año (por encima del 100%), pero debido al rápido crecimiento inicial disminuyen fuertemente durante los primeros años, llegando a valores del 20% a partir de los 6-7 años, y quedando situados en torno al 15% a los 15 años. El alargamiento del turno hasta los 20 años hace que estos valores bajen nuevamente hasta situarse cerca del 10%. Los volúmenes unitarios de la especie aprovechados en turno clásico (14-15 años), se sitúan en torno a 0,1-0,2 m3/pie, que para una densidad media de 1250 pies/ha supone un aprovechamiento medio de 125-250 m3/ha, si bien esta cifra puede variar, como es lógico, según calidad y edad de corta, habiendo rodales que permiten obtener hasta 400- 500 m3/ha en edades superiores a los 15 años y buenas estaciones. En una referencia a datos obtenidos en masas que habían sobrepasado ampliamente el turno de máxima renta en especie, en concreto al volumen de un rodal con una edad de 40 años medido en este trabajo, se obtuvo un volumen aproximado de 830 m3/ha, donde los árboles dominantes cubicaban ente 2 y 3 m3/pie, con una densidad de 700 pies/ha. En el caso de Galicia, y tras entrevistar sobre la posterior renovación del vuelo a 11 propietarios privados (que aportaron el 25% de las parcelas de inventario de este trabajo), se constató que la mayoría optan actualmente por destoconar y repoblar de nuevo con plántulas de vivero, debido fundamentalmente a los avances en la mejora genética y a la simplicidad en los cuidados de un brinzal frente a un chirpial, todo ello aun existiendo la posibilidad de rebrote en el eucalipto. También son muchas las parcelas de gestión empresarial que son de primer turno, tanto de semilla como de variedades clonales. Por tanto, todas las herramientas desarrolladas en este trabajo se han centrado en las masas de primer turno de la especie, y además se debe tener en cuenta que el modelo existente (García y Ruiz, 2003) ya cubre las masas de rebrote en la región. Juan Daniel García Villabrille 10 2.3. Modelización del crecimiento y gestión forestal Los sistemas de gestión forestal pueden clasificarse de forma sencilla en dos grupos (Gadow, 2001): Rotation Forest Management (RFM) o sistemas de gestión por turno y Continuous Cover Forestry (CCF) o sistemas de cubierta forestal continua. Los eucaliptales del N y NO de España, aunque en su ciclo productivo no se aplican claras y generalmente tampoco se realizan podas, encajan en el grupo de los RFM. En estos sistemas se dan ciclos repetitivos de crecimiento divididos en fases más o menos comunes en todos los casos: instalación inicial de la masa (plantación, rebrote o regeneración), posibles tratamientos selvícolas durante la vida de la masa (podas, desbroces, clareos y claras), y una corta final en la que se extraen todos los pies de la masa a una edad (o turno) determinada. Tras la corta final de todos los pies de la masa, ésta se renueva por rebrote o plantación para iniciar un nuevo ciclo. Un caso particular de los sistemas RFM son los sistemas de gestión Short Rotation Forestry Management (SRFM), que se caracterizan por ciclos de corta muy reducidos y por densidades de plantación muy elevadas (a veces incluso con filas dobles), dado que suelen aplicarse a especies de rápido crecimiento que tratan de cubrir la creciente demanda de fuentes de energía renovables. Aunque existen muchas referencias al término SRFM desde principios del siglo XX, en McKay (2011) se puede encontrar una amplia revisión sobre los métodos aplicados en este tipo de sistemas. Las plantaciones de eucalipto en Galicia podrían tratarse de esta manera, aunque de momento son casi inexistentes las plantaciones diseñadas específicamente como SRFM. 2.3.1. Modelos forestales de crecimiento Los modelos forestales son abstracciones simplificadas de la realidad (Sterba, 2001), o dicho de otro modo, expresiones matemáticas que sirven para predecir el crecimiento y la producción de sistemas forestales (árboles individuales, rodales, montes, etc.) y, por tanto, útiles para la gestión a cualquier nivel (García, 1988). Los modelos forestales de crecimiento suelen comprender los cambios en el crecimiento, la mortalidad y la estructura de la masa, elementos que pueden expresarse, como se ha mencionado, mediante expresiones matemáticas, pero también mediante la lógica para relacionarlas o como el código de programación que sirve para implementarlos en un simulador informático (Diéguez-Aranda et al., 2009, p. 12). El origen de los modelos forestales se remonta a las primeras tablas de producción, elaboradas por Paulsen en 1795 (en: Assmann, 1970), y que han evolucionado hasta los actuales modelos dinámicos. Dichas tablas son cuadros numéricos que tratan de mostrar la evolución en el tiempo de las variables de un rodal forestal regular o coetáneo, en un Introducción 11 ámbito geográfico determinado, para las distintas clases de calidad de estación y las diferentes intervenciones (Madrigal, 1991). Si bien en las masas de eucalipto no se realizan intervenciones sobre la densidad y, por tanto, unas tablas de producción de existencias normales o de selvicultura media (modelos estáticos) podrían describir de una forma aceptable su evolución, los modelos dinámicos son más indicativos de la relación causa-efecto entre la densidad del rodal en un momento concreto y la evolución de las variables que definen dicho rodal (área basimétrica, volumen o mortalidad, entre otros) (Diéguez-Aranda et al., 2009, p. 16). 2.3.2. Clasificación del modelo de crecimiento desarrollado Son varias las clasificaciones de los modelos de crecimiento que se pueden encontrar en la bibliografía existente (p. ej.: García, 1988; Erviti, 1991; Vanclay, 1994, p. 5; Davis et al., 2001, p. 186), y aunque éstas no son excluyentes entre sí, atendiendo a la división que se puede encontrar en Diéguez-Aranda et al. (2009, p. 13), el modelo desarrollado en esta tesis se puede definir como un modelo de rodal, dinámico, determinístico y empírico. A continuación se hace una breve descripción de cada una de las categorías mencionadas. Los modelos de rodal son los de menor nivel de detalle, ya que emplean variables de masa. Son también los que menos información de partida requieren para estimar el probable desarrollo de las variables dasométricas mediante las funciones de transición, proporcionan buenas estimaciones y son los más robustos a largo plazo (Burkhart, 2003). Los modelos dinámicos requieren datos reales de crecimiento, obtenidos al hacer al menos dos mediciones de las mismas parcelas (Diéguez-Aranda et al., 2009, p. 16). Estos modelos permiten predecir los cambios en las variables de estado a partir de un momento cualquiera, ofreciendo resultados consistentes (Cieszewski y Bailey, 2000). Por su parte, los modelos determinísticos se caracterizan porque generan siempre el mismo resultado ante el mismo estado inicial, ya que no incluyen componentes aleatorios (Erviti, 1991). Por último, los modelos empíricos están construidos a partir de los datos experimentales recogidos en campo, y que representan la realidad en un lugar y momento determinados. Se caracterizan porque las variables son fáciles de obtener y por las buenas estimaciones cuantitativas (Diéguez-Aranda et al., 2009, p. 17). Su principal desventaja es que sólo producen buenas estimaciones si el rango de datos del gestor está dentro del correspondiente a la muestra de ajuste del modelo a utilizar. 2.3.3. Importancia de los modelos para la gestión forestal La importancia de los modelos forestales reside en la capacidad para predecir el desarrollo de un sistema forestal (árbol, grupo de árboles o rodal) y, por tanto, facilitar la Juan Daniel García Villabrille 12 selección de la mejor alternativa de gestión en función de los objetivos fijados previamente (Diéguez-Aranda et al., 2009, p. 11). Dichos objetivos, que pueden cambiar a lo largo del tiempo, deben cubrir las necesidades o demandas tanto de los propietarios como de la sociedad (Davis et al., 2001, p. 3). Por su parte, la gestión forestal incluye una fase de planificación y otra de toma de decisiones a distintos niveles (región, empresa, monte, rodal o árbol) dependiendo del alcance o el nivel de información disponible (Valsta, 1993). Todos estos niveles deben ser considerados con el fin de satisfacer las necesidades que se platean en cada uno de ellos (Diéguez-Aranda et al., 2009, p. 11). En el caso de E. globulus en Galicia, uno de los principales objetivos de la gestión es la de disponer de modelos de predicción de la producción en volumen y/o biomasa a nivel de rodal que ofrezcan información verosímil. 2.4. Modelos existentes para Eucalyptus globulus Labill. Según Wang y Baker (2007), a pesar de la importancia económica a nivel mundial de la especie, el desarrollo de modelos de crecimiento y producción para la misma no ha tenido tanta atención como ha ocurrido con otras especies, debido a la relativamente reciente expansión del mercado de la celulosa y a que la gestión está repartida entre numerosos agentes, lo cual no conduce a compartir información. Otra razón, como se deduce de la bibliografía existente, es el adelanto en recursos y conocimiento generado en los países del Hemisferio Norte, donde predominan las coníferas y frondosas autóctonas, entre las cuales no se encuentra el eucalipto, aunque fuese introducido en Europa a mediados del siglo XIX (García y Ruiz, 2003), y a pesar de su rápida expansión mundial desde mediados del siglo XX (Turnbull, 1999). En cualquier caso, los modelos disponibles para Eucalyptus globulus son muy diversos respecto al nivel de resolución, habiéndolos de árbol individual y de rodal. Una característica común en todos ellos es la gran variedad de tratamientos selvícolas que se aplican a las plantaciones y, por tanto, a las parcelas de toma de datos, que surgen como combinación de tres elementos: el desbroce, la fertilización y el marco de plantación. Por ello, se debe tener en cuenta que los modelos reflejan todos estos factores, propios de la selvicultura aplicada en cada lugar, y que son difíciles de esquematizar o estructurar, pero que se deben entender como incluidos. En el continente australiano Goodwin y Candy (1986) utilizaron datos procedentes de una única plantación para elaborar un modelo de árbol individual con un índice independiente de la distancia para dos salidas, una de altura y otra de crecimiento diametral. La distribución diamétrica fue ajustada mediante una distribución Beta y la mortalidad fue establecida a partir de valores próximos a los marcados por la regla de los 3/2 del autoaclareo (Reineke, 1933). Inions (1992) desarrolló un modelo de predicción de Introducción 13 área basimétrica y volumen de rodal con datos procedentes de parcelas temporales, y un módulo de predicción de altura dominante e índice de sitio a partir del análisis de tronco de 87 árboles. Un trabajo de mayor cobertura regional es el de Wong et al. (2000), quienes ajustaron un sistema de tres componentes para la predicción de altura dominante, área basimétrica y volumen en rodales del sureste de Australia y Gippsland en seis especies de eucaliptos, de los cuales E. globulus, E. nitens y E. viminalis fueron modelizados de forma separada. El más reciente es el modelo de clases dimensionales de Wang y Baker (2007), para plantaciones de primer turno en las zonas de Central Victoria, Gippsland y Green Triangle. Este trabajo incluye funciones de transición en altura dominante, área basimétrica y mortalidad, así como un módulo de desagregación de clases diamétricas a nivel de rodal, y a nivel de árbol individual una función de perfil, una ecuación de espesor de corteza y una relación altura-diámetro. Fuera de Australia algunos de los primeros trabajos aparecen en Portugal, siendo los modelos GLOBUS y EUSOP desarrollados en 1990 y 1993 respectivamente (Tomé et al., 1995). El modelo GLOBUS fue elaborado a partir de datos procedentes de la región central costera, y puede ser empleado para predecir el crecimiento en altura dominante y área basimétrica en plantaciones de primer turno a partir de datos de inventarios o a partir de datos de rodal (área basimétrica, densidad e índice de sitio), sin incluir ningún modelo de mortalidad. El modelo EUSOP realiza las mismas funciones, pero solamente a partir de datos de inventario, aunque puede ser empleado para primer y segundo turno. Este modelo también incluye una tasa de mortalidad del 2% por defecto. Posteriormente, se han elaborado nuevas y mejoradas versiones denominadas genéricamente GLOBULUS. Así, se desarrollaron los modelos GLOBULUS V 1.0.0 en 1998 (Tomé et al., 1998), y las versiones 2.0 (Tomé et al., 2000), 2.1 (Tomé et al., 2001) y finalmente 3.0 (Tomé et al., 2006). Las sucesivas mejoras se pueden resumir a partir de las características de la última versión, que fue desarrollada a partir de una base de datos que cubre todo el país, dividido en regiones productivas a nivel biogeográfico donde se tiene en cuenta el número de días de lluvia y la altitud, para plantaciones de primer turno así como para masas de rebrote. Este modelo está dividido en dos módulos, uno de inicialización con seis ecuaciones y otro de proyección con once ecuaciones. Como características más reseñables, GLOBULUS 3.0 incorpora la estimación de carbono, permite simular la transición entre turnos y simular el comportamiento de la masa antes de la selección de brotes. Este modelo está disponible para su uso de forma gratuita a través de la dirección web http://home.isa.utl.pt/~joaopalma/modelos/webglobulus/. También se han llevado a cabo en Portugal otros trabajos, como el Modelo SOP de Amaro (1997), de predicción de altura dominante y área basimétrica, a partir de una red de datos de plantaciones de primer turno de la empresa de celulosa portuguesa SOPORCEL, o el modelo de árbol individual GLOBTREE (Soares y Tomé, 2003), elaborado a partir de datos de parcelas permanentes, de ensayos de espaciamiento y fertilización en plantaciones de primer turno de la zona costera Centro y Norte de Material y métodos 21 parcelas instaladas inicialmente, debido a diversas causas: incendios, ataques severos de plagas o enfermedades o porque habían desaparecido parcial o totalmente tras una corta. Además, y tras un análisis exploratorio, en cada módulo (funciones de transición y/o funciones de salida) se empleó un número de parcelas diferente, ya que se excluyeron algunas de ellas por los motivos expuestos más adelante. 3.2. Árboles tipo 3.2.1. Criterios de selección En algunas de las parcelas de crecimiento se procedió a seleccionar y apear, siempre que fue posible, tres árboles tipo. Estos árboles se eligieron fuera de los límites de la parcela pero próximos a ella, en una banda de unos 20 m alrededor de su perímetro. Dichos pies se seleccionaron por su asignación a tres categorías, apeando un árbol dominante, uno intermedio y uno sumergido. Dado que se ha trabajado con tres materiales de reproducción diferentes, un objetivo fue obtener una muestra equilibrada entre ellos, de tal forma que se muestreasen, al igual que en el inventario, la mayor variedad posible respecto a calidades y edades, intentando llegar a una cifra de 90-100 árboles. También se tuvo en cuenta para su asignación que no presentasen ningún defecto de forma (bifurcaciones, deformaciones del fuste, excesiva ramosidad, guía terminal no vigorosa, etc.). Finalmente se derribaron 87 árboles tipo, 29 de cada material de reproducción (brinzal, clones F0 y clones F1), correspondientes a 13 parcelas de brinzal, 9 parcelas de clones F0 y 10 parcelas de clones F1. Dichos árboles sirvieron para la construcción de las tarifas de cubicación y ecuaciones de biomasa de árbol individual. En la Figura 3.2 se muestran los pares de datos altura-diámetro empleados para el ajuste de ecuaciones de árbol individual. 0 10 20 30 40 50 0 10 20 30 40 50 Diámetro normal (cm) Altura total (m) Inventario Brinzal Clon F0 Clon F1 0 5 10 15 20 25 30 0 5 10 15 20 25 30 Diámetro normal (cm) CD - CH Dominante Inermedio Sumergido Figura 3.2. Gráfico de alturas frente a diámetros de los árboles tipo inventariados. Izquierda: separación por materiales de reproducción sobre datos de inventario; derecha: separación por clases sociológicas sobre las clases diamétricas y de altura muestreadas (CD-CH). Juan Daniel García Villabrille 22 3.2.2. Toma de datos En monte, una vez apeados y desramados los árboles seleccionados, se midió la altura del tocón (hst) con cinta métrica desde el suelo hasta la sección de corte (con apreciación a 0,001 m), la longitud del tronco apeado desde la sección de corte hasta la punta (con apreciación a 0,1 m), se marcaron y midieron diámetros con corteza (di) a 0 m (hst), 0,3 m, 0,6 m, 1,3 m y 2 m de altura desde el tocón, y posteriormente cada metro hasta la troza de un metro completa cuyo diámetro en punta delgada (dpd) estuviese más próximo a 7 cm (tamaño comercial para celulosa) con forcípula (con precisión de 0,001 m). Una vez efectuadas y anotadas las mediciones anteriores, se procedió a trocear el árbol en trozas. Luego se midió el peso de cada troza con báscula de precisión (con apreciación de 0,02 kg), el número y peso de ramas insertas en cada troza (con apreciación de 0,02 kg) considerando dos tipos, muertas y/o vivas, y dividiendo en dos los diámetros de inserción (d): ramas con d < 2,5 cm y ramas con 2,5 cm ≤ d < 7 cm. Tras el pesaje, se extrajo una rodaja de la base de cada troza de unos 5-6 cm de altura, además de la rodaja de la parte superior de la última troza, siendo el total de rodajas de cada árbol el número de trozas más uno. Por último, se pesó toda la copa con una báscula de precisión (con apreciación de 0,02 kg), tomando un 20% en peso como submuestras, tanto de las ramas a lo largo del fuste como de la copa (Pérez-Cruzado et al., 2011), que se trasladó a laboratorio para su procesado. En laboratorio se procedió a la separación de la corteza de la madera en las rodajas, y se midió el espesor de corteza con calibre (con apreciación de 0,001m). También se subdividió la muestra de ramas con diámetro menor de 2,5 cm en tres fracciones menores: hojas, ramillos (d < 0,6 cm) y ramas finas (0,6 cm ≤ d < 2,5 cm). Todas las fracciones fueron pesadas en báscula de precisión (con apreciación de 0,001 kg). Posteriormente fueron secadas en estufa hasta pesada constante a 105°. Finalmente, se tomó el peso de cada una de las fracciones procesadas, en báscula de precisión (con apreciación de 0,001 kg). 3.3. Elaboración de los datos En el Anexo I figura la notación empleada, para lo que se siguieron las pautas establecidas en Diéguez-Aranda et al. (2009, p. 249). 3.3.1. Variables de rodal y resumen de inventarios A partir de los datos de las tres campañas de inventario se determinaron las siguientes variables dasométricas: edad (t, años), densidad (N, pies/ha), área basimétrica (G, m2/ha), diámetro medio cuadrático (dg, cm), diámetro medio aritmético ( d , cm), diámetro Material y métodos 23 dominante (d0, cm), altura media (h, m), altura dominante (H0, m) e índice de Hart- Becking (IH,%). Cuando la obtención de una variable requiere tener en cuenta el número de pies, éstos fueron siempre árboles vivos y de tamaño inventariable (ht > 1,3 m). En la Tabla 3.1 se muestra un cuadro resumen del número total de parcelas medidas, separadas por materiales de reproducción e inventario; y en la Tabla 3.2 se muestra un resumen de las principales variables de rodal para todas las parcelas medidas en cada inventario. Tabla 3.1. Número de parcelas medidas separadas por materiales parentales e inventario. Brinzal Clon F0 (Anselmo) Clon F0 (C17) Clon F1 (Odiel) Clon F1 (Tinto) Clon F1 (Piedras) Clon F1 (C6) Chirpial Total Inventario 1 73 25 1 10 4 1 1 13 128 Inventario 2 66 25 1 10 4 1 1 10 118 Inventario 3 52 25 1 10 4 0 1 10 103 Inventario 4 0 0 0 0 4 0 0 0 4 Tabla 3.2. Estadísticos descriptivos por inventario de la edad, número de pies por hectárea, área basimétrica y altura dominante para todas las parcelas en cada inventario. Inventario 1 (n = 128) Inventario 2 (n = 118) Inventario 3 (n = 99) Media Máximo Mínimo Desv. est Media Máximo Mínimo Desv. est Media Máximo Mínimo Desv. est t 10 40 1 5 11 41 2 5 11 27 3 4 N 1124 2000 595 282 1082 1680 594 261 1061 1631 595 252 G 16,2 61,1 0,2 9,0 17,1 63,0 0,9 9,0 17,6 40,4 1,1 7,5 H0 18,9 39,5 3,1 6,9 20,0 40,3 5,8 6,6 20,6 41,6 5,7 5,9 Tabla 3.2. (continuación). Sólo clon Tinto (F1). Inventario 4 (n = 4) Media Máximo Mínimo Desv. est t 6 6 6 0 N 1422 1508 1296 90 G 13,1 14,9 11,5 1,4 H0 16,2 17,9 14,3 1,5 Nota: t es la edad (años), N es el número de pies por hectárea (pies/ha), G es el área asimétrica del rodal (m2/ha) y H0 es la altura dominante (m). En adelante se hace referencia a todas estas variables dasométricas con la correspondiente notación, añadiendo el subíndice 1 ó 2, según se trate de una variable considerada de un estado inicial o final respectivamente. Así, por ejemplo, H01 y G1 se Juan Daniel García Villabrille 24 considerarían las variables altura dominante y área basimétrica en el instante inicial t1, mientras H02 y G2 son las mismas variables en el instante final t2. Respecto a la calidad de los datos, hay varios factores que influyen en la misma relacionados con la medición de las parcelas permanentes, destacando entre ellos el período de tiempo entre mediciones, la tasa de crecimiento de la especie estudiada o las influencias climáticas (García, 1988). En la práctica, la duración del período entre mediciones depende en gran medida de la tasa de crecimiento de la especie (Gadow et al., 1999). En esta tesis, el período ha sido de un año, aunque no es lo más recomendable, pero de esa manera se ha pretendido adquirir la máxima información posible y obtener al menos tres inventarios. Por otra parte, la uniformidad en los períodos de medición tiene influencia directa en los errores de estimación, siendo el caso más sencillo frente a la posibilidad de adquirir datos provenientes de intervalos de medición irregulares (García, 1988). En este sentido, las mediciones llevadas a cabo en este trabajo han sido realizadas en la misma fecha, de años sucesivos, considerando un margen de ±15 días respecto a la fecha de la primera medición. Por último, otro factor cuya influencia está relacionada con la duración del período entre mediciones son las prácticas selvícolas, fundamentalmente desbroces, al afectar éstos en gran medida a la reducción de competencia y posibilitar un mayor crecimiento (González-Río et al., 1997). Si bien la mejor forma de controlar este aspecto es realizar la medición de las parcelas en el momento de la intervención, la escasa capacidad para controlar este factor imposibilitó esta labor. 3.3.2. Variables de árbol individual A partir de los datos de campo se determinaron las variables volumen individual y biomasa individual por fracciones para todos los árboles apeados. Para el desarrollo de tarifas de cubicación individual se calculó el volumen de las trozas mediante la fórmula de Smalian. Dado que este método sobrestima el volumen real (García, 1995, p. 12), se optó por tomar cuatro diámetros (a 0 m, 0,3 m, 0,6 m y 1,3 m desde la sección de corte) en la troza basal, donde se produce una concavidad notoria. Además, resulta más fácil de medir el espesor de corteza en las extremidades de las trozas que en su centro. El volumen de cada árbol se calculó como la suma de los volúmenes de todas las trozas que los conformaban. El peso de las trozas en verde fue transformado a biomasa asignando una proporción de materia seca (MS) en cada troza igual a la media del contenido en MS de las rodajas que la delimitan mediante la siguiente expresión: h s P P MS = [3.1] Material y métodos 25 donde Ps es la media del peso seco de las rodajas superior e inferior a la troza considerada y Ph es la media del peso húmedo de las rodajas superior e inferior a la troza considerada. Para la estimación del peso seco de las trozas se evaluaron dos métodos inicialmente. Empleando trozas de dos metros con las tres rodajas correspondientes (inferior, central y superior), y tomando como dato observado la MS de la rodaja central, se estimaron los valores de MS de esta rodaja como: i) la media aritmética de la MS de las rodajas inferior y superior, y ii) la media de la MS de las rodajas inferior y superior ponderada por el volumen de las trozas inferior y superior respectivamente. La mejor estimación correspondió a la media aritmética, siendo, por tanto, el método empleado para la estimación de la MS de las trozas. Por otra parte, el análisis gráfico de la humedad a lo largo del fuste, representada por la humedad de las rodajas en altura, no reveló ninguna tendencia concreta. Para el análisis de las diferencias entre los diferentes materiales de reproducción se obtuvieron datos de 29 árboles de cada tipo (brinzales, clon F0 y clon F1), siendo las muestras equilibradas. En el caso de la roca madre, 34 árboles fueron apeados en suelos de esquistos/pizarras, mientras que 53 lo fueron en suelos graníticos. En la Tabla 3.3 se muestran los estadísticos descriptivos de todas las variables medidas o estimadas en los árboles tipo. Tabla 3.3. Estadísticos descriptivos de la muestra de árboles tipo. dn ht t phtot wtot vicc visc Brinzal (n = 29) Media 14,7 16,7 10,0 242,5 105,4 0,180 0,154 Máximo 27,3 28,1 17,0 748,0 313,0 0,582 0,486 Mínimo 3,6 5,4 2,0 4,5 1,8 0,002 0,001 Desv. est. 7,2 6,5 4,2 225,3 99,3 0,177 0,149 Clon F0 (n = 29) Media 12,3 14,1 6,9 149,4 64,6 0,106 0,091 Máximo 25,2 24,6 9,9 636,8 299,9 0,502 0,444 Mínimo 3,6 4,9 5,0 6,9 3,0 0,003 0,002 Desv. est. 5,6 5,1 1,3 161,7 72,4 0,131 0,114 Clon F1 (n = 29) Media 11,7 13,3 7,8 137,1 57,3 0,099 0,084 Máximo 23,9 22,7 9,8 519,1 256,0 0,398 0,340 Mínimo 2,4 4,0 5,0 3,6 1,5 0,001 0,001 Desv. est. 6,2 5,4 1,3 137,2 61,2 0,107 0,092 Nota: dn: diámetro normal (cm); ht: altura total (m): t: edad (años): phtot: peso verde en kg; wtot: biomasa o materia seca (MS, kg): vicc: volumen individual del fuste con corteza (m3) hasta un diámetro en punta delgada de 7 cm; visc: volumen individual del fuste sin corteza (m3) hasta un diámetro en punta delgada (dpd) de 7 cm considerando la corteza. Juan Daniel García Villabrille 26 Además del volumen y la biomasa, existen otras características (o variables) consideradas de interés, que fueron sometidas a diferentes test estadísticos. Dichas variables fueron: − Biomasa de madera relativa (respecto de la biomasa total), denotada por ww.rel. − Biomasa de corteza relativa (respecto de la biomasa total), denotada por wb.rel. − Biomasa de fracciones relativa (respecto de la biomasa total), incluyendo todas las ramas y hojas, denotada por wbr.rel. − Longitud de copa relativa (respecto de la altura total, ht), denotada por cl.rel. − Densidad estimada en verde (o aparente) de madera, denotada por dap.ma (en kg/m3). − Densidad estimada en verde (o aparente) de troza, incluyendo madera y corteza, denotada por dap.tr (en kg/m3). − Densidad básica de madera, denotada por dbas.ma (en kg/m3). − Densidad básica de troza, incluyendo madera y corteza, denotada por dbas.tr (en kg/m3). La finalidad de estos análisis tuvo una doble intención: − Aportar información de cara al momento del aprovechamiento y la compra-venta, ya que la forma habitual de cuantificación es el peso en destino. − Identificar las posibles diferencias en biomasa de las fracciones del árbol por materiales de reproducción. Del peso total en verde del árbol, la parte “pulpable” se extrae de la madera del fuste, que en el aprovechamiento es dividido en trozas de longitud variable entre 2,5 a 5 m, a veces con corteza y otras sin corteza. Dichas trozas son dispuestas en pilas, durante un período de tiempo variable hasta un máximo que queda fijado en los contratos, generalmente 48 horas, hasta que son cargadas y transportadas al centro de recepción. Puesto que el muestreo de los árboles tipo supuso un trabajo continuo durante un año y medio para toda Galicia, y las condiciones climáticas momentáneas no son homogéneas en la misma época del año en lugares distintos, el análisis del contenido en humedad no se pudo analizar según el momento de la corta. Para ello se debería desarrollar un experimento concreto, muestreando en las mismas parcelas durante todo el año a intervalos regulares. Sin embargo, y aunque no fue posible dentro de esta tesis, sí se dispuso de información procedente de una experiencia previa (Rojo-Alboreca y Rodríguez-Da Costa, 2010), consistente en muestreos mensuales de parcelas de chirpiales de E. globulus en seis localidades distintas de 1 a 6 años de edad durante un año y cuatro meses en San Sadurniño (A Coruña). En la Figura 3.3 se presentan los gráficos obtenidos de ese trabajo, donde se muestra la humedad mensual, estimada como Material y métodos 27 la media para todas las parcelas muestreadas (con edades diferentes), y las variaciones de la humedad (en base seca) para cada parcela y edad inicial. En dicho trabajo se observa una falta de tendencias en la humedad cuando se consideran distintas parcelas, aunque se intuyen variaciones estacionales para el mismo punto de muestreo. 0 10 20 30 40 50 60 70 80 90 1 2 3 4 5 6 7 8 9 10 11 12 Promedio mensual humedad 0 50 100 150 200 250 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 Figura 3.3. izquierda) eje X: meses; eje Y: Variación de la humedad mensual media; derecha) eje X: meses; eje Y: variación mensual de la humedad en base seca para la misma edad inicial y estación. Fuente: Rojo-Alboreca y Rodríguez-Da Costa (2010). A tenor de estos datos, en esta tesis se consideró oportuno inicialmente ajustar una ecuación de peso verde en el momento del apeo, asumiendo que no existe variación estacional, pero sí se deben considerar los plazos operativos en el aprovechamiento, para lo que se deberían haber pesado todas las rodajas a 24 y 48 horas. Durante este tiempo, el secado debe ser lo más parecido respecto a la temperatura de secado al aire libre en cualquier época del año. Otro factor que dificulta este ajuste es que tanto el descortezado como el transporte de la madera suponen una ventilación forzada, y son factores difíciles de controlar. Otra opción de aproximación al peso en verde en el momento del pesaje es la de transformar el volumen en verde a peso, y estimar las pérdidas de humedad a 24 o 48 h. Para ello se necesita una ecuación de volumen, desarrollada en el Capítulo 6 (Funciones de salida), la densidad aparente en verde para el fuste con o sin corteza, y las variaciones de la humedad de 0 a 48 horas. De las 1134 rodajas obtenidas en el muestreo de árboles tipo, se tomaron 52 discos de madera descortezados y 52 discos de corteza. Dichas muestras fueron dispuestas apiladas, simulando su disposición en monte, y pesadas a 24 y 48 h para estimar la pérdida de humedad durante este período de tiempo, siendo el secado a temperatura ambiente seca de 27,5ºC y temperatura ambiente húmeda de 19,5ºC, con ligera ventilación artificial. En la Tabla 3.4 se exponen los estadísticos descriptivos de la densidad por materiales de reproducción, empleando solamente la estimación hecha para las trozas que no contenían ninguna rama inserta. Juan Daniel García Villabrille 28 Esto fue así debido a que se observa una tendencia creciente de la densidad estimada cuando aumenta el número de ramas con d ≥ 2,5 cm insertas en la troza (Figura 3.4). A su vez, las inserciones de las ramas suponen un aumento de volumen cuyo peso queda registrado en la báscula pero no en las mediciones para volumen, por lo que se buscó la troza “ideal” en la medida de lo posible. Además, todas las trozas pertenecen a alturas mayores de 1,3 m y menores que la altura de la base de la troza que contiene el diámetro en punta delgada (con dpd = 7 cm), ya que la mayor parte de las trozas pertenecen a esta parte del fuste (83% aprox.). También los datos en los extremos mostraron valores muy diferentes de humedad y densidad estimada, por lo que fueron eliminados, quedando finalmente 350 trozas comprendidas entre 1,3 m y la troza anterior a la troza superior, sin ninguna rama inserta. La muestra de brinzal, con mayor edad máxima para la muestra completa de trozas, también fue reducida hasta igualar la edad máxima de las trozas en clones (10 años), equilibrando el número de muestras y las edades, para poder realizar un análisis más completo. Tabla 3.4. Descriptivos de la densidad aparente de las trozas con y sin corteza por materiales de reproducción. Variable Mat. reprod. Media Máximo Mínimo Desv. est. Densidad aparente verde troza sc (kg/m3) Brinzal1 (n = 201) 1122,1 1506,8 919,6 84,5 Brinzal2 (n = 55) 1110,0 1291,9 919,6 73,2 Clon F0 (n = 72) 1144,9 1422,5 959,2 97,9 Clon F1 (n = 77) 1163,5 1591,7 978,3 127,4 Densidad aparente verde troza cc (kg/m3) Brinzal1 (n = 201) 1111,5 1470,2 940,2 81,1 Brinzal2 (n = 55) 1099,6 1270,0 959,7 60,7 Clon F0 (n = 72) 1130,8 1394,6 952,4 92,7 Clon F1 (n = 77) 1146,2 1563,7 978,6 115,4 Densidad básica troza sc (kg/m3) Brinzal1 (n = 201) 557,2 768,6 421,1 57,8 Brinzal2 (n = 55) 513,1 620,3 421,1 42,7 Clon F0 (n = 72) 557,1 728,5 432,2 58,3 Clon F1 (n = 77) 545,3 728,2 450,6 68,6 Densidad básica troza cc (kg/m3) Brinzal1 (n = 201) 527,7 716,3 406,2 51,0 Brinzal2 (n = 55) 489,2 581,2 406,2 36,9 Clon F0 (n = 72) 529,1 675,4 414,6 52,1 Clon F1 (n = 77) 514,6 666,4 422,0 58,8 1 muestra de brinzal completa (edad máxima 17 años); 2 muestra de brinzal parcial (edad máxima 10 años) Material y métodos 29 012345 900 1000 1100 1200 1300 1400 1500 Nº ramas insertas Dens. ap. verde madera (kg/ m3) Brinzal clon F0 clon F1 1000 1100 1200 1300 1400 Material de reproducción Dens. ap madera verde (kg/ m3) 012345 400 500 600 700 Nº ramas insertas Dens. básica madera (kg/ m3) Brinzal clon F0 clon F1 400 500 600 700 800 Material de reproducción Dens. básica madera (kg/ m3) Figura 3.4. Variación de la densidad aparente y la densidad básica (kg/m3) según el número de ramas con d ≥ 2,5 cm insertas en las trozas (izquierda) y por materiales de reproducción (derecha). En la Tabla 3.5 se muestran los estadísticos descriptivos para las pérdidas de humedad en los discos. A partir de los valores medios en pérdidas de humedad, los valores medios de densidad verde y las tarifas de cubicación se puede obtener una aproximación al peso fresco en el momento del pesaje. Tabla 3.5: Pérdida de humedad en corteza y madera por horas, secado natural. Pérdida de humedad (%) n = 52 Media Máximo Mínimo Desv.est. Corteza 0-24 h 10,76 29,00 1,77 7,12 Corteza 0-48 h 22,93 45,59 5,78 11,37 Madera 0-24 h 3,53 9,91 0,95 1,92 Madera 0-48 h 8,11 19,39 3,39 3,54 Por otra parte, con el fin de completar el muestreo y atender a uno de los principales problemas presentes en las masas de esta especie, se estableció una metodología para la medición (estimada) de la defoliación causada por el gorgojo del eucalipto (Gonipterus scutellatus Gyll.) a nivel de árbol individual. Dicha metodología está basada en una escala de clasificación de la defoliación establecida por la Estación Fitopatológica do Areeiro (EFA), que establece cinco grupos de defoliación para la copa. La clasificación dentro de un grupo concreto de defoliación (def) es el resultado de una operación Juan Daniel García Villabrille 30 matemática que multiplica el grado de defoliación (grdef) de cada una de las hojas seleccionadas en los ramillos por el porcentaje de hojas afectadas (af) en dichos ramillos. La aplicación de esta metodología en campo fue la siguiente: una vez apeado el árbol se tomaron del tercio superior de la copa cuatro ramillos de unas 15-30 hojas de cuatro puntos opuestos de la misma. Se procedió de esa manera ya que el tercio superior de la copa es la zona con más acceso a la luz y, por tanto, donde con mayor probabilidad se puede encontrar el gorgojo. En cada ramillo se evaluaron dos características: defoliación de cada hoja (grdef) y afección (af). La variable grdef es un valor en porcentaje asignado de forma visual, y por tanto subjetiva, para la falta de limbo foliar en los ramillos seleccionados. Dicho valor fue asignado tras una doble evaluación (de las personas que recogieron las muestras), y su estimación es mejorable con la práctica. La variable af es también un porcentaje, cuyo valor fue asignado de forma objetiva ya que consiste en el recuento de hojas afectadas sobre el total de hojas en el ramillo. La escala final para establecer la variable def es la que se muestra en la Tabla 3.6, y en la Tabla 3.7 los valores de defoliación obtenidos para los árboles tipo. Tabla 3.6: Niveles de defoliación (Fuente: ENCE, procedente de EFA Areeiro). Clase def Intervalo Nula 0 0-10% Baja 1 11-25% Moderada 2 26-66% Alta 3 67-99% Árbol seco 4 100% Tabla 3.7: Estadísticos descriptivos para defoliación en la muestra de árboles tipo. Mat. reprod. Media Máximo Mínimo Desv. est. Brinzal 0,55 2 0 0,73 Clon F0 (Anselmo) 1,2 3 0 0,97 Clon F1 (Odiel) 1 3 0 1,2 Durante la toma de muestras en árboles tipo y en posteriores comprobaciones sobre los datos se constató que la defoliación es un fenómeno que afecta principalmente a las copas de los árboles en luz, es decir de los de mayores dimensiones, mientras los sumergidos estuvieron menos afectados de forma sistemática. Por último, esta variable se introdujo sólo en los modelos de biomasa foliar, por descartes previos sobre otras fracciones y modelos, de dos formas principales: por un lado como una variable explicativa más, y por otro como variable asociada a la expansión de un parámetro. Juan Daniel García Villabrille 38 − Unas funciones de transición que definen la tasa de cambio del sistema y que dependen del estado inicial. − Unas funciones de control que regulan las modificaciones en las variables de estado, debidas a cambios instantáneos originados por tratamientos selvícolas. En este caso se ha optado por desarrollar un modelo de masa o rodal completo, complementado por una función de desagregación de clases diamétricas. Las razones para la elección de este tipo de modelo fueron: − El tipo de información disponible: dos, tres y cuatro inventarios de parcelas permanentes espaciados un período de un año (lo que corresponde a uno, dos o tres años de crecimiento respectivamente), en los que se han tomado datos con el fin de estimar variables de rodal. − Los modelos de rodal son los más adecuados para la planificación de la gestión de plantaciones forestales, ya que representan un buen compromiso entre generalidad (entendida como su posible aplicación a un amplio rango de situaciones) y precisión de las estimaciones (García, 1988). Además, teniendo en cuenta las necesidades actuales de gestión de la especie en Galicia (ver apartados 2.2 y 2.4), este nivel de detalle se puede considerar suficiente. − Un modelo de rodal es, generalmente, el paso inicial a llevar a cabo cuando se programa el desarrollo de modelos de crecimiento para cualquier especie. − Al estar formada parte de la muestra por datos procedentes de masas de brinzales y clones, se considera oportuno analizar la diferente distribución de la producción por clases diamétricas, aunque el único destino de las plantaciones (a efectos prácticos) sea la celulosa, que no requiere una clasificación de los productos. − El mejor comportamiento estadístico de las variables de rodal frente a las de árbol individual es otro argumento más de cara a la construcción de un modelo de rodal completo. 4.2. Estructura y funciones del modelo Como ya se ha comentado, la elaboración del modelo de crecimiento de la presente tesis tiene como objetivo la proyección de determinadas variables dasométricas, de manera que pueda ayudar al gestor en la planificación y gestión de las masas de E. globulus. Para conseguir este objetivo es necesario desarrollar una serie de módulos (ecuaciones) que en conjunto constituyen el modelo global de crecimiento y producción. En el modelo desarrollado, el estado de un rodal (monoespecífico y regular) de la especie en un determinado momento (a la edad t) queda definido por tres variables, dos de ellas medidas o conocidas: la altura dominante (H0) y la densidad (N); y una tercera Metodología del modelo de crecimiento 39 que puede ser medida o predicha: el área basimétrica (G). Las funciones de transición de altura dominante, área basimétrica y mortalidad (como habitualmente se denomina a la que indica la disminución del número de pies/ha) proyectan estas variables a partir de sus valores iniciales, o vector de estado (H01, G1 y N1), y de las edades inicial y final de proyección (t1 y t2), utilizando para ello ecuaciones en diferencias algebraicas, que definen tasas de cambio en el sistema. El modelo incorpora además otras funciones, denominadas de salida, como la de desagregación diamétrica, y las tarifas de cubicación de volumen hasta un diámetro en punta delgada de 7 cm (Vcc y Vsc) y de biomasa total o por fracciones (Wtot, Ww, Wb, Wbr, Wt y Wl). En la Figura 4.1 se muestra la estructura básica del modelo global de crecimiento. En ella se ha empleado la simbología descrita en Diéguez-Aranda et al. (2009, p. 125) para cada componente del modelo. Las variables de entrada se han denotado por el símbolo . Los procesos que utilizan funciones matemáticas para transformar variables de entrada en variables intermedias o en variables de salida se han indicado mediante el símbolo . Por último, las variables resultado o de salida se han indicado mediante el símbolo . Estas variables pueden ser a su vez variables de entrada en otras funciones, en cuyo caso se han denotado mediante el símbolo . Se asume que la altura dominante de la masa permanece inalterada por los tratamientos selvícolas, sobre todo teniendo en cuenta la ausencia de claras en la especie. Por tanto, la altura dominante depende tan sólo de la edad correspondiente y de la calidad de la estación, medida por el índice de sitio. Éste a su vez se determina mediante una función de crecimiento en altura dominante, a partir del valor de la altura dominante y la edad del rodal. La ecuación de crecimiento en área basimétrica se utiliza para proyectar esta variable en el tiempo. Al no realizarse claras, se simula a partir del valor de dicha variable tomado en campo, y si no se dispone de los datos de un inventario, la proyección se hace a partir del valor predicho mediante una función de inicialización de área basimétrica. Esta función no alcanza el grado de certeza del inventario diamétrico (el más sencillo), siendo éste, por tanto, recomendable. El área basimétrica futura se estima a partir de la existente en un determinado momento, por lo que la función lleva implícita la historia de la masa. Aunque se probaron otras formas adicionales para los modelos de altura dominante y área basimétrica (parámetros diferentes según tipos de suelo y/o material de reproducción y variables climáticas), los modelos más sencillos resultaron mejores. No obstante, existen modelos para la especie que proponen parámetros para estas funciones basándose en diferencias regionales (Wang y Baker, 2007) o bien modelos que incluyen variables relacionadas con la estación, como la cota o la precipitación anual (Tomé et al., 2006). La evolución de la densidad (expresada mediante el número de pies/ha) puede plantearse a partir del valor de dicha variable en un instante cualquiera. Es decir, a cualquier edad se puede simular una densidad y comprobar los efectos sobre el resto de Juan Daniel García Villabrille 40 variables del rodal. En las parcelas medidas no se efectuaron claras, y no se consideró la mortalidad debida a fenómenos extraordinarios puntuales, por lo que la evolución de la densidad está determinada por la mortalidad natural existente (por competencia por agua, luz y nutrientes), que queda definida por una función denominada de mortalidad. Figura 4.1. Estructura básica del modelo dinámico de crecimiento desarrollado. Nota: t1, t2 = edad del rodal al inicio y al final de la proyección, respectivamente; H01, H02, N1, N2, G1, G2, V1, V2, W1, W2, C1, C2 = altura dominante, número de pies/ha, área basimétrica, volumen/ha, biomasa/ha y carbono/ha a las edades inicial t1 y final t2 de la proyección, respectivamente; dpd = diámetro en punta delgada. Fuente: Diéguez-Aranda et al. (2009). La función de desagregación permite conocer el número de árboles de cada clase diamétrica en un momento determinado. Ello posibilita una valoración económica adaptativa (Yang et al., 2004), ya que, en este caso, además de la madera destinada a celulosa se podría evaluar también la cantidad de madera destinada a biomasa (por diámetros reducidos) o a madera estructural (por diámetros muy gruesos), ambas con demanda creciente. Se debe tener en cuenta que la toma de datos se centró en aquellas V2 hasta dpd = 7cm W2 total o por fracciones C2 total o por fracciones Crecimiento = V2 - V1 Crecimiento = W2 - W1 Crecimiento = C2 - C1 V1 hasta dpd = 7cm W1 total o por fracciones C1 total o por fracciones Funciones de salida Funciones de salida Variables de estado (t2) N2 G2 Variables de estado (t1) H01 N1 G1 Funciones de transición H02 Metodología del modelo de crecimiento 41 masas que representaban el grueso del mercado, es decir, las destinadas a la producción de celulosa. El volumen y la biomasa de rodal en un momento determinado se obtienen mediante ecuaciones estáticas que incluyen el valor de la altura dominante, el área basimétrica y/o la densidad en dicho momento. Una vez que se dispone del valor de la producción total (volumen y/o biomasa) se puede estimar el crecimiento medio o corriente del rodal, tanto en volumen como en biomasa, sin más que dividir esos valores entre la edad, o bien por la diferencia entre los mismos entre dos edades consecutivas. Con esa información es posible comparar diferentes alternativas selvícolas y elegir la más adecuada en función de los objetivos de gestión propuestos. También se incluyen ecuaciones de volumen y biomasa de árbol individual, que permiten realizar cálculos con mayor grado de detalle. Dos son los casos en los que se pueden aplicar. Por un lado, cuando se dispone de un inventario pie a pie, permiten estimar el volumen o la biomasa del rodal por agregación de volúmenes individuales, siendo los valores por hectárea obtenidos más aproximados a la realidad que los obtenidos con las ecuaciones de rodal. Por otro lado, permiten estimar el volumen o la biomasa correspondiente al árbol medio de cada una de las clases diamétricas estimadas con la función de desagregación. Es necesario mencionar que la función de desagregación sólo proporciona el número de árboles de un diámetro determinado, siendo desconocida la altura de dichos árboles. Dado que esta variable es necesaria para la estimación de volumen y/o biomasa, se debe recurrir a una relación altura-diámetro que proporcione la altura total de los árboles para cada rodal. Crecente-Campo et al. (2010) desarrollaron un modelo (ver apartado 2.4. Modelos existentes para Eucalyptus globulus Labill.) que fue evaluado sobre los datos de partida de esta tesis, obteniendo valores del REMC = 1,8 m sobre la muestra de árboles tipo y REMC = 1,66 m sobre los datos de inventario. El mismo modelo ajustado a ambas muestras no obtuvo convergencia, y otras formas (lineales y no lineales) produjeron errores superiores, por lo que se recomienda como el más apropiado para estimar alturas de árboles individuales. Para utilizar este modelo se necesita conocer el diámetro del árbol, así como la altura dominante y el diámetro dominante del rodal. Funciones de transición 43 5. Funciones de transición En este capítulo se aborda el ajuste de las tres funciones de transición que componen el modelo dinámico, es decir, aquellas que predicen los valores futuros para las variables de estado del rodal: altura dominante (H0, m), área basimétrica (G, m2/ha) y número de pies por hectárea (N, pies/ha). 5.1. Introducción 5.1.1. Altura dominante La calidad de estación representa las características inherentes de un determinado lugar, caracterizado por unos factores ambientales homogéneos (bióticos, edáficos y climáticos). A su vez, la productividad de la estación, que es reflejo de la calidad, es la capacidad de dicha estación para producir volumen o biomasa para un determinado genotipo y un régimen selvícola determinado, estando determinada por las condiciones ambientales y afectada por la propia gestión (Skovsgaard y Vanclay, 2008). Por otra parte, resulta necesario hacer estimaciones precisas de la productividad que permitan determinar el turno óptimo o evaluar una especie antes de la plantación (Álvarez- González et al, 2004a), o también modelizar la producción y el crecimiento actual y futuro, definir estratos de cara al inventario o un aprovechamiento bajo el principio de la sostenibilidad (García, 1983). Debido a que no es fácil hacer estimaciones directas del volumen o la biomasa, es necesario establecer un indicador (o índice) de la calidad de estación que, en un marco de gestión práctica, debe ser cuantitativo, objetivo, fácil de obtener e independiente de la densidad del rodal (Burkhart y Tomé, 2012, p. 131). Este índice puede estar basado en factores intrínsecos a la masa (evolución de alguna variable de rodal con la edad) o en factores extrínsecos (clima, suelo, vegetación, etc.), siendo la evolución de la altura dominante uno de los más adecuados por su baja dependencia de la densidad del rodal y de los tratamientos selvícolas habituales (entre los que no se incluyen las claras por lo alto), así como por estar estrechamente relacionado con la producción en volumen. Por tal motivo, la calidad de estación se suele expresar mediante el índice de sitio (IS), definido como la altura dominante a una edad determinada llamada edad de referencia. Además de lo expuesto anteriormente, y como ya se comentó en el Capítulo 1, se debe tener en cuenta que el actual modelo de calidad de estación para E. globulus en Juan Daniel García Villabrille 44 Galicia (García y Ruiz, 2003) presenta limitaciones, por lo que se ha desarrollado una función de transición de la altura dominante. 5.1.2. Área basimétrica Otra variable de interés para la gestión forestal es el área basimétrica, la cual es una medida fundamental de la densidad, que tiene en cuenta a la vez el tamaño del árbol medio y el número de pies por unidad de superficie (Gadow y Hui, 1999), está directamente relacionada con el volumen y es una herramienta básica para planificar las intervenciones selvícolas (Gadow et al, 2001). El área basimétrica, junto con el número de pies por unidad de superficie, son las variables que se emplean con mayor frecuencia para medir la densidad de una masa forestal cuando el objetivo es predecir su crecimiento y su producción en madera, siendo muy útiles cuando se trata con rodales no aclarados de una determinada edad y calidad de estación (Clutter et al., 1983, p. 70). En el modelo dinámico de crecimiento, por tanto, se ha ajustado una función de transición de área basimétrica del rodal. También se incluye una función que permite estimar el área basimétrica inicial del rodal y tener así un punto de partida para comenzar la simulación cuando no se conoce el valor real de dicha variable, empleando para ello otras variables de estado o bien relaciones entre éstas. Esta función de inicialización proporciona estimaciones sujetas a un error mayor de lo deseable, por lo que únicamente debe utilizarse cuando el rodal no esté todavía establecido o cuando no se disponga de datos de un inventario que permitan determinar el área basimétrica real del rodal. 5.1.3. Mortalidad Por último, es necesario un modelo de mortalidad que estime la disminución del número de pies a lo largo del turno si se pretenden obtener proyecciones realistas. Esto es especialmente importante en masas no aclaradas, como en el caso del eucalipto, donde la competencia intra-específica es continua (Amateis et al., 1997); sin embargo, la mortalidad es uno de los componentes menos comprendidos de los procesos naturales de crecimiento (Álvarez-González et al., 2004b). Vanclay (1994, p. 173) distingue dos tipos de mortalidad natural: regular y catastrófica. La mortalidad regular, o auto-aclareo, es debida a la competencia por luz, agua y nutrientes en un rodal (Peet y Christensen, 1987), pero también por la normal incidencia de plagas, enfermedades y condiciones climatológicas adversas. La mortalidad catastrófica es el resultado de perturbaciones aleatorias o riesgos tales como fuego, viento, nieve o ataques por insectos fuera de lo normal. La función de transición que permite estimar la disminución del número de pies, denominada función de mortalidad, sólo tiene en cuenta la mortalidad regular. Funciones de transición 45 Si bien existe la posibilidad de emplear modelos que estiman la probabilidad de supervivencia, o mortalidad, de los árboles de una masa (Woollons, 1998), teniendo en cuenta que el período de medición de esta tesis fue de sólo dos años, no se considera apropiado el desarrollo de este tipo de modelos. Además, trabajos previos (Diéguez- Aranda et al., 2005a; Tomé et al., 2006; Álvarez-González et al., 2010; Gómez-García et al., 2014 y 2015) han demostrado que las ecuaciones dinámicas ofrecen un buen resultado, obviando la naturaleza binomial de la mortalidad. 5.2. Material y Métodos 5.2.1. Datos Para desarrollar el modelo de crecimiento y producción se consideraron finalmente todas las parcelas de brinzales, clon Anselmo (F0), clon Odiel (F1) y clon Tinto (F1) que se inventariaron al menos dos veces, con edades inferiores a 21 años y que no presentaron evoluciones anómalas en las variables altura dominante, área basimétrica o número de pies. Así, se emplearon 16 parcelas de crecimiento medidas en dos ocasiones, 81 parcelas medidas tres veces y 4 parcelas medidas en cuatro ocasiones, lo que se tradujo en 279 intervalos de crecimiento. A partir de las mediciones realizadas en los diferentes inventarios de las parcelas se calcularon las variables de rodal necesarias para realizar los diferentes ajustes (Tabla 5.1). Del total de parcelas, 59 corresponden a suelos de granitos y 38 a suelos de pizarras/esquistos. Como consideración, los datos, y por tanto los modelos ajustados, reflejan la influencia de daños por defoliación (provocados por Mycosphaerella spp. y/o Gonipterus scutellatus Gyll.), fenómeno presente en mayor o menor grado en todas las masas muestreadas y que puede reducir el crecimiento (Merino et al., 2003). Dicho fenómeno se localiza, como se pudo comprobar durante la toma de datos, principalmente en la parte superior de las copas, situadas en luz, por lo que en muchas ocasiones reduce, al menos, la altura total de los árboles según la época y la intensidad del ataque. Juan Daniel García Villabrille 46 Tabla 5.1. Estadísticos descriptivos por inventario de las variables de estado de las parcelas empleadas en el ajuste de las funciones de transición. Brinzal Inventario 1 (n = 59) Inventario 2 (n = 59) Inventario 3 (n = 47) Media Máximo Mínimo Desv. est. Media Máximo Mínimo Desv. est. Media Máximo Mínimo Desv. est. t 12 20 1 4 13 21 2 4 13 20 3 3 N 1068 1825 661 253 1032 1670 661 240 1007 1631 625 246 G 18,5 38,6 0,2 8,8 19,5 40,8 0,9 8,8 19,7 40,4 1,1 9,0 H0 20,7 39,1 3,1 7,2 21,8 40,3 5,9 6,9 21,9 41,6 5,7 6,9 Clon F0 (Anselmo) Inventario 1 (n = 24) Inventario 2 (n = 24) Inventario 3 (n = 22) Media Máximo Mínimo Desv. est. Media Máximo Mínimo Desv. est. Media Máximo Mínimo Desv. est. t 6,8 11 4 2 7,8 12 5 2 9,1 13 6 2 N 1145 1540 780 206 1115 1540 772 202 1076 1540 731 203 G 12,2 19,5 3,4 4,7 13,3 20,7 4,7 4,7 14,6 21,3 6,7 4,4 H0 15,1 23,9 7,7 4,0 16,6 24,3 8,6 3,9 17,9 25,1 11,2 3,6 Clon F1 (Odiel y Tinto) Inventario 1 (n = 14) Inventario 2 (n = 14) Inventario 3 (n = 12) Media Máximo Mínimo Desv. est. Media Máximo Mínimo Desv. est. Media Máximo Mínimo Desv. est. t 5 8 3 2 6 9 4 2 7,17 10 5 2 N 1154 1534 800 268 1135 1508 760 265 1148 1508 720 276 G 8,2 17,0 3,0 4,3 9,8 18,7 4,6 4,3 13,3 21,7 8,8 3,5 H0 12,5 19,7 7,2 4,1 14,8 20,4 9,0 3,8 17,9 23,4 12,4 3,4 Tabla 5.1. (continuación). Clon F1 (Tinto) Inventario 4 (n = 4) Media Máximo Mínimo Desv. Est. t 6 6 6 0 N 1422 1508 1296 90 G 13,1 14,9 11,5 1,4 H0 16,2 17,9 14,3 1,5 Nota: t es la edad (años), N es el número de pies por hectárea (pies/ha), G es el área asimétrica del rodal (m2/ha) y H0 es la altura dominante (m). 5.2.2. Desarrollo de los modelos La alternativa más empleada para el desarrollo de modelos de crecimiento es utilizar ecuaciones dinámicas de la forma Y2 = f (Y1, t1, t2), donde Y2 es el valor de la función a la Funciones de transición 47 edad t2 e Y1 es el valor de la misma función a la edad t1. Estas ecuaciones dinámicas cumplen dos propiedades esenciales: (i) son invariantes con respecto a la edad de referencia (base-age invariance, BAI), es decir, cualquier par de puntos (t1, Y1) de una curva concreta puede emplearse para hacer predicciones a la edad t2 y las estimaciones obtenidas serán siempre las mismas; y (ii) son invariantes con respecto al intervalo de proyección (path invariance), de modo que el resultado de proyectar primero de t0 a t1 y luego de t1 a t t2 es el mismo que el obtenido en un solo paso de t0 a t2. Además, en el caso de que la variable de estado sea la altura dominante, una única ecuación dinámica define conjuntamente la curva de crecimiento en esta variable y el modelo de índice de sitio. Bailey y Clutter (1974) fueron los primeros en formalizar la propiedad de invarianza con respecto a la edad de referencia, aunque en un sentido más amplio, puesto que según dichos autores esta propiedad implica también que la estimación de los parámetros de la ecuación dinámica se realice empleando una metodología de ajuste que sea a su vez invariante en la edad de referencia, es decir, que las estimaciones de los parámetros no dependan de un punto arbitrariamente seleccionado como referencia y que permita el empleo de todos los datos disponibles. Estos autores propusieron una metodología para obtener ecuaciones dinámicas con estas características que se conoce en la literatura forestal como Algebraic Difference Approach (ADA). Dicha metodología esencialmente se basa en la premisa de que es posible modelizar el crecimiento de una determinada variable con una familia de curvas originadas por un modelo de crecimiento base con todos los parámetros comunes excepto uno, que es un parámetro específico de la estación forestal. Así, se podrían obtener, por tanto, para cada modelo de crecimiento base, tantas ecuaciones dinámicas como parámetros tenga dicho modelo. Cieszewski y Bailey (2000) propusieron una generalización de la metodología ADA, denominada Generalized Algebraic Difference Approach (GADA) cuya principal ventaja radica en que se pueden obtener ecuaciones dinámicas a partir de un modelo de crecimiento base permitiendo que más de un parámetro varíe con las condiciones específicas de la estación forestal. Para muchos modelos base, esta es la única posibilidad de generar familias de curvas polimórficas con múltiples asíntotas, una propiedad importante de los sistemas de calidad de estación (Cieszewski, 2002). El primer paso de la metodología GADA es seleccionar una ecuación base e identificar en ella los parámetros que se desea que sean específicos de la estación. Posteriormente, debe definirse explícitamente cómo cambian dichos parámetros entre las diferentes estaciones reemplazándolos con funciones explícitas de X (una variable independiente oculta que describe la productividad de la estación como resultado de las prácticas de gestión, las condiciones del suelo y los factores ecológicos y climáticos) y nuevos parámetros. De este modo, la ecuación base bidimensional inicialmente seleccionada (Y = f (t)) se expande en una ecuación tridimensional Y = f (t, X) que describe tanto cambios transversales como longitudinales con dos variables independientes t y X. Debido a que X Juan Daniel García Villabrille 48 no se puede medir de una forma fiable e incluso no se puede definir funcionalmente, el último paso de la metodología GADA consiste en sustituir X por condiciones iniciales equivalentes que representan observaciones puntuales de las dos variables observables t e Y (p. ej., (t1, Y1)), de manera que el modelo se pueda definir implícitamente (p. ej., Y2 = f (t2, t1, Y1)) y sea aplicable (Cieszewski y Bailey, 2000; Cieszewski, 2002). Durante este proceso a menudo se eliminan parámetros redundantes, obteniéndose un modelo que tiene un número de parámetros menor o igual que la ecuación base original. 5.2.3. Ajuste de los modelos Las variables de estado que se modelizan con ecuaciones dinámicas (altura dominante, número de pies por hectárea o área basimétrica) son medidas en campo y están sujetas a errores de medición. Al ajustar la ecuación dinámica se asume que la variable Y no tiene error cuando está a la derecha de la igualdad pero sí tiene error cuando está en la parte izquierda de la igualdad, lo que plantea un problema puesto que, en esta situación, las técnicas de estimación convencionales aplicadas a modelos lineales dan lugar a estimaciones de los parámetros sesgadas (Myers, 1986, p. 358). Por extensión, el problema es similar en modelos no lineales como las ecuaciones dinámicas, por lo que los valores de la variable de estado que aparecen en la parte derecha de la igualdad deben ser valores predichos que no pueden ser conocidos hasta que los parámetros hayan sido estimados (Krumland y Eng, 2004). Aunque existen diferentes metodologías para resolver este problema (p. ej., Gregoire et al., 1995; Cieszewski et al., 2000), uno de los métodos más sencillos de aplicar es el propuesto por Cieszewski et al. (2000) basado en el empleo de variables dummy, que en realidad es una simplificación de los modelos de efectos mixtos propuestos por Lindstrom y Bates (1990). Información detallada de esta metodología y de sus posibles aplicaciones se puede encontrar en Cieszewski (2003). En la metodología de variables dummy los valores iniciales de la ecuación dinámica son los mismos para todos los datos de una unidad de muestreo determinada (árbol o parcela). La edad inicial puede elegirse arbitrariamente para cada individuo, dentro de unos límites (p. ej., no puede ser cero). El valor de la variable de estado correspondiente a la edad inicial para cada unidad de muestro se estima simultáneamente con el resto de los parámetros globales de la ecuación dinámica. A diferencia de las metodologías de estimación tradicionales que requieren una elección de la edad base previa al ajuste para forzar a la ecuación a pasar por el punto variable de estado/edad, el método de variables dummy tiene en cuenta que cada medición tiene error y, por tanto, no fuerza a que la ecuación pase por ningún punto en concreto. En lugar de eso, la curva se ajusta a las tendencias individuales de los datos originales. Para aplicar este método se necesita contar con, al menos, dos mediciones para cada unidad de muestreo (árbol o parcela), y Funciones de transición 55 5.2.5. Comparación y selección de modelos Según Weiskittel et al. (2011, p. 279), una de las consideraciones importantes en el desarrollo de un modelo son los criterios empleados para su selección. Así, la comparación de modelos se basó en los estadísticos de ajuste R2, indicativo de la variabilidad de la muestra explicada por el modelo, REMC y E, indicativos del error en las unidades de la variable dependiente, y el AIC, indicativo de la mejor forma para el modelo a partir de los datos de partida, aunque no asegura una extrapolación adecuada (Weiskittel et al., 2011, p. 281). Las expresiones de estos cuatro estadísticos quedaron expuestas en el apartado 3.5.1. También se realizaron análisis gráficos de los residuos (ei) frente a los valores predichos ( i y ˆ) y de éstos frente a los valores observados ( i y), así como del error relativo (ER) por clases de edad (t) (ver apartado 3.5). 5.3. Resultados 5.3.1. Altura dominante En la Tabla 5.4 se muestran los resultados del ajuste para los modelos de altura dominante (M1 a M9). Todos los parámetros en todos los modelos fueron significativamente distintos de cero a un nivel de confianza del 95% (α = 0,05), aunque fue necesario fijar en 0 el parámetro b1 de los modelos de M3 y M6, así como el b2 del modelo de M9, que se corresponden con las formulaciones GADA de los modelos de Korf, Hossfeld y Bertalanffy-Richards respectivamente. Los estadísticos de ajuste fueron similares en todos los modelos, con valores de R2 en torno a 0,992 sin sesgo. El valor de REMC varió entre 0,673 y 0,733 m. Respecto al valor del AIC, los menores valores correspondieron a modelos derivados de la función de Hossfeld, en concreto a los modelos M4 (Hossfeld ADA, resuelto por a1, con AIC = -141) y M6 (Hossfeld GADA, con b1 = 0 y AIC = -143). Tanto para las funciones de Korf como de Hossfeld, las formas GADA presentaron menores valores del AIC que las correspondientes formulaciones ADA, y entre éstas fueron mejores los modelos en los que el parámetro específico de la estación (X) fue a1. En el caso de la función de Bertalanffy-Richards el menor AIC correspondió a la forma ADA resuelta también por a1 (M7). Si bien los valores de R2, REMC y E son similares entre sí para todos los modelos, a partir de los valores del AIC se pueden hacer diferencias, teniendo en cuenta que una reducción en 10 unidades de este estadístico indica diferencias significativas (Burnham y Anderson, 2002). No obstante, dada la similitud de los estadísticos mencionados, la Juan Daniel García Villabrille 56 selección del mejor modelo se basó en consideraciones gráficas y del comportamiento deseado para el modelo. Tabla 5.4. Parámetros estimados y estadísticos obtenidos en el ajuste de modelos de altura dominante. Modelo Parámetros Ajuste independiente b1 b2 b3 R2 REMC E AIC Korf M1 3,559 (0,000) 0,5811 (0,000) 0,9930 0,6807 0,001 -136,8 M2 59,67 (0,000) 0,5315 (0,000) 0,9927 0,6947 0,000 -125,5 M3 0* 14,11 (0,000) 0,6037 (0,000) 0,9931 0,6784 0,001 -138,8 Hossfeld M4 - 14,91 (0,000) 1,377 (0,000) 0,9932 0,6747 -0,000 -141,2 M5 41,05 (0,000) 1,187 (0,000) 0,9926 0,7016 -0,001 -120,0 M6 0* 514,8 (0,000) 1,39 (0,000) 0,9932 0,6728 0,000 -143,4 Bertalanffy- Richards M7 - 0,1523 (0,000) 1,277 (0,000) 0,9931 0,6793 -0,000 -138,1 M8 41,04 (0,000) 0,8993 (0,000) 0,9919 0,7325 -0,003 -95,94 M9 0,1552 (0,000) 0* 4,357 (0,000) 0,9930 0,6807 -0,000 -136,9 *Indica parámetro fijado en dicho valor (0); el valor entre paréntesis indica el nivel de significación. Tras un análisis gráfico previo de las proyecciones de todos los modelos, se presenta exclusivamente el correspondiente a los modelos M3, M4 y M6 (Figura 5.1). Todos ellos representaron adecuadamente los datos de partida, siendo muy similares hasta los 10 años. En edades avanzadas (> 15 años) y calidades bajas, el modelo M4 predice menos altura que los otros dos (M3 y M6), mientras en calidades altas es el modelo M6 el que predice los menores valores. En general, el modelo M3 es el que mejor representa las tendencias observadas. Respecto a la selección de una edad de referencia (tref), y aunque Álvarez-González et al. (2004) y otros autores (Diéguez-Aranda et al., 2006a; Barrio-Anta et al., 2008) recomiendan emplear como criterio la edad a la que se obtiene el menor error relativo (ER), en este caso se ha seleccionado la edad de referencia de 7 años por los mismos motivos que lo hicieron García y Ruiz (2003), es decir, a petición de la empresa ENCE (colaborador principal en el aporte de datos). El gráfico de error relativo (Figura 5.2) muestra que todos los modelos se comportaron de forma similar a partir de los cuatro años, y el ER descendió paulatinamente a pesar de la disminución en el número de datos en edades elevadas, llegando a su mínimo a los 21 años. Sin embargo, teniendo en cuenta que no existen datos suficiente a partir de los 15 años, que el turno de la especie puede fijarse a partir de los 13 años (Álvarez-González et al., 2005a) como posteriores simulaciones con el modelo desarrollado en esta tesis demostraron, y que la edad de referencia debe situarse antes del turno (Goelz y Burk, 1992), los valores del ER (%) Funciones de transición 57 entre los 6 y los 9 años son muy similares, por lo que no se rechaza la elección arbitraria realizada para tref = 7 años con ninguno de los modelos. Figura 5.1. Comparación de las proyecciones de los modelos M3, M4 y M6 sobre los datos observados de altura dominante frente a la edad. Las proyecciones corresponden a índices de sitio de 7, 13, 19 y 25 m a la edad de referencia de 7 años. Figura 5.2. Error relativo para los modelos de altura dominante M3, M4 y M6 y número de datos observados (n) frente a la edad. Finalmente, tras seleccionar el modelo M3 (Korf GADA) por su mejor reflejo de las tendencias en las clases extremas, éste se ajustó empleando variables dummy (o categóricas) para tratar de identificar diferencias según los materiales de reproducción y el tipo de suelo (granitos frente a esquistos/pizarras). Los parámetros que multiplican a las variables dummy de estos ajustes no resultaron significativos, por lo que se descartó el Edad (años) Altura dominante (m) 10 20 30 40 5 10 15 20 M3 M4 M6 Observado 4 6 8 10 12 14 0 1 2 3 4 5 6 Edad (años) ER(%) 0 5 10 15 20 25 30 Nº datos M3 M4 M6 n Juan Daniel García Villabrille 58 empleo de un modelo con parámetros específicos del tipo de planta empleado (brinzales frente a clones, y éstos a su vez agrupados y en categorías diferentes) o del tipo de suelo (granitos frente a esquistos). El último análisis correspondió a las proyecciones del modelo seleccionado (M3) sobre los datos de partida por materiales de reproducción (Figura 5.3). Las predicciones se adaptan a los datos y, como se puede comprobar, los escasos datos pertenecientes a clones F1 (Odiel y Tinto) se sitúan en calidades media-altas (IS = 19 aprox.), y los máximos son similares a los del clon F0 (Anselmo), no habiendo sido posible obtener datos de masas clonales en la máxima calidad simulada (IS = 25 m), la cual quedó establecida a partir de los valores observados para los brinzales al principio y final del intervalo de proyección. Figura 5.3. Curvas de calidad seleccionadas (modelo M3) superpuestas sobre los datos observados separados por materiales de reproducción frente a la edad para índices de sitio de 7, 13, 19 y 25 m a la edad de referencia de 7 años. 5.3.2. Área basimétrica En la Tabla 5.5 se muestran los resultados del ajuste para los modelos de área basimétrica. Todos los parámetros fueron significativamente diferentes de cero a un nivel de confianza del 95%, aunque al igual que para los ajustes de altura dominante, fue necesario fijar b1 a 0 en los modelos de M3, M6 y el parámetro b2 a 0 en el modelo M9, correspondientes a los modelos base de Korf, Hossfeld y Bertalanffy-Richards en sus Edad (años) Altura dominante (m) 10 20 30 40 5 10 15 20 Brinzal Clon Anselmo Clon Odiel 5 10 15 20 10 20 30 40 Clon Tinto Observado Predicho Funciones de transición 59 formulaciones GADA respectivamente. Los estadísticos de todos los modelos fueron similares, con valores de R2 en torno a 0,994 y E despreciable. En los valores de REMC se identifican ligeras diferencias, siendo el mejor modelo el M3 (Korf GADA, con b1 = 0) con REMC = 0,71 m2/ha. Los valores del AIC obtenidos por los modelos derivados de la función de Korf fueron mejores que los de las otras dos funciones (Hossfeld y Bertalanffy-Richards), siendo el mejor el modelo M3 (AIC = -115). Además, al igual que para la función de altura dominante, los mejores valores de AIC de las formulaciones GADA parecen indicar que las curvas de área basimétrica deberían ser polimórficas con asíntotas diferentes, aunque el segundo mejor valor del AIC (-106) correspondió a la función de Korf resuelta por a2 que genera curvas polimórficas pero con una asíntota común. En general, las formas GADA fueron mejores que las formas ADA, de forma similar a lo encontrado por Barrio et al. (2006, 2008) y Castedo-Dorado et al. (2007). Tabla 5.5. Parámetros estimados y estadísticos obtenidos en el ajuste de modelos de área basimétrica. Modelo Parámetros Ajuste independiente b1 b2 b3 R2 REMC E AIC Korf M1 6,513 (<0,001) 1,077 (<0,001) 0,9948 0,7261 -0,003 -100,8 M2 44,76 (<0,001) 0,7705 (<0,001) 0,9950 0,7193 -0,000 -106,1 M3 0* 22,74 (<0,001) 1,074 (<0,001) 0,9951 0,7071 -0,001 -115,6 Hossfeld M4 - 39,24 (<0,001) 1,955 (<0,001) 0,9945 0,7505 -0,008 -82,41 M5 37,51 (<0,001) 1,48 (<0,001) 0,9942 0,7691 -0,009 -68,76 M6 0* 514,8 (<0,001) 1,39 (<0,001) 0,9949 0,7262 -0,005 -100,8 Bertalanffy- Richards M7 - 0,2093 (<0,001) 2,130 (<0,001) 0,9943 0,7627 -0,008 -73,42 M8 39,86 (<0,001) 1,218 (<0,001) 0,9935 0,8131 -0,011 -37,71 M9 0,2066 (<0,001) 0* 6,802 (<0,001) 0,9946 0,7421 -0,005 -88,72 *Indica parámetro fijado en dicho valor (0); el valor entre paréntesis indica el nivel de significación. Tras un análisis gráfico previo de las proyecciones de todos los modelos, se decidió comparar finalmente las proyecciones de los modelos M3, M4 y M9 (Figura 5.4) que fueron los que mejor representaron los datos de partida entre los 5 y los 14 años, especialmente en las calidades intermedias. El modelo M2 (curvas polimórficas de asíntota común) y el M6 (curvas polimórficas con asíntotas diferentes) quedaron descartados por no adaptarse a las trayectorias observadas a pesar de los meojres valores del AIC que obtuvieron en comparación con los modelos M4 y M9. Juan Daniel García Villabrille 60 Figura 5.4. Comparación de las proyecciones de los modelos M3, M4 y M9 sobre los datos observados de área basimétrica. Las proyecciones corresponden a valores de área basimétrica de 4, 11, 18 y 25 m2/ha a la edad de 7 años. Por otra parte, algunas trayectorias observadas sufren modificaciones del ritmo normal de crecimiento provocadas por unas determinadas condiciones de clima, suelo y topografía (Fernández-López, 1982, en referencia a calidad de estación), aunque también es conocido el efecto positivo de la eliminación de la competencia y la fertilización (Ruiz et al., 2008). En el caso de los crecimientos negativos observados en edades avanzadas, y tras descartar errores de medición, la mortalidad es el único factor responsable. También la densidad de plantación es un factor que puede alterar la percepción inicial de la verdadera calidad. Así, de nuevo las parcelas más jóvenes (< 5 años) con buenos crecimientos en grosor pero no en las mayores densidades, no quedaron situadas en las curvas superiores, siendo deseable un modelo que no haga muchas diferencias inicialmente. En la Figura 5.5 se muestra el gráfico del error relativo para los modelos M3, M4 y M9. Los tres modelos produjeron valores similares partir de los 4 años, siendo ligeramente mejores los del modelo M3. Edad (años) Área basimétrica (m2/ha) 10 20 30 40 5 10 15 20 M3 M4 M9 Observado Funciones de transición 61 Figura 5.5. Error relativo para los modelos de área basimétrica M3, M4 y M9 y número de datos observados (n) frente a la edad. Finalmente se seleccionó el modelo M3 (Korf GADA) para la predicción de área basimétrica, y al igual que para la función de atura dominante se realizó un ajuste con variables dummy para tratar de identificar diferencias en función del material de reproducción y del tipo de suelo de la plantación. Los parámetros para identificar diferencias entre materiales de reproducción no fueron significativos, mientras que para el tipo de suelo sí resultaron significativos (Tabla 5.6), cuyo REMC disminuyó 0,02 m2/ha frente al modelo reducido (o único). Por otra parte, sería deseable tener mejor representación en cada suelo de todas las calidades y edades, ya que las plantaciones en suelos graníticos se concentran en clases de área basimétrica intermedias y sólo en edades de 5 a 15 años, lo que no justifica los cambios de forma de las curvas, sobre todo en edades jóvenes, donde mayores diferencias hay y menos datos se tienen (Figura 5.6). De igual forma, Wang y Baker (2007) rechazaron un modelo con parámetros específicos según la región para el caso de la altura dominante por falta de datos en todas las calidades, a pesar de obtener parámetros significativos. Tampoco la mejora obtenida sobre el REMC (2,8%) se considera suficiente para proponer el empleo de un modelo completo (pizarras y granitos por separado). 4 6 8 10 12 14 0 5 10 15 20 Edad (años) ER(%) 0 5 10 15 20 25 30 Nº datos M3 M4 M9 n Juan Daniel García Villabrille 62 Tabla 5.6. Parámetros estimados en los ajustes del modelo M3 (Korf GADA, con b1 = 0) de área basimétrica según el tipo de suelo, con b2 = b20 + b21·I y b3 = b30 + b31·I, donde I = 0 para esquistos y 1 para granitos. Modelo Parámetros Ajuste b20 b21 b30 b31 REMC M3 27,08 (<0,001) -12,13 (<0,001) 1,115 (<0,001) -0,5092 (0,032) 0,68 El valor entre paréntesis indica el nivel de significación. Figura 5.6. Comparación de curvas de área basimétrica generadas por los modelos reducido (único) y completo (pizarras y granitos por separado) para valores de área basimétrica de 4, 11, 18 y 25 m2/ha a los 7 años. Respecto a las variables climáticas recopiladas (precipitación total en un año, en l/m2; días de precipitación > 1 mm; días de precipitación > 10 mm; días de helada), todas fueron incluidas en el análisis de variables independientes, con el fin de establecer la existencia o no de una relación entre el crecimiento en área basimétrica y la precipitación o las temperaturas extremadamente bajas, intentando dotar de mayor sentido biológico a las predicciones. En la Figura 5.7 se muestra el incremento anual en área basimétrica frente a las variables climáticas consideradas. En todos los casos se observó una falta de relación clara entre los crecimientos en área basimétrica y dichas variables, y posteriores ajustes, lineales y no lineales, confirmaron nula o escasa capacidad explicativa de dichas variables. Finalmente se decidió no incluir tampoco dichas variables en los ajustes de los modelos. El hecho de que finalmente no fuese útil ninguna de las variables climáticas recopiladas podría ser indicativo de que no se producen situaciones limitantes para la Edad (años) Área basimétrica (m2/ha) 10 20 30 40 5 10 15 20 Modelo reducido Modelo completo-pizarras Modelo completo-granitos Funciones de transición 63 especie respecto a la precipitación anual o por heladas en las zonas de Galicia donde se muestrearon las parcelas, que teóricamente es la idónea para este eucalipto. Wang y Baker (2007) señalan que las precipitaciones en el sureste de Australia se sitúan entre 600-1000 mm/año, produciendo situaciones limitantes para el crecimiento. Esto facilitaría la aparición de diferencias por regiones, lo que puede haber sido un factor importante en el caso del modelo de área basimétrica elaborado por dichos autores, así como para el de Tomé et al. (2006) desarrollado para Portugal, donde la precipitación media anual para todas las regiones es de 810 mm/año (máx. = 2623 mm/año, mín. = 372 mm/año) como se extracta de Tomé et al. (2001). Figura 5.7. Crecimiento en área basimétrica anual (m2/ha) frente a diversas variables climáticas. Como se comentó, son muchos los factores que dificultan el ajuste de una función de área basimétrica que mejore los errores obtenidos, como el desconocimiento de los tratamientos aplicados (p. ej., desbroces y fertilización fundamentalmente) o la cantidad de datos según edades y calidades, siendo también recomendable aumentar el período entre mediciones para esta especie (Gadow et al., 1999) con el fin de obtener trayectorias a medio-largo plazo. Por último, las proyecciones del modelo M3 seleccionado (Korf GADA, con b1 = 0) sobre los datos de partida separados por materiales de reproducción (Figura 5.8) no muestran evidencias de un comportamiento erróneo de las curvas en ningún caso. 1000 1500 2000 2500 3000 3500 0 1 2 3 4 Precipitación total 1 año (L/ m 2 ) incremento de G ( m2/ ha) 100 150 200 0 1 2 3 4 Días precipitación > 1 mm 20 40 60 80 100 0 1 2 3 4 Días precipitación > 10 mm incremento de G ( m2/ ha) 0 10 20 30 40 50 60 0 1 2 3 4 Días de helada Juan Daniel García Villabrille 64 Figura 5.8. Curvas de área basimétrica generadas por el modelo M3 superpuestas sobre los datos observados separados por materiales de reproducción frente a la edad para valores de área basimétrica de 4, 11, 18 y 25 m2/ha a la edad de 7 años. 5.3.2.1. Función de inicialización de área basimétrica Tras ajustes previos que incluyeron un mayor número de modelos, los mejores resultados correspondieron a seis modelos en los cuales todos los parámetros fueron significativos al nivel del 5%. El modelo M1 correspondió a un modelo de la forma 321 00 a aa tNHaG ⋅⋅⋅= , el modelo M2 de la forma NaHaG a200 1−⋅= y los modelos M3 a M6 fueron modelos derivados de la función seleccionada para la proyección de área basimétrica. A su vez, en los modelos M3 y M4 el parámetro específico de la estación se sustituyó por un modelo alométrico ( 1 0a ISaX ⋅= ), siendo el M4 compatible con la función de proyección, en el modelo M5 se empleó un modelo de variable combinada ( NaISaX a20 1−⋅= ) y en el modelo M6 se empleó una expresión similar al M5 pero sustituyendo la densidad por la densidad de plantación ( pl aNaISaX 20 1−⋅= ) ya que este dato no cambia a lo largo del turno. Los mejores resultados correspondieron al modelo M2, con R2 = 0,874, REMC = 2,9 m2/ha y E = 0, si bien presenta el inconveniente de poder predecir valores negativos, por lo que se descartó inicialmente. Del resto de modelos, todos indicaron ligera sobrestimación, pero los que menor sesgo (E = -0,08) produjeron fueron el M3 y M5 con REMC = 3,1 m2/ha y REMC = 2,9 m2/ha respectivamente, ambos derivados de la función de Korf, si bien no aseguran la compatibilidad con el modelo ajustado para la Edad (años) Área basimétrica (m2/ha) 10 20 30 40 5 10 15 20 Brinzal Clon Anselmo Clon Odiel 5 10 15 20 10 20 30 40 Clon Tinto Observado Predicho Funciones de transición 71 Tabla 5.8. Resultados del test de correlación de Pearson entre residuos del ajuste independiente. Residuos ρxy p-valor H0 frente a G 0,538 <0,001 H0 frente a N 0,056 0,346 G frente a N 0,219 <0,001 Si bien no se encontró correlación ente los residuos del modelo de altura dominante y los del modelo de mortalidad (p-valor = 0,346), en los otros dos casos sí se detectó correlación (p-valor <0,001), por lo que se procedió al ajuste simultáneo de los tres componentes del modelo dinámico, cuyos estadísticos se muestran en la Tabla 5.9. Todos los parámetros de los tres modelos fueron significativos al nivel de significación del 5%. Tabla 5.9. Resultados del ajuste simultáneo (NSUR) para las funciones de transición de altura dominante, área basimétrica y mortalidad. Función Parámetros Ajuste (n = 279) R2 REMC E Altura dominante a2 =13,90 (<0,001) a3 = 0,5989 (<0,001) 0,9931 0,68 -0,000 Área basimétrica b2 =21,16 (<0,001) b3 = 0,9906 (<0,001) 0,9951 0,71 -0,006 Mortalidad c0 = 0,1995 (<0,001) 0,9897 31 -0,091 El valor entre paréntesis indica el nivel de significación. Los errores obtenidos son similares a los de Tomé et al. (2001) para Eucalyptus globulus en Portugal, con muestras de tamaño n > 11 000 para altura dominante (REMC = 0,83 m), área basimétrica (REMC = 0,69 m2/ha) y mortalidad (REMC = 28 pies/ha). Con respecto al modelo de Wang y Baker (2007), también para E. globulus, con muestras n > 2000 los valores de REMC obtenidos fueron de 1,76 m, 1,63 m2/ha y 24 pies/ha para altura dominante, área basimétrica y mortalidad respectivamente, siendo ligeramente inferiores. El análisis gráfico de los residuos de los modelos para las tres funciones de transición (altura dominante, área basimétrica y mortalidad) no mostró evidencias claras de heterocedasticidad, por lo que no fue necesario realizar el ajuste empleando regresión ponderada (Figura 5.14). Juan Daniel García Villabrille 72 Figura 5.14. Gráfico de residuos frente a valores predichos para el modelo de altura dominante (m) (arriba), área basimétrica (m2/ha) (centro) y mortalidad (pies/ha) (abajo) en el ajuste simultáneo. 5.4. Resumen de funciones de transición del modelo dinámico En la Tabla 5.10 se muestra un resumen de las funciones que componen el modelo dinámico de crecimiento propuesto para rodales regulares de Eucalyptus globulus en Galicia, cuyo destino principal es la celulosa. En caso de que no se conozca el valor del área basimétrica basta con sustituir G1 en la función de transición de dicha variable por el valor predicho por la ecuación de inicialización Gt, que depende del índice de sitio y la edad (IS, t). Para calcular IS basta con sustituir la edad t2 en la función de transición de la altura dominante por 7 años. 10 20 30 40 -2 0 2 Predichos Residuos 0 10 20 30 40 -3 -2 -1 0 1 2 3 Predichos Residuos 800 1000 1200 1400 1600 -2 0 2 Predichos Residuos Funciones de transición 73 Tabla 5.10. Resumen de las funciones de transición que componen el modelo dinámico de crecimiento para rodales regulares de Eucalyptus globulus en Galicia. Función de transición para la altura dominante (H0) ( )        ⋅         −⋅= −5989,0 2 1 102 90,13 exp exp t X XH ( ) ( ) ( )      ⋅−+⋅⋅+⋅⋅= 2 01 5989,0 1 5989,0 101 5989,0 1 5989,0 11 90,134 2 1HLnttHLnttX Función de transición para el área basimétrica (G) ( )        ⋅         −⋅= −9906,0 2 1 12 16,21 exp exp t X XG ( ) ( ) ( )      ⋅−+⋅⋅+⋅⋅= 2 1 9906,0 1 9906,0 11 9906,0 1 9906,0 11 16,214 2 1GLnttGLnttX Función de inicialización para el área basimétrica (G) ( )        ⋅         −⋅= −8689,0 1 1 23,78 exp exp t X XG t 2564,0 1819,1 ISX ⋅= Función de transición de mortalidad (N) 2 2 1 2 2 5,0 12 100100 0,1995 − −                      −       ⋅+= tt NN En la Figura 5.15 se muestran las curvas de evolución resultantes de aplicar las funciones de transición seleccionadas sobre las trayectorias observadas en cada parcela. Juan Daniel García Villabrille 74 Figura 5.15. Curvas de evolución con la edad de las variables de estado altura dominante (H0, m), área basimétrica (G, m2/ha) y número de pies por hectárea (N, pies/ha) generadas con las funciones de transición que constituyen el modelo dinámico de crecimiento para Eucalyptus globulus en Galicia, superpuestas a las trayectorias observas. Arriba: curvas para índices de sitio de 7, 13, 19 y 25 m a la edad de referencia de 7 años. Centro: curvas para valores de área basimétrica de 4, 11, 18 y 25 m2/ha a la edad de 7 años. Abajo: curvas de mortalidad para densidades de plantación de 800, 1200, 1600 y 2000 pies/ha. Edad (años) 500 1000 1500 2000 5 10 15 20 Número de pies (pies/ha) 500 1000 1500 2000 N 0 10 20 30 40 Área basimétrica (m2/ha) 0 10 20 30 40 G 0 10 20 30 40 Altura dominante (m) 0 10 20 30 40 H0 Observado Predicho Funciones de salida 75 6. Funciones de salida En este capítulo se aborda el ajuste de las funciones de salida que completan el modelo de crecimiento, es decir, aquellas que predicen los valores de la producción, actuales o futuros, de volumen y biomasa de árbol individual y de rodal. También se incluye una función de desagregación de rodal que permite conocer el reparto de la producción por clases diamétricas. 6.1. Introducción 6.1.1. Tarifas de cubicación La estimación de las existencias en volumen de las masas forestales es uno de los objetivos más importantes que se le presentan al gestor forestal. Tradicionalmente este problema se ha solventado estimando el volumen de los árboles individuales mediante tarifas o tablas de cubicación de árbol individual, que son funciones que dependen de variables sencillas de medir (generalmente el diámetro normal y la altura total). Otra de las posibles soluciones para estimar el volumen de los árboles es mediante el empleo de funciones de perfil, justificadas en el caso de que existan diferentes destinos comerciales en función del grosor de la madera. En el caso del eucalipto en Galicia, cuyas masas están destinadas principalmente a la producción de celulosa, y por tanto a la trituración, no existe una primera necesidad de desarrollar funciones de perfil, y además ya existe un modelo (Diéguez-Aranda et al, 2009, p. 208). Este tipo de ecuaciones se consideran importantes para árboles o masas cuyo destino sea la madera estructural o piezas especiales como las bateas. Además de las tarifas de cubicación de árbol individual, también se consideran importantes las tarifas de rodal, que sirven para estimar el volumen de una masa, generalmente a partir de las variables de estado H0, G y/o N, y que son consideradas funciones de salida del modelo (García, 1988). En la literatura forestal han surgido diferentes aproximaciones para determinar las producciones futuras en volumen de una masa (Clutter et al., 1983, p. 94), siendo recomendable estimar las variables necesarias para su cálculo en un momento y aplicar la ecuación de cubicación en dicho momento. Además, en esta tesis se proveen todas las herramientas para aplicar esta vía. Aunque E. globulus es una de las especies más importantes en Galicia, y de que el producto para la obtención de la celulosa es la madera del árbol exclusivamente, aún no Juan Daniel García Villabrille 76 existen herramientas que estimen el volumen de este producto. Por otra parte, es de interés para la industria la estimación del porcentaje de corteza en volumen (y en peso) ya que el aprovechamiento no siempre se realiza de la misma forma, eliminando o no la corteza, por lo que se debe tener cuenta que la logística se verá afectada. Así, se han desarrollado tarifas de cubicación de árbol individual y de rodal, con y sin corteza para madera procedente de plantaciones de primer turno destinadas a celulosa. En una primera aproximación, las tarifas de cubicación se pueden dividir en tarifas simples o tarifas compuestas, según sean adecuadas respectivamente para una única especie o para más de una (Diéguez-Aranda et al., 2003, p. 196). En este trabajo se ajustan tarifas simples al provenir los datos de una única especie. Además, se testará la necesidad de ajustar modelos diferentes según los diferentes materiales de reproducción y también en función del tipo de suelo. Las tarifas de cubicación de árbol individual son ecuaciones que relacionan el volumen de un árbol con variables dendrométricas, o variables de entrada, generalmente sencillas de medir. Dependiendo del número de variables que intervengan en la ecuación se pueden clasificar en tarifas de una, dos o tres entradas, según el volumen esté en función de una, dos o tres variables de árbol, respectivamente. El número de variables que incluya la tarifa debe ser reducido y, además, dichas variables deben ser sencillas de medir en campo, con el fin de facilitar su uso, y estar fuertemente correlacionadas con el volumen y débilmente correlacionadas entre sí. En general, la precisión en la estimación del volumen obtenida con una tarifa de cubicación aumenta a medida que aumenta el número de entradas (variables), aunque el exceso de variables explicativas puede suponer un esfuerzo de muestreo mayor (López et al., 2003) y añadir colinealidad entre las variables (Peña, 2002, p. 557). 6.1.2. Ecuaciones de biomasa Se conoce como biomasa a la cantidad de materia que existe en un determinado ecosistema por unidad de superficie o de volumen (Newbould, 1970; Pardé, 1980). La mayoría de los estudios realizados hasta la fecha tienen como objetivo conocer cómo influye la acumulación y distribución de la biomasa en el equilibrio del ecosistema (flujos de energía y de nutrientes), evaluar las posibilidades energéticas e industriales en determinadas regiones o determinar el potencial de almacenamiento del carbono atmosférico en los sistemas forestales a medio y largo plazo, en relación a la contribución de las masas arboladas al ciclo global del carbono (Parresol, 1999; van Laar y Akça, 1997, p. 183). Por otro lado, diferentes intereses industriales, energéticos y ambientales han dirigido su atención a la biomasa forestal, especialmente las industrias de trituración y de pasta de papel, dada su importancia en Galicia, que demandan cada vez en mayor medida biomasa forestal no maderable (Bermúdez y Piñeiro, 2001). Estudios previos revelan que Galicia, Funciones de salida 77 amparándose en su elevado rendimiento forestal y en su desarrollado sector de transformación de la madera, podría ser capaz de asimilar buena parte de estos recursos arbóreos no maderables para el aprovechamiento industrial o energético (Bermúdez y Touza, 2000; Bermúdez y Piñeiro, 2001). Entre las estrategias propuestas se encuentra el aprovechamiento de fracciones arbóreas no maderables tras el aprovechamiento forestal, o el establecimiento de plantaciones y programas de manejo especialmente diseñados para la producción de biomasa (Álvarez González et al., 2005a). Otro aspecto que se debe tener en cuenta es que la biomasa representa el peso seco, por lo que no es un dato relevante a la hora de realizar una compra-venta de madera. Su utilidad radica, sobre todo, en conocer las proporciones o cantidades de cada uno de los elementos del árbol y así poder cuantificar la materia seca de madera pulpable o de las fracciones no pulpables, los nutrientes y el carbono acumulado en una masa forestal. Este carbono, denominado Bono o Crédito de Carbono, es uno de los sistemas establecidos en el Protocolo de Kioto (1997) como medio de reducción de los gases de efecto invernadero, existiendo un mercado internacional. La estimación de biomasa, al igual que el volumen, se ha solventado mediante ecuaciones que dependen de variables sencillas de medir. El principal problema de la estimación de la biomasa es que se trata de un sistema compuesto por diferentes partes o fracciones, y por tanto debe coincidir la estimación hecha para todo el árbol con la suma de la biomasa por fracciones. Así, Zellner (1962) desarrolló una solución aplicable a este tipo de sistemas, denominado Seemingly Unrelated Regression (SUR) (ver Anexo II, apartado II.3), que asegura la aditividad en este tipo de sistemas. Por otra parte, el desarrollo de ecuaciones de biomasa requiere de una serie de operaciones sobre los datos que pueden ser empleadas, a su vez, para conocer otras características de la madera tales como la densidad aparente (en peso verde), la densidad básica (en peso seco) o la humedad. De todas estas características, se consideran de utilidad directa en la compra-venta tanto la densidad aparente en verde como las variaciones de la humedad de la madera, ya que permiten estimar el peso de las trozas, siendo ésta la forma habitual de enajenación. El hecho de que se reduzca el número de variables en comparación con las tarifas de cubicación se debe a que los sistemas aditivos generalmente presentan limitaciones de convergencia en el ajuste, y más cuando el número de ecuaciones es elevado como en el caso de la biomasa, para la que se pueden definir al menos cuatro fracciones diferentes: madera, corteza, ramas y hojas. A su vez, las ramas se suelen subdividir en tres fracciones: ramas gruesas, ramas finas y ramillos, lo que puede aumentar hasta seis el número de componentes finales de los modelos, existiendo otras posibilidades como la separación entre ramas vivas y muertas, por ejemplo. Aunque ya existen modelos de biomasa de árbol individual en Galicia para E. globulus, como los de Brañas et al. (2000a) o Álvarez-González et al. (2005a), dichos sistemas presentan carencias respecto al muestreo considerando las mejoras aportadas Juan Daniel García Villabrille 78 por Pérez-Cruzado y Rodríguez-Soalleiro (2011), y tampoco existen tarifas de rodal. En esta tesis, además del análisis de diferentes características (o variables) para los diferentes materiales de reproducción y de las mejores ecuaciones de forma independiente, se han desarrollado sistemas de ecuaciones aditivas. A nivel de árbol individual se ha hecho una comparación entre un sistema “completo”, que considera todas las fracciones posibles según tamaños, y uno más simple denominado “reducido”, que sólo diferencia tres fracciones (ramas, corteza y madera), añadiendo posteriormente la biomasa de hojas mediante una ecuación independiente. También se ha ajustado un sistema completo para la estimación de la biomasa a nivel de rodal. Al igual que el modelo dinámico y las tarifas de cubicación, los modelos están indicados para madera procedente de masas de eucalipto destinadas a celulosa. Por otra parte, para considerar los evidentes daños por goniptero (Gonipterus scutellatus, Gyll.) presentes en las masas de eucalipto en Galicia, se ha evaluado por primera vez el comportamiento de una variable de estimación subjetiva, denominada “defoliación”, para introducir el efecto de los daños a nivel de biomasa foliar. 6.1.3. Función de desagregación Una función de desagregación es una herramienta que permite conocer o estimar, en un momento determinado, el número de árboles que hay en cada clase diamétrica, a partir de la información de una variable de rodal (variable agregada). Estas funciones informan sobre la estructura del rodal (Gorgoso et al., 2012) o del paisaje (Rubin et al., 2006), permitiendo planificar con mayor detalle las intervenciones y el destino de los productos (Scolforo et al., 2003). En las masas de eucalipto de Galicia la división de los productos es quizá más sencilla que en otras especies, como los pinos por ejemplo, ya que no se suelen considerar diferentes destinos para un mismo fuste. Tanto la celulosa como la energía requieren un producto triturado, y solamente en el mucho menos común caso de la madera estructural (en forma de madera laminada) se requiere obtener un elevado número de piezas por aserrado. A su vez, para cada destino, se espera extraer el máximo provecho de cada árbol (excluyendo las hojas), de forma que la producción de energía admite todo el árbol, mientras la celulosa o la madera estructural admiten toda la madera hasta un diámetro en punta delgada, permitiendo que la parte no aprovechada sea materia prima para bioenergía. Por otra parte, el reciente empleo de clones de Eucalyptus globulus en repoblación suscita interés por conocer el tipo de estructura de masa que se genera en las diferentes variedades clonales, donde se asume una gran semejanza de los árboles que conforman un rodal. Por esa razón, se considera importante el ajuste de una función de distribución que permita estimar el número de pies por clases diamétricas, así como describir la forma de la distribución en cada material de reproducción. Funciones de salida 79 Como función de distribución se seleccionó la función de Weibull, y en la obtención de parámetros se seleccionó el método de los momentos (ver Anexo II Análisis estadístico). Esta metodología ha sido empleada, entre otros autores, por Hyink (1980), Burk y Burkhart (1984), Páscoa (1987), Condés (1997), Del Río (1999), Scolforo et al. (2003), Palahí et al. (2007), Fonseca et al. (2009) o Gorgoso et al. (2012). Esta técnica brinda generalmente mejores estimaciones, incluso en proyecciones a largo plazo (Cao et al., 1982; Reynolds et al., 1988; Borders y Patterson, 1990; Torres-Rojo et al., 2000). Además, presenta la ventaja de que puede garantizar la compatibilidad entre las predicciones del área basimétrica estimada explícitamente (a través de un modelo de crecimiento en área basimétrica) y la predicha a partir de la distribución diamétrica (Hyink, 1980; Knoebel et al., 1986). En algunos casos (Baldwin y Feduccia, 1987; Cao, 2004) se ha utilizado una metodología híbrida, recuperando un parámetro a partir de un momento de la distribución y otro a partir de un percentil. 6.2. Material y Métodos 6.2.1. Datos 6.2.1.1. Tarifas de cubicación y ecuaciones de biomasa Los datos empleados para el ajuste de las tarifas de cubicación de árbol individual corresponden a los 87 árboles tipo descritos en el apartado 3.3.2 (Variables de árbol individual). Estos se pudieron dividir según materiales de reproducción en muestras equilibradas de tamaño n = 29, y según el tipo de suelo 34 árboles se apearon en suelos de esquistos/pizarras y 53 en suelos graníticos. En el caso de las tarifas de cubicación y ecuaciones de biomasa de rodal, los datos para el ajuste provienen de la suma de los volúmenes individuales y pesos estimados para los árboles de las parcelas de inventario, posteriormente transformados a la unidad superficial de referencia, la hectárea. Para el desarrollo de las tarifas de volumen y ecuaciones de biomasa de rodal el número de parcelas diferentes empleadas en el análisis fue de 114, correspondientes a 307 inventarios, incluyendo todas las pertenecientes a brinzales, clones F0 (Anselmo) y clones F1 (Odiel y Tinto) que fueron medidas al menos en una ocasión. En la Figura 6.1 se muestra el gráfico de dispersión altura dominante-área basimétrica para los datos de rodal. En la Tabla 6.1 se exponen los estadísticos descriptivos de las variables de rodal para las parcelas empleadas en el ajuste de tarifas de volumen y ecuaciones de biomasa de rodal. Juan Daniel García Villabrille 80 Figura 6.1. Gráfico de dispersión altura dominante-área basimétrica de los datos empleados en el ajuste de tarifas de cubicación de rodal. Tabla 6.1. Estadísticos descriptivos por inventario y materiales de reproducción de las variables elaboradas para las parcelas empleadas en el ajuste de las funciones de salida de volumen y biomasa de rodal. Brinzal Inventario 1 (n = 71) Inventario 2 (n = 64) Inventario 3 (n = 51) Media Máximo Mínimo Desv. est Media Máximo Mínimo Desv. est Media Máximo Mínimo Desv. est t 12 20 1 3 12 21 2 3 13 20 3 3 N 1111 1897 661 292 1058 1680 661 266 1028 1631 625 251 G 18,7 38,6 0,2 8,7 19,6 40,8 0,9 8,6 20,0 40,4 1,1 8,7 dg 14,4 25,3 1,4 4,5 15,1 25,7 3,5 4,3 15,4 26,3 4,2 4,5 d 13,2 23,6 1,3 4,2 14,0 24,0 3,2 4,1 14,4 24,7 3,8 4,2 d0 22,3 37,3 2,4 6,9 23,2 38,2 5,8 6,7 23,4 39,3 5,8 6,8 h 15,6 28,7 2,2 5,1 16,5 29,1 4,2 4,9 17,0 29,6 4,6 5,1 H0 20,9 39,1 3,1 6,9 21,9 40,3 5,9 6,7 22,1 41,6 5,7 6,7 Clon F0 (Anselmo) Inventario 1 (n = 25) Inventario 2 (n = 25) Inventario 3 (n = 25) Media Máximo Mínimo Desv. est Media Máximo Mínimo Desv. est Media Máximo Mínimo Desv. est t 7 11 4 1,84 8 12 5 1,84 9 13 6 1,84 N 1154 1540 780 207 1122 1540 772 201 1101 1540 731 205 G 11,9 19,5 3,4 4,8 13,0 20,7 4,7 4,9 13,8 21,3 5,4 4,9 dg 11,4 16,9 5,7 3,1 12,1 17,8 6,8 3,0 12,6 18,7 7,1 2,9 d 10,9 16,0 5,5 2,9 11,7 16,9 6,5 2,8 12,2 17,8 6,8 2,8 d0 15,5 25,0 8,3 3,9 16,4 27,0 10,0 4,0 17,1 27,8 10,8 3,9 h 12,2 17,5 6,0 3,4 13,5 18,5 6,7 3,4 14,3 19,7 7,2 3,5 H0 14,9 23,9 7,7 4,1 16,4 24,3 8,6 4,0 17,1 24,9 9,4 4,0 Clon F1 (Odiel y Tinto) Inventario 1 (n = 14) Inventario 2 (n = 14) Inventario 3 (n = 14) Media Máximo Mínimo Desv. est Media Máximo Mínimo Desv. est Media Máximo Mínimo Desv. est t 5 8 3 1,62 6 9 4 1,62 7 10 5 1,62 N 1154 1534 800 268 1135 1508 760 265 1125 1508 720 268 G 8,2 17,0 3,0 4,3 9,8 18,7 4,6 4,3 12,6 21,7 5,9 3,9 dg 9,5 14,8 5,1 3,2 10,5 15,9 6,4 3,0 12,0 16,8 8,8 2,5 d 8,9 13,7 5,0 2,9 9,9 15,0 6,1 2,8 11,4 16,1 8,5 2,4 d0 14,2 21,5 7,3 4,6 15,6 22,4 9,3 4,0 18,0 23,1 13,0 3,1 h 9,6 16,3 5,4 3,4 11,5 17,8 7,1 3,5 13,4 18,8 9,4 3,3 H0 12,4 19,7 7,2 4,1 14,9 21,1 9,0 3,9 17,1 23,4 11,2 3,7 0 10 20 30 40 10 20 30 40 Área basimétrica ( m 2 / ha) Altura dominante (m) Funciones de salida 87      − − − ⋅      − ⋅       =b ax c e b ax b c xf 1 )( [6.4] donde x es la variable aleatoria, a es un parámetro de situación que define el origen de la función, b un parámetro de escala y c un parámetro de forma que es el responsable de la asimetría de la función. Todos estos parámetros toman valores mayores o iguales que cero y según su variación la función de Weibull presenta formas muy diversas. Si c = 1, la función de Weibull coincide con la distribución exponencial; si 1 < c < 3,6 la distribución presenta asimetría positiva; si c = 3,6 se aproxima a la distribución normal; y si c > 3,6 presenta asimetría negativa. En el caso de que se asigne un valor cero al parámetro de situación a, se obtiene la denominada función de densidad de Weibull biparamétrica (de dos parámetros), que toma la siguiente expresión:       − − ⋅       ⋅       =b x c e b x b c xf 1 )( [6.5] La ecuación que se ha empleado para realizar el ajuste de las distribuciones de frecuencias relativas acumuladas se obtiene integrando la expresión [6.4], cuya expresión matemática es la siguiente: cc b ax xb ax c edxe b ax b c xF      − −      − − − −=⋅⋅      − ⋅       = ∫ 1 )( 0 1 [6.6] siendo F(x) la función de distribución de Weibull con la que se obtiene la frecuencia relativa acumulada de árboles con un diámetro menor o igual a x. Como ya se comentó, antes de obtener los parámetros de una función de distribución, es necesario ajustar modelos de diámetro medio aritmético (d) a partir de variables de estado del rodal. Los modelos seleccionados para predicción del diámetro medio aritmético se muestran en la Tabla 6.6, tomados de Diéguez-Aranda et al. (2009) todos ellos y de la forma mostrada en la siguiente expresión: ( ) β ,exp Xdd g−= [6.7] Juan Daniel García Villabrille 88 Tabla 6.6. Modelos ajustados para la predicción de diámetro medio aritmético. Modelo Expresión M1 ( ) tbHbbdd g⋅+⋅+−= 2010 exp M2 ( ) tbHbbdd g/exp 2010 +⋅+−= M3 ( ) NbHbbdd g⋅+⋅+−= 2010 exp M4 ( ) tbISbbdd g⋅+⋅+−= 210 exp M5 ( ) NbtbISbbdd g⋅+⋅+⋅+−= 3210 exp M6 ( ) 010 exp Hbbdd g⋅+−= Nota: d es el diámetro medio aritmético (cm), dg es el diámetro medio cuadrático (cm), b0, b1, b2,.., son los parámetros a determinar, t es la edad (años), H0 es la altura dominante (m), N es el número de pies por hectárea (pies/ha) e IS es el índice de sitio (m). 6.2.3. Ajuste de modelos En el Anexo II (Análisis estadístico) se exponen los fundamentos básicos de cada tipo de ajuste, abordando ahora las particularidades de cada uno de ellos. 6.2.3.1. Tarifas de cubicación En el ajuste de tarifas de cubicación es común que la dispersión de los residuos de las regresiones por mínimos cuadrados (OLS) tienda a aumentar con los valores predichos, apareciendo un caso típico de heterocedasticidad (falta de homogeneidad en la varianza del error). Para corregir este efecto se puede optar por transformar las variables (tomando logaritmos en ambos términos de la expresión), o mediante regresión ponderada (Parresol, 1999). En este caso se optó por emplear regresión ponderada, para lo cual se asoció a cada observación un peso igual a la inversa de la varianza de su error, 1/σ2, durante el proceso de ajuste. Siguiendo a Neter et al. (1996), se empleó una función potencial del tipo: k ii X= 2 σ [6.8] donde se debe determinar el valor del exponente k. Para ello se utilizó la metodología propuesta por Park (1966), que consiste en emplear los errores del modelo ajustado sin pesos, êi, como variable dependiente en el modelo potencial de varianza del error, es decir: Funciones de salida 89 k ii Xe ⋅= γ 2 ˆ [6.9] o bien: ii LnXkLneLn ⋅+= γ 2 ˆ [6.10] donde γ es una constante y la variable X toma diferentes valores en función de las variables del modelo ajustado, siendo los factores de ponderación resultantes de la forma 1/Xk con X = d, h, d·h, etc., para tarifas de árbol individual, o X = G, H0, G·H0, etc., para tarifas de rodal. Para tratar de identificar diferencias entre materiales de reproducción (brinzales, clones F0 y clones F1), suelos (granitos o esquistos/pizarras) y promotor (individual o empresa) se ajustaron modelos usando la metodología de variables dummy. Esta metodología consiste en la expansión de los parámetros ai del modelo de la siguiente forma general: r+jiij Ia+a ⋅ 1 [6.11] donde aij son los parámetros del modelo e Ir identifica las diferencias de una muestra dividida en dos grupos, como es en el caso de tratar de identificar diferencias entre tipos de suelos o entre brinzales y clones (éstos en un único grupo). En el caso de los materiales de reproducción, en el que la muestra se puede dividir en tres grupos, la expresión general de cada uno de los parámetros ai expandidos es la siguiente: sjir+jiij IaIa+a ⋅+⋅ +2 1 [6.12] donde aij son los parámetros del modelo e Ir e Is identifican las diferencias de una muestra dividida en tres grupos. Asociado a cada ajuste con variables dummy se emplearon dos test generalmente utilizados para caracterizar diferencias entre regiones geográficas (Álvarez-González et al., 2005b; Calama et al., 2003; Huang et al., 2000) y así poder evaluar la necesidad de desarrollar modelos específicos. Esos test fueron el test L de Lakkis-Jones (Khattree y Naik, 1999) y el estadístico F* de la “suma de cuadrados no lineal” (Bates y Watts, 1988). Para testar los resultados de los ajustes iniciales de las tarifas de cubicación individual y asegurar la validez del resto del procedimiento se aplicó el test F (Montgomery y Peck, 1992). Juan Daniel García Villabrille 90 El ajuste de tarifas de cubicación se llevó a cabo por mínimos cuadrados (OLS) para la selección de modelos, y por mínimos cuadrados ponderados (WLS) cuando la varianza del error de los modelos seleccionados presentó heterocedasticidad, empleando la función nls del software R (R Core Team, 2014). 6.2.3.2. Ecuaciones de biomasa El ajuste de sistemas de ecuaciones de biomasa se basó en la técnica de ajuste simultáneo de ecuaciones aparentemente no relacionadas (SUR), aunque existen otras metodologías que también permiten cumplir con la aditividad y obtener resultados similares (Tomé et al., 2001; Dong et al., 2015). No obstante, antes de ajustar un determinado sistema de ecuaciones, y siguiendo a Baskerville (1972), se realizó un análisis exploratorio de la varianza de las diferentes fracciones de biomasa según las diferentes clases diamétricas. En dicho análisis se observó que dicha varianza no era constante, por lo que el autor recomienda el ajuste de sistemas de ecuaciones linealizadas (empleando logaritmos). En estos casos, además, se debe considerar que la conversión de las estimaciones logarítmicas insesgadas de la media y la varianza de nuevo a unidades aritméticas no es directa. Esto resulta del hecho de que si la distribución del Ln (Y) (observado) para una Xi (variable independiente) dada es normal, la distribución de Y puede no ser normal y además será sesgada. De hecho, si la distribución es normal en forma logarítmica, la solución de la ecuación linealizada para una Xi y la determinación del antilogaritmo de Ln (Yi) proporcionan la mediana de la distribución aritmética sesgada en lugar de la media (Finney, 1941; Brownlee, 1967). Para la corrección del sesgo de las estimaciones, Brownlee (1967) propuso emplear la siguiente expresión:        += 2 ˆ exp ˆ 2 σ µ Y [6.13] o lo que es lo mismo: ( ) ( )( )         ⋅=         ⋅= 2 expexp 2 exp ˆ exp ˆ 22 σσ µ est YLnY [6.14] donde Y ˆ es el valor predicho deshaciendo la transformación logarítmica, µ ˆ es el valor estimado para Ln (Y) denotado por Ln (Y)est, y  es la varianza del error. Por otra parte, aunque no existieron diferencias significativas entre materiales de reproducción para el volumen, sí se analizaron las posibles diferencias que pudieran existir en las variables ww.rel, wb.rel, wbr.rel, cl.rel, dap.ma, dap.tr, dbas.ma y dbas.tr, empleando test estadísticos de normalidad (Shapiro-Wilk, 1965) y test de comparación de medias (test t, Funciones de salida 91 test de Wilcoxon y test de Kruskal-Wallis) (ver Anexo I. Notación, y Anexo II. Análisis estadístico). Para los modelos independientes de biomasa de corteza y madera de árbol individual, además, se aplicó el proceso de cross validation. Para ello se realizaron divisiones de la muestra original, manteniendo un 70% de los datos para el ajuste y un 30% para la validación, y se estableció como criterio de selección el mejor REMC promedio para todas las iteraciones de cada modelo. En ambas fracciones se realizaron 1000 divisiones de la muestra, o iteraciones, sobre los mejores modelos seleccionados en fases previas de ajuste. El ajuste de sistemas aditivos de biomasa se llevó a cabo por regresión aparentemente no relacionada (SUR), empleado la función nlsystemfit del software R (R Core Team, 2014). 6.2.3.3. Función de desagregación En este trabajo se empleó el valor de a = 0 (función Weibull biparamétrica), ya que ajustes previos para predicción del diámetro mínimo no resultaron en predicciones aceptables (R2 = 0,25). Wang y Baker (2007) obtuvieron ajustes del diámetro mínimo superiores (R2 = 0,55) para Eucalyptus globulus. García-Guemes et al. (2002) o Palahí et al. (2007) emplearon a = 0, mientras Gorgoso et al. (2012) testaron a = 0, a = dmin, a = 0,5·dmin y el estimador de Zanakis (1979), obteniendo los mejores resultados con la función Weibull biparamétrica (a = 0). Por otra parte, Condés (1997) señala que los ajustes con la función Weibull biparamétrica no alcanzan la precisión adecuada cuando se trata de modelizar la distribución diamétrica de parcelas con diámetros medios elevados, pero esta situación raramente se presenta en Eucalyptus globulus en turnos de celulosa, siendo el diámetro medio para los datos empleados entre 12 y 16 años de 14,99 cm. 6.2.4. Comparación de distribuciones por materiales de reproducción Con el fin de analizar la estructura de las distribuciones diamétricas en cada material de reproducción se calculó un índice que indica la forma de la curva en cada parcela. Dicho índice, denominado genéricamente índice de forma de la curva (ifc), elaborado para analizar la amplitud diamétrica en cada material de reproducción, es la relación entre la máxima frecuencia relativa observada en cada parcela y el rango de su distribución, cuya expresión es la siguiente: i i dd xF ifc )( ))(max( minmax − = [6.15] Juan Daniel García Villabrille 92 donde max (F (x))i es la máxima frecuencia relativa observada en la parcela i y dmax y dmin son los diámetros extremos de la distribución de la parcela i. Este índice no tiene en cuenta una distribución de referencia; simplemente muestra, para un rango de valores, su grado de concentración. Puesto que no se instalaron parcelas de todos los materiales de reproducción en todas las edades, la comparación de este índice se limitó al rango de edad en el que están presentes todos ellos, hasta un máximo de 10 años. 6.2.5. Comparación de modelos Como se expuso en el apartado 3.5 (Evaluación de modelos), la comparación de modelos se basó en los estadísticos de ajuste R2, REMC y E, y en el análisis gráfico de los residuos (ei) frente a los valores predichos ( i y ˆ) y de éstos frente a los observados (yi). 6.3. Resultados 6.3.1. Tarifas de cubicación de árbol individual En la Tabla 6.7 se muestra el resultado de la estimación de los parámetros y los estadísticos del ajuste de los modelos analizados, para la selección de la tarifa de cubicación de árbol individual con corteza. En los modelos que contienen término independiente (M4 y M7) dicho parámetro no resultó significativo, quedando descartados aquellos modelos que carecen de sentido biológico. A excepción del modelo M8, que incluye el diámetro normal (dn) y la altura de copa (hblc), el resto obtuvieron estadísticos muy similares, con pequeñas diferencias, con R2 > 0,993. El REMC, que fue menor de 0,012 m3 para estos modelos, correspondió al 6,4 y el 4,8% del volumen medio de todos los árboles inventariados entre 14 a 16 y 14 a 23 años respectivamente, siendo estos intervalos en los que se sitúa el momento de la corta en Galicia (Álvarez-González et al, 2005a). El modelo con mejores estadísticos y con todos los parámetros significativos fue el M1, el mismo que se utilizó en las tarifas existentes hasta la fecha en Galicia para Ecualyptus globulus (Diéguez-Aranda et al., 2009), también empleado por Tomé et al. (2007) en Portugal, incluso para tarifas sin corteza, y uno de los mejores en el ajuste de diferentes modelos realizados por Idoiagabeitia y Fernández-Manso (2001) para la especie en Vizcaya. También Bi (1994) empleó este mismo modelo para seis especies diferentes de eucaliptos (E. cypellocarpa, E. fastigiata, E. nitens, E. obliqua, E. radiata y E. viminalis). A su vez, es el modelo que contiene las variables más frecuentemente tomadas en cualquier inventario (diámetro normal dn y altura total ht), lo que supone una ventaja frente a resultados que hubiesen requerido otras variables más difíciles de obtener Funciones de salida 93 como el diámetro normal sin corteza (dsc), o bien de localización subjetiva como la altura de copa (hblc). Tabla 6.7. Resultados del ajuste para los modelos de tarifas de cubicación con corteza de árbol individual. Modelo Parámetros Ajuste b0 b1 b2 b3 b4 R2 REMC E M1 3,236·10 - 5 (<0,001) 1,896 (<0,001) 1,1106 (<0,001) 0,9936 0,0116 -0,001 M2 3,019·10 - 5 (<0,001) 1,030 (<0,001) 0,9932 0,0119 0,000 M3 3,309·10 - 5 (<0,001) 0,9932 0,0119 -0,000 M4 -2,425·10 - 3 (0,352) 3,697·10 - 5 (<0,001) 1,873 (<0,001) 1,094 (<0,001) 0,9936 0,0116 0,000 M5 3,26·10 - 5 (<0,001) 1,899 (<0,001) 1,098 (<0,001) 8,957·10 - 3 (0,825) 0,9935 0,0116 -0,001 M6 3,107·10 - 5 (<0,001) 9,952·10 - 1 (<0,001) 0,02915 (0,474) 0,9932 0,0119 0,000 M7 -2,405·10 - 3 (0,359) 3,716·10 - 5 (<0,001) 1,875 (<0,001) 1,084 (<0,001) 7,3·10 - 3 (0,856) 0,9935 0,0117 0,000 M8 9,524·10 - 5 (<0,001) 2,255 (<0,001) 0,4608 (<0,001) 0,9772 0,0219 -0,001 M9 3,554·10 - 5 (<0,001) 1,87597 (<0,001) 1,159 (<0,001) 0,9924 0,0125 -0,001 El valor entre paréntesis indica el nivel de significación. Tras seleccionar el modelo M1 como el mejor de la fase de ajuste, se realizaron posteriores ajustes con variables dummy para comprobar si existían diferencias en las muestras en función de los diferentes materiales de reproducción y en función de los distintos suelos de procedencia, y por tanto proponer un modelo completo (considerando esas diferencias) o bien uno reducido. Entre ambas opciones resultó mejor el modelo que crea diferencias según el tipo de suelo, con R2 = 0,9942 y REMC = 0,0109, que mejoran los valores obtenidos en el ajuste inicial. Respecto a los parámetros, para el caso de los materiales de reproducción únicamente el parámetro que acompañó al diámetro normal (dn), y sólo para los clones F0, fue significativamente diferente de cero al nivel de α = 5%. En el caso del modelo en función del suelo sólo fueron significativamente diferentes de cero los parámetros que acompañaron a la variable de la altura (ht), con valores 1,204 (p-valor < 0,0001) para pizarras y 0,917 (p-valor < 0,001) para granitos. Los posteriores resultados de los test F* (suma extra de cuadrados no lineal) y L (Lakkis-Jones) para cada uno de los supuestos iniciales, habiendo fijado el nivel de significación al α = 5%, indicaron que no sería necesario elaborar modelos diferentes según los diferentes materiales de reproducción (F* con p-valor = 0,228; L con p-valor = 0,182). En el caso del tipo de suelo, pareció necesario proseguir con el análisis de un modelo completo, habiendo una mejora del 14,3% en la suma de cuadrados del error (SSE). El Juan Daniel García Villabrille 94 posterior análisis gráfico de los residuos de este modelo reveló una marcada heterocedasticidad, por lo que fue ajustado mediante regresión ponderada y expandiendo sólo el parámetro que había resultado significativo en el ajuste inicial (b2). En esta ocasión, el parámetro expandido dejó de ser significativamente diferente de cero (pvalor > α = 0,05), por lo que finalmente se retomó el modelo reducido como la mejor opción. Al igual que para el caso del modelo completo, el análisis gráfico de residuos también indicó que éstos eran heterocedásticos, por lo que finalmente se ajustó un modelo con peso = 1/(dn·ht) 2,1113. En la Tabla 6.8 se muestran los resultados de este ajuste, con R2 = 0,993 y REMC = 0,012, habiendo desaparecido la ligera sobrestimación del ajuste inicial, siendo los valores predichos muy similares a los valores observados en todo el intervalos de predicción (Figura 6.3). Tabla 6.8. Resultados del ajuste con pesos para volumen individual con corteza. Modelo Parámetros Ajuste b 0 b 1 b 2 R 2 REMC E M1 2,993·10 - 5 (<0,001) 1,973 (<0,001) 1,060 (<0,001) 0,9933 0,012 0,000 El valor entre paréntesis indica el nivel de significación. Dichos estadísticos fueron comparados con los que generó el modelo disponible hasta la fecha (Diéguez-Aranda et al., 2009) sobre los datos de partida, con REMC = 0,015 y E = -0,001. Los valores de los estadísticos de este modelo, aunque inferiores al ajustado en este trabajo, se consideran igualmente válidos teniendo en cuenta que se desarrolló a partir de datos de árboles apeados en turno, sólo procedentes de semilla, de dos localizaciones distintas y con menor volumen de datos por árbol, si bien la muestra era de 165 árboles (mayor que los 87 empleados en esta ocasión). Además, este modelo permite calcular volúmenes por destinos comerciales, ya que consta también de una función de perfil compatible. Figura 6.3. Gráficos de valores observados frente a valores predichos para el modelo de volumen individual con corteza (m3/pie). 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.0 0.1 0.2 0.3 0.4 0.5 0.6 Predichos Observados Funciones de salida 95 En el caso del volumen individual sin corteza, el mejor modelo en la fase de selección también resultó ser el M1. Al igual que para el caso del volumen con corteza, se realizaron ajustes de modelos completos y reducidos para los mismos casos de estudio (materiales de reproducción y suelos). Si bien en el caso de los materiales de reproducción ningún parámetro del ajuste fue significativamente diferente de cero, en el caso de los tipos de suelo los dos parámetros b1 y b2, que acompañan a dn y ht respectivamente, sí lo fueron. Igualmente fue necesario un ajuste con regresión ponderada, ya que los residuos del modelo completo fueron heterocedásticos, pero al contrario que en el caso del volumen con corteza, los parámetros expandidos siguieron siendo significativamente diferentes de cero para un nivel de probabilidad del 95%. El resultado de los test F* y L realizados fue diferente según el test aplicado, no siendo necesario desarrollar modelos según el test F* (p-valor = 0,054) pero sí según el test L (p-valor = 0,043), por lo que se consideró que sería necesario otro contraste. Así, se compararon los estadísticos del modelo completo para suelos con los del modelo reducido, ajustados ambos con pesos, siendo el factor de ponderación empleado para el modelo reducido 1/(dn·ht) 2,2575. En la Tabla 6.9 se muestra el resultado del ajuste de ambos modelos. Tabla 6.9. Resultados del ajuste con pesos para volumen individual sin corteza. Modelo Parámetros Ajuste b0 b10 b11 b20 b21 R2 REMC E M1 reducido 2,0944·10-5 (<0,001) 1,997 (<0,001) 1,101 (<0,001) 0,9902 0,0122 -0,000 M1 completo 2,018·10-5 (<0,001) 1,761 (<0,001) 0,3031256 (0,012) 1,343 (<0,001) -0,2931 (0,016) 0,9899 0,0123 -0,001 El valor entre paréntesis indica el nivel de significación. El modelo reducido, o independiente de la muestra, ofreció mejores resultados que un modelo que tiene en cuenta el tipo de suelo. Este resultado coincidió con otro contraste empleado, que fue el ajuste de modelos de volumen de corteza donde la expansión de los parámetros no resultó significativa para ninguno de los factores estudiados (materiales de reproducción y suelo). Tampoco el ajuste de un modelo de volumen de corteza a partir de una relación directa de esta misma variable con el volumen con corteza se consideró adecuado, ya que supondría asumir que el porcentaje de corteza es fijo para todas las edades, cuando se pudo comprobar que éste disminuye con el aumento de las dimensiones del árbol. En la Figura 6.4 se muestra el gráfico de valores observados frente a valores predichos, y al igual que en caso del volumen con corteza, se considera que hay buena distribución de los valores presentados. Juan Daniel García Villabrille 96 Figura 6.4. Gráficos de valores observados frente a valores predichos para el modelo de volumen individual sin corteza (m3/pie). 6.3.2. Tarifas de cubicación de rodal Como se ha indicado en el apartado anterior, en el caso de las tarifas de árbol individual no se encontraron diferencias entre los materiales de reproducción, y por esa razón el análisis para las tarifas de rodal se centró únicamente en la selección del mejor modelo. En la Tabla 6.10 se muestra el resultado de los ajustes para volumen de rodal con corteza. A excepción de los modelos M7 y M11, el resto de los modelos obtuvieron todos sus parámetros significativos al nivel de probabilidad del 95%. De forma general, los modelos con un nivel de información medio ofrecieron los mejores estadísticos, y las mejores variables fueron la densidad (N), el área basimétrica (G), la altura dominante (H0) y el diámetro dominante (d0), lo que coincide con los resultados de la selección de variables por stepwise mediante la función stepAIC del software R (R Core Team, 2014). Cabe destacar que la edad (t), dato que habitualmente es conocido a priori por el promotor, siempre fue una variable acompañada por parámetros significativamente diferentes de cero, pero sin embargo ninguno de los modelos en los que se incluyó fue el mejor. De todos los modelos ajustados, el mejor R2 correspondió al modelo M3 (0,9973), y también el REMC fue el más bajo (5,56 m3/ha), aunque respecto al E = -0,3098, tres modelos obtuvieron menores valores absolutos (M1, M5 y M6), pero sólo uno de ellos (M1) estuvo del lado de la subestimación. Para dicho modelo (M3), se realizó un análisis de correlación entre las variables, con la función cor.test del software R (R Core Team, 2014), y se comprobó que ésta era muy alta entre d0 y el resto de variables (> 0,9). Por otra parte, la mejora que se produce por introducir la variable d0 es muy pequeña. Finalmente se seleccionó como mejor modelo el M5, que sólo depende de las variables que se proyectan en el modelo dinámico (N, G y H0) y cuyos estadísticos de ajuste fueron muy similares a los del modelo M3. 0.0 0.1 0.2 0.3 0.4 0.5 0.0 0.1 0.2 0.3 0.4 0.5 Predichos Observados Funciones de salida 103 muestra: por un lado la muestra de brinzal y por otro la muestra compuesta por ambos clones, siendo el modelo testado de la forma 21 0b t b nbr hdbw ⋅⋅= . El primer ajuste indicó que la expansión del parámetro b0 no era significativamente diferente de cero para α = 0,05, por lo que se reajustó el modelo, expandiendo b1 y b2. El resultado indicó que la expansión del parámetro b2 dejó de ser significativa, y así sucesivamente hasta que finalmente se concluyó que el modelo reducido (sin identificar grupos) era el único con todos los parámetros significativos. Por otra parte, y ya que los test indicaron que la densidad básica en la madera de clones es mayor, entonces se deberían obtener mayores valores de biomasa en clones que en brinzales para un mismo diámetro. No obstante, se considera que existen limitaciones en este trabajo, ya que el muestreo destructivo no fue diseñado para la estimación de la densidad, sino para la estimación de biomasa. En un análisis, se realizó un ajuste de modelos alométricos sencillos, relacionando la biomasa de madera con el diámetro, sin tener en cuenta otras variables. El primer ajuste (ajuste1) consistió en un modelo para cada uno de los tres materiales de reproducción con todos los datos. El segundo ajuste (ajuste2) consistió también en un ajuste para cada material de reproducción, pero limitando la edad máxima observada de brinzales a la de los clones (muestra más joven), que fue de 10 años. Por último se agruparon los clones, ajustando un modelo para clones y otro para brinzales (ajuste3), además de mantener esta muestra en la misma edad máxima (10 años). En la Figura 6.7 se muestran las predicciones de los modelos ajustados, y las correspondientes densidades básicas de las muestras del ajuste. En el primer ajuste se observa que, si bien es esperable que las curvas generadas para clones estén por encima, dada la mayor densidad básica de la madera, al ser la muestra más joven y, por tanto, menos densa, aun no se ha producido dicho efecto. En el segundo y tercer ajuste se igualan las predicciones, pero se debe tener en cuenta que la edad media de brinzales era ligeramente inferior a la edad media en clones, por lo que tampoco se puede asegurar, a partir de estos ajustes, que los clones presenten mayor densidad básica cuando se emplean todas las trozas para el ajuste. Analizando las densidades para la misma edad máxima, se comprobó que se muestrearon diámetros mayores en las masas clonales, siendo masas de mayores densidades medias (N = 1140 pies/ha en clones frente a N = 1070 pies/ha en brinzales). Si bien no existieron diferencias respecto de la biomasa de madera a nivel de árbol individual, tanto las diferencias en densidad básica como estos ajustes indican que es esperable, pero no demostrable con total seguridad a partir de los datos, una mayor producción en biomasa, y posiblemente un mayor crecimiento en diámetro, para masas clonales frente a brinzales de igual edad. Juan Daniel García Villabrille 104 Figura 6.7. Arriba: Predicciones de los modelos de biomasa de madera (ww, kg/pie) por materiales de reproducción; Debajo: Densidad básica de madera (dbas.ma, kg/m3) para las muestras de ajuste. De izquierda a derecha: ajuste1, ajuste2 y ajuste3. A partir de estos resultados, y aunque las curvas elaboradas empleando modelos por separado son similares, se ajustó una ecuación de biomasa con variables dummy para diferenciar madera de brinzal de madera clonal dentro del sistema aditivo completo. Los errores de este modelo fueron superiores a los de los modelos presentados más adelante en la Tabla 6.18, por lo que no se consideró oportuno, para los datos de partida, emplear ecuaciones diferentes. Se asume que estos resultados están condicionados por la falta de datos de árboles tipo clonales entre los 10 y los 17 años (máximo muestreado para brinzal), a pesar de tener estados dimensionales similares y una cantidad de datos bastante equilibrada entre las muestras. Por otra parte, y como validación de las tarifas de cubicación de árbol individual y de los resultados de densidad en verde, se probó el método de estimación de pesos en verde en dos aprovechamientos de montes privados para los que se tuvo la oportunidad de realizar predicciones: 1.- El primer aprovechamiento correspondió a un monte de brinzal de 0,63 ha y 40 años, en el cual se montaron dos parcelas de 350 y 500 m2 durante el primer inventario. La estimación final, obtenida como el volumen con corteza multiplicado por una densidad de 1,150 t/m3 (estimada a partir de la muestra completa de trozas) y por unas pérdidas estimadas de humedad del 10%, sólo pudo ser comparada con la cifra de venta final, pagada al precio de 32 €/t. El resultado final estimado fue de 18 870 €, siendo pagados al vendedor 18 900 €, lo cual resulta en una diferencia de 30 € por defecto (error del 0,16%). 2.- En el segundo caso, se realizó una estimación para un aprovechamiento de 25 ha de brinzal, en el que se montaron seis parcelas de muestreo de 200 m2, dos para cada una de las tres diferentes calidades identificadas en la visita de campo. Para las estimaciones 5 10 15 20 25 0 50 100 150 200 250 Diámetro normal (cm) Predicho ww (kg/pie) Brinzal Clon F0 Clon F1 0 5 10 15 20 25 30 0 50 100 150 200 250 Diámetro normal (cm) Brinzal Clon F0 Clon F1 0 5 10 15 20 25 30 0 50 100 150 200 250 Diámetro normal (cm) Brinzal Clon Brinzal clon F0 clon F1 400 500 600 700 800 Material de reproducción dbas.ma (kg/ m3) Brinzal clon F0 clon F1 400 500 600 700 800 Material de reproducción Brinzal Clon 400 500 600 700 800 Material de reproducción Funciones de salida 105 se emplearon la misma densidad verde estimada y las mismas pérdidas de humedad que en el caso anterior. La cifra obtenida, 5871 toneladas, pudo ser comparadas con los pesajes en báscula del comprador, que fueron de 3293 toneladas de madera descortezada y 2154 toneladas con corteza. Se estimó el porcentaje de corteza para la masa a partir de las tarifas con y sin corteza, y su volumen fue sumado al volumen de la madera extraída sin corteza, resultado finalmente 5842 toneladas con corteza, habiendo una diferencia final de 29 toneladas por exceso (error del 0,49%). En el momento de realización de la primera prueba aún no se había ajustado la tarifa de rodal, por lo que necesariamente se calcularon los volúmenes a partir de cubicaciones individuales transformadas a la hectárea posteriormente. En el segundo caso se siguió el mismo procedimiento, pero además se emplearon las tarifas de rodal (con muestra no definitiva) para comparar el resultado. En este caso, el resultado estimado con las tarifas de rodal fue de 5933 toneladas, lo que supondría un error por exceso de 91 toneladas, lo que triplica el error de la estimación pie a pie, si bien sigue siendo un error bajo. Así, se considera que las tarifas de rodal son realistas, siendo aconsejable el inventario pie a pie en los momentos de mayor trascendencia, como el período de establecimiento de la masa (1 a 3 años), la edad de referencia (7 años), el momento aproximado del turno según máxima renta en volumen o biomasa (12-15 años) o el momento de la corta más allá del turno. Parte de la importancia de estas pruebas radica no sólo en las predicciones del modelo, sino que facilita la localización de rodales de propiedad privada de edades superiores a los 15 años, lo cual también es indicativo del tipo de inversión que supone un eucaliptal. El propietario privado no sólo busca una especie que de ingresos altos, además alarga el turno bajo el criterio de la oferta que satisface sus estimaciones o necesidades momentáneas. Esto está alejado del máximo rendimiento financiero, principalmente por falta de gestión (Díaz-Balteiro, 2008). 6.3.4. Ecuaciones de biomasa de árbol individual En la Tabla 6.16 se muestra el resultado del ajuste independiente (sin aditividad) de las mejores ecuaciones para cada una de las fracciones del árbol. Todos los parámetros fueron significativamente diferentes de cero para α = 0,05, obteniendo los mejores valores de los estadísticos con modelos que incluyeron variables poco frecuentes. Así, se obtuvieron como variables útiles la longitud de copa (cl) para predicción de biomasa de hojas y el espesor de corteza a la altura normal (ecn) para predicción de biomasa en madera y corteza. En esta fase se realizó una validación cruzada, o cross validation, para las fracciones de corteza y madera, centrado en el análisis de pocas variables combinadas de distintas maneras. El resultado indicó que los modelos expuestos en la Tabla 6.16 fueron seleccionados en más del 40% de las ocasiones frente a otros modelos más sencillos (Tabla 6.17). Juan Daniel García Villabrille 106 Para la predicción de biomasa de hojas se presenta, además de la mejor ecuación (denotada como Hojas1), el modelo que incluyó la variable defoliación (def) (denotado por Hojas2). A pesar de ser una variable de medición subjetiva, mejoró el modelo que sólo incluía dn y ht. Sin embargo, la variable cl parece ser mejor para cuantificar la biomasa foliar. Esta variable, además de ofrecer una medida que varía con las dimensiones del árbol y la edad, está condicionada por la propia defoliación. De cualquier forma, para un inventario diamétrico en el que no se disponga de los datos de longitudes de copa se podría emplear el modelo que incluye la defoliación, usando un valor de def establecido visualmente. Como se comentó, la defoliación es mayor en los árboles en luz, por lo que la asignación de la defoliación debe considerar al menos dos estados diferentes en el rodal, diferenciando el grado de defoliación de los pies sumergidos por un lado y el de los intermedios y dominantes por otro. En el caso de las masas gallegas, la defoliación podría ser el factor de mayor influencia en la dificultad para mejorar las estimaciones a nivel de biomasa en las fracciones terminales. Tabla 6.16. Resultados del ajuste de modelos de biomasa de árbol individual sin aditividad. Fracción Modelo Parámetros Ajuste b0 b1 b2 b3 R2 REMC Hojas1 210 b t cbb nl hdw l ⋅= ⋅− 2,395 (<0,001) -0,031 (<0,001) -2,032 (<0,001) 0,6893 1,89 Hojas2 210 b t defbb nl hdw ⋅= ⋅− 2,914 (<0,001) 0,04186 (0,019) -2,306 (<0,001) 0,6467 2,01 Ramillos 21 0b t b sct hdbw ⋅⋅= 0,0613 (0,046) 3,181 (<0,001) -1,823 (<0,001) 0,7572 0,68 R. gruesas + finas 1 0b nbr dbw ⋅= 0,0838 (0,00) 1,601 (<0,001) 0,7810 2,15 Madera 221 0b t ebb nw hdbw cn ⋅⋅= ⋅− 0,0089 (<0,001) 1,982 (<0,001) 0,0339 (0,00) 1,222 (<0,001) 0,9863 7,71 Corteza 321 0b cn b t b nb ehdbw ⋅⋅⋅= 0,0296 (0,00) 1,438 (<0,001) 0,6524 (<0,001) 0,6567 (<0,001) 0,9728 1,32 Total 321 0b b b t b ntot dhdbw ⋅⋅⋅= 0,0268 (0,00) 1,640 (<0,001) 0,8938 (<0,001) 0,3311 (0,010) 0,9887 8,64 El valor entre paréntesis indica el nivel de significación. Funciones de salida 107 Tabla 6.17. Resultados de cross validation para modelos de biomasa de corteza y biomasa de madera de árboles individuales. Modelo biomasa corteza veces mejor REMC Modelo biomasa madera veces mejor REMC 321 0b cn b t b nb ehdbw ⋅⋅⋅= 435 221 0b t ebb nw hdbw cn ⋅⋅= ⋅− 436 21 0b cn b nb edbw ⋅⋅= 230 321 0b b b t b nw dhdbw ⋅⋅⋅= 257 321 0b t ebb nb hdbw cn ⋅⋅= ⋅− 213 21 0b t b nw hdbw ⋅⋅= 219 21 0b t b nb hdbw ⋅⋅= 122 321 0b blc b t b nw hhdbw ⋅⋅⋅= 88 António et al. (2007) observaron que la inclusión de variables de copa para predicción de la fracción de hojas en E. globulus implicaba una mejora del 29% de la suma de cuadrados del error (SSE). En este trabajo, la inclusión de cl en el modelo de hojas supuso una mejora del 17,83% del SSE, pero no así para las ramas, para las que António et al. (2007) obtuvieron una mejora el 19% en dicha fracción. En este sentido, los resultados obtenidos por António et al. (2007) no concuerdan para las fracciones menores (ramas y hojas), a pesar de que cl informa sobre las dimensiones relativas a la copa aportando mejoras sobre los modelos (Carvalho y Parresol, 2003). La alta variabilidad de las fracciones menores se ve afectada por dos circunstancias principales: a veces no existen las ramas gruesas y los daños por defoliación afectan de forma importante a la cantidad de ramillos y, sobre todo, de hojas. Para la estimación de biomasa de madera y de corteza, la inclusión de la variable espesor de corteza (ecn) supuso una mejora de las estimaciones, sobre todo, y como era de esperar, en el caso de la corteza. A partir de los resultados de estos ajustes (biomasa de madera y corteza) se comprobó mediante el test de Wilcoxon si existían diferencias en la relación ecn/dn para cada material de reproducción. El resultado indicó que los clones mantienen similitud (p-valor = 0,9191) mientras ambos son diferentes a los brinzales (pvalor = 0,02 brinzal frente a clon F0; p-valor = 0,03 brinzal frente a clon F1). De esta forma, cuando los parámetros estimados por variables dummy para separar los materiales de reproducción no resultaron significativos, un modelo que incluya ecn como variable explicativa podría ofrecer estimaciones más específicas, sin tener que separar las muestras para el ajuste. Wei y Borralho (1997) encontraron en un estudio sobre Eucalyptus urophylla que el grosor de corteza y el grosor de corteza relativo tienen alta heredabilidad1 (h2 = 0,45 y h2 = 0,40 respectivamente), mientras López et al. (2002) indican un valor medio (h2 = 0,32) para la corteza de Eucalyptus globulus. Se debe tener 1 Heredabilidad: definida según la FAO como la proporción de la variación fenotípica para un carácter cuantitativo en un medio ambiente determinado que es debida a la variación genética. Se considera que una heredabilidad h2 > 0,3 es elevada. Juan Daniel García Villabrille 108 en cuenta que la medición de esta variable (ecn) no resulta tan sencilla como la obtención de los diámetros o las alturas. En la Tabla 6.18 se muestran los resultados del ajuste de los dos sistemas de ecuaciones de biomasa de árbol individual planteados inicialmente. Todos los parámetros fueron significativamente diferentes de cero al nivel de significación del 5%. El sistema completo de árbol individual fue el empleado para las posteriores estimaciones de todas las fracciones por separado a nivel de rodal, ya que los porcentajes de carbono descritos para la especie (Brañas et al., 2000a) están definidos para cada fracción y no por agrupaciones más simplistas, pudiendo realizar estimaciones más específicas de carbono para cada fracción a nivel de rodal. Además, los valores de los estadísticos de ajuste indican que el modelo completo puede ser empleado sin incrementar el error en biomasa total de forma importante (REMC = 9,74 kg frente a REMC = 9,30 kg del modelo reducido). A su vez, el valor de REMC para la biomasa total de elementos lignificados, excluyendo las hojas, fue de 9 kg. La comparación de estas ecuaciones con los estadísticos del ajuste de sistemas aditivos (Tabla 6.19) indica que la mejora en cada fracción al emplear modelos más complejos puede ser muy importante, si bien no se asegura la minimización del error para la biomasa total. Aunque es visible el aumento en la variabilidad explicada por las ecuaciones independientes respecto del sistema aditivo completo, en el caso de la biomasa de madera, donde se considera más importante la estimación, las diferencias son menores del 5%. Tabla 6.18. Resultados del ajuste de sistemas aditivos de biomasa de árbol individual. Fact. correc. Sesgo Parámetros Ajuste Modelo Fracción K b0 b1 b2 R2 REMC Sistema completo Hojas 1,332 exp (-2,616) (<0,001) 1,356 (<0,001) 0,5992 2,14 Ramillos 1,163 exp (-3,375) (<0,001) 1,379 (<0,001) 0,6147 0,87 R. gruesas + finas 1,139 exp (-3,433) (<0,001) 1,928 (<0,001) 0,6946 2,55 Madera 1,013 exp (-4,365) (<0,001) 1,801 (<0,001) 1,244 (<0,001) 0,9851 8,06 Corteza 1,028 exp (-4,977) (<0,001) 2,037 (<0,001) 0,4996 (<0,001) 0,9241 2,20 Total 0,9856 9,74 Sistema reducido Hojas 2,398 (<0,001) -0,03099 (<0,001) -2,032 (<0,001) 0,6893 1,89 Ramas 1,097 exp (-2,757) (<0,001) 1,75 (<0,001) 0,7944 2,64 Madera 1,013 exp (-4,302) (<0,001) 1,899 (<0,001) 1,126 (<0,001) 0,9849 8,09 Corteza 1,028 exp (-4,879) (<0,001) 2,081 (<0,001) 0,4225 (0,004) 0,9223 2,23 Total 0,9869 9,30 El valor entre paréntesis indica el nivel de significación. Funciones de salida 109 Tabla 6.19. Mejoras en porcentaje para los ajustes de biomasa individual con ecuaciones independientes. Fracción Hojas Ramillos Ramas Madera Corteza Total Mejora REMC (%) 11,68 21,83 15,68 4,30 40,00 11,29 Atendiendo a las ramas gruesas, Ruiz-Peinado et al. (2011) hicieron distintas agrupaciones y/o restricciones para esta fracción en diferentes especies de frondosas, siendo incluida en una fracción junto con el fuste para E. globulus, manteniendo niveles de MEF = 0,97. Pérez-Cruzado y Rodríguez-Soalleiro (2011) testaron además de la longitud de copa (cl) otras variables de copa en sistemas de ecuaciones aditivas para E. nitens en Galicia, y no se justificaron las mejoras para las fracciones menores ya que se incrementaron los errores en las fracciones de mayor representación e importancia (biomasa de madera y biomasa total). Álvarez-González et al. (2005a) emplearon variables de rodal para ajustar un sistema aditivo en E. globulus, pero sólo se presentó mejor ajuste para la fracción de ramas finas (R2 = 0,80), aunque en esta tesis se agruparon con la fracción de ramas gruesas. En un ajuste previo por separado se obtuvo el mismo valor (R2 = 0,80), empleando como variables explicativas dn y cl, pero el valor de REMC se pudo reducir hasta un 60% al considerarlas como una única fracción. Los resultados mejoran los errores obtenidos por Brañas et al. (2000a) para todas las fracciones, pero la diferencia respecto a los errores obtenidos por Álvarez-González et al. (2005a) es mucho mayor. Respecto a estos trabajos, el aumento en la intensidad del muestreo destructivo realizado en esta ocasión supuso la principal mejora, aunque se debe tener en cuenta que la variabilidad explicada por el modelo en cada fracción no es muy diferente. En un trabajo de biomasa individual para 13 especies de eucaliptos en Australia, Bi et al. (2004) desarrollaron sistemas aditivos, haciendo una reducción de la fracción de ramas similar a la que se presenta en el modelo reducido en este estudio. Para nueve de esas especies se obtuvieron coeficientes de determinación R2 > 0,9 en la fracción conjunta de ramas. Otra opción de agrupación de las fracciones es la que se presenta en Ruiz-Peinado et al. (2012), agregando la fracción de ramas gruesas al fuste, ajustando las ramas finas por separado y agregando las hojas a los ramillos, consiguiendo valores de R2 de 0,97, 0,86 y 0,83 respectivamente, si bien los errores fueron del orden de los obtenidos por Álvarez-González et al. (2005a). En la Figura 6.8 se muestran los gráficos de residuos frente a valores predichos para el sistema completo de biomasa de árbol individual. Los residuos no muestran tendencias, por lo que no se consideró necesario el empleo de regresión ponderada. Juan Daniel García Villabrille 110 Figura 6.8. Gráficos de residuos frente a valores predichos en forma logarítmica para el sistema de biomasa de árbol individual, W = biomasa en kg/pie. 6.3.5. Ecuaciones de biomasa de rodal En la Tabla 6.20 se muestran los resultados del ajuste del sistema de biomasa de rodal. Todos los parámetros fueron significativos al nivel de confianza del 5%, obteniendo coeficientes de determinación superiores a 0,94 en todas las fracciones, siendo buenos resultados para las fracciones menores (hojas y ramas de todos los tamaños), con mayor variabilidad. En la Figura 6.9 se muestra el gráfico de valores observados frente a valores predichos para el sistema de ecuaciones de biomasa de rodal. A pesar de las mejoras en el ajuste de las fracciones terminales (hojas y ramillos) en comparación con el ajuste de árbol individual, se sigue detectando una gran dispersión de las predicciones en torno a la recta 1:1. Si bien esto es reflejo de la gran variabilidad de estas fracciones, el sesgo a nivel de rodal para todas las predicciones fue prácticamente inexistente. Las fracciones con mayor peso e importancia (madera y corteza) presentaron una baja dispersión en torno a la recta 1:1, pero sin embargo se detectó subestimación para valores elevados, siendo una consecuencia de menor importancia que la sobrestimación (Satoo y Madgwick, 1982, p. 24). -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 -4 -2 0 2 4 Log (Whojas) Residuos estandarizados -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 -2 -1 0 1 2 Log (Wramillos) -1 0 1 2 3 -2 -1 0 1 2 Log (Wramas) -2 -1 0 1 2 3 -2 -1 0 1 2 Log ( W corteza ) Predichos Residuos estandarizados -1 0 1 2 3 4 5 -2 -1 0 1 2 Log ( W madera ) Predichos 0 1 2 3 4 5 -3 -2 -1 0 1 2 3 Log ( W total ) Predichos Funciones de salida 111 Tabla 6.20. Resultados del ajuste de los sistemas de ecuaciones aditivas para la estimación de biomasa de rodal. Factor correcc. sesgo Parámetros Ajuste Modelo Fracción K b0 b1 b2 R2 REMC Sistema completo Rodal Hojas 1,003 exp (6,288) (<0,001) 0,6724 (<0,001) 0,9480 284,4 Ramillos 1,003 exp (5,431) (<0,001) 0,6809 (<0,001) 0,9529 119,4 R. gruesas + finas 1,0001 exp (6,068) (<0,001) 0,9581 (<0,001) 0,9995 63,71 Madera 1,001 exp (5,060) (<0,001) 1,001 (<0,001) 1,018 (<0,001) 0,9947 3589 Corteza 1,001 exp (4,607) (<0,001) 1,012 (<0,001) 0,4641 (<0,001) 0,9972 240,0 Total 0,9955 3942 El valor entre paréntesis indica el nivel de significación. Figura 6.9. Gráfico de valores observados frente a valores predichos para los modelos de biomasa de rodal (t/ha). En referencia a modelos de rodal, Tomé et al. (2001) realizaron un único ajuste de biomasa para el modelo GLOBULUS 2.1, que fue una ecuación de biomasa total. Para el resto de fracciones ajustaron ecuaciones de proporción para el fuste en la biomasa total, proporción de la corteza en el fuste y proporción de las hojas en la copa. De esa manera, la biomasa de madera y de ramas la obtuvieron por simple diferencia entre salidas de las 0 1 2 3 4 5 6 0 1 2 3 4 5 6 Whojas Observados 0.0 0.5 1.0 1.5 2.0 2.5 3.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Wramillos 0 5 10 15 0 5 10 15 W ramas 0 5 10 15 20 0 5 10 15 20 W corteza Predichos Observados 0 50 100 150 200 250 0 50 100 150 200 250 W madera Predichos 0 50 100 150 200 250 300 0 50 100 150 200 250 300 W total Predichos Juan Daniel García Villabrille 112 ecuaciones ajustadas, obteniendo un error medio cuarático o EMC = 121,01 (REMC = 11), lo que significa casi triplicar el error del sistema ajustado en este trabajo. Posteriormente, Tomé et al. (2006) ajustaron un sistema aditivo de rodal para el modelo GLOBULUS 3.0, como mejora frente al sistema anterior, si bien no indican sus estadísticos de ajuste. De este modelo se seleccionó la ecuación de biomasa de madera, y ésta fue aplicada sobre los datos de partida de esta tesis, obteniendo valores de R2 = 0,97, REMC = 8,6 t/ha y E = 5,6 t/ha, considerándose aceptables. Se debe tener en cuenta que este modelo de biomasa emplea todas las variables de estado del rodal (H0, N, G, t e IS), estando definida tref en 10 años. Por tanto, se cambió el dato IS obtenido con el modelo desarrollado por el resultante de aplicar la clasificación del modelo GLOBULUS 3.0. Para utilizar en Galicia este modelo portugués sería necesario que existiesen datos fiables de estaciones meteorológicas cercanas a las parcelas, pero como se comprobó durante la recopilación de las variables climáticas para el análisis inicial, no se asegura la cercanía. Por su parte, el uso a través de internet de la aplicación pública elaborada en Portugal con el modelo GLOBULUS 3.0 revela una gran variabilidad en los días de lluvia a considerar en cada caso, aun cuando las parcelas estuviesen muy próximas (http://home.isa.utl.pt/~joaopalma/modelos/webglobulus/). De esta forma se resta fiabilidad para su empleo a partir de los datos menos específicos que se encontraron para el caso de Galicia a partir de las fuentes disponibles: http://www2.meteogalicia.es/galego/observacion/estacions/estacions.asp. Finalmente se hizo una estimación para la muestra de ajuste empleando un dato fijo de 180 días, obteniendo R2 = 0,95, REMC = 10,7 t/ha y E= 7,6 t/ha, lo que empeoró las estimaciones iniciales, siendo mejor el empleo de un IS estimado con el modelo de altura dominante ajustado. Esta prueba se consideró importante, ya que el modelo GLOBULUS 3.0 podría ser una referencia para utilizar en el caso de no existir ecuaciones de rodal, indicando que el desarrollo de modelos para Galicia era necesario para obtener estimaciones más aproximadas, además de simplificar el número de variables independientes en los modelos. Por otra parte, la transformación de biomasa a carbono se puede realizar mediante los porcentajes de carbono descritos en Brañas et al. (2000a), establecidos a partir de muestras tomadas en Galicia, para cada una de las fracciones en que se han dividido los sistemas denominados completos, tanto de árbol individual como de rodal. En el caso de emplear el modelo reducido se deberá establecer un valor único para toda la fracción de ramas a partir de los valores de cada uno sus componentes, aunque también se podrían emplear los porcentajes descritos en Tomé et al. (2006), cuya división del sistema de biomasa es igual a la del sistema reducido de esta tesis. Los porcentajes de las fracciones similares de uno y otro trabajo no son iguales, pero tampoco existe una razón para descartar ninguno de ellos. Los porcentajes de carbono según Brañas et al. (2000) son del 52% en hojas, 46,4% en ramillos, 45,3% en ramas, 42,5% en corteza y del 45,2% en madera. Por su parte, los porcentajes de carbono descritos en Tomé et al. (2006) son del Evaluación del modelo 119 7. Evaluación del modelo Este capítulo se ha divido en tres partes. Primero se realiza una evaluación global del modelo, analizando los errores de las funciones de transición y de salida a partir de datos reales. La segunda parte se centra en la comparación del modelo de calidad de estación y las predicciones de volumen y biomasa con otros modelos disponibles para la especie utilizando los mismos datos sobre los que se evalúa el modelo ajustado. La tercera parte se centra en la simulación de diferentes alternativas respecto a densidades de plantación y calidades, para la predicción de los turnos de máxima renta en volumen y las producciones correspondientes. 7.1. Evaluación global del modelo ajustado 7.1.1. Introducción La validación es una manera de ayudar a decidir al usuario si un determinado modelo es o no adecuado para su uso. La forma ideal de realizar la validación de un modelo es a partir de datos independientes, pero no siempre es posible conseguirlos. Como ya se comentó, y aunque se realizaron dos pruebas en aprovechamientos reales descritas en el capítulo dedicado a las funciones de salida, no fue posible la validación sobre muestras independientes. Según algunos autores (Huang, 2002; Kozak y Kozak, 2003), si los datos usados en la parametrización han sido tomados según un diseño de muestreo adecuado y si el modelo se ha construido sobre una base teórica con significado biológico, la calidad de los ajustes probablemente reflejará la calidad de las predicciones dentro de la misma población. De forma inicial, los modelos de esta tesis se han elaborado a partir de los datos de una red de parcelas que cubren la combinación de diferentes edades, densidades, calidades, materiales de reproducción y tratamientos selvícolas de las plantaciones de Eucalyptus globulus en Galicia, por lo que la realidad actual se considera bastante bien representada, salvo por la falta de datos de crecimiento correspondientes a edades superiores a los 13 años en clones F0 y a los 10 años en clones F1. Por otra parte, se cuenta con los estadísticos de ajuste de cada uno de los modelos sobre la muestra de ajuste que ya se mostraron en cada capítulo, y que a priori indican buenos resultados en todos los módulos desarrollados. Juan Daniel García Villabrille 120 Así, se ha realizado una evaluación de los modelos ajustados por medio de las variables de estado del rodal y las variables de salida de interés para las masas de eucalipto, es decir, el volumen y la biomasa del rodal, sobre una muestra aleatoria de 270 intervalos de crecimiento. Según Huang (2002), basándose en la experiencia y en el tratamiento de diferentes modelos de crecimiento y producción, pueden utilizarse los siguientes umbrales de aceptación o rechazo cuando se necesitan usar estadísticos cuantitativos en la toma de decisiones: − Si el error medio de la predicción está dentro de ± 10% del valor de la media observada al 95% de nivel de significación, el modelo resulta aceptable, siempre y cuando la validación gráfica muestre tendencias razonables. − Si el error medio de la predicción es mayor que ± 10%, pero menor que ± 20% del valor de la media observada, la aceptación o rechazo del modelo depende de las circunstancias específicas. Este es, por tanto, un intervalo de indecisión, requiriéndose nuevos datos para poder dar una respuesta definitiva. − Si el error medio de la predicción es mayor que ± 20% del valor de la media observada, el modelo puede ser rechazado. Sin embargo, aún cuando se puede decir que los modelos que producen más de ± 10% de error son menos deseables, el nivel de aceptación o rechazo es algo subjetivo. 7.1.2. Estadísticos de evaluación Además de los valores de los estadísticos obtenidos en los ajustes para cada una de las funciones, y que se mostraron en los capítulos precedentes, se emplearon los siguientes dos estadísticos: por un lado, la raíz del error medio cuadrático (REMC), cuya expresión se recoge en el apartado 3.5.1 (Estadísticos de ajuste), y por otro lado el error crítico (Ecrit, expresado como un porcentaje de la media), reordenando el estadístico de Freese (1960): ( ) y yy E crit n iii crit 2 1 2 2ˆ χτ ∑ = −⋅ = [7.1] donde n es el número total de observaciones del conjunto de datos utilizados, i y es el valor observado, i y ˆes su predicción con el modelo ajustado, y es la media de los valores observados, τ es el valor de distribución normal estándar para un determinado nivel de probabilidad ( τ =1,960 para α = 0,05 ), y 2 crit χ se obtiene para α = 0,05 y n grados de libertad. Si el error máximo especificado expresado como un porcentaje de la Evaluación del modelo 121 media observada está dentro del límite del error crítico, el test 2 n χ indicará que el modelo produce estimaciones no satisfactorias; en caso contrario, indicará que las predicciones son aceptables. 7.1.3. Resultados En la Tabla 7.1 se presentan los resultados de la evaluación de los modelos ajustados para las funciones de transición del modelo. Para la altura dominante, la REMC y el Ecrit correspondientes fueron 1,04 m y 9,67% respectivamente. Respecto al área basimétrica, el modelo no descendió del umbral de Ecrit = 10%, obteniendo una REMC = 1,2 m2/ha y un Ecrit = 12,58%. El modelo de mortalidad fue el que menos error cometió de las tres funciones de transición ajustadas, con una REMC = 42 pies/ha y un Ecrit = 7,21%. Por tanto, según Huang (2002), los valores obtenidos para altura dominante y mortalidad son aceptables para su uso práctico, mientras el modelo de área basimétrica se encuentra en el intervalo de indecisión (Ecrit > 10% pero < 20%), aunque más cerca del límite de aceptación que del de rechazo. Estos resultados quedaron corroborados por las pruebas F realizadas, en las que se rechaza la existencia de tendencias sistemáticas en las tres funciones de transición (Figura 7.1), por lo que pueden considerarse aceptables los tres modelos. Tabla 7.1. Estadísticos obtenidos en la validación del modelo para las variables de rodal altura dominante (H0, m), área basimétrica (G, m2/ha) y número de pies por hectárea (N, pies/ha). Altura dominante (m) Área basimétrica (m2/ha) Mortalidad (pies/ha) REMC Ecrit REMC Ecrit REMC Ecrit 1,04 9,67 1,16 12,58 42 7,21 Figura 7.1. Gráfico de valores observados frente a valores predichos para altura dominante (izquierda), área basimétrica (centro) y mortalidad (derecha) en la validación del modelo. Línea continua: ajuste; línea punteada: recta 1:1. 5 10 15 20 25 30 35 40 5 10 15 20 25 30 35 40 Predichos Observados 0.4342 + 0.9786 x R2= 0.971 F valor 0.958 , Prob.> F 0.723 0 10 20 30 40 0 10 20 30 40 Predichos 0.4353 + 0.9791 x R2= 0.978 F valor 0.959 , Prob.> F 0.73 600 800 1000 1200 1400 1600 1800 600 800 1000 1200 1400 1600 1800 Predichos -27.9881 + 1.0242 x R2= 0.971 F valor 1.049 , Prob.> F 0.696 Juan Daniel García Villabrille 122 En la Tabla 7.2 se muestran los valores de la REMC y el Ecrit para volumen de rodal con y sin corteza en la validación del modelo. Los valores de REMC obtenidos fueron 11,58 m3/ha para Vcc y 9,97 m3/ha para Vsc. Si bien el Ecrit se situó por encima del 10% para ambos modelos, la prueba F rechazó la existencia de tendencias sistemáticas (Figura 7.2), siendo aptos para su empleo. Tabla 7.2. Estadísticos obtenidos en la validación del modelo para las variables de salida volumen con corteza (Vcc, m3/ha) y volumen sin corteza (Vsc, m3/ha). Vcc (m3/ha) Vsc (m3/ha) REMC Ecrit REMC Ecrit 11,58 14,92 9,97 15,15 Figura 7.2. Gráfico de valores observados frente a valores predichos para volumen de rodal con corteza (izquierda) y volumen de rodal sin corteza (derecha) en la validación del modelo. Línea continua: ajuste; línea punteada: recta 1:1. En la Tabla 7.3 se muestran los valores de los estadísticos obtenidos para biomasa total (Wtot, t/ha) y biomasa de madera de rodal (Wma, t/ha) en la validación del modelo. Dado el alto porcentaje que supone la biomasa de madera frente a la biomasa total, en ambos casos se obtuvieron estadísticos muy similares, como es lógico, sin que la variabilidad de las fracciones menores afecte sensiblemente al modelo. Así, la REMC fue de 6,8 t/ha para Wtot y 5,8 t/ha para Wma.. Tabla 7.3. Estadísticos obtenidos en la validación del modelo para las variables de salida biomasa total (Wtot, t/ha) y biomasa de madera (Wma, t/ha). Wtot (t/ha) Wmad (t/ha) REMC Ecrit REMC Ecrit 6,8 14,99 5,79 16,43 0 100 200 300 400 500 600 0 100 200 300 400 500 600 Predichos Observados 2.0429 + 0.9935 x R2= 0.988 F valor 0.987 , Prob.> F 0.915 0 100 200 300 400 500 0 100 200 300 400 500 Predichos 1.5143 + 0.9954 x R2= 0.988 F valor 0.991 , Prob.> F 0.94 Evaluación del modelo 123 En la Figura 7.3 se muestran las rectas de regresión entre valores reales y predichos para biomasa de madera y biomasa total de rodal (t/ha) en la validación del modelo. Los valores del coeficiente de determinación resultantes fueron altos, R2 = 0,986, tanto en biomasa de madera como total. Aunque el Ecrit fue superior al 10% en ambos casos, la prueba F indicó la inexistencia de sesgos, por lo que los modelos son válidos para la predicción de biomasa. Figura 7.3. Gráfico de valores observados frente a predichos para biomasa de rodal (t/ha) de madera (izquierda) y biomasa total (derecha) en la validación del modelo. Línea continua: ajuste; línea punteada: recta 1:1. Las predicciones de diámetro medio aritmético (d) se emplearon para recuperar los parámetros de la función de Weibull por el método de los momentos, y se hizo una prueba bajo un supuesto posible, en el que se establecieron dos destinos principales de la masa: energía y celulosa. Las clases diamétricas límite consideradas fueron de 10 cm para energía, y entre 10 y 40 cm para celulosa. Como se vio en el aparado correspondiente a la función de desagregación (6.3.6), el ajuste de distribuciones no se considera el más deseable debido al alto porcentaje de rechazos, por lo que se recomienda emplear este procedimiento sólo a nivel informativo, pero no como base de decisión. En la Tabla 7.4 se muestran los resultados de la validación para el modelo de d por materiales de reproducción. Los errores se situaron en el entorno de 1 cm o menos, siendo los brinzales los que obtuvieron peores resultados frente a las variedades clonales, con un Ecrit que se sitúa en la barrera del 10%, mientras el clon F0 arrojó un valor de Ecrit = 3,98% y el clon F1 un valor de Ecrit = 5,35%. En la Figura 7.4 se muestra la recta de ajuste entre los valores reales y predichos, siendo clara la subestimación, aunque el R2 fue de 0,946 y la prueba F rechazó la existencia de tendencias como también ocurrió para los tres materiales de reproducción por separado (F = 1,03, Prob. > F = 0,78). Así, para brinzal se obtuvo un R2 = 0,938 (F = 1,027, Prob. > F = 0,87), para clon F0 un valor R2 = 0,954 (F = 0,914, Prob. > F = 0,714) y para clon F1 un R2 = 0,9 (F = 0,853, Prob. > F = 0,579). 0 50 100 150 200 250 300 0 50 100 150 200 250 300 Predichos Observados 0.5245 + 1.0088 x R 2 = 0.986 F valor 1.018 , Prob.> F 0.886 0 50 100 150 200 250 300 350 0 50 100 150 200 250 300 350 Predichos 0.7309 + 1.0068 x R2= 0.986 F valor 1.014 , Prob.> F 0.912 Juan Daniel García Villabrille 124 Tabla 7.4. Estadísticos obtenidos en la validación del modelo para la variable de salida diámetro medio aritmético para la muestra conjunta y por materiales de reproducción. Muestra conjunta Brinzal Clon F0 Clon F1 REMC Ecrit REMC Ecrit REMC Ecrit REMC Ecrit 0,9 12,35 1,1 10,04 0,55 3,98 0,76 5,35 Figura 7.4. Gráfico de valores observados frente a valores predichos en la validación del modelo para diámetro medio aritmético (d, cm) empleando H0 y N del modelo ajustado para predicción de d, y los parámetros c y b de la función de Weibull biparamétrica recuperados por el método de los momentos. Línea continua: ajuste; línea punteada: recta 1:1. En la Tabla 7.5 se muestran los resultados de la validación del modelo para la predicción del número de pies en función del destino posible (energía o celulosa), aunque se debe tener en cuenta que las masas estudiadas están principalmente destinadas a celulosa, y también el diseño del muestreo. Los menores valores del Ecrit para la predicción de pies destinados a energía correspondieron al clon F1, si bien las parcelas pertenecen a las edades más tempranas de media, y algunas poseen las mayores densidades de todas las muestreadas, por lo que se asemejan más a una plantación energética que los brinzales, por ejemplo. El clon F0 es el único material de reproducción para el que se consiguió descender del 20% (Ecrit = 17,32%) y sólo en el caso de predicción de pies destinados a celulosa, aunque los brinzales también se situaron en un valor límite con Ecrit = 20,29%. Por tanto, no se recomienda la evaluación para diferentes destinos a partir de los modelos ajustados, si bien éstos mostraron buena capacidad de predicción para el diámetro medio aritmético. 5 10 15 20 25 5 10 15 20 25 Predichos Observados 0.2808 + 1.0166 x R2= 0.946 F valor 1.033 , Prob.> F 0.787 Evaluación del modelo 125 Tabla 7.5. Estadísticos obtenidos en la validación del modelo para la predicción del número de pies en función del destino posible (energía o celulosa), para la muestra conjunta y por materiales de reproducción. Muestra conjunta Brinzal Clon F0 Clon F1 Predicción REMC Ecrit REMC Ecrit REMC Ecrit REMC Ecrit Energía (pies/ha) 152 61,35 112 45,18 135 28,77 257 23,47 Celulosa (pies/ha) 144 42,7 100 20,29 129 17,32 251 51,89 7.2. Comparaciones con otros modelos Una vez evaluado el comportamiento del modelo, se realizaron diferentes comparativas con el modelo dinámico disponible hasta la fecha para Eucalyptus globulus en Galicia (García y Ruiz, 2003), así como con los modelos de Tomé et al. (2006) y de Wang y Baker (2007). Como se indicó en la introducción, obtener un modelo de calidad de estación que representase adecuadamente las masas de primer turno en Galicia era una de las principales necesidades. Así, se abordó la comparación gráfica del modelo de calidad de estación ajustado frente a los modelos comentados. Posteriormente se realizó una evaluación de las funciones de transición de esos otros modelos sobre la muestra empleada para la validación en esta tesis. Dado que la muestra empleada es independiente de la que sirvió para el ajuste del modelo elaborado en la presente tesis, este proceso, al igual que ocurrió con la relación altura-diámetro elaborada por Crecente-Campo et al. (2010), puede ser considerado una validación de dichos modelos, teniendo en cuenta que dos de ellos no fueron desarrollados para Galicia (Tomé et al., 2006 y Wang y Baker, 2007). 7.2.1. Modelo de calidad de estación La evaluación de los modelos se hizo sobre la misma muestra aleatoria de 270 intervalos de crecimiento que se emplearon para la evaluación global del modelo (apartado 7.1. Evaluación global del modelo). La obtención de las curvas de calidad de estación a partir del modelo de Tomé et al. (2006) requirió la recopilación de datos referentes a la cantidad de días de lluvia con precipitación superior a 0,1 mm. Para ello se emplearon los datos obtenidos de www.meteogalicia.es referentes a los últimos 8 años, limitación debida a la falta de datos anteriores a esas fechas. En la Tabla 7.6 se muestra el número de días promedio anuales con precipitación superior a 0,1 mm empleados para cada provincia. Juan Daniel García Villabrille 126 Tabla 7.6. Días anuales con precipitación > 0,1 mm por provincias, indicando la estación meteorológica A Coruña (Santiago EOAS) Lugo (Burela) Pontevedra (Mondariz) 158 173 140 En la Tabla 7.7 se muestran el error (REMC) de los referidos modelos. A simple vista, el modelo de Tomé et al. (2006) produce errores similares al modelo ajustado, y tanto este modelo como el de Wang y Baker (2007) parecen más apropiados que el modelo elaborado por García y Ruiz (2003), que produjo las peores estimaciones, además de no poder predecir en siete de los casos por ser la altura al inicio del intervalo mayor de 33,8 m (asíntota del modelo para plantaciones). Tabla 7.7. Estadísticos de ajuste sobre los datos de partida de diferentes modelos de calidad disponibles para E. globulus (García y Ruiz, 2003; Tomé et al., 2006; Wang y Baker,2007). Modelo REMC Modelo ajustado 1,04 García y Ruiz (2003) 1,27 Tomé et al, (2006) 1,03 Wang y Baker (2007) 1,16 El modelo de Wang y Baker (2007) produjo mejores estimaciones que el modelo de García y Ruiz (2003), disminuyendo 11 cm la REMC frente a éste, aunque el modelo de Tomé et al. (2006) fue el mejor para la muestra seleccionada (REMC = 1 m). Respecto a los valores máximos, Wang y Baker (2007), empleando el modelo de Bertalanffy- Richards resuelto por a2, fijaron una asíntota con un valor de 50 m. Tomé et al. (2006) que ajustaron el modelo de Korf resuelto por a2, expandieron el parámetro que fija la asíntota en función de los días de precipitación, permitiendo así que un modelo que genera curvas polimórficas de asíntota común se comporte como un modelo de curvas polimórficas con diferentes asíntotas, entre un mínimo de 29,07 m (término independiente del modelo que rige el parámetro de la asíntota) y un máximo, para estaciones con 173 días de lluvia, por ejemplo, de 78,9. Por otra parte, el modelo aquí ajustado proporciona asíntotas de 69,5 m y 26,3 para IS = 25 e IS = 7 respectivamente, siendo similares a los otros modelos. Como referencia a máximos registrados de los que se tiene constancia (Souto da Retorta, Chavín, provincia de Lugo), entre 55 y 65 m y existiendo masas de 40 m a los 20 años incluidas en la muestra de ajuste, un máximo valor posible de 78,9 m no parece fácilmente alcanzable, siendo más aproximada la asíntota del modelo ajustado, si bien esto es poco relevante en masas destinadas a celulosa. Respecto al empleo del modelo de Tomé et al. (2006), en el caso de estaciones con calidades diferentes (pero próximas entre sí), la de menor calidad pero con más días de Evaluación del modelo 127 lluvia podría alcanzar mayores alturas. Otro caso que genera dudas es el de dos estaciones en una región climáticamente homogénea, en referencia a épocas de lluvias y de sequías, en las que habiendo la misma precipitación anual (l/m2), ésta se concentre en diferente número de días. Tampoco resulta sencillo el empleo de este modelo en Galicia, ya que la posibilidad de obtener el número de días de lluvia en cualquier punto para un punto concreto es difícil, al menos con la sensibilidad que muestra el modelo de Tomé et al. (2006), hecho que se puede comprobar a través de su aplicación web. Respecto al tipo de curvas generadas por los modelos de García y Ruiz (2003) y Wang y Baker (2007), con asíntota común, no se consideran adecuadas, al menos para representar la situación de las masas gallegas donde la falta de selvicultura a partir de edades avanzadas desemboca en estados de decaimiento, acompañados de ataques por goniptero, presente en toda la región, lo que limita el crecimiento en altura. En la Figura 7.5 se muestra una comparación de las curvas de calidad de estación ajustadas en esta tesis con las correspondientes a los modelos de García y Ruiz (2003), Tomé et al. (2006) y Wang y Baker (2007), con IS = 7, 13, 19 y 25 m en tref = 7 años. Como se puede comprobar en dicha figura, existen diferencias notorias en el crecimiento de las calidades más bajas entre los modelos para Galicia y los de Portugal (Tomé et al., 2006) y Australia (Wang y Baker, 2007), aun cuando se seleccionó el modelo que tenía mayor tasa de crecimiento en dichas calidades (apartado 5.3.1, Resultados de altura dominante). Esto pone de manifiesto la necesidad de recopilar más datos de parcelas de baja calidad, si bien en ocasiones las plantaciones que no muestran aptitud para el crecimiento son sustituidas antes de llegar al turno, teniendo constancia de cambios por E. nitens y E. regnans (éste en un único caso). No obstante, el modelo de García y Ruiz (2003) para plantaciones también predice el mismo crecimiento en altura dominante en calidades bajas, que a su vez es inferior al de las curvas para el grueso de su estudio, que fueron los chirpiales. Por otra parte, en la calidad máxima seleccionada (IS = 25 m) y para edades avanzadas (> 15 años), el modelo de Tomé et al. (2006) es muy similar al modelo ajustado, mientras el modelo de Wang y Baker (2007) presenta mayor tasa de crecimiento. Juan Daniel García Villabrille 128 Figura 7.5. Curvas de altura dominante con IS = 7, 13, 19 y 25 m en tref = 7 años frente a la edad para el modelo de calidad de estación ajustado en esta tesis, superpuesto a las curvas predichas por los modelos de García y Ruiz (2003), Tomé et al. (2006) y Wang y Baker (2007) y los datos observados,. 7.2.2. Evaluación de otros modelos En la Tabla 7.8 se presentan los resultados de la validación de las funciones de transición de los otros modelos existentes para la especie, utilizando la misma muestra de datos empleada para la evaluación del modelo desarrollado en esta tesis. El modelo de Tomé et al. (2006), con REMC = 1,03 m y Ecrit = 9,52%, parece igualmente válido para la estimación de crecimiento en altura dominante que el modelo ajustado. Por otra parte, el modelo de García y Ruiz (2003) se mostró como el menos apropiado, siendo su REMC el más elevado entre los modelos testados (REMC = 1,27 m). Respecto al área basimétrica, ninguno de los modelos estuvo por debajo del Ecrit = 10%, siendo el mejor el de Wang y Baker (2007), con REMC = 1,1 m2/ha y Ecrit = 11,84%. El modelo ajustado, con REMC = 1,16 m2/ha y Ecrit = 12,58% (Tabla 7.1), mejoró las estimaciones del modelo de García y Ruiz (2003) con amplia diferencia. El modelo de área basimétrica de Tomé et al. (2006), quizá debido al elevado número de variables intervinientes relacionadas con las características de la estación, dificulta las predicciones sobre muestras independientes fuera del rango geográfico. Tanto éste como el modelo de García y Ruiz (2003) fueron descartados para una fase posterior en la evaluación de salidas (volumen y biomasa), ya que la variable G interviene en ambas, dados los altos valores tanto de REMC como de Ecrit obtenidos. Edad (años) Altura dominante (m) 10 20 30 40 5 10 15 20 Modelo ajustado García y Ruiz (2003) Wang y Baker (2007) Tomé et al. (2006) Evaluación del modelo 135 Como se puede comprobar, para la calidad inferior se obtienen valores mayores de volumen con el modelo ajustado en esta tesis, mientras que en las otras dos, y especialmente en la mejor calidad, se obtienen volúmenes menores. También se observa que prácticamente no hay diferencias entre los volúmenes en la peor calidad si se comparan las opciones con y sin mortalidad, mientras que sí hay grandes diferencias en las otras dos calidades, siendo mayores las estimaciones que sí consideran mortalidad, aun a pesar de que en las tablas de Fernández-López no se consideró esa disminución del número de pies, como ya se ha comentado. Por otra parte, hay que volver a incidir en el hecho de que se ha utilizado en las estimaciones del modelo ajustado en esta tesis la función de inicialización de área basimétrica, con los comentados errores que eso supone. Aceptando que los valores aportados por Fernández-López (1985) son representativos de la realidad de su época, la reducción de la producción que se observa en las calidades altas utilizando el modelo ajustado en esta tesis podría ser debida a los problemas surgidos en la mayoría de las plantaciones desde los años 90 por enfermedades (Mycosphaerella spp.) y/o sobre todo por la plaga defoliadora (Gonipterus scutellatus), ya que pueden reducir de forma importante la producción (Merino et al., 2003). Aunque el goniptero ha sido tratado por control biológico con el parasitoide Anaphes nitens (Pérez et al., 2003), y en las entrevistas con productores se indicaron mejoras puntuales en el vigor asociadas a campañas sanitarias con este enemigo del gorgojo, se constató la presencia de daños por defoliación en la totalidad de las parcelas instaladas para esta tesis. Además, las continuas rotaciones de la especie, que suponen la extracción a través de la biomasa de una importante proporción de los nutrientes del sistema (Álvarez- González et al. 2005a), derivan en una menor disponibilidad de nutrientes para los sucesivos turnos. También las prácticas de acordonar los restos y los propios aprovechamientos suponen una elevada extracción de nutrientes, en suelos que comparativamente tienen pocos (Merino et al., 2003). Por otra parte, el empleo de planta procedente de mejora genética, o reproducción clonal, además de los programas de fertilización y una mejor evaluación del sitio antes de plantación pueden ser los causantes del aumento que se refleja en la comparativa con las producciones de la calidad más baja. En este sentido, y tras diversas entrevistas, se sabe que muchos propietarios llevan a cabo programas, más o menos sofisticados en su elaboración, que pueden incluir análisis de suelos, fertilización en el momento de la plantación (aunque no en todas las plantaciones, según Merino et al., 2003) y búsqueda de planta de calidad. También se realizan controles de la vegetación arbustiva por trituración, lo que facilita el aporte de nutrientes a la plantación, si bien no se aplican con la misma intensidad y/o frecuencia en todos los casos. Como se acaba de indicar, existen justificaciones para tratar de explicar las diferencias entre los dos modelos en las calidades extremas, pero sin embargo no se puede asegurar que tales diferencias sean reales, pues el modelo de esta tesis se ha ajustado con poca información en las calidades extremas. Sin embargo, los resultados Juan Daniel García Villabrille 136 para las calidades intermedias, donde mayor cantidad de datos se consiguió, son muy similares a los indicados por Fernández-López (1985). En cualquier caso, el empleo del modelo de esta tesis para la predicción de volúmenes en un momento determinado o en intervalos cortos, así como las pruebas realizadas en aprovechamientos reales, como se ha discutido en capítulos anteriores, indicaron que sus estimaciones son realistas. EucaTool®: aplicación web para estimar el crecimiento y la producción de las plantaciones de Eucalyptus globulus de Galicia 137 8. EucaTool®: aplicación web para estimar el crecimiento y la producción de las plantaciones de Eucalyptus globulus de Galicia Como se ha podido comprobar en los apartados anteriores, el modelo de crecimiento y producción para las plantaciones de Eucalyptus globulus de Galicia desarrollado en esta tesis está compuesto por diferentes ecuaciones, módulos o funciones (de transición y de salida), así como por las relaciones entre las mismas. De tal manera, la aplicación práctica del modelo resulta muy compleja para sus potenciales usuarios, que son los selvicultores y gestores forestales de Galicia y de las regiones próximas, así como los investigadores, docentes y estudiantes de las diferentes facultades y centros de investigación forestales. Para solucionar este importante inconveniente, y con el objeto de cumplir con el último e imprescindible paso de cualquier investigación (y muy especialmente de la realizada desde el sector público), que consiste en difundir y transferir los resultados obtenidos al sector de la sociedad interesado en ellos, lo que por otra parte da sentido a todo el esfuerzo y recursos dedicados a la investigación, se ha diseñado y desarrollado un simulador informático denominado EucaTool® para facilitar la aplicación del modelo elaborado en esta tesis doctoral. EucaTool® es una aplicación web o de software “en la nube”, de acceso e utilización gratuita, que ha sido desarrollada para estimar el crecimiento y la producción de plantaciones de semilla y clonales de eucalipto (Eucalyptus globulus Labill.) en Galicia. Para ello, EucaTool® implementa el modelo dinámico de crecimiento y producción desarrollado en esta tesis, válido para plantaciones clonales y no clonales de la especie en primera rotación de Galicia y que, como ya se ha mencionado, integra diferentes módulos o funciones de transición para la altura dominante (curvas de calidad de la estación), número de pies por hectárea (función de mortalidad) y área basimétrica, junto con ecuaciones o tarifas de cubicación, de biomasa y de carbono. Concretamente, la aplicación permite estimar de forma muy sencilla el volumen, la biomasa y el contenido de carbono de dichas plantaciones a partir de datos de árboles individuales, de clases diamétricas o de valores medios de rodales, así como predecir también su crecimiento y producción futura (indicando la edad óptima de corta por máxima renta en especie) a partir de la medición tres variables medias del rodal: número de árboles por hectárea, altura dominante y área basimétrica; y la edad correspondiente. Juan Daniel García Villabrille 138 EucaTool® permite generar informes de las simulaciones realizadas, que se exponen en vista preliminar para ser impresos. Con los informes también se genera un código QR para poder acceder a las simulaciones digitales desde cualquier lector QR. Puede accederse a EucaTool® desde cualquier dispositivo con conexión a internet desde la dirección: http://app.eucatool.com. Además, toda la información relativa a su uso está publicada en una web ligada al aplicativo: http://www.eucatool.com En el Anexo III se incluye el manual de uso de EucaTool®, en el que se detallan todas las utilidades de la aplicación así como la manera de utilizarla. Conclusiones 139 9. Conclusiones A continuación se exponen las principales conclusiones derivadas de los diferentes análisis llevados a cabo: − Se ha desarrollado un modelo dinámico de crecimiento y producción para plantaciones clonales en primer turno y de brinzales de Eucalyptus globulus destinadas a celulosa en Galicia. Los datos disponibles para el ajuste de las funciones de transición del modelo dinámico, aunque de buena calidad, son mejorables, ya que no se pudieron obtener datos de plantaciones clonales mayores de 11 y 13 años para clones F1 y F0 respectivamente, ni se cubren las calidades extremas en el intervalo completo desde 1 a 21 años; además, los períodos entre los inventarios fueron cortos (1 año). Por tanto, se recomienda la continuidad y ampliación de la red de parcelas, recogiendo la información de los tratamientos aplicados en cada caso. − Los mejores modelos para las funciones de altura dominante (e índice de sitio) y área basimétrica fueron funciones polimórficas con diferentes asíntotas (Korf GADA), con errores críticos del 9,7% y 12,6% respectivamente. El mejor modelo para la función de mortalidad fue la ecuación diferencial propuesta por Woollons (1998) expresada en forma de diferencias algebraicas (ADA), que genera curvas sigmoides con mayores tasas de mortalidad para mayores densidades iniciales, y que obtuvo un error crítico del 7,2%. − El rango de índices de sitio observados fue de 7-25 m a la edad de referencia de 7 años. El rango de áreas basimétricas observadas fue de 4-25 m2/ha a la edad de 7 años. La selección de la edad de referencia fue evaluada en ambos casos mediante el error relativo, resultando un valor del 3,8% para el índice de sitio y del 5,6% para el área basimétrica. − Las funciones de salida de volumen y biomasa proporcionaron estimaciones por debajo del límite de rechazo del 20% en proyecciones a uno, dos y/o tres años, considerándose las tarifas de cubicación muy buenas para la estimación de la producción a partir de datos de campo, como quedó probado en cubicaciones reales. El modelo para rodales en los que el volumen es superior a 450 m3/ha es Juan Daniel García Villabrille 140 muy útil, sobre todo si te tiene en cuenta que existe un elevado número de rodales con grandes volúmenes y en edades superiores a la del turno, inherente al tipo de propiedad existente en Galicia y al modelo de gestión que ésta conlleva. − Los diferentes análisis llevados a cabo para tratar de identificar diferencias en las funciones de crecimiento y de salida según los materiales de reproducción (brinzales y clones), tipos de suelo (granitos/esquistos) y variables climáticas (relativas a precipitación y heladas) indicaron que la mejor solución era un modelo único (o reducido) en todos los casos. − Los test realizados sobre determinadas características de la madera, a partir de muestras de trozas comparables, indicaron que la densidad básica (kg/m3) es mayor para los clones que para brinzales, lo que implicaría una mayor producción de pasta a partir de árboles clonales que de semilla a igualdad en las variables de estado del rodal. No obstante, los resultados de densidad básica no deben tomarse como concluyentes, ya que las muestras no se recogieron ni se trataron específicamente para este propósito. − El análisis de las partes del árbol indicó que la longitud de copa relativa (respecto de la altura total), la biomasa relativa de ramas de todos los tamaños y hojas (respecto de la biomasa total) y el espesor de corteza relativo (respecto del diámetro normal) son mayores en los clones que en los brinzales. Por tanto, cuando el ajuste de modelos por separado no fue significativo, los modelos reducidos que incluyeron las variables longitud de copa (para biomasa de hojas), o espesor de corteza (para biomasa de corteza y madera), que además supusieron mejoras en las estimaciones, permitirían identificar las diferencias entre brinzales y clones. La variable "defoliación", aunque subjetiva, mejoró los resultados del modelo de biomasa de hojas en comparación con un modelo que no la incluyó. − La recuperación de los parámetros de la función Weibull por el método de los momentos no proporcionó los resultados deseables a nivel de parcela y, por tanto, se deberían probar otras metodologías, aunque la predicción a nivel de muestras grandes parece ser adecuada. Asumiendo los altos errores críticos en la etapa de evaluación (de diferentes destinos comerciales), la predicción del diámetro medio aritmético indicó errores críticos muy bajos, con un valor del 3,9% para clones F0, y del 5,4% para clones F1. Para el caso de los brinzales, las predicciones situaron los errores en el valor límite del 10%. En el análisis de las distribuciones Conclusiones 141 diamétricas se encontraron diferencias, siendo las de clones F1 y brinzales similares y menos apuntadas que las de clones F0. − Para facilitar la aplicación práctica del modelo desarrollado en esta tesis doctoral se desarrolló una aplicación web llamada EucaTool®, de acceso y uso gratuitos a través de la dirección http://app.eucatool.com. Esta aplicación permite estimar el volumen, la biomasa y el carbono, así como el crecimiento y la producción futura de rodales de Eucalyptus globulus destinadas a celulosa, utilizando datos de árboles individuales, de clases diamétricas o variables de rodal. The main conclusions derived from the study are now presented in English: − A dynamic growth and yield model was developed for first rotation clonal and seedling Eucalyptus globulus stands growing in Galicia and destined for pulp production. Data available for this dynamic model, although of good quality, could be improved, as no data from clonal plantations beyond 11 and 13 years were available for clones F1 and F0 respectively, and the extreme site qualities for the entire age range (1 to 21 years) were not sampled. Moreover, the interval between inventories was short (1 year). New inventories of the established plots and new plots are recommended, in order to collect more information and treatments applied in each case. − The best models for both dominant height (and site index) and basal area were polymorphic functions with different asymptotes (Korf GADA), with critical errors of respectively 9.7% and 12.6%. The best model for the mortality function was that proposed by Woollons (1998), using the algebraic difference approach (ADA), which generates a sigmoidal curve model, which implies higher mortality rates for higher initial stand stocking, with a critical error of 7.2%. − The observed range of site quality was 7-25 m at a reference age of 7 years. The observed range of basal area classes was 4-25 m2/ha at the same age (7 years). Selection of the reference age was assessed by the relative error, yielding a value of 3.8% for dominant height, and 5.6% for basal area. − Output functions for volume and biomass over time (1 to 3 years) provided critical errors below the 20% cut-off, and were considered very good for Juan Daniel García Villabrille 142 predicting production rates from field data, as proven for real predictions. The stand volume equation for cases in which volume is over 450 m3/ha is very useful, particularly considering the large number of high volume stands with ages beyond the optimal rotation, which is inherent in the ownership model for Galicia. − Assessment of different growth and output functions depending on the reproductive material (seedlings/clones), soil (granite/schist) and climatic variables (related to rainfall and frosts) indicated that the best solution was a single (or reduced) model in all cases. − Tests for wood characteristics, using comparable samples of logs, indicated that the basic density (kg/m3) decreased in the order Anselmo (F0 clone) > Odiel (F1 clone) > seedlings, which implies greater pulp production from clonal than from seedling trees for a similar stand state vector. However, the basic density results should not be considered conclusive, as samples were not collected or specifically treated for this purpose. − Analysis of different tree parts indicated that the relative crown length (in relation to the total height), relative biomass of all sizes branches and leaves (in relation to total biomass) and relative bark thickness (in relation to diameter at breast height) are greater in clones than in seedlings. Therefore, after rejection of a complete model, a reduced model including crown length (for leaf biomass), or bark thickness (for bark and wood biomass) as variables proved that model improved and could help to identify the differences between seedlings and clones. The variable “defoliation”, although subjective, and no better than the crown length, improved estimates of leaf biomass relative to a model that did not include the variable. − The moment-based parameter recovery method applied to the Weibull function did not provide the desired results at plot level and therefore other approaches should be tried. However, prediction for large samples seem to be adequate. Despite the large critical errors obtained when evaluating estimates of tree number in each diameter class (for different commercial purposes), prediction of mean diameter indicated very low critical errors, with values of 3.9% for F0 clones and 5.4% for F1 clones. In the case of seedlings, critical error was close to Conclusiones 143 the limit of 10%. In addition, the diameter distributions of F0 clones were more peaked than for F1 clones and seedlings, which had similar distributions. − To enable practical application of the developed model, a web application called EucaTool®, which is cost free, accessible and usable at http://app.eucatool.com, was developed. This application enables prediction of volume, biomass and carbon content and of growth and future production of Eucalyptus globulus pulp stands, using data from individual trees, diameter classes or stand variables. Referencias 145 10. Referencias Alder, D. 1980. Estimación del volumen forestal y predicción del rendimiento. Estudios FAO Montes 22/2, Roma, 198 pp. Álvarez, J. G., Schröder, J., Rodríguez R., Ruiz, A. D. 2002. Modelling the effect of thinnings on the diameter distributions of even-aged Maritime pine stands. Forest Ecology and Management. 165: 57-65. Álvarez-González, J.G., Barrio Anta, M., Diéguez Aranda, U., Rojo Alboreca, A. 2004a. Metodología para la construcción de curvas de calidad de estación. Cuad. Soc. Esp. Cienc. For. 18: 303-309. Álvarez-González, J. G., Castedo-Dorado, F., Ruiz-González, A. D., López-Sánchez, C.A., von Gadow K. 2004b. A two-step mortality model for even-aged stands of Pinus radiata D. Don in Galicia (Northwestern Spain). Annals of Forest Science. 61: 439-448. doi:10.1051/forest:2004037. Álvarez-González, J. G., Balboa, M. A., Merino, A., Rodríguez-Soalleiro, R. 2005a. Estimación de la biomasa arbórea de Eucalyptus globulus y Pinus pinaster en Galicia. Recursos Rurais. 1 (1): 21-30. Álvarez González, J. G., Ruiz González, A. D., Rodríguez Soalleiro, R., Barrio Anta, M. 2005b. Ecoregional site index models for Pinus pinaster in Galicia (northwestern Spain). Annals of Forest Science. 62 (2): 115-127. Álvarez-González J. G., Zingg, A., Gadow, Kv. 2010. Estimating growth in beech forests: a study based on long term experiments in Switzerland. Annals of Forest Science. 67 (3): 307. doi:10.1051/forest/2009113. Amaro, A. 1997. Modelaçao do crescimento de povoamentos de Eucalyptus globulus Labill. em 1ª rotaçao. O património florestal de uma empresa. Ph. D. Thesis. Inv. Op. Eng. Sist. Instituto Superior Técnico, Universidade Técnica de Lisboa. (Inédito). Amaro, A., Tomé, M., Themido, I. 1997. ALBA, un modelo de crescimento em área basal para povoamentos de eucalipto. En: Actas del "II Congreso Forestal Español- Irati 97". Puertas, F. y Rivas, M. (Eds.). Pamplona, 23-27 junio 2001. Mesa IV, Tomo IV, pp. 39-44. Amaro, A., Reed, D., Soares, P. (Eds.) 2003. Modelling forest systems. CABI publishing. Amateis, R. L., Burkhart, H. E., Jiping, L. 1997. Modeling survival in juvenile and mature loblolly pine plantations. For. Ecol. Manage. 90: 51-58. António, N., Tomé, M., Tomé, J., Soares, P., Fontes, L. 2007. Effect of tree, stand, and site variables on the allometry of Eucalyptus globulus tree biomass. Canadian Journal of Forest Research, 37: 895-906.