# Modélisation de la croissance et de la production du liège en Tunisie

Tahar SGHAIER¹
Salah GARCHI¹
Thouraya AZIZI¹

¹ Institut national de recherches en génie rural, eaux et forêts (INRGREF)
Rue Hédi El Karray El Menzah IV
BP 10, 2080 Ariana
Tunisie

Auteur correspondant /
Corresponding author :

Tahar SGHAIER -
sghaier.tahar@iresa.agrinet.tn

<figure><img src="image_1.png" /><figcaption>**Photo 1.**
Les chênes-lièges, *Quercus suber* L., sont démasclés pour favoriser la formation du liège femelle et en faire par la suite la récolte de liège de bonne qualité.
Photo T. Sghaier.</figcaption></figure>

Doi : 10.19182/bft2020.346.a31805 - Droit d'auteur © 2020, Bois et Forêts des Tropiques - © Cirad - Date de soumission : 7 novembre 2019 ;
date d'acceptation : 15 juillet 2020 ; date de publication : 25 décembre 2020.

<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

cirad

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

Licence Creative Commons :
Attribution - 4.0 International.
Attribution-4.0 International (CC BY 4.0)

**Citer l'article / To cite the article**
Sghaier T., Garchi S., Azizi T., 2020. Modélisation de la croissance et de la production du liège en Tunisie. Bois et Forêts des Tropiques, 346 : 3-20.
Doi : https://doi.org/10.19182/bft2020.346.a31805
# RÉSUMÉ

## Modélisation de la croissance et de la production du liège en Tunisie

L'absence d'outils nécessaires pour l'estimation des productions, la gestion et l'aménagement des peuplements de chêne-liège (Quercus suber L.) en Tunisie constitue un handicap pour le développement et la conservation de la suberaie dans le pays. Les objectifs de ce travail étaient de préciser l'impact de l'intensité de démasclage des arbres de chêne-liège sur la croissance et la production du liège, d'élaborer des modèles d'estimation de la croissance et la biomasse du liège à l'échelle de l'arbre, et enfin de vérifier des coefficients de démasclage appliqués en Tunisie et leur impact sur la production nationale de liège. Les données collectées proviennent d'un essai multisites de démasclage et de deux inventaires réalisés dans différents chantiers de démasclage durant deux campagnes successives de récolte de liège. Trois sites expérimentaux représentatifs des forêts productives de chêne-liège en Tunisie ont fait l'objet d'un même essai où trois coefficients de démasclage ont été appliqués sur des arbres appartenant à trois classes de diamètre. Des mesures de l'épaisseur du liège à hauteur de poitrine ont été effectuées annuellement durant une rotation de 12 ans, à la fin de laquelle le liège produit sur chaque arbre a été récolté et pesé. Les résultats obtenus ont montré un effet « coefficients de démasclage » non significatif et un effet « classes de diamètre » très hautement significatif sur la croissance en épaisseur du liège. Le poids de liège récolté est affecté significativement par les coefficients de démasclage et les classes de diamètre. Deux modèles de croissance et de biomasse de liège ont été développés à l'échelle de l'arbre avec des précisions respectives de 96 % et 93 %. Les coefficients de démasclage réellement appliqués en Tunisie varient de 0,6 à 3,6 avec une valeur moyenne de 1,7. Les pertes dans les quantités de liège récolté issues des deux inventaires ont été respectivement estimées, en cas d'application des coefficients de démasclage étudiés (2, 2,5 et 3), à 14 %, 32 % et 50 %.

**Mots-clés** : chêne-liège, coefficient de démasclage, croissance, biomasse, modèle, Tunisie.

# ABSTRACT

## Modelling cork oak growth and cork production in Tunisia

The lack of necessary tools to estimate cork production and assess the management of cork oak forests (Quercus suber L.) in Tunisia is handicapping their development and conservation. The aims of this study were: (i) to specify the impact of first-layer bark removal on tree growth and cork production, (ii) to develop models for estimating cork oak growth and biomass on the scale of individual trees, (iii) to check the bark removal coefficients applied in Tunisia and their impact on the country's cork production. The data were collected from bark removal experiments on several sites and from two inventories conducted at different sites during two successive cork harvesting campaigns. The same experiment was conducted at three trial sites that were representative of productive cork oak forests in Tunisia, with three bark removal coefficients applied to trees of three different diameter classes. Measurements of cork thickness at breast height were taken yearly over a 12-year rotation, at the end of which the cork produced by each tree was harvested and weighed. The results obtained show that the effects of the bark removal coefficients on cork growth and thickness were non-significant, while the effects of diameter class were highly significant. The weight of the harvested cork was significantly affected both by the bark removal coefficients and the diameter class. Two models were developed for the individual tree scale for cork oak growth and biomass, with an accuracy of 96% and 93% respectively. The bark removal coefficients actually applied in Tunisia vary from 0.6 to 3.6 with an average of 1.7. Losses of cork quantities shown by the two inventories were estimated at 14%, 32% and 50% respectively for the three bark removal coefficients applied and analysed (2, 2.5 and 3).

**Keywords:** cork oak, bark removal coefficient, growth, biomass, model, Tunisia.

# RESUMEN

## Modelización del crecimiento y la producción de corcho en Túnez

La ausencia de las herramientas necesarias para la estimación de las producciones, la gestión y la ordenación de la masa de alcornoque (Quercus suber L.) en Túnez constituye un reto para el desarrollo y la conservación del alcornocal en el país. Los objetivos de este estudio fueron (i) precisar el impacto de la intensidad de descorche de los alcornoques en el crecimiento y la producción de corcho, (ii) elaborar modelos de estimación del crecimiento y la biomasa del corcho a escala del árbol y (iii) verificar los coeficientes de descorche aplicados en Túnez y su impacto en la producción nacional de corcho. Los datos recogidos provienen de un ensayo en diferentes localizaciones de descorche y de dos inventarios realizados en diferentes zonas de descorche durante dos campañas sucesivas de recogida de corcho. Tres zonas experimentales representativas de los bosques productivos de alcornoque en Túnez fueron objeto de un mismo ensayo, en el que se aplicaron tres coeficientes de descorche en árboles pertenecientes a tres clases de diámetro. Anualmente se efectuaron medidas de espesor del corcho a la altura del pecho durante una rotación de doce años, al final de la cual el corcho producido por cada árbol fue recogido y pesado. Los resultados obtenidos mostraron un efecto «coeficientes de descorche» no significativo y un efecto «clases de diámetro» muy altamente significativo sobre el crecimiento en grosor del corcho. El peso del corcho recogido resulta significativamente afectado por los coeficientes de descorche y las clases de diámetro. Se desarrollaron dos modelos de crecimiento y de biomasa de corcho a escala del árbol con precisiones respectivas del 96 % y del 93 %. Los coeficientes de descorche realmente aplicados en Túnez varían de 0,6 a 3,6, con un valor medio del 1,7. Las pérdidas en las cantidades de corcho recogido provenientes de los dos inventarios se estimaron, en caso de aplicación de los coeficientes de descorche estudiados (2, 2,5 y 3), en 14%, 32% y 50% respectivamente.

**Palabras clave:** alcornoque, coeficiente de descorche, crecimiento, biomasa, modelo, Túnez.
Introduction

Le chêne-liège, Quercus suber L., est une essence endémique du domaine méditerranéo-atlantique du bassin méditerranéen. Il appartient au plan systématique à l'ordre des Fagales, à la famille des Fagacées, à la sous-famille des Quercineae ou Quercoidées et au genre Quercus (Bouchafra et Fraval, 1991). Ce genre comporte 200 à 500 espèces (Natividade, 1956) dont Quercus suber L. qui est la seule espèce de chêne produisant du liège (Bonnet, 2019). Le chêne-liège est considéré parmi les espèces forestières les plus importantes des forêts méditerranéennes pour ses multiples rôles écologiques, socio-économiques et environnementaux. La superficie mondiale des forêts de chêne-liège est de l'ordre de 2 306 000 ha dont 37,5 % se trouvent en Afrique du Nord. Le chêne-liège occupe 376 000 ha au Maroc, 400 000 ha en Algérie (Bendaanoun, 1996, in Nouri, 2009) et 90 423 ha en Tunisie dont 70 208 ha à l'état pur et 20 215 ha en mélange principalement avec le chêne zeen (DGF, 2005).

Le liège est une matière première naturelle très versatile et porteuse de caractéristiques uniques relatives à sa légèreté, son élasticité et sa compressibilité (Bonnet, 2019). Il est imperméable aux liquides et aux gaz et présente des capacités d'isolation thermique et acoustique excellentes. Il s'agit d'un tissu mort qui joue le rôle d'un revêtement protecteur de l'arbre (Veillon, 1998). Lorsque son épaisseur dépasse 10 mm, le liège protège l'arbre contre le feu (Ghalem, 2006).

L'exploitation de l'arbre de chêne-liège se fait en plusieurs étapes et donne des qualités de liège différentes tout au long de sa vie. Le premier écorçage se fait lorsque l'arbre atteint 25-30 ans, avec une circonférence du tronc à hauteur de poitrine de plus de 70 cm, et donne un liège de structure irrégulière et très dur appelé « liège mâle » ou « liège vierge », ou bien encore « liège blanc » (Bonnet, 2019). Il est davantage utilisé dans les revêtements de sols ou comme isolant. Le deuxième écorçage est effectué 9 à 12 ans plus tard et donne le liège de première production qui est une matière encore irrégulière, mais moins dure et de qualité moyenne. Ce n'est qu'avec le troisième écorçage qu'on obtient un liège plus noble avec des qualités requises pour la fabrication de bouchons de liège. À partir de cette phase et à des intervalles réguliers de 9 à 12 ans, le chêne-liège fournit du liège de qualité pendant environ un siècle et demi (Bonnet, 2019). En effet, d'après le même auteur, l'opération d'écorçage peut être renouvelée jusqu'à 16 fois si l'arbre est correctement exploité.

L'un des principaux soucis des gestionnaires forestiers et des chercheurs concernant la gestion du chêne-liège est l'impact de l'intensité ou du coefficient de démasclage appliqué sur l'arbre et sa production de liège durant la prochaine rotation ou période de croissance du liège (Paulo et Tomé, 2017). Pour déterminer la hauteur de la partie du tronc à écorcer ou à démascler, les gestionnaires forestiers appliquent un coefficient de démasclage qui dépend de la

circonférence du tronc à hauteur de poitrine (CBH). Le coefficient de démasclage à appliquer est variable et dépend surtout de l'état de l'arbre et de la qualité de la station. Quant à la durée de la rotation de récolte, celle-ci est déterminée selon l'épaisseur du liège marchand et elle est variable d'un pays à un autre et diffère aussi dans un même pays en fonction de la fertilité des stations (Kouba, 2005).

Les modèles de prédiction de la croissance et la production du liège constituent de nos jours des outils indispensables pour la gestion et l'aménagement des forêts de chêne-liège, ceci pour deux raisons. D'une part, pour l'aménagement des forêts, ces modèles sont des outils aussi importants que les modèles d'indices de sites, de distribution de tiges ou de qualité du liège. D'autre part, pour l'économie, l'évaluation de la production de liège à l'échelle locale, régionale et nationale permet une meilleure programmation de l'approvisionnement de l'industrie en matière première et l'exportation des produits manufacturés (Vázquez et Pereira, 2005). Les études de la croissance, de la production et de la qualité du liège sont nombreuses dans le bassin méditerranéen, parmi lesquelles on peut citer celles développées au Portugal (Ferreira et al., 1991 ; Ribeiro et Tomé, 2002 ; Pereira et Tomé, 2004 ; Vázquez et Pereira, 2005 ; Paulo et Tomé, 2010 ; Pereira, 2011 ; Paulo et Tomé, 2017 ; Costa et al., 2020), en Espagne (Montero, 1987 ; Montero et Cañellas, 2003 ; Sánchez-González et al., 2007, 2008 ; Pasalodos-Tato et al., 2018a, b), en Tunisie (Sebei et al., 2004 ; Garchi et Sghaier, 2008 ; Sghaier et al., 2011 ; Stiti et al., 2014, 2019), en Algérie (Messaoudene et Metna, 2010 ; Dehane et Ghefar, 2017 ; Naggar et al., 2019) et au Maroc (El Alami et al., 2019).

Nous nous sommes proposé dans ce travail, en nous basant sur le traitement de données collectées à partir d'un essai multisites et à long terme de démasclage et de deux inventaires réalisés dans différents chantiers de récolte de liège :

* d'étudier l'impact de l'intensité de démasclage du chêne-liège sur la croissance et la production du liège ;

* d'élaborer des modèles d'estimation de la croissance et la biomasse du liège à l'échelle de l'arbre ;

* de vérifier les coefficients de démasclage réellement appliqués en Tunisie et mesurer leur impact sur la production nationale de liège.

<figure><img src="image_6.png" /><figcaption>**Figure 1.** Répartition du chêne-liège en Tunisie (Garchi et Sghaier, 2008).</figcaption></figure>
# Matériel et méthodes

## Matériel expérimental et données collectées

Un essai multisites de démasclage de liège a été installé dans trois sites représentatifs de la forêt productive de chêne-liège en Tunisie, Ain-Draham, Chihia et Tébénia (figure 1).

Les trois sites expérimentaux (Ain-Draham, Chihia et Tébénia) sont situés dans la région de Kroumirie au nord-ouest de la Tunisie (Garchi et Sghaier, 2008).

* Le site de Ain-Draham est situé dans un bioclimat humide à hiver tempéré. Il est caractérisé par un relief accidenté et un sol brun forestier très humifère. Sur le plan floristique, le groupement à *Quercus suber*, *Arbutus unedo*, *Pistacia lentiscus*, *Erica arborea* est individualisé au sein de ce site.

* Le site de Chihia appartient au bioclimat humide supérieur à hiver tempéré avec un relief accidenté dans son ensemble et un sol brun lessivé. Le groupement floristique de ce site est composé de *Quercus suber* et *Cytisus triflorus* et un maquis où domine *Erica scoparia*.

* Le site de Tébénia se trouve dans un bioclimat humide supérieur à hiver tempéré. Ce site est caractérisé par des sols bruns lessivés à bruns forestiers. Sur le plan floristique, le groupement climacique à *Quercus suber* et *Cytisus triflorus* est défini dans ce site.

Le dispositif expérimental adopté pour l'ensemble des trois sites est un dispositif complètement aléatoire à deux facteurs étudiés. Il s'agit d'un premier facteur relatif aux dimensions des tiges de chêne-liège et d'un second facteur relatif aux coefficients de démasclage appliqués. Dans chaque site, les arbres qui constituent le dispositif expérimental ont été répartis en trois classes de circonférence désignées par $C_1$, $C_2$ et $C_3$. Trois coefficients de démasclage ont par la suite été appliqués pour chacune des trois classes de circonférence. En plus du coefficient de démasclage théoriquement utilisé par les gestionnaires forestiers en Tunisie, correspondant à deux fois la circonférence mesurée à 1,30 m et désigné par $K_1$, deux autres coefficients de démasclage $K_2$ et $K_3$ ont été testés, correspondant respectivement à 2, 2,5 et 3 fois la circonférence du tronc de l'arbre à 1,30 m de hauteur ou à hauteur de poitrine (CBH). Le nombre de répétitions pour chaque combinaison « classes de circonférence x coefficients de démasclage » est de neuf arbres. Les essais ont été installés en juillet 1996 pour les deux sites de Ain Draham et Chihia, et en juillet 1997 pour le site de Tébénia.

Des mesures annuelles de l'épaisseur du liège à hauteur de poitrine ont été effectuées durant toute la durée de rotation fixée à 12 ans. À la fin de la rotation de 12 ans, le liège de chaque arbre des différents dispositifs expérimentaux a été récolté et les paramètres suivants ont été mesurés :

* diamètre à 1,30 m sur écorce DBH (m) ;

* épaisseur du liège à 1,30 m EP<sub>BH</sub> (mm) ;

* épaisseur du liège à la limite supérieure de la zone démasclée EP<sub>LS</sub> (mm) ;

* hauteur de la zone démasclée H<sub>DEM</sub> (m) ;

* biomasse du liège récolté (kg).

Parallèlement aux essais de démasclage, deux inventaires visant à vérifier les coefficients de démasclage pratiqués réellement par l'administration forestière tunisienne et l'estimation de la production de liège ont été effectués durant les campagnes de récolte du liège de 2016 et 2017. En 2016, 401 arbres de différentes dimensions ont été réprés dans huit chantiers de démasclage (Ain-Boulehya, Ain-Hamraya, Ain-Serj, Ain-Zena, El-Mzara, Oued Zene, Rouii, Tébénia) et ont fait l'objet des mesures du diamètre à hauteur de poitrine (DBH), de l'épaisseur du liège à hauteur de poitrine (EP<sub>BH</sub>), de la hauteur de démasclage (H<sub>DEM</sub>) et de la biomasse du liège récolté. Les mesures de la biomasse du liège réalisées en 2016 et celles effectuées dans les dispositifs expérimentaux de démasclage ont servi à l'évaluation d'un modèle de prédiction de la biomasse de liège à l'échelle de l'arbre.

En 2017, 40 placettes de forme circulaire de 10 arbres chacune ont été choisies aléatoirement dans six chantiers de démasclage (Adissa, Ain-Hamraya, Ain-Snoussi, Ain-Sobah, Babouche, Oued Zene). Le rayon de chaque placette, dont le centre est matérialisé par un arbre, correspond à la distance qui sépare l'arbre le plus éloigné de la placette de celui du centre. Dans chacune de ces placettes, tous les arbres renfermant du liège femelle ont fait l'objet des mesures du diamètre à hauteur de poitrine (DBH), de l'épaisseur du liège à hauteur de poitrine (EP<sub>BH</sub>) et de la hauteur de démasclage (H<sub>DEM</sub>). En plus de la vérification des coefficients de démasclage pratiqués par les gestionnaires forestiers en parallèle aux mesures de 2016, ces placettes de mesures ont servi à l'estimation du rendement du liège à l'échelle des triages inventoriés par le biais du modèle développé, d'une part, pour les coefficients de démasclage réellement pratiqués et, d'autre part, pour les coefficients 2, 2,5 et 3 fois la circonférence à hauteur de poitrine testés dans le cadre de cette étude.

## Méthodes statistiques adoptées

### Effet de la grosseur des tiges et de l'intensité de démasclage sur la croissance et la production du liège

Pour étudier l'effet de la grosseur des tiges et de l'intensité de démasclage, et leur interaction sur la croissance en épaisseur et la biomasse du liège, des analyses individuelles intrasites et des analyses globales intersites ont été effectuées.

### Analyses intrasite

Une analyse de la variance à deux critères de classification selon un modèle croisé fixe, réalisée sur les données mesurées à la fin de la rotation (après 12 ans de croissance), a été effectuée pour chaque station. L'égalité des variances a été vérifiée préalablement à l'aide du test de Bartlett (Dagnelie, 2011).
## Analyses intersites

Pour regrouper les trois dispositifs en une seule analyse et tester les effets sites, coefficients de démasclage et leur interaction, une analyse de la variance à trois critères de classification selon un modèle partiellement hiérarchisé fixe a été utilisée.

En cas d'absence d'interaction (coefficients de démasclage x classes de circonférence) dans le cas des analyses intrasite et (coefficients de démasclage x sites expérimentaux) dans le cas des analyses intersites, des comparaisons multiples de moyennes des modalités des facteurs (coefficients de démasclage, classes de circonférence et sites expérimentaux) dont les effets sont significatifs ont été effectuées en recourant à la méthode de Newman-Keuls (Dagnelie, 2011).

## Prédiction de la biomasse du liège à l'échelle de l'arbre

Afin d'établir une équation qui permette d'estimer la biomasse du liège récoltée à l'échelle de l'arbre, différentes relations linéaires et non linéaires ont été ajustées et comparées en considérant comme variables exogènes les caractéristiques dendrométriques de l'arbre et du liège produit. L'équation retenue devait avoir tous les coefficients significatifs et assurer la meilleure précision (R²<sub>aj</sub> maximum), la normalité de la distribution des résidus et la constance de ces derniers en fonction des valeurs prédites (homosédasticité).

## Modélisation de la croissance en épaisseur du liège

Les fonctions de croissance décrivent des variations dans la taille globale d'un organisme ou d'une population avec l'âge ; elles peuvent également décrire les changements d'une variable particulière d'un arbre ou d'un peu- plement avec l'âge, dans ce cas l'épaisseur du liège. La croissance en épaisseur du liège d'un arbre, représentée par des valeurs observées non négatives, constitue un phénomène limité et non décroissant. Le modèle de croissance utilisé doit donc présenter une asymptote horizontale pour les valeurs du temps tendant vers l'infini et ne pourra pas prendre des valeurs négatives. Différents modèles présentant ces caractéristiques sont cités dans la littérature. Il s'agit de fonctions continues, non négatives, non décroissantes et bornées (Debouche, 1979).

Dans la plupart des cas, un modèle de croissance est élaboré en vue de pouvoir prédire la taille moyenne ou le taux de croissance probable à un âge donné et/ou de faciliter la comparaison de la croissance estimée à d'autres estimations publiées. Ce qui, mathématiquement, se traduit par les notions (Abouba, 2014) suivantes :

* vitesse de croissance absolue (Va) : dérivée de la taille par rapport au temps ;

* vitesse de croissance relative (Vt) : vitesse de croissance absolue rapportée à la taille de l'individu, ce qui permet de définir la forme des courbes ;

* vitesse moyenne (Vm) : hauteur moyenne de la courbe des vitesses de croissance sur l'axe de la taille ;

* temps de croissance (T) : rapport entre la taille maximale et la vitesse moyenne.

<table><caption>Tableau I. Fonctions de croissance étudiées.</caption><thead><tr><td>Fonction</td><td>Nom</td><td>Équation</td><td>V<sub>m</sub> (mm/an)</td><td>T (années)</td></tr></thead><tbody><tr><td>F1</td><td>Johnson-Schumacher</td><td>Y = a<sub>1</sub> x exp(- (a<sub>3</sub> / (t - a<sub>2</sub>))</td><td>a<sub>1</sub> / (4a<sub>3</sub>)</td><td>4a<sub>3</sub></td></tr><tr><td>F2</td><td>Gompertz</td><td>Y = a<sub>1</sub> x exp(-exp(- (t - a<sub>2</sub>) / a<sub>3</sub>))</td><td>a<sub>1</sub> / (4a<sub>3</sub>)</td><td>4a<sub>3</sub></td></tr><tr><td>F3</td><td>Gauss modifié</td><td>Y = a<sub>1</sub> x (1 - exp(- ((t - a<sub>2</sub>)<sup>2</sup>) / a<sub>3</sub>))</td><td>√(π / 2) a<sub>1</sub> / (2a<sub>3</sub>)</td><td>2√(2/π) a<sub>3</sub></td></tr><tr><td>F4</td><td>Logistique</td><td>Y = a<sub>1</sub> / (1 + exp(- (t - a<sub>2</sub>) / a<sub>3</sub>))</td><td>a<sub>1</sub> / (6a<sub>3</sub>)</td><td>6a<sub>3</sub></td></tr></tbody></table>

Le tableau I présente les fonctions de croissance étudiées à trois paramètres (Debouche, 1979) avec les valeurs de la vitesse moyenne V<sub>m</sub> (mm/an) et le temps de croissance T (années). Ces deux dernières caractéristiques, liées au réalisme biologique des fonctions de croissance étudiées, sont considérées parmi les critères de comparaison de ces fonctions.

Les trois paramètres de ces différentes fonctions désignent :

* Y : la variable étudiée (épaisseur du liège en mm) ;

* t : le temps (âge du liège en années) ;

* a₁ : l'asymptote horizontale, correspondant à la valeur maximale vers laquelle tend Y ;

* a₂ : la valeur de t qui situe la courbe sur l'axe des x par son point d'inflexion ou son origine ;

* a₃ : paramètre exprimé dans les mêmes unités que t et qui mesure l'étalement du phénomène de croissance sur l'axe des x et donc la vitesse de croissance.

Étant donné le caractère longitudinal des données utilisées pour la modélisation de la croissance en épaisseur du liège des arbres, c'est-à-dire des données répétées sur les mêmes individus (arbres), il y a forcément une autocorrélation entre les résidus successifs d'un même arbre. Pour contrôler cette autocorrélation et augmenter la précision des estimations des fonctions ajustées, un processus auto-
régressif sur les erreurs issues des observations en provenance d'un même individu a été incorporé dans le modèle d'ajustement des paramètres des différentes fonctions étudiées (Diéguez-Aranda et al., 2006). Il s'agit de la fonction « structure autorégressive d'erreurs successives » CAR(x) exprimée par :

$$e_{ij} = \sum_{n=1}^{x} d_n \rho_n^{t_{i-j}(n)} e_{i(j-n)} + \varepsilon_{ij} \quad (\text{équation 1})$$

où :

• $e_{ij}$ est la différence entre l'épaisseur du liège observée et celle estimée de l'arbre i à l'âge j;

• $d_n$ est une variable « dummy » qui prend la valeur 1 si les mesures correspondent à des âges différents pour un même arbre et 0 sinon ;

• $\rho_n$ est le paramètre du nème ordre autorégressif à estimer ;

• $t_{ij} - t_{i(j-n)}$ est la durée (en années) qui sépare la jème de la (j - n)ème observation.

Pour corriger l'autocorrélation des résidus et déterminer l'ordre de la fonction CAR(x) à utiliser, des graphiques présentant les résidus de la variable explicative contre les résidus de l'observation précédente de chaque arbre ont été présentés et examinés.

Toutes les fonctions de croissance ajustées ont été comparées sur la base de leur capacité d'ajustement, leur capacité de prédiction, leur réalisme biologique ainsi que la normalité et l'homoscédasticité de la distribution de leurs résidus.

## Capacité d'ajustement

Sont alors calculés les éléments suivants :

• biais (moyenne des résidus) :

$$biais = \frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_i) \quad (\text{équation 2})$$

• écart-type résiduel (Root Mean Square Errors) :

$$RMSE = \sqrt{\frac{\sum_{i=1}^{n} (y_i - \hat{y}_i)^2}{n-p}} \quad (\text{équation 3})$$

• coefficient de détermination ajusté :

$$R_{aj}^2 = 1 - \frac{(n-1) \sum_{i=1}^{n} (y_i - \hat{y}_i)^2}{(n-p) \sum_{i=1}^{n} (y_i - \bar{y})^2} \quad (\text{équation 4})$$

Ce coefficient représente la part de variabilité de la variable dépendante qui est expliquée par la régression en tenant compte du nombre de paramètres p dans le modèle (n étant le nombre d'observations).

## Capacité de prédiction

La capacité de prédiction des différentes fonctions étudiées est basée sur le calcul des résidus PRESS (Prediction Sum of Squares) ou erreurs de prédiction (Sánchez-González et al., 2005). Ces résidus sont équivalents aux résidus obtenus en omettant chaque fois une observation. Il s'agit d'ajuster les paramètres de l'équation avec les autres observations restantes, estimer la réponse de l'observation omise et com-

parer la valeur prédite avec celle observée : $y_i - \hat{y}_{i,i} = e_{i,i} (i = 1, 2, ..., n)$ où $y_i$ est le volume observé, $\hat{y}_{i,i}$ est la valeur estimée de l'observation i (lorsque cette dernière est absente des données utilisées pour l'ajustement des paramètres de l'équation de régression) et $n$ est le nombre d'observations. Chaque équation ou fonction possède dans ce cas n résidus PRESS et la quantité PRESS (somme des carrés des écarts de prédiction) est ainsi définie :

$$PRESS = \sum_{i=1}^{n} (y_i - \hat{y}_{i,i})^2 \quad (\text{équation 5})$$

Plus la valeur de PRESS est proche de la somme des carrés des écarts résiduels, meilleure est la capacité de prédiction de l'équation de régression. Les résidus PRESS ont été également utilisés pour calculer les critères de comparaison suivants :

• biais de prédiction : $biais_p = \frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_{i,i}) \quad (\text{équation 6})$;

• écart-type résiduel de prédiction :

$$RMSE_p = \sqrt{\frac{\sum_{i=1}^{n} (y_i - \hat{y}_{i,i})^2}{n-p}} \quad (\text{équation 7})$$

• coefficient de détermination de prédiction :

$$R_p^2 = 1 - \frac{\sum_{i=1}^{n} (y_i - \hat{y}_{i,i})^2}{\sum_{i=1}^{n} (y_i - \bar{y})^2} \quad (\text{équation 8})$$

## Réalisme biologique

Le réalisme biologique des fonctions étudiées est basé sur le calcul de la vitesse moyenne annuelle de la croissance et le temps sur lequel s'étale la croissance du liège.

## Normalité de la distribution des résidus

La violation possible de la normalité de la distribution des erreurs a été examinée à l'aide du test de normalité de Ryan-Joiner (Ryan et Joiner, 1976) et de la présentation graphique des quantiles normaux (QQ-probability plots) ou droite de Henry.

Toutes les analyses statistiques ont été réalisées à l'aide des procédures GLM de SAS/STAT (SAS, 2004) et MODEL de SAS/ETS (SAS, 2004).

## Résultats

### Effet du coefficient de démasclage et de la grosseur des tiges sur l'épaisseur du liège (essais de démasclage)

### Analyses intrasites

Le tableau II présente, pour chacun des trois sites expérimentaux et chacune des variables mesurées, les valeurs et la signification du test F de Fisher issues de l'analyse de la variance à deux critères de classification relatives aux tests des effets « coefficients de démasclage », « classes de circonférence » et leur interaction sur la croissance en épaisseur du liège.

D'après les résultats du tableau II, seule la dimension du tronc des arbres (classes de circonférence) affecte significativement ($\alpha \le 0,001$) la croissance en épaisseur du liège mesurée à hauteur de poitrine (EP<sub>BH</sub>) et à la limite
supérieure de la zone démasclée (EP<sub>LS</sub>), ceci pour les trois sites expérimentaux. Le même tableau montre également l'absence d'effet significatif des coefficients de démasclage pratiqués ainsi que celui de l'interaction « classes de circonférence » x « coefficients de démasclage ». Concernant la biomasse du liège récolté à l'échelle de l'arbre dans les trois sites expérimentaux, l'analyse de la variance (tableau II) a montré, comme attendu, un effet très hautement significatif (α ≤ 0,001) de la grosseur du tronc, un effet significatif (α ≤ 0,05) à hautement significatif (α ≤ 0,01) des coefficients de démasclage et l'absence de l'effet interaction entre les deux facteurs étudiés.

L'absence d'interaction signifie que les deux facteurs étudiés agissent indépendamment l'un de l'autre et qu'il est possible de comparer les modalités de chaque facteur significatif à part. Le tableau III présente, par variable et par site expérimental, les valeurs moyennes des différentes modalités des deux facteurs étudiés et les résultats de la comparaison multiple de ces moyennes par le biais du test de Newman-Keuls (Dagnelie, 2011).

Le tableau III montre que l'épaisseur du liège aux deux niveaux de mesure augmente avec la taille de l'arbre. L'écart enregistré entre les deux classes de circonférence extrêmes varie en fonction du site entre environ 10 et 12 mm pour l'épaisseur moyenne du liège à hauteur de poitrine, et est de 48 à 57 kg/arbre pour la biomasse moyenne du liège, ceci en faveur des arbres les plus gros. Les trois coefficients de démasclage testés donnent en moyenne des valeurs semblables du point de vue de l'épaisseur du liège mais, au contraire, des écarts élevés du point de vue de la biomasse, surtout entre les deux coefficients extrêmes (K<sub>1</sub> = 2 et K<sub>3</sub> = 3) qui varient en moyenne d'environ 8 à 19 kg/arbre selon le site expérimental.

## Analyse intersites (analyse globale)

Le tableau IV présente les résultats de l'analyse de la variance intersites (regroupement des essais) réalisée sur les trois variables mesurées. Il s'agit d'une analyse de la variance à trois critères de classification partiellement hiérarchisés. L'objet de cette analyse est de tester principalement l'effet « sites » (considéré comme un facteur fixe) et l'interaction « sites x coefficients de démasclage ».

Les résultats du tableau IV montrent un effet

## Tableau II.

Valeurs et signification du test F de Fisher relatives aux analyses de la variance par site expérimental et par variable mesurée.

<table><thead><tr><th rowspan="2">Variables mesurées</th><th rowspan="2">Sites</th><th colspan="3">Valeurs du test F de Fisher (F observé)</th><th rowspan="2">RMSE</th></tr><tr><th>Coefficients de démasclage</th><th>Classes de circonférence</th><th>Interaction</th></tr></thead><tbody><tr><td rowspan="3">EP<sub>BH</sub></td><td>Ain-Draham</td><td>1,04</td><td>8,61***</td><td>1,07</td><td>6,9171</td></tr><tr><td>Chihia</td><td>0,94</td><td>20,63***</td><td>0,42</td><td>5,7453</td></tr><tr><td>Tébénia</td><td>0,52</td><td>23,21***</td><td>0,40</td><td>5,4433</td></tr><tr><td rowspan="3">EP<sub>LS</sub></td><td>Ain-Draham</td><td>0,74</td><td>9,59***</td><td>1,31</td><td>5,4947</td></tr><tr><td>Chihia</td><td>0,47</td><td>18,14***</td><td>0,37</td><td>5,6861</td></tr><tr><td>Tébénia</td><td>1,15</td><td>12,30***</td><td>0,33</td><td>4,9275</td></tr><tr><td rowspan="3">Biomasse</td><td>Ain-Draham</td><td>2,81*</td><td>32,50***</td><td>1,61</td><td>17,7441</td></tr><tr><td>Chihia</td><td>2,37*</td><td>50,14***</td><td>0,19</td><td>14,5269</td></tr><tr><td>Tébénia</td><td>6,87**</td><td>51,58***</td><td>1,71</td><td>15,3041</td></tr></tbody><tfoot><tr><td colspan="6">* : effet significatif pour α = 0,05 ; ** : effet significatif pour α = 0,01 ; *** : effet significatif pour α = 0,001 ; EP<sub>BH</sub> : épaisseur du liège en mm à hauteur de poitrine (DBH) ; EP<sub>LS</sub> : épaisseur du liège en mm à la limite supérieure de la surface démasclée ; biomasse : poids du liège en kg par arbre.</td></tr></tfoot></table>

## Tableau III.

Comparaison de moyennes par site expérimental et par variable mesurée.

<table><thead><tr><th rowspan="2">Variables mesurées</th><th rowspan="2">Sites</th><th colspan="3">Classes de circonférence</th><th colspan="3">Coefficients de démasclage</th></tr><tr><th>C<sub>1</sub></th><th>C<sub>2</sub></th><th>C<sub>3</sub></th><th>K<sub>1</sub></th><th>K<sub>2</sub></th><th>K<sub>3</sub></th></tr></thead><tbody><tr><td rowspan="3">EP<sub>BH</sub></td><td>Ain-Draham</td><td>39,20<sup>b</sup></td><td>48,47<sup>a</sup></td><td>49,12<sup>a</sup></td><td>43,53<sup>a</sup></td><td>45,87<sup>a</sup></td><td>47,39<sup>a</sup></td></tr><tr><td>Chihia</td><td>33,28<sup>c</sup></td><td>37,04<sup>b</sup></td><td>45,37<sup>a</sup></td><td>39,52<sup>a</sup></td><td>37,30<sup>a</sup></td><td>38,88<sup>a</sup></td></tr><tr><td>Tébénia</td><td>33,89<sup>b</sup></td><td>44,34<sup>a</sup></td><td>46,12<sup>a</sup></td><td>40,74<sup>a</sup></td><td>40,82<sup>a</sup></td><td>42,08<sup>a</sup></td></tr><tr><td rowspan="3">EP<sub>LS</sub></td><td>Ain-Draham</td><td>35,85<sup>b</sup></td><td>44,40<sup>a</sup></td><td>42,53<sup>a</sup></td><td>39,47<sup>a</sup></td><td>41,66<sup>a</sup></td><td>41,65<sup>a</sup></td></tr><tr><td>Chihia</td><td>29,09<sup>c</sup></td><td>33,85<sup>b</sup></td><td>40,51<sup>a</sup></td><td>35,55<sup>a</sup></td><td>34,27<sup>a</sup></td><td>33,62<sup>a</sup></td></tr><tr><td>Tébénia</td><td>28,25<sup>b</sup></td><td>34,99<sup>a</sup></td><td>36,39<sup>a</sup></td><td>32,45<sup>a</sup></td><td>32,19<sup>a</sup></td><td>35,00<sup>a</sup></td></tr><tr><td rowspan="3">Biomasse</td><td>Ain-Draham</td><td>21,92<sup>c</sup></td><td>40,98<sup>b</sup></td><td>79,44<sup>a</sup></td><td>36,19<sup>b</sup></td><td>50,67<sup>a</sup></td><td>55,47<sup>a</sup></td></tr><tr><td>Chihia</td><td>9,52<sup>c</sup></td><td>27,27<sup>b</sup></td><td>57,87<sup>a</sup></td><td>28,20<sup>b</sup></td><td>30,48<sup>b</sup></td><td>35,98<sup>a</sup></td></tr><tr><td>Tébénia</td><td>20,73<sup>c</sup></td><td>37,63<sup>b</sup></td><td>73,70<sup>a</sup></td><td>34,80<sup>b</sup></td><td>49,15<sup>a</sup></td><td>48,11<sup>a</sup></td></tr></tbody><tfoot><tr><td colspan="8">Les moyennes de même lettre en exposant ne présentent pas de différences significatives, selon le test de Newman-Keuls.</td></tr></tfoot></table>

## Tableau IV.

Résultats de l'analyse de la variance : essai multisites.

<table><thead><tr><th rowspan="2">Sources de variation</th><th colspan="3">Valeurs du test F de Fisher (F observé)</th></tr><tr><th>EP<sub>BH</sub></th><th>EP<sub>LS</sub></th><th>Biomasse du liège</th></tr></thead><tbody><tr><td>Coefficients</td><td>1,11</td><td>0,33</td><td>7,95<sup>**</sup></td></tr><tr><td>Sites</td><td>17,24<sup>***</sup></td><td>28,80<sup>***</sup></td><td>13,03<sup>***</sup></td></tr><tr><td>Classes / Sites</td><td>16,35<sup>***</sup></td><td>13,53<sup>***</sup></td><td>44,18<sup>***</sup></td></tr><tr><td>Coefficients x Sites</td><td>0,78</td><td>1,04</td><td>1,13</td></tr><tr><td>(Coefficients x classes) / Sites</td><td>0,70</td><td>0,69</td><td>1,27</td></tr><tr><td>Écart-type résiduel (RMSE)</td><td>6,0253</td><td>5,3791</td><td>15,7544</td></tr></tbody></table>
<caption>Tableau V.</caption>

Résultats des comparaisons multiples de moyennes : essai multisites.

<table><thead><tr><th rowspan="2">Variables mesurées</th><th colspan="3">Sites</th><th colspan="3">Coefficients de démasclage</th></tr><tr><th>Ain-Draham</th><th>Chihia</th><th>Tébénia</th><th>K<sub>1</sub></th><th>K<sub>2</sub></th><th>K<sub>3</sub></th></tr></thead><tbody><tr><td>EP<sub>BH</sub></td><td>45,60<sup>a</sup></td><td>38,57<sup>c</sup></td><td>41,45<sup>b</sup></td><td>41,26<sup>a</sup></td><td>41,33<sup>a</sup></td><td>43,02<sup>a</sup></td></tr><tr><td>EP<sub>LS</sub></td><td>41,92<sup>a</sup></td><td>34,48<sup>b</sup></td><td>33,21<sup>b</sup></td><td>35,82<sup>a</sup></td><td>36,04<sup>a</sup></td><td>36,76<sup>a</sup></td></tr><tr><td>Biomasse</td><td>47,45<sup>a</sup></td><td>31,55<sup>b</sup></td><td>44,02<sup>a</sup></td><td>35,72<sup>b</sup></td><td>39,96<sup>ab</sup></td><td>43,97<sup>a</sup></td></tr></tbody><tfoot><tr><td colspan="7">Les moyennes de même lettre en exposant ne présentent pas de différences significatives, selon le test de Newman-Keuls.</td></tr></tfoot></table>

« sites » très hautement significatif de l'épaisseur du liège mesurée aux deux niveaux du tronc et l'absence des effets « coefficients de démasclage » et de l'interaction « coefficients de démasclage x sites ». Le même tableau confirme les résultats obtenus lors des analyses intrasites, c'est-à-dire un effet hautement significatif « classes/sites » et un effet non significatif « (coefficients x classes) / sites ». L'absence d'interaction « coefficients de démasclage x sites » et celle des « coefficients de démasclage » signifie que les trois coefficients de démasclage ne présentent pas de différences significatives entre leurs moyennes entre tous les sites confondus, et que la performance des sites n'est pas affectée positivement ou négativement par un coefficient de démasclage particulier. Ce dernier résultat prouve que les deux coefficients de démasclage les plus élevés (K<sub>2</sub> = 2,5 et K<sub>3</sub> = 3) peuvent être appliqués dans les trois sites étudiés et dans d'autres sites similaires sans réduire la croissance en épaisseur du liège par comparaison au coefficient de

démasclage le plus faible (K<sub>1</sub> = 2).

Pour ce qui est de la biomasse du liège, comme attendu d'après les résultats des analyses intrasites également, le tableau IV montre un effet « coefficients de démasclage » hautement significatif, un effet « sites » très hautement significatif et une absence d'interaction des deux facteurs. Le tableau V présente les résultats des comparaisons multiples de moyennes (épaisseurs et biomasse du liège) des sites pour tous les coefficients de démasclage et de celles des coefficients de démasclage dans tous les sites.

Le tableau V montre que les trois sites se distinguent l'un de l'autre pour la croissance en épaisseur du liège mesurée à hauteur de poitrine (EP<sub>BH</sub>), avec le site de Ain-Draham en tête, suivi du site de Tébénia puis celui de Chihia. Pour l'épaisseur du liège mesurée à la limite supérieure de la zone démasclée (EP<sub>LS</sub>), le site de Ain-Draham occupe toujours la première position et les deux autres sites la deuxième position avec des épaisseurs moyennes similaires. En ce qui concerne la biomasse du liège, les deux sites de Ain-Draham et Tébénia occupent la première place et dépassent celui de Chihia avec respectivement 16 et 12 kg/arbre environ. Les trois coefficients de démasclage se répartissent en trois classes avec le coefficient 3 en tête, le coefficient 2 en troisième position et le coefficient 2,5 en position intermédiaire et chevauchant entre la première et la troisième position. Le gain de poids de liège par arbre entre les deux coefficients extrêmes (K<sub>1</sub> = 2 et K<sub>3</sub> = 3) est en moyenne de 8 kg environ.

## Durée de rotation

La durée de rotation de récolte du liège pratiquée en Tunisie est fixée par l'administration forestière à 12 ans pour toutes les stations et peut dépasser ce délai dans certains cas telle l'absence de moyens financiers pour la récolte. La figure 2 présente les valeurs moyennes de l'épaisseur du liège mesurée à hauteur de poitrine (EP<sub>BH</sub>) pour chacun des trois sites expérimentaux ainsi que les limites de l'épaisseur du liège optimale pour la production de bouchons naturels qui représentent le produit le plus important du liège. Ces limites se situent entre 34 et 44 mm. Les segments aux sommets des barres des histogrammes représentent l'intervalle de confiance autour de l'épaisseur moyenne pour un degré de confiance de 95 %. La normalité de la distribution des mesures de l'épaisseur du liège pour chaque site expérimental a été vérifiée préalablement à l'aide du test de normalité de Ryan-Joiner (Ryan et Joiner, 1976).

En situant l'épaisseur moyenne du liège produit après 12 ans de croissance et ses limites de confiance obtenues au niveau de chaque site expérimental par rapport aux deux épaisseurs limites définissant la classe optimale d'épaisseur permettant la production de bouchons naturels, on constate ce qui suit.

<figure><img src="image_7.png" /><figcaption>Figure 2.</figcaption></figure>

Épaisseur moyenne du liège à hauteur de poitrine (EP<sub>BH</sub>) par site expérimental et classe d'épaisseurs optimales pour la production de bouchons naturels.
* Sur le site de Ain-Draham, plus de 95 % des arbres produisent du liège dont l'épaisseur dépasse largement la zone d'épaisseur optimale. Il est à signaler qu'après prélèvement des bouchons naturels le reste des plaques de liège va se transformer en déchets et servir à la trituration pour produire des bouchons agglomérés. Ainsi, pour optimiser la valorisation des plaques de liège, il ne faut pas que l'épaisseur de ces dernières dépasse la limite supérieure souhaitée par l'industrie du liège. Dès lors et suite aux résultats obtenus dans le site de Ain-Draham, il est possible de réduire la durée de rotation à moins de 12 ans étant donné la qualité de cette station du point de vue de la croissance du liège.

* Sur les sites de Chihia et Tébénia, les épaisseurs moyennes du liège et leurs limites de confiance sont situées en pleine zone d'épaisseur optimale de production de bouchons naturels avec une légère supériorité du site de Tébénia. La durée de rotation de 12 ans peut donc être adoptée pour ces deux sites, voire légèrement réduite pour le site de Tébénia.

## Estimation du poids du liège à l'échelle de l'arbre (essais de démasclage et inventaire)

Afin d'établir un modèle qui permette de prédire la biomasse (poids) du liège récolté à l'échelle de l'arbre, ont été utilisées les mesures d'arbres individuels réalisées en 2016 dans huit chantiers de démasclage et celles effectuées après 12 ans de croissance dans les trois dispositifs expérimentaux des essais de démasclage.

Les trois caractéristiques de l'arbre de chêne-liège les plus corrélées avec la biomasse du liège récolté sont le DBH ($r = 0,83$ ; $p < 0,0001$), la hauteur de la surface démasclée ou $H_{Dem}$ ($r = 0,84$ ; $p < 0,0001$) et l'épaisseur du liège à hauteur de poitrine ou $EP_{BH}$ ($r = 0,64$ ; $p < 0,0001$). Une équation entre la biomasse et chacune de ces caractéristiques a été établie (équations 9, 10 et 11 ; figures 3, 4 et 5). Une relation plus précise entre la biomasse du liège récolté et l'ensemble des trois caractéristiques a été également é laborée (équation 12). Comme pour les tarifs de cubage des arbres forestiers, l'une ou l'autre de ces équations peut être utilisée selon la facilité de l'acquisition et/ou la disponibilité des données (caractéristiques mesurées) et la précision souhaitée des estimations à obtenir.

## Relation entre la biomasse du liège et le DBH

La relation qui lie la biomasse du liège au DBH est une fonction puissance (équation 9). Le DBH explique à lui seul environ 75 % de la variabilité de la biomasse du liège observée au niveau des arbres de chêne-liège considérés dans cette étude :

$$ \text{Biomasse} = 3,518 \times 10^{-3} \text{DBH}^{2,382} \text{ (équation 9)} \\ \text{avec } R^{2}_{aj} = 0,748 \text{ et RMSE} = 8,652 \text{ kg.} $$

La figure 3 présente, sur le même graphique, la projection des biomasses observées en fonction du DBH, supposées à la courbe issue des données simulées à partir de l'équation 9.

<figure><img src="image_8.png" /><figcaption>**Figure 3.**
Données observées et courbe de prédiction de la biomasse du liège en fonction du DBH.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 4.**
Données observées et droite de prédiction de la biomasse du liège en fonction de la hauteur de la zone démasclée.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 5.**
Données observées et courbe de prédiction de la biomasse du liège en fonction de son épaisseur à hauteur de poitrine.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 6.**

QQ-probability plots des résidus (a) et projection des résidus en fonction des biomasses estimées (b).</figcaption></figure>

### Relation entre la biomasse du liège et la hauteur de démasclage ($H_{DEM}$)

L'équation 10 décrit la relation liant la biomasse du liège à la hauteur de la zone démasclée :

$$ \text{Biomasse} = -14,868 + 17,122 \times H_{\text{DEM}} \text{ (équation 10)} $$

avec $R_{aj}^2 = 0,706$ et $\text{RMSE} = 9,359$ kg.

La figure 4 présente, sur le même graphique, la projection des biomasses observées en fonction de la hauteur de démasclage ($H_{DEM}$), superposées à la droite issue des données simulées à partir de l'équation 10.

### Relation entre la biomasse et l'épaisseur du liège à hauteur de poitrine ($EP_{BH}$)

La prédiction de la biomasse du liège en fonction de son épaisseur à hauteur de poitrine peut se faire à l'aide de l'équation 11. Par comparaison aux deux autres caractéristiques de l'arbre de chêne-liège (DBH et $H_{DEM}$), l'épaisseur du liège à hauteur de poitrine ($EP_{BH}$) semble être la variable la moins précise comme prédicteur de la biomasse de liège produit :

$$ \text{Biomasse} = 2,737 \times 10^{-2} \times e^{(5,577 \times 10^{-2} EP_{BH})} \text{ (équation 11)} $$

Avec $R_{aj}^2 = 0,468$ et $\text{RMSE} = 12,468$ kg.

La figure 5 schématise l'évolution de la biomasse du liège en fonction de son épaisseur, d'une part, à partir des données observées et, d'autre part, à partir des données simulées (courbe) à l'aide de l'équation 11.

### Relation entre la biomasse du liège et DBH, $H_{DEM}$ et $EP_{BH}$

L'équation 12 regroupe les trois variables ou caractéristiques du liège et de l'arbre utilisées comme variables exogènes pour prédire la biomasse du liège récoltée. En se basant sur les valeurs des coefficients de détermination ajustés ($R_{aj}^2$) et celles des écarts-types résiduels (RMSE), cette dernière équation se présente comme la plus précise par rapport au trois précédentes (équations 9, 10 et 11).

$$ \text{Biomasse} = 1,408 \times 10^{-2} \text{DBH}^{1,192} H_{\text{DEM}}^{0,688} \text{EP}_{\text{BH}}^{0,682} \text{ (équation 12)} $$

avec $R_{aj}^2 = 0,932$ et $\text{RMSE} = 4,495$ kg.

Dans cette équation, la biomasse est exprimée en kilogrammes par arbre (kg/arbre), le DBH en centimètres (cm), la hauteur de démasclage ou $H_{DEM}$ en mètres (m) et l'épaisseur du liège mesurée à hauteur de poitrine ou $EP_{BH}$ en millimètres (mm).

La figure 6a présente les quantiles normaux (QQ-probability plots) des résidus issus de l'équation 12 et la figure 6b la projection des mêmes résidus en fonction des biomasses estimées. Ces deux graphiques montrent que la normalité de la distribution des résidus peut être acceptée et que leur répartition autour de l'axe de valeurs nulles peut être considérée comme constante et aléatoire.

La figure 7 montre la forte liaison entre les biomasses prédites et celles observées ($r = 0,9656$ ; $p < 0,0001$).

<figure><img src="image_12.png" /><figcaption>**Figure 7.**

Biomasses estimées en fonction des biomasses observées.</figcaption></figure>

### Prédiction de la croissance du liège en fonction de son âge et du diamètre de l'arbre (essais de démasclage)

#### Choix de l'ordre du modèle autorégressif CAR(x) pour corriger l'autocorrélation entre résidus successifs issus d'un même arbre

La figure 8 (première colonne) montre que, en absence de l'incorporation de la structure autorégressive CAR(x) dans l'ajustement des paramètres de la fonction F1, une tendance linéaire de la corrélation entre les résidus successifs est observée. Après incorporation de la fonction CAR(x) avec un premier ordre (deuxième colonne), puis un deuxième ordre (troisième colonne), la tendance linéaire est considérablement réduite avec CAR(1) puis disparaît complètement avec CAR(2). Ainsi, la structure autorégressive de deuxième ordre (CAR(2)), permettant le contrôle de l'autocorrélation entre
<figure><img src="image_13.png" /><figcaption>**Figure 8.** Lag1-résidus et Lag2-résidus de l'épaisseur du liège en fonction des résidus issus de la fonction F1 ajustée sans tenir compte de la structure autorégressive des erreurs (première colonne), et en utilisant la structure autorégressive des erreurs du premier et du second ordre (deuxième et troisième colonnes, respectivement).</figcaption></figure>

résidus successifs d'un même arbre, a été incorporée dans le processus d'estimation des paramètres des quatre fonctions de croissance étudiées (tableau I).

## Estimation des paramètres et comparaison des fonctions de croissance étudiées

Les valeurs estimées des paramètres de chacune des quatre fonctions étudiées et les caractéristiques statistiques correspondantes, relatives aux différents critères d'évaluation de ces fonctions, figurent au tableau VI. La figure 9 présente une projection des résidus issus de chaque fonction de croissance étudiée en fonction des scores normaux ou des quantiles normaux (QQ-probability plots). La même figure 9 qui présente également les résultats du test de normalité de la distribution de ces résidus selon le test de Ryan-Joiner montre que la distribution de ces dernières ne s'écarte pas trop de celle d'une distribution normale. La normalité de la distribution des résidus issus des quatre fonctions ajustées peut donc être acceptée.

<figure><img src="image_14.png" /><figcaption>**Figure 9.** QQ-probability plots des résidus issus des quatre fonctions étudiées et test de normalité de Ryan-Joiner.</figcaption></figure>
Tableau VI.

Estimation des paramètres et critères de comparaison : ajustement avec CAR(2).

<table><thead><tr><th rowspan="2">Fonction</th><th rowspan="2">Par.</th><th rowspan="2">Estim.</th><th rowspan="2">p > |t|</th><th colspan="3">Capacité d'ajustement</th><th colspan="3">Capacité de prédiction</th><th colspan="2">Réalisme biologique</th></tr><tr><th>Biais</th><th>RMSE</th><th>R<sup>2</sup><sub>aj</sub></th><th>Biais<sub>p</sub></th><th>RMSE<sub>p</sub></th><th>R<sup>2</sup><sub>p</sub></th><th>V<sub>m</sub></th><th>T</th></tr></thead><tbody><tr><td rowspan="5">F1</td><td>a<sub>1</sub></td><td>227,050</td><td>&lt; 0,0001</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>a<sub>2</sub></td><td>-9,836</td><td>&lt; 0,0001</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>a<sub>3</sub></td><td>36,918</td><td>&lt; 0,0001</td><td>-0,0385</td><td>2,888</td><td>0,9438</td><td>-0,0388</td><td>2,895</td><td>0,9437</td><td>1,54</td><td>146,67</td></tr><tr><td>ρ<sub>1</sub></td><td>1,123</td><td>&lt; 0,0001</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>ρ<sub>2</sub></td><td>0,856</td><td>&lt; 0,0001</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="5">F2</td><td>a<sub>1</sub></td><td>74,649</td><td>&lt; 0,0001</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>a<sub>2</sub></td><td>7,631</td><td>&lt; 0,0001</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>a<sub>3</sub></td><td>8,002</td><td>&lt; 0,0001</td><td>-0,0605</td><td>2,893</td><td>0,9436</td><td>-0,0607</td><td>2,900</td><td>0,9435</td><td>2,33</td><td>32,01</td></tr><tr><td>ρ<sub>1</sub></td><td>1,122</td><td>&lt; 0,0001</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>ρ<sub>2</sub></td><td>0,856</td><td>&lt; 0,0001</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="5">F3</td><td>a<sub>1</sub></td><td>75,954</td><td>&lt; 0,0001</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>a<sub>2</sub></td><td>-5,222</td><td>&lt; 0,0001</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>a<sub>3</sub></td><td>19,226</td><td>&lt; 0,0001</td><td>-0,0507</td><td>2,890</td><td>0,9438</td><td>-0,0510</td><td>2,897</td><td>0,9436</td><td>2,46</td><td>30,68</td></tr><tr><td>ρ<sub>1</sub></td><td>1,123</td><td>&lt; 0,0001</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>ρ<sub>2</sub></td><td>0,856</td><td>&lt; 0,0001</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="5">F4</td><td>a<sub>1</sub></td><td>53,419</td><td>&lt; 0,0001</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>a<sub>2</sub></td><td>7,424</td><td>&lt; 0,0001</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>a<sub>3</sub></td><td>3,604</td><td>&lt; 0,0001</td><td>-0,1061</td><td>2,910</td><td>0,9430</td><td>-0,1063</td><td>2,917</td><td>0,9428</td><td>2,47</td><td>21,62</td></tr><tr><td>ρ<sub>1</sub></td><td>1,121</td><td>&lt; 0,0001</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>ρ<sub>2</sub></td><td>0,857</td><td>&lt; 0,0001</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody><tfoot><tr><td colspan="12">Par. : paramètre ; Estim. : estimation ; RMSE : écart-type résiduel ; R<sup>2</sup><sub>aj</sub> : coefficient de détermination ajusté ; Biais<sub>p</sub> : biais PRESS ; RMSE<sub>p</sub> : écart-type résiduel PRESS ; R<sup>2</sup><sub>p</sub> : coefficient de détermination de prédiction ; V<sub>m</sub> : vitesse de croissance moyenne (mm/an) ; T : temps de croissance (années).</td></tr></tfoot></table>

D'après les résultats obtenus, basés sur la signification des paramètres estimés, les valeurs des différents critères statistiques caractérisant la capacité d'ajustement, la capacité de prédiction et le réalisme biologique de chaque fonction de croissance étudiée (tableau VI), la fonction de Gompertz (F2) semble être la mieux indiquée pour prédire la croissance en épaisseur du liège en fonction de son âge. Ce choix est basé essentiellement sur la valeur de l'asymptote horizontale (valeur du paramètre $a_1$) et le réalisme biologique caractérisé par la vitesse moyenne annuelle de croissance et le temps sur lequel s'étale la période de cette croissance.

## Modèle dynamique de croissance en épaisseur du liège

La fonction retenue pour prédire la croissance en épaisseur du liège en fonction de son âge (fonction F2) est une fonction statique puisqu'elle ne tient compte que de l'âge du liège. Afin de rendre cette fonction dynamique (dépendant des caractéristiques dendrométriques de l'arbre de chêne-liège), le premier paramètre ($a_1$) relatif à l'asymptote horizontale, correspondant à la valeur maximale vers laquelle tend l'épaisseur du liège quand l'âge tend vers l'infini, a été exprimé à l'aide d'une fonction linéaire du diamètre à hauteur de poitrine (DBH). Le choix d'exprimer le paramètre $a_1$ à l'aide d'une relation linéaire en fonction du

DBH est basé, d'une part, sur la tendance linéaire observée entre l'épaisseur et le DBH (figure 10) et, d'autre part, sur le fait que l'asymptote horizontale ne doit pas être fixe pour tous les arbres et doit dépendre des caractéristiques dendrométriques de chaque arbre qui varient en fonction de la fertilité des stations. En plus, les résultats des essais de démasclage ont montré que le diamètre de l'arbre influe significativement sur la croissance en épaisseur du liège.

<figure><img src="image_15.png" /><figcaption>Figure 10.</figcaption></figure>

Relation entre l'épaisseur de liège et le DBH (données des essais de démasclage et celles des deux inventaires de 2016 et 2018).
## Tableau VII.

Estimation des paramètres et caractéristiques statistiques de l'équation 13.

<table><thead><tr><th rowspan="2">Fonction</th><th rowspan="2">Par.</th><th rowspan="2">Estim.</th><th rowspan="2">p &gt; |t|</th><th colspan="3">Capacité d'ajustement</th><th colspan="3">Capacité de prédiction</th></tr><tr><th>Biais</th><th>RMSE</th><th>R<sup>2</sup><sub>aj</sub></th><th>Biais<sub>p</sub></th><th>RMSE<sub>p</sub></th><th>R<sup>2</sup><sub>p</sub></th></tr></thead><tbody><tr><td rowspan="5">F2</td><td>a<sub>2</sub></td><td>7,495</td><td>&lt; 0,0001</td><td rowspan="5">-0,0569</td><td rowspan="5">2,529</td><td rowspan="5">0,9569</td><td rowspan="5">-0,0570</td><td rowspan="5">2,536</td><td rowspan="5">0,9568</td></tr><tr><td>a<sub>3</sub></td><td>7,915</td><td>&lt; 0,0001</td></tr><tr><td>b<sub>0</sub></td><td>42,283</td><td>&lt; 0,0001</td></tr><tr><td>b<sub>1</sub></td><td>0,703</td><td>&lt; 0,0001</td></tr><tr><td>&rho;<sub>1</sub></td><td>1,015</td><td>&lt; 0,0001</td></tr><tr><td>&rho;<sub>2</sub></td><td>0,815</td><td>&lt; 0,0001</td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody><tfoot><tr><td colspan="10">Par. : paramètre ; Estim. : estimation ; RMSE : écart-type résiduel ; R<sup>2</sup><sub>aj</sub> : coefficient de détermination ajusté ; Biais<sub>p</sub> : biais PRESS ; RMSE<sub>p</sub> : écart-type résiduel PRESS ; R<sup>2</sup><sub>p</sub> : coefficient de détermination de prédiction ; V<sub>m</sub> : vitesse de croissance moyenne (mm/an) ; T : temps de croissance (années).</td></tr></tfoot></table>

Après avoir exprimé le paramètre $a_1$ de la fonction F2 à l'aide d'une relation linéaire en fonction du DBH (équation 13), cette dernière a été réajustée à nouveau sur les mêmes données utilisées pour les ajustements du tableau VI. Les résultats de ce nouvel ajustement figurent au tableau VII :

$$EP_{BH} = (b_0 + b_1 \times DBH) \times \exp\left(-\exp\left(-\frac{t - a_2}{a_3}\right)\right) \quad (\text{équation 13}).$$

Ainsi, pour un arbre donné, le modèle dynamique qui permet de prédire la croissance de son épaisseur du liège en tenant compte de son diamètre (DBH) et de l'âge du liège est le suivant :

$$EP_{BH} = (42,283 + 0,703 \times DBH) \times \exp\left(-\exp\left(-\frac{t - 7,495}{7,915}\right)\right) \quad (\text{équation 14})$$

avec $R_{aj}^2 = 0,957$ et $RMSE = 2,529$ mm.

## Caractérisation des coefficients de démasclage réellement appliqués par l'administration forestière tunisienne (inventaires 2016 et 2017)

En profitant des deux inventaires de démasclage réalisés durant les deux campagnes de déliégeage de 2016 et 2017 répartis sur 14 chantiers de démasclage, les coefficients de démasclage réellement pratiqués ont été calculés (tableau VIII).

Les résultats obtenus de ces deux inventaires (tableau VIII et figure 11) montrent ce qui suit :

* les coefficients de démasclage minimaux appliqués sont presque constants et ne dépendent pas de la dimension des tiges ;

* les coefficients maximaux et moyens diminuent avec l'augmentation de la dimension des tiges ;

* la hauteur moyenne de la partie démasclée évolue très lentement en fonction de la dimension des tiges et semble se stabiliser à partir d'un diamètre supérieur à 50 cm.

Il ressort de ces résultats qu'en moyenne les ouvriers chargés de l'opération de démasclage pratiquent le déliégeage à des hauteurs qui ne varient pas beaucoup d'un arbre à l'autre quelle que soit la dimension de ce dernier.

## Tableau VIII.

Statistiques descriptives des coefficients de démasclage appliqués par l'administration forestière tunisienne (résultats des inventaires réalisés en 2016 et 2017).

<table><thead><tr><th rowspan="2">Classes de diamètre (cm)</th><th rowspan="2">Nombre d'arbres</th><th colspan="4">Hauteurs de démasclage appliquées (m)</th><th colspan="4">Coefficients de démasclage appliqués</th></tr><tr><th>Minimum</th><th>Maximum</th><th>Moyenne</th><th>CV%</th><th>Minimum</th><th>Maximum</th><th>Moyenne</th><th>CV%</th></tr></thead><tbody><tr><td>25</td><td>200</td><td>0,70</td><td>2,60</td><td>1,68</td><td>23,43</td><td>0,89</td><td>3,55</td><td>2,04</td><td>23,17</td></tr><tr><td>35</td><td>288</td><td>0,70</td><td>3,60</td><td>2,01</td><td>22,85</td><td>0,70</td><td>3,47</td><td>1,81</td><td>23,15</td></tr><tr><td>45</td><td>168</td><td>0,90</td><td>3,90</td><td>2,14</td><td>24,81</td><td>0,57</td><td>2,70</td><td>1,53</td><td>25,24</td></tr><tr><td>55</td><td>68</td><td>1,20</td><td>3,70</td><td>2,31</td><td>20,75</td><td>0,76</td><td>2,31</td><td>1,35</td><td>22,37</td></tr><tr><td>65 et plus</td><td>28</td><td>1,40</td><td>3,45</td><td>2,34</td><td>18,88</td><td>0,70</td><td>1,57</td><td>1,09</td><td>19,73</td></tr></tbody><tfoot><tr><td>Total</td><td>752</td><td>0,70</td><td>3,9</td><td>1,99</td><td>25,43</td><td>0,57</td><td>3,55</td><td>1,74</td><td>27,85</td></tr></tfoot></table><p>CV% : coefficient de variation.</p>
Tableau IX.

Épaisseur du liège observée par classe de grosseur des tiges
(résultats des inventaires réalisés en 2016 et 2017).

<table><thead><tr><th rowspan="2">Classes de diamètre (cm)</th><th rowspan="2">Nombre d'arbres</th><th colspan="3">Épaisseur du liège (mm)</th><th rowspan="2">CV%</th></tr><tr><th>Minimum</th><th>Maximum</th><th>Moyenne</th></tr></thead><tbody><tr><td>25</td><td>200</td><td>10,5</td><td>49,0</td><td>27,28</td><td>26,08</td></tr><tr><td>35</td><td>288</td><td>15,0</td><td>55,0</td><td>30,56</td><td>25,24</td></tr><tr><td>45</td><td>168</td><td>12,5</td><td>72,5</td><td>33,10</td><td>27,92</td></tr><tr><td>55</td><td>68</td><td>17,5</td><td>65,0</td><td>34,45</td><td>27,97</td></tr><tr><td>65 et plus</td><td>28</td><td>20,0</td><td>60,0</td><td>35,46</td><td>28,28</td></tr></tbody><tfoot><tr><th>Total</th><td>752</td><td>10,5</td><td>72,5</td><td>30,95</td><td>27,63</td></tr></tfoot></table>

CV% : coefficient de variation.

<figure><img src="image_16.png" /><figcaption>Figure 11.</figcaption></figure>

Valeurs observées des coefficients de démasclage (min, max et moyenne) et des hauteurs moyennes de démasclage appliquées (résultats des inventaires réalisés en 2016 et 2017 sur 752 arbres).

<figure><img src="image_17.png" /><figcaption>Figure 12.</figcaption></figure>

Biomasse de liège estimée (kg/arbre) en fonction des classes de diamètres observés dans les 40 placettes inventoriées en 2017.

Cette pratique a comme conséquence de surexploiter les arbres de petites dimensions où le coefficient de démasclage peut dépasser 3,5 fois la circonférence du tronc, et de sous-exploiter les gros arbres où ce même coefficient de démasclage peut devenir inférieur à l'unité. Il est donc clair, d'après ces résultats, qu'il y a un manque à gagner en ce qui concerne la récolte du liège. Celui-ci peut être très important, surtout pour les arbres de grandes dimensions qui sont capables de donner plus de liège.

Le tableau IX montre comment évolue l'épaisseur du liège à la récolte (rotation de 12 ans) en fonction du diamètre du tronc dans les différents chantiers de démasclage des deux inventaires de 2016 et 2017. Ces résultats corroborent les résultats similaires obtenus dans les essais de démasclage.

Pour estimer le rendement du liège récolté en se basant sur les coefficients réellement appliqués et ceux testés dans le cadre des essais de démasclage, les données issues des mesures des placettes de l'inventaire de 2017 ont été utilisées. Le tableau X présente les caractéristiques moyennes des 40 placettes mesurées en 2017 et les rendements du liège estimés (équation 12) dans le cas des coefficients appliqués réellement et ceux relatifs aux trois coefficients de démasclage testés, soit $K_1 = 2$, $K_2 = 2,5$ et $K_3 = 3$ fois la circonférence à hauteur de poitrine. Le rendement en liège est calculé sur la base des arbres démasclés par placette, rapporté à l'hectare.

En se basant sur les valeurs des rendements moyens estimés par hectare (tableau X), nous constatons que les rendements des trois coefficients $K_1$ (2 x CBH), $K_2$ (2,5 x CBH) et $K_3$ (3 x CBH) dépassent en moyenne ceux des coefficients réellement pratiqués de 14 %, 32 % et 50 % respectivement.

La figure 12, qui présente l'évolution de la biomasse moyenne de liège estimée par arbre en fonction des classes de diamètres observés dans les 40 placettes inventoriées en 2017, montre l'importance des écarts entre les coefficients de démasclage réellement pratiqués et ceux correspondant à 2, 2,5 et 3 fois la circonférence du tronc à hauteur de poitrine, surtout pour les arbres de grandes dimensions. Ces différences présentent un manque à gagner non négligeable pour l'administration forestière tunisienne.
# Discussion

## Effet de l'intensité de démasclage et de la grosseur des tiges sur la croissance en épaisseur et la production du liège

Les résultats obtenus dans le cadre de cet essai mulsites d'écorçage de liège ont montré que les deux coefficients de démasclage les plus élevés appliqués ($K_2 = 2,5$ et $K_3 = 3$) n'ont pas affecté négativement et d'une façon significative la croissance en épaisseur du liège quels que soient le site expérimental et la classe de grosseur des tiges et ceci par comparaison au coefficient classique ($K_1 = 2$) théoriquement appliqué par les gestionnaires forestiers en Tunisie. La grosseur des tiges affecte positivement et significativement la croissance du liège quels que soient le site expérimental et le coefficient de démasclage appliqué. Dans une étude similaire au Portugal, Paulo et Tomé (2017) ont étudié l'effet de trois coefficients de démasclage ($K_1 = 1,5$, $K_2 = 2$, $K_3 = 2,5$) appliqués durant la première extraction du liège sur la croissance future de l'épaisseur du liège dans deux régions de différents potentiels de productivité (faible et haute productivité). Les résultats obtenus ont montré que l'augmentation de l'intensité de démasclage a eu un petit effet négatif significatif sur la croissance en épaisseur du liège dans la région de faible productivité, tandis qu'aucun effet significatif n'a été détecté dans la région de haute productivité. Dans une autre étude visant à mesurer l'impact de l'intensité de démasclage sur la qualité du liège produit, El Alami et al. (2019) ont testé au Maroc l'effet de quatre intensités de démasclage ($K_1 = 1$, $K_2 = 1,5$, $K_3 = 2$, $K_4 = 2,5$) sur les propriétés physiques et mécaniques du liège formé (la densité, la compression et le module d'élasticité). Les résultats obtenus ont montré qu'aucune différence significative n'a été détectée entre les quatre coefficients de démasclage appliqués et ceci pour les trois caractéristiques mesurées.

Contrairement à l'épaisseur, le poids du liège récolté est affecté positivement et significativement par le coefficient de démasclage et par la dimension du tronc de l'arbre. Les écarts enregistrés entre les deux classes de circonférence extrêmes varient selon les sites entre environ 10 et 12 mm pour l'épaisseur moyenne du liège et entre 48 et 57 kg/arbre pour la biomasse moyenne du liège en faveur des arbres les plus gros. Les trois coefficients de démasclage testés donnent en moyenne des valeurs semblables du point de vue de l'épaisseur du liège mais présentent au contraire des écarts importants du point de vue de la biomasse, surtout entre les deux coefficients extrêmes ($K_1 = 2$ et $K_3 = 3$) qui varient en moyenne de 8 à 19 kg/arbre selon le site expérimental.

En ce qui concerne les sites expérimentaux, le site de Ain-Draham et celui de Chihia représentent les deux sites extrêmes du point de vue de la croissance

et de la production du liège. Les écarts enregistrés entre ces deux sites sont en moyenne de l'ordre de 7 mm pour l'épaisseur du liège et de 16 kg/arbre pour la production de liège, en faveur du site de Ain-Draham.

## Durée de rotation

En se basant sur les valeurs moyennes de l'épaisseur du liège obtenues après une rotation de 12 ans dans les trois sites expérimentaux et leurs limites de confiance avec un degré de confiance de 95 %, et en les comparant aux limites optimales de l'épaisseur du liège demandées par l'industrie du liège pour la production de bouchons naturels, la durée de rotation peut être réduite à 9-10 ans pour le site de Ain-Draham, 10-11 ans pour le site de Tébénia, et maintenue à 12 ans pour le site de Chihia. Les durées de rotation déterminées dans le cadre de cette étude corroborent celles adoptées en Espagne et au Portugal qui varient de 9 à 14 ans en fonction de la fertilité du site et la densité du peuplement (Montero et Cañellas, 2003 ; Pereira et Tomé, 2004). Il est à signaler que, bien que la durée minimum de rotation acceptée selon les recommandations officielles soit de 9 ans, des durées plus courtes sont autorisées dans certaines régions d'Espagne (Pasalodos-Tato et al., 2018a).

## Modèle dynamique de croissance en épaisseur du liège

Le modèle de croissance en épaisseur du liège d'arbre individuel a été développé à partir des mesures répétées effectuées annuellement durant les 12 années de croissance sur les arbres des trois sites expérimentaux. Il s'agit d'un modèle non linéaire ajusté sur les données de 217 arbres avec une précision d'environ 96 % et qui permet de prédire la croissance en épaisseur du liège en tenant compte du diamètre du tronc de l'arbre de chêne-liège et de l'âge du liège. Sánchez-González et al. (2008) ont développé un modèle similaire selon l'approche des équations en différence algébrique généralisée (GADA). Ajusté sur des données de 432 arbres provenant de la forêt de chêne-liège du parc naturel « Los Alcornocales » situé dans l'extrême Sud-Ouest

<caption>Tableau X.</caption>

Caractéristiques moyennes des 40 placettes inventoriées en 2017.

<table><thead><tr><th>Variables</th><th></th><th>Minimum</th><th>Maximum</th><th>Moyenne</th><th>CV%</th></tr></thead><tbody><tr><td>Densité observée (arbres/ha)</td><td></td><td>76</td><td>631</td><td>271</td><td>46,14</td></tr><tr><td>Densité démasclée (arbres/ha)</td><td></td><td>68</td><td>568</td><td>235</td><td>42,83</td></tr><tr><td>DBH (cm)</td><td></td><td>24,92</td><td>61,00</td><td>41,04</td><td>17,740</td></tr><tr><td>Hauteur démasclée (m)</td><td></td><td>1,53</td><td>3,060</td><td>2,16</td><td>14,75</td></tr><tr><td>Épaisseur du liège (mm)</td><td></td><td>23,10</td><td>46,29</td><td>29,51</td><td>16,46</td></tr><tr><td>Coefficients de démasclage appliqués</td><td></td><td>1,25</td><td>2,18</td><td>1,77</td><td>13,05</td></tr><tr><td>Coefficients appliqués</td><td></td><td>15,560</td><td>113,617</td><td>46,368</td><td>46,22</td></tr><tr><td>Rendement</td><td>K1 : (2 x CBH)</td><td>15,809</td><td>135,166</td><td>52,666</td><td>48,58</td></tr><tr><td>du liège</td><td>K2 : (2,5 x CBH)</td><td>18,433</td><td>157,594</td><td>61,405</td><td>48,58</td></tr><tr><td>(Qx/ha)</td><td>K3 : (3 x CBH)</td><td>20,896</td><td>178,656</td><td>69,612</td><td>48,58</td></tr></tbody><tfoot><tr><td colspan="6">CBH : circonférence à hauteur de poitrine ou à 1,30 m<br>de hauteur ; CV% : coefficient de variation.</td></tr></tfoot></table>
de l'Espagne, ce modèle qui explique 98 % de la variabilité totale permet de prédire, à l'échelle de l'arbre, l'épaisseur du liège à n'importe quel âge connaissant la même épaisseur mesurée à n'importe quel autre âge. Avec une précision de 96 %, comparable à celle obtenue par le modèle développé en Espagne, le modèle développé dans le cadre de notre étude, qui nécessite comme variables exogènes le diamètre du tronc de l'arbre de chêne-liège et l'âge du liège au lieu de l'épaisseur du liège mesurée à un âge donné, est plus pratique du point de vue de l'utilisation. Vu la faible variation du diamètre des arbres d'une année à une autre, notre modèle peut être utilisé non seulement sur des données issues des mesures directes dans des peuplements de chêne-liège mais aussi sur celles issues des inventaires forestiers.

## Modèle dynamique de production du liège

Le modèle de prédiction de la biomasse du liège à l'échelle de l'arbre a été développé à partir des données de 589 arbres. Avec une précision de 93 %, le modèle é laboré permet de prédire la biomasse du liège d'un arbre en fonction de son diamètre du tronc sur écorce (DBH), de l'épaisseur du liège mesuré à hauteur d'homme (EP<sub>BH</sub>) et de la hauteur de la partie démasclée du tronc (H<sub>DEM</sub>). La précision relativement élevée (93 %) de notre modèle corrobore les résultats obtenus par Vázquez et Pereira (2005) qui ont démontré que les caractéristiques dendrométriques à l'échelle de l'arbre permettent d'expliquer à elles seules 80 à 90 % de la variabilité totale du poids du liège récolté.

Différentes études ont utilisé, sous diverses formes, les mêmes variables explicatives pour prédire le poids du liège à l'échelle de l'arbre. Montero (1987) a ajusté un modèle linéaire sans ordonnée à l'origine qui fait intervenir le produit de la circonférence du tronc sur écorce et la hauteur de la partie démasclée (CBH x H<sub>DEM</sub>) comme variable explicative. Sánchez-González et al. (2007), puis Pasalodos-Tato et al. (2018b) ont proposé à leur tour une formule simple basée sur le calcul du volume de liège (CBH x EP<sub>BH</sub> x H<sub>DEM</sub>) puis sa multiplication par la densité du liège. Dans une autre étude, Ferreira et Carvalho Oliveira (1991) ont opté pour un simple modèle linéaire pour prédire le poids du liège en utilisant séparément les deux variables explicatives CBH et H<sub>DEM</sub>. Finalement, dans le cadre d'une étude sur le potentiel de production d'une petite forêt locale de chêne-liège en Tunisie, Stiti et al. (2019) ont ajusté sur des données issues de 44 arbres localisés dans la forêt de Ain-Snoussi un modèle non linéaire de prédiction du liège à l'échelle de l'arbre avec une précision de l'ordre de 77 %, qui fait intervenir comme variables explicatives le produit (CBH x H<sub>DEM</sub>) et l'âge du liège à la place de son épaisseur.

## Coefficients de démasclage réellement appliqués en Tunisie et leur impact sur la production nationale de liège

D'après les résultats des deux inventaires réalisés durant les deux campagnes de démasclage 2016 et 2017 dans différents chantiers de récolte de liège, les coefficients de démasclage appliqués par les ouvriers ont varié de 0,6 à

3,6 fois la circonférence du tronc à hauteur de poitrine avec une valeur moyenne de 1,7 et un coefficient de variation de 28 %. Une étude plus approfondie, basée sur la répartition des arbres démasclés par classes de grosseur des tiges, a montré que les coefficients de démasclage minimaux appliqués sont presque toujours les mêmes quelle que soit la taille de l'arbre, que les coefficients maximaux et moyens diminuent avec l'augmentation de la dimension des arbres et que la hauteur moyenne de la partie démasclée évolue très lentement en fonction de la dimension des arbres et semble se stabiliser à partir de certains diamètres des troncs, au-delà de 50 cm. Ces résultats montrent que l'opération de démasclage en Tunisie s'effectue en général à des hauteurs plus ou moins constantes qui ne dépendent pas de la dimension des arbres démasclés. Le fait de démascler les arbres de chêne-liège jusqu'à une hauteur constante sans tenir compte de leur dimension engendre une surexploitation des arbres de petites dimensions où le coefficient de démasclage peut dépasser 3,5 fois la circonférence du tronc, et une sous-exploitation des gros arbres où le coefficient de démasclage peut descendre sous l'unité, ce qui représente un manque à gagner qui peut être très important dans le cas des gros arbres. Par comparaison aux coefficients de démasclage réellement pratiqués, le manque à gagner du point de vue de la récolte de liège a été respectivement estimé dans le cadre de cette étude à 14 %, 32 % et 50 %, pour les trois coefficients étudiés, soit 2, 2,5 et 3 fois la circonférence du tronc à hauteur de poitrine.

## Conclusion

Les résultats obtenus et les modèles développés dans le cadre de cette étude améliorent nos connaissances sur la croissance et la production du liège et constituent, de nos jours, des outils indispensables pour la gestion et la conservation des forêts de chêne-liège en Tunisie.

L'absence d'effet négatif significatif de l'augmentation de l'intensité de démasclage sur la croissance en épaisseur du liège démontre qu'il est possible d'augmenter la production de liège sans nuire à la vigueur de l'arbre de chêne-liège, du moins dans les stations à haut potentiel de production. Afin de vérifier sur le long terme l'effet de ces coefficients de démasclage plus sévères sur la croissance du liège et la vigueur des arbres démasclés, les mêmes essais sont suivis durant la deuxième rotation et les résultats seront disponibles dès la fin de 2020 et en 2021. Toutefois et en attendant les résultats de la deuxième rotation, même si aucune différence significative du point de vue de la croissance en épaisseur du liège n'a été détectée durant la première rotation entre les trois coefficients de démasclage appliqués dans les trois sites expérimentaux de notre étude, le coefficient de démasclage le plus sévère testé (K<sub>3</sub> = 3) doit être utilisé avec prudence dans les stations les moins fertiles. D'après les résultats des deux inventaires réalisés dans différents chantiers de démasclage en Tunisie, où le coefficient de démasclage réellement pratiqué était en moyenne de l'ordre de 1,7 et diminue pour les gros arbres, les pertes dans la production de liège sont respectivement estimées à 14 %, 32 % et 50 % pour les trois coefficients
étudiés (2, 2,5 et 3). Ces pertes peuvent être considérables surtout dans les stations à forte productivité, caractérisées par des arbres de gros diamètres.

Le modèle de prédiction de la croissance en épaisseur du liège qui utilise comme prédicteurs (variables exogènes) le diamètre du tronc de l'arbre de chêne-liège et l'âge du liège est un modèle dynamique de forte précision ($R_{aj}^2$ = 96%) et facile à utiliser. Connaissant le diamètre d'un arbre, ce modèle permet aux gestionnaires forestiers d'estimer l'évolution de la croissance en épaisseur du liège année par année et de déterminer par conséquent en avance la durée nécessaire de rotation pour atteindre une épaisseur commerciale donnée. Ce modèle peut être utilisé également pour classer les différents peuplements de chêne-liège selon leur potentiel de croissance ou de production du liège.

Le modèle d'estimation de la biomasse du liège développé est également un modèle de forte précision ($R_{aj}^2$ = 93%) et utilise des prédicteurs fiables et faciles à mesurer (diamètre du tronc, épaisseur du liège à hauteur de poitrine, hauteur de la partie démasclée). La connaissance préalable de la production attendue d'une forêt donnée de chêne-liège représente sans doute un atout fort intéressant pour les gestionnaires forestiers pour bien réussir la planification et la programmation de la campagne de récolte du liège (main-d'œuvre, transport, stockage, etc.). Comme pour les autres productions forestières vendues en adjudications (cas de la vente du bois sur pied en Tunisie), le modèle d'estimation de la biomasse du liège sur pied permet aux gestionnaires forestiers d'estimer préalablement la valeur de la récolte sur des bases scientifiques.

L'utilisation conjointe des deux modèles de prédiction de la croissance et la production du liège élaborés dans le cadre de cette étude permettra aux gestionnaires forestiers de déterminer le moment opportun de la récolte du liège et de fixer par conséquent la durée de rotation (à l'échelle de la parcelle, du peuplement ou de la forêt), d'estimer la biomasse du liège à récolter et de répartir cette dernière par classes d'épaisseur.

Finalement et comme perspectives de cette étude, trois autres modèles liés à la gestion et l'aménagement des forêts de chêne-liège en Tunisie méritent d'être développés. Il s'agit d'un modèle d'indices de sites, un modèle de distribution des tiges par classes de grosseur et un modèle d'estimation de la qualité du liège sur pied. En plus des deux modèles déjà élaborés, une fois les trois nouveaux modèles développés, les gestionnaires forestiers disposeront de tous les outils nécessaires à l'optimisation de la gestion et la conservation des forêts de chêne-liège en Tunisie.

## Références

Abouba M., 2014. Modélisation de la croissance en diamètre des espèces de pin du site de Haouaria en Tunisie. PFE master en biostatistique, Gembloux Agro-Bio Tech, Université de Liège, Belgique, 81 p.

Bendanoun M., 1996. La Mamora, un patrimoine national vital pour l'avenir de Salé et de la région. Journée d'étude « Sauvons la Mamora », Salé, 25 mars 1996, ENFI-MEABR, Maroc, 31 p.

Bonnet S., 2019. Le chêne-liège, un arbre généreux. The forest time. https://www.the-forest-time.com/fr/guides-des-pays-et-regions/portugal/le-chene-liège-un-arbre-genereux-090976213

Bouchafra A., Fraval A., 1991. Présentation du chêne-liège et de la subéraie. In : Villemant C., Fraval A. La faune du chêne-liège. Rabat, Maroc, Actes Éditions, 1-26.

Costa A., Barbosa I., Pestana M., Miguel C., 2020. Modelling bark thickness variation in stems of cork oak in south western Portugal. European Journal of Forest Research, 139: 611-625. https://doi.org/10.1007/s10342-020-01273-9

Dagnelie P., 2011. Statistique théorique et appliquée. Tome 2. Inférence statistique à une et à deux dimensions. Bruxelles, Belgique, De Boeck, 736 p.

Debouche C., 1979. Présentation coordonnée de différents modèles de croissance. Revue de Statistique Appliquée, 27 (4) : 5-22.

Dehane B., Ghefar M., 2017. Étude de la variabilité de croissance du liège dans le Nord Algérien. Geo-Eco-Trop, 41 (1) : 39-54.

DGF, 2005. Deuxième Inventaire national forestier et pastoral : Résultats de l'inventaire forestier sur l'échelon de Jendouba (Gouvernorats de Jendouba, Béja et Bizerte). Direction générale des forêts, Tunisie, 129 p.

Diéguez-Aranda U., Grandas-Arias J. A., Alvarez-Gonzalez J. G., Gadow K., 2006. Site quality curves for birch stands in north-western Spain. Silva Fennica, 40 (4): 631-644. https://doi.org/10.14214/sf.319

El Alami A., Hakam A., Kabouchi B., Famiri A., Bouchti A., Zianai M., 2019. Effect of the stripping intensity on physical and mechanical properties of cork. ARPN Journal of Engineering and Applied Sciences, 14 (4): 921-926.

Ferreira M. C., Carvalho Oliveira A. M., 1991. Modelling cork oak production in Portugal. Agroforestry Systems, 16: 41-54. https://doi.org/10.1007/bf00053196

Garchi S., Sghaier T., 2002. Effet de l'intensité de démasclage sur la formation du liège. Annales de la Recherche Forestière au Maroc, 35 : 88-95.

Garchi S., Sghaier T., 2008. Étude de l'effet du coefficient de démasclage et de la grosseur des tiges sur la croissance en épaisseur du liège : résultats d'un essai multi-sites après 10 ans de croissance. Revue Internationale de Géologie, de Géographie et d'Écologie Tropicales, 32 : 21-28.

Ghalem A., 2006. Étude typologique, stratégie de réhabilitation et réaction du milieu après incendie : cas de la subéraie de Hafir et Zarieffet (Wilaya de Tlemcen). PFE cycle ingénieur, Faculté des sciences, Université de Tlemcen, Algérie, 94 p.

Kouba Y., 2005. Contribution à l'étude de la filière liège en Algérie. PFE cycle ingénieur, Université Mouloud Memmari de Tizi-Ouzou, Algérie, 71 p.

Messaoudene M., Metna B., 2010. Étude de la qualité du liège des subéraies orientales de Tizi-Ouzou. Integrated Protection in Oak Forests/Wet Palearctic Regional Section (IOBC/WPRS) Bulletin, 57 : 33-40.

Montero G., 1987. Modelos para cuantificar la producción de corcho en Alcornocales en función de la calidad de la estación y de los tratamientos selvícolas, Ph.D. Thesis, INIA, Madrid,
España, 277 p.

Montero G., Cañellas I., 2003. The silviculture of cork oak woodlands in Spain. Silva Lusitana, 11: 1-19.

Naggar O., Bouhraoua R. T., Torres E., Zedek M., 2019. Étude de l'influence des facteurs du milieu sur la croissance et la production du liège de la subéraie du Parc national de Theniet El Had (Tissemsilt, Algérie). Bois et Forêts des Tropiques, 342 (4): 41-54. https://doi.org/10.19182/bft2019.342.a31793

Natividade V. J., 1956. Subériculture. Nancy, France, École nationale des eaux et des forêts, 281 p.

Nouri M., 2009. Facteurs pédo-climatiques et évolution de la subéraie tunisienne : propriétés physico-chimiques et hydro-dynamiques des sols dans les forêts de chêne-liège (Quercus suber L.). Thèse de doctorat en sciences agronomiques, INAT, Tunisie, 173 p.

Pasalodos-Tato M., Pukkala T., Cañellas I., Sánchez-González M., 2018a. Optimizing the debarking and cutting schedule of cork oak stands. Annals of Forest Science, 75 (61). https://doi.org/10.1007/s13595-018-0732-8

Pasalodos-Tato M., Alberdi I., Cañellas I., Sánchez-González M., 2018b. Towards assessment of cork production through National Forest Inventories. Forestry, 91: 110-120. https://doi.org/10.1093/forestry/cpx036

Paulo J. A., Tomé M., 2010. Predicting mature cork biomass with *t* years of growth from one measurement taken at any other age. Forest Ecology and Management, 259 (10): 1993-2005. https://doi.org/10.1016/j.foreco.2010.02.010

Paulo J. A., Tomé M., 2017. Does debarking intensity during the first cork extraction affect future cork thickness? Annals of Forest Science, 74 (66). https://doi.org/10.1007/s13595-017-0662-x

Pereira, H., 2011 Cork: Biology, Production and Uses. Elsevier, 336 p.

Pereira H., Tomé M., 2004. Non-wood products: cork oak. In: Burley J., Evans J., Youngquist J. A. (eds). Encyclopedia of forest sciences. Oxford, United Kingdom, Elsevier, 613-620. https://doi.org/10.1016/b0-12-145160-7/00278-7

Ribeiro F., Tomé M., 2002. Cork weight prediction at tree level. Forest Ecology and Management, 171 (3): 231-241. https://doi.org/10.1016/s0378-1127(01)00780-0

Ryan T. A., Joiner B. L., 1976. Normal probability plots and tests for normality. State College, PA, USA, Pennsylvania State University, 12 p.

SAS Institute Inc., 2004b. SAS/ETS 9.1 User's Guide. Cary, NC, USA, SAS Institute Inc, 5136 p. https://support.sas.com/documentation/onlinedoc/91pdf/sasdoc_91/stat_ug_7313.pdf

Sánchez-González M., Calama R., Cañellas I., Montero G., 2007. Variables influencing cork thickness in Spanish cork oak forests: A modelling approach. Annals of Forest Science, 64: 301-312. https://doi.org/10.1051/forest:2007007

Sánchez-González M., Cañellas I., Montero G., 2008. Base-age invariant cork growth model for Spanish cork oak (Quercus suber L.) forests. European Journal of Forestry Research, 127: 173-182. https://doi.org/10.1007/s10342-007-0192-4

Sánchez-González M., Tomé M., Montero G., 2005. Modelling height and diameter growth of dominant cork oak trees in Spain. Annals of Forest Science, 62: 633-643. https://doi.org/10.1051/forest:2005065

Sebei H., Albouchi A., Rapp M., El Aouni M. H., 2004. Productivité en biomasse du chêne-liège dans une séquence de dégradation de la subéraie à Cytise de Kroumirie (Tunisie). Annals of Forest Science, 61: 347-361. https://doi.org/10.1051/forest:2004028

Sghaier T., Garchi S., Ammari Y., Santiago R., Rosello M. E., 2011. Estimation de la qualité du liège en pile dans le parc à liège de Bellif (Tunisie Septentrionale) selon la méthode de l'IPROCOR. Revue Internationale de Géologie, de Géographie et d'Écologie Tropicales, 35: 51-68.

Stiti B., Hasni H., Ben Rjab S., Khaldi A., 2014. Évaluation de l'épaisseur et de la qualité du liège de reproduction dans la forêt de Bellif, nord-ouest de la Tunisie. Integrated Protection in Oak Forests. IOBC-WPRS Bulletin, 101: 69-73.

Stiti B., Chaar H., Khaldi K., Henchi B., 2019. Analyse et modélisation de la production du liège en vue d'une gestion optimale : cas de la forêt de Ain Snoussi, Nord-Ouest de la Tunisie. Journal of New Sciences, Agriculture and Biotechnology, 65 (4): 4088-4095.

Vázquez J., Pereira H., 2005. Mixed models to estimate tree oven-dried cork weight in Central and Southern Portugal. Forest Ecology and Management, 213 (1-3): 117-132. https://doi.org/10.1016/j.foreco.2005.03.047

Veillon S., 1998. Guide technique de subériculture dans les Pyrénées-Orientales. Typologie de peuplements et étude préliminaire. Mémoire stage fin d'études FIF, ENGREF, France, 73 p.

### Sghaier et al. – Contribution des auteurs

<table><thead><tr><td>Rôle du contributeur</td><td>Noms des auteurs</td></tr></thead><tbody><tr><td>Conceptualisation</td><td>T. Sghaier, S. Garchi</td></tr><tr><td>Gestion des données</td><td>T. Sghaier, S. Garchi</td></tr><tr><td>Analyse formelle</td><td>T. Sghaier</td></tr><tr><td>Acquisition du financement</td><td>T. Sghaier, S. Garchi</td></tr><tr><td>Enquête et investigation</td><td>T. Sghaier, S. Garchi</td></tr><tr><td>Méthodologie</td><td>T. Sghaier</td></tr><tr><td>Gestion de projet</td><td>T. Sghaier, S. Garchi</td></tr><tr><td>Ressources</td><td>INRGREF</td></tr><tr><td>Logiciels</td><td>T. Sghaier</td></tr><tr><td>Supervision</td><td>T. Sghaier, S. Garchi</td></tr><tr><td>Validation</td><td>T. Sghaier, S. Garchi</td></tr><tr><td>Visualisation</td><td>T. Sghaier, T. Azizi</td></tr><tr><td>Écriture – Préparation de l'ébauche originale</td><td>T. Sghaier, T. Azizi</td></tr><tr><td>Écriture – Révision et édition</td><td>T. Sghaier, T. Azizi</td></tr></tbody></table>