Gauthier Ligot
email: gligot@uliege.be

Thibaut Gheysen

Jérôme Perin

Romain Candaele

Francois De Coligny

Alain Licoppe

Philippe Lejeune

J Ér Ôme Perin

Franc ¸ois De Coligny

From the simulation of forest plantation dynamics to the quantification of bark-stripping damage by ungulates

Introduction

Excessive ungulate densities are known to harm agricultural and forestry production [START_REF] Ward | Modelling the cost of roe deer browsing damage to forestry[END_REF]. In the case of timber production, ungulates are known to browse seedlings, rub their antlers on tree trunks or strip the bark from young trees. Among such damage, bark-stripping by red deer (Cervus elaphus L.) is especially detrimental to certain tree species such as Norway spruce (Picea abies(L.) Karst.) whose bark is particularly sought after by red deer and whose damaged stems are especially prone to fungal infection and decay [START_REF] Burn ¸eviča | Fungi inhabiting bark stripping wounds made by large game on stems of Picea abies (L.) Karst. in Latvia[END_REF][START_REF] Cukor | Scots pine (Pinus sylvestris L.) demonstrates a high resistance against bark stripping damage[END_REF]. Given that the density of large wild ungulates has substantially increased over the last century in temperate forests [START_REF] Milner | Temporal and spatial development of red deer harvesting in Europe : biological and cultural factors: red deer harvesting in Europe[END_REF][START_REF] Apollonio | European ungulates and their management in the 21st century[END_REF], the extent of this damage is of concern and deserves a rigorous quantification.

To date, the drivers and characteristics of bark-stripping damage on spruce trees have been relatively well described (Gill, 1992b;[START_REF] Vospernik | Probability of bark stripping damage by red deer (Cervus elaphus) in Austria[END_REF][START_REF] Candaele | Mitigation of bark stripping on spruce: the need for red deer population control[END_REF]Cukor et al., 2019a;[START_REF] Widén | Nutrient balancing or spring flush -What determines spruce bark stripping level by red deer? For[END_REF]. We also have valuable information on the subsequent decay in the stems and welldetailed models of forest stand dynamics and forest management [START_REF] Burn ¸eviča | Fungi inhabiting bark stripping wounds made by large game on stems of Picea abies (L.) Karst. in Latvia[END_REF][START_REF] Čermák | Stem decay by Stereum sanguinolentum after red deer damage in the Českomoravská vrchovina Highlands[END_REF][START_REF] Vacek | Bark stripping, the crucial factor affecting stem rot development and timber production of Norway spruce forests in Central Europe[END_REF]. However, combining all this knowledge to evaluate the impact of bark-stripping from planting through to final harvest at stand scale has rarely been undertaken and described ( Čermák et al., 2004a). Therefore, it remains difficult to assess how to adapt forest management to different levels of bark-stripping damage.

The prevalence of bark-stripping damage is known to increase with ungulate density and depends on tree size, hunting management and environmental characteristics (Gill, 1992b;[START_REF] Candaele | Mitigation of bark stripping on spruce: the need for red deer population control[END_REF][START_REF] Vospernik | Probability of bark stripping damage by red deer (Cervus elaphus) in Austria[END_REF][START_REF] Widén | Nutrient balancing or spring flush -What determines spruce bark stripping level by red deer? For[END_REF]. Red deer mostly strip the bark of relatively young spruce trees as soon as the stem is stiff and accessible while the bark is still thin [START_REF] Vospernik | Probability of bark stripping damage by red deer (Cervus elaphus) in Austria[END_REF][START_REF] Ligot | Modeling recent bark stripping by red deer (Cervus elaphus) in south Belgium coniferous stands[END_REF]. The prevalence of bark-stripping has been monitored using forest inventories in sensitive stands. In relatively large forest areas (over 1000 ha), typical values of annual bark-stripping rates range between 0 and 12% [START_REF] Candaele | Mitigation of bark stripping on spruce: the need for red deer population control[END_REF]. This rate is correlated to deer abundance but the relationship saturates beyond a certain ungulate density [START_REF] Ligot | Modeling recent bark stripping by red deer (Cervus elaphus) in south Belgium coniferous stands[END_REF][START_REF] Candaele | Mitigation of bark stripping on spruce: the need for red deer population control[END_REF]. The other factors that can affect the frequency of bark-stripping damage are winter harshness, forest composition, distance from roads and urban areas, canopy opening, availability of other food resources such as beechnuts and acorns, and topographic position [START_REF] Ligot | Modeling recent bark stripping by red deer (Cervus elaphus) in south Belgium coniferous stands[END_REF][START_REF] Candaele | Mitigation of bark stripping on spruce: the need for red deer population control[END_REF][START_REF] Vospernik | Probability of bark stripping damage by red deer (Cervus elaphus) in Austria[END_REF][START_REF] Widén | Nutrient balancing or spring flush -What determines spruce bark stripping level by red deer? For[END_REF][START_REF] Konôpka | Influence of tree species and size on bark browsing by large wild herbivores[END_REF]. The damage might occur both in winter and summer (Gill, 1992b;[START_REF] Widén | Nutrient balancing or spring flush -What determines spruce bark stripping level by red deer? For[END_REF]. Nevertheless, in summer, the bark is more loosely attached to the stem, and so deer typically remove longer strips of bark than in winter. According to [START_REF] Gheysen | A regional inventory and monitoring setup to evaluate bark peeling damage by red deer (Cervus elaphus) in coniferous plantations in Southern Belgium[END_REF], the strip can be twice as long in summer (30.2 cm in average) as in winter (15.6 cm).

After bark is damaged on spruce trees, wounded tissues are often infected by fungi, causing subsequent wood decay and discoloration [START_REF] Löffler | Zur ausbreitung der wundfäule in der fichte[END_REF]Čermák et al., 2004a;[START_REF] Mäkinen | Increment and decay in Norway spruce and Scots pine after artificial logging damage[END_REF][START_REF] Burn ¸eviča | Fungi inhabiting bark stripping wounds made by large game on stems of Picea abies (L.) Karst. in Latvia[END_REF]. Various decaycausing fungi have been identified; of which the most important one is Stereum sanguinolentum. After infection by S. sanguinolentum, a decay column forms around the wound [START_REF] Burn ¸eviča | Fungi inhabiting bark stripping wounds made by large game on stems of Picea abies (L.) Karst. in Latvia[END_REF]. The vertical spread of the decay depends on wound age, wound size, and various environmental factors. It varies greatly across studies [START_REF] Löffler | Zur ausbreitung der wundfäule in der fichte[END_REF]Čermák et al., 2004a;[START_REF] Čermák | Stem decay by Stereum sanguinolentum after red deer damage in the Českomoravská vrchovina Highlands[END_REF][START_REF] Mäkinen | Increment and decay in Norway spruce and Scots pine after artificial logging damage[END_REF]. The spread is generally fast during the first years after infection but subsides after 20-30 years ( Čermák et al., 2004a;[START_REF] Čermák | Stem decay by Stereum sanguinolentum after red deer damage in the Českomoravská vrchovina Highlands[END_REF]. The resulting height of the decay column can reach at most 4-6 m depending on the rate of decay spread and so on environmental conditions ( Čermák et al., 2004a;[START_REF] Čermák | Stem decay by Stereum sanguinolentum after red deer damage in the Českomoravská vrchovina Highlands[END_REF][START_REF] Vacek | Bark stripping, the crucial factor affecting stem rot development and timber production of Norway spruce forests in Central Europe[END_REF]. For example, [START_REF] Vacek | Bark stripping, the crucial factor affecting stem rot development and timber production of Norway spruce forests in Central Europe[END_REF] observed that stem rot reached maximum 4.5 m (mean = 1.9 m) at one site and maximum 6.0 m (mean = 3.1 m) at another site, in the Czech Republic.

A significant proportion of bark-stripped stems may contain decayed wood of low economic value. According to [START_REF] Vacek | Bark stripping, the crucial factor affecting stem rot development and timber production of Norway spruce forests in Central Europe[END_REF], the average proportion of decayed wood in damaged trees can range between 30% and 39%. Similar or even higher values were found in the Czech Republic [START_REF] Čermák | Stem decay by Stereum sanguinolentum after red deer damage in the Českomoravská vrchovina Highlands[END_REF]Čermák et al., 2004b) and in Belgium [START_REF] Heyninck | Pertes en volume et en valeur d'épicéas adultes touchés par des dégâts de cervidés[END_REF]. Often, the trunk sections that contain decayed wood cannot be sawn and can be used only by pulp and paper industries. The value of these sections containing decayed wood is thus very limited [START_REF] Heyninck | Pertes en volume et en valeur d'épicéas adultes touchés par des dégâts de cervidés[END_REF].

These losses can be attenuated by forest management. In western Europe, spruce plantations are usually regularly thinned before the final cut (clear-cut). The damaged trees can then be thinned preferentially over the healthy trees [START_REF] Vacek | Bark stripping, the crucial factor affecting stem rot development and timber production of Norway spruce forests in Central Europe[END_REF][START_REF] Heyninck | Pertes en volume et en valeur d'épicéas adultes touchés par des dégâts de cervidés[END_REF]. However, in most cases, it is unrealistic to fell every damaged tree in one or a few thinnings because bark damage are often aggregated within stands [START_REF] Hahn | Position, size, and spatial patterns of bark stripping wounds inflicted by red deer (Cervus elavus L.) on Norway spruce using generalized additive models in Austria[END_REF]. Forest managers can protect stands against ungulate damage: plantations can be fenced, or individual trees can be protected by scarifying their trunks to stimulate resin flows (Fig. 3), fitting plastic sleeves or nets around their trunks, wrapping branches around their trunks, or using chemicals [START_REF] Trout | Protection of trees from mammal damage[END_REF]. Individual protection can be installed, at about the time of the first thinning, on all trees or only on a selection of the most valuable trees that are expected to be harvested last (crop trees) [START_REF] Perin | De nouvelles normes sylvicoles pour les futaies pures équiennes d'épicéa et de douglas en appui à la gestion de la forêt publique en Wallonie[END_REF].

Assessing the financial impact of bark-stripping damage on timber production is complex because it requires taking into account multiple variables and processes on a long time scale (> 50 years). Consequently, most studies have focused on certain processes but rarely make a comprehensive assessments of financial losses [START_REF] Gill | The economic implications of deer damage: a review of current evidence[END_REF] . [START_REF] Gill | The economic implications of deer damage: a review of current evidence[END_REF] illustrated how such an assessment could be made with elementary models, but acknowledged that their approach was limited and that further models had to be developed to estimate yield losses, particularly for Norway spruce. In line with what has been proposed to assess the economic impacts of browsing by roe deer [START_REF] Ward | Modelling the cost of roe deer browsing damage to forestry[END_REF][START_REF] Rakotoarison | Analyse et modélisation de la gestion du grand gibier: cas de la région aquitaine[END_REF] or moose [START_REF] Wam | Taking timber browsing damage into account: A density dependant matrix model for the optimal harvest of moose in Scandinavia[END_REF], an assessment can be made by modeling forest dynamics and interactions with ungulate populations [START_REF] Weisberg | Predicting understorey vegetation cover from overstorey attributes in two temperate mountain forests[END_REF]. Once models are established, they can be used to predict the costs and revenues of different management strategies in ranging treatment and/or environmental conditions (e.g., ungulate density).

To fill this gap, we coupled a model of the forest dynamics and management of even-aged coniferous stands with models of bark-stripping damage. We then conducted a virtual experiment to address several questions: Does bark-stripping damage substantially harm the profitability of Norway spruce plantations? Should rotation be shortened in the stands that are highly impacted by bark-stripping damage? Is it more cost-effective to protect plantations with fences or to protect individual trees against bark-stripping at the time of the first thinning (i.e., after some damage has already occurred)?

Materials and methods

Study area

Different data sets were collected in Southern Belgium (Wallonia, approx. 50 • N, 5 • E, Fig. 1), mostly in the Ardenne region. The climate in Southern Belgium is warm temperate, with no dry season and a mild summer (classified as Cfb in the Köppen-Geiger scheme). Elevation varies between 20 m and 700 m and is generally highest in Ardenne. Annual rainfall and mean annual temperature range respectively between 800 mm.year -1 and 10.5 • C and 1400 mm.year -1 and 7.5 • C (averages computed for 1981-2010 by the Belgian Royal Meteorology Institute). Annual rainfall are positively correlated to elevation whereas mean annual temperature are negatively correlated to elevation. The climate in Ardenne is thus colder and wetter than in the other areas of Wallonia. About 60% of the forest area of Wallonia is located in Ardenne. Norway spruce plantations cover around 26% of the total productive forest area. In 2016, they constituted 40% of the standing stock and 50% of the total wood production. Most of these plantations of Norway spruce (92%) are located in Ardenne [START_REF] Lejeune | L'épicéa wallon: état de la ressource en 2021[END_REF]. Spruce plantations are mostly found on brown acidic soil with sometimes poor drainage. They are normally thinned every 5-10 years once the top height reaches 13-20 m and are clearccut when the top height reaches 30 m [START_REF] Perin | De nouvelles normes sylvicoles pour les futaies pures équiennes d'épicéa et de douglas en appui à la gestion de la forêt publique en Wallonie[END_REF]. Red deer, the roe deer (Capreolus capreolus) and the wild boar (Sus scrofa) are the three native ungulate species living in the study area. Within the study area, the estimated density of red deer populations ranged between 0 and 16.5 animals/km 2 and the reported number of shot red deer ranged between 0 and 6.7 animals/km 2 /year [START_REF] Candaele | Mitigation of bark stripping on spruce: the need for red deer population control[END_REF]. Roe deer or wild boar are also relatively abundant but their density cannot be accurately estimated [START_REF] Candaele | Mitigation of bark stripping on spruce: the need for red deer population control[END_REF].

Bark-stripping inventory data

Since 2003, regional monitoring of bark-stripping damage has been carried out in the study area. This monitoring system has been described in detail by [START_REF] Gheysen | A regional inventory and monitoring setup to evaluate bark peeling damage by red deer (Cervus elaphus) in coniferous plantations in Southern Belgium[END_REF] and [START_REF] Candaele | Mitigation of bark stripping on spruce: the need for red deer population control[END_REF]. The sampling units were installed on the nodes of a 200 × 200 m sampling grid that fall into spruce plantations of 8-36 years old (Fig. 2a). In each samplig unit, three circular plots of up to 10 m radius were laid out on a north-south transect (Fig. 2b). In each plot, the six tree closest to the plot center were measured. The sampling units had been measured annually between mid-April and mid-May since 2003. Between 2004 and 2015, over 4,000 plots and 60,000 trees were measured each year [START_REF] Candaele | Mitigation of bark stripping on spruce: the need for red deer population control[END_REF].

Models of forest dynamics and management

A model of forest stand dynamics was developed to simulate the evolution of mono-specific plantations of Norway spruce, Douglas-fir (Pseudotsuga Menziesii (Mirb.) Franco.) and larch (Larix decidua Mill. and L. kaempferii (Lamb.) Carrière) in Southern Belgium. This model, named GYMNOS, is an empirical distance-independent tree model of stand dynamics with an annual time step. It is implemented as a JAVA module in the open-source Capsis4 platform [START_REF] Dufour-Kowalski | Capsis: an open software framework and community for forest growth modelling[END_REF]. It includes equations and algorithms used to simulate stand initialization, top-height growth [START_REF] Perin | Modelling the top-height growth and site index of Norway spruce in Southern Belgium[END_REF], individual tree growth [START_REF] Perin | Distance-independent tree basal area growth models for Norway spruce, Douglas-fir and Japanese larch in Southern Belgium[END_REF], individual tree mortality and thinning (Table 1).

The model had previously been thoroughly tested and used to define reference sylvicultural scenarios such as local yield tables [START_REF] Perin | De nouvelles normes sylvicoles pour les futaies pures équiennes d'épicéa et de douglas en appui à la gestion de la forêt publique en Wallonie[END_REF]. The model is open-source and can be downloaded free (hdl.handle.net/2268/260468), for example to run simulations with different parameters (e.g., site index, plantation density or thinning date and intensity) from those we chose in this study.

Top-height growth

The growth of stand top-height, defined as the mean height of the 100 largest trees per hectare, is an important variable in our stand dynamics model. The evolution of stand top-height is predicted using the model of Perin et al. (2013) (Eq. 1). The top-height in a stand (Hdom y ) that is y years old is predicted knowing the top-height (Hdom 50 ) at 50 years. The top-height reached at 50 years is an indicator of site productivity called the site index (SI), and is used to define yield classes [START_REF] Perin | Modelling the top-height growth and site index of Norway spruce in Southern Belgium[END_REF]. For instance, in sites of average fertility, the top-height of a spruce stand reaches approximately 27 m at 50 years of age. In this model, tree age is defined as the biological age of the trees (plantation age + planted seedling age assumed to be 4 years).

Stand initialization

At stand initialization, the girth at breast height of tree (gbh i ) is drawn from a log-normal distribution (Eq. 2) whose parameters depend on plantation density (N) and top-height (Hdom). To fit this model, we collected data about tree diameter in 332 additionnal plots of 50-1850 m 2 in even-aged stands of Norway spruce, Douglas-fir or larch ranging in age from 9 to 93 years.

The height of a tree (h i ) was estimated using a model depending on the mean girth of the 100 largest trees (Gdom), tree girth (gbh i ) and stand top-height (Hdom, Eq. 3). This model was fitted with another dataset containing records of tree diameter and height of 13,236 trees in 720 circular plots of 1018 m 2 located in even-aged stands of Norway spruce, Douglas-fir or larch ranging also in age from 9 to 93 years.

Growth

Tree growth is predicted using the model of [START_REF] Perin | Distance-independent tree basal area growth models for Norway spruce, Douglas-fir and Japanese larch in Southern Belgium[END_REF]. Tree annual girth increment (∆gbh i ) depends on initial tree girth (gbh i,y-1 ), initial stand top-height (Hdom y-1 ), estimated top-height reached at the end of the year (Hdom y ) and initial stand basal area (BA y-1 ). This model was calibrated for stands taller than 10 m in top-height. Under 10 m top-height, the annual girth increment (between year y and y -1) is derived from Eq. 2 computing the difference between µ values (Eq. 5). Once tree girth is updated, tree height increment is computed with Eq. 6 derived from Eq. 3.

Timber volume

The volume of trunk, from stump to the height of the section 7 cm in diameter, is predicted using local species-specific allometric relationships [START_REF] Dagnelie | Cubage des arbres et des peuplements forestiers : Tables et équations[END_REF] that depend on tree girth (gbh i ) and stand top-height (Hdom, Eq. 7).

i ∼ Log-N (µ, σ 2 ); µ = 0.5 • log(cg 2 ) -σ 2 ; σ = 0.0245 • N 0.317 (2) cg = 40000 • π • 13.3 • rdi N 1/1.65 ; rdi = 1 -exp(-2.61 • 10 -6 • N • Hdom 2 ) H-D allometry h i = 0.799 • Hdom • arctan tan( 1 0.799 ) • gbh i Gdom tan( 1 0.799 ) • gbh i Gdom (3) Girth growth (h i ≥ 10) ∆gbh i = (gbh 2 i,y-1 + ∆g i • 4 • π) 0.5 -gbh i,y-1 (4
)

∆g i = P 2 • gbh i,y-1 -M • A + (M • A + gbh i,y-1 ) 2 ) -4 • A • gbh i,y-1 0.5 A = 3.98 • Hdom 0.780 y P = 0.216 + 0.801 • (Hdom y -Hdom y-1 ) M = 1 + exp(0.135 • Hdom y -0.185 • BA y-1 )
Girth growth (h i < 10)

∆gbh i = exp(log(gbh i,y-1 ) + µ i,y -µ i,y-1 ) + gbh i,y-1 (5) Height growth ∆h i = 0.799 • (Hdom y -Hdom y-1 ) • artan tan( 1 0.799 ) • gbh i,y Gdom y-1 (6)
Timber volume

V i = 0.0135 -0.00128 • gbh i + 0.0000457 • gbh 2 i -7.70 • 10 -8 .gbh 3 i (7) -0.00114 • Hdom + 2.58 • 10 -6 • gbh 2 i • Hdom Taper function g10 i = 5.36 + 1.07 • gbh i -0.00194 • gbh 2 i + 7.47 • 10 -7 • gbh 3 i (8) -0.416 • h i + 2.86 • 10 -5 • gbh 2 i • h i g h,i = a h + b h • g10 i + c h g10 2 i (9) Tree mortality rdi = N N max ; N max = 40000 • π • 13.286 • cg -1.65 (10) rdi max = rdi y-1 + (1 -min(1, rdi 8.5 y-1 ))(rdi -rdi y-1 ) -0.5 max(0, rdi y-1 -1) (11)
s i = u i • gbh i -min(gbh) max(gbh) -min(gbh) u i ∼ U[0, 1] (12)
Thinning The volume of decayed wood resulting from bark-stripping damage is computed using local allometric equations [START_REF] Dagnelie | Tables de cubage des arbres et des peuplements forestiers[END_REF]. First, the girth at a reference height level (10% of tree height, g10 i ) is computed (Eq. 8).

s i = S • u i + (1 -S) • gbh i -gbh * m u i ∼ U[0, 1] (13) gbh * = T • (max(gbh) -min(gbh)) + min(gbh) m = max(gbh * -min(gbh), max(gbh) -gbh * ) + 1 cut damaged /N damaged cut healthy /N healthy = 1.5 (14) Bark-stripping rate τ summer,s = f summer (age s ; µ = 3.08, θ = 0.194) • BSR • (36 -8) 0.854 • 0.2 (15) τ winter,s = f winter (age s ; µ = 2.68, θ = 0.440) • BSR • (36 -8) 0.950 • 0.8 (16) Prob. of damage P(being damaged, i) = 0.0144 • dbh i • exp(-0.125 • dbh i ) + ε i ( 17
)
ε i ∼ N (0, 5.27 • 10 -3 ) Decay spread dl = exp(0.769 + 0.336 • log(Ks) + 0.150 • log(w) + 0.605 • log(l) (18) +0.336 • log(p) + 0.545 • log(rw))
Then, the girth at different height levels (g h ) is computed with equation 9. The chosen height levels were 10, 50, 130, 250, 350 cm in height and 30% of tree height. In this equation, the parameters a h , b h and c h depends on the chosen height level and they were calibrated by [START_REF] Dagnelie | Tables de cubage des arbres et des peuplements forestiers[END_REF].

Using these estimates, the girth at any height can be estimated by linear interpolation. The volume of the base of a trunk is estimated by summing the volume of successive trunk sections and assuming that their shape is a truncated cone [START_REF] Rondeux | La mesure des arbres et des peuplements forestiers[END_REF].

Mortality

Mortality is modeled using a self-thinning curve defining the maximum tree density (N max ) observed in a stand of a given quadratic mean girth (cg) and yields a relative density index (rdi, Eq. 10).

The conventional approach consists in simulating mortality once rdi > 1. However, this approach led to simulations with unrealistically large rdi values and abrupt changes in stand density. A smooth function, depending on the initial stand rdi (rdi y-1 ), was therefore developed to simulate mortality in a more gradual way (Eq. 11). Using this approach, tree mortality is simulated so long as stand rdi is greater than rdi max .

To determine which trees die, a score s i is computed for every tree (Eq. 12). Trees are sorted according to this score and the algorithm sequentially removes the trees with the lowest scores so long as rdi > rdi max . The score includes a stochastic term u i , whose values are drawn from a uniform distribution bounded between 0 and 1.

Thinning algorithm

Thinnings are defined by values of thinning intensity, thinning types and a stochastic coefficient. The thinning intensity corresponds to a target value of stand density, standing volume or relative density index (rdi). The thinning type (T ) ranges between 0 and 1 and translates a forester's preference for harvesting small or large trees. Thinning whose types is 0, 1, and 0.5 harvest respectively the smallest trees, the largest trees and trees of intermediate girth (gbh * ). The stochastic coefficient (S) also ranges between 0 and 1 and indicates the amount of added random noise. If T = 0 and S = 0, then only the smallest trees are harvested whereas if S = 1 then the trees are harvested independently of tree size. The algorithm computes a score s i for each tree, and sequentially harvests the tree with the lowest s i (Eq. 13) until thinning intensity is reached.

This algorithm was further improved to account for the presence of trees with bark-stripping damage because foresters may preferentially thin damaged trees over healthy ones. Using additionnal data measured in 59 thinned stands (10,882 trees) where the proportion of bark-stripping damage ranged between 8% and 99%, we found that the forester's preference for harvesting damaged trees was best modeled by computing the mean odds ratio between damaged and healthy trees (Eq. 14). To ensure that this odds ratio is maintained in the simulation, the thinning algorithm is applied first on damaged trees and secondly on healthy trees with adapted thinning intensities.

Crop trees can be selected and protected against barkstripping damage at the time of the first thinning. These trees are, if possible, not harvested before the final cut. The crop trees were selected preferentially among healthy large trees using the same algorithm as that used to simulate thinnings with S = 0.3 and T = 1.

Models of bark-stripping damage

The abundance of bark-stripping damage is modulated by one variable: the average bark-stripping rate observed at the landscape scale (BSR). This rate is defined as the average proportion of healthy or damaged trees presenting new (maximum 1 year-old) damage, in stands 8-36 years old. It is the main indicator that can be computed from bark-stripping inventories (Section 2.2). Assuming the forest landscape to be composed of pure even-aged spruce stands of ranging age, it is reasonable to assume that ungulates preferentially gather in stands of a certain age that provide adequate food resources and shelter. The bark-stripping rate at the stand level will then vary with stand age. We assumed that bark-stripping rate was nil at age 0, peaked at an intermediate age and was very low in old stands. Like [START_REF] Candaele | Mitigation of bark stripping on spruce: the need for red deer population control[END_REF] using the same bark-stripping inventory data, we assumed that the proportion of the bark-stripping rate of winter and summer damage peaked, respectively, in stands 21 and 12 years old. In addition, we assumed that only 5% of total bark stripping damage occurred in stands younger than 8 years and older than 36 years. This assumption seemed reasonable because most damage has been observed, in different independent studies, in stands 8-36 years old [START_REF] Vospernik | Probability of bark stripping damage by red deer (Cervus elaphus) in Austria[END_REF][START_REF] Jerina | Red deer (Cervus elaphus) bark stripping on spruce with regard to spatial distribution of supplemental feeding places[END_REF][START_REF] Gill | The economic implications of deer damage: a review of current evidence[END_REF][START_REF] Konôpka | Influence of tree species and size on bark browsing by large wild herbivores[END_REF]. With these assumptions (distribution modes, proportion of damage in stands 8-36 years old), we computed the parameters of the two corresponding probability density functions ( f summer and f winter in Eq. 15 and 16). These distributions are used to distribute the damage among stands depending on their age. To simulate the amount of damage corresponding to a chosen BSR, the damage probability functions are multiplied by a first coefficient (C = BSR/ 36 8 f (x) • dx). The result is then multiplied by a second coefficient to ensure that the proportions of summer and winter damage are 20% and 80%, respectively, and to compute the rates of summer and winter new bark-stripping damage (τ summer,s and τ winter,s ).

Next, to assign this damage among the trees of one stand, we used the data of bark-stripping inventories to model the probability of a tree being damaged according to its diameter (dbh, Eq. 17). This probability was used to compute a score for every tree. The damage was assigned to the trees with the highest score. We note that it contains a stochastic term ε i drawn from a normal distribution.

From 2003 to 2005, the positions and sizes of 1820 summer damage instances and 9955 winter damage instances were recorded. Based on these records, the height, width, and length of winter and summer damage wounds were estimated. The damage width and height were not significantly different between summer and winter. However, the summer damage wounds were generally longer [START_REF] Gheysen | A regional inventory and monitoring setup to evaluate bark peeling damage by red deer (Cervus elaphus) in coniferous plantations in Southern Belgium[END_REF]. The distribution of damage height, width, and length was modeled with Weibull and log-normal distributions (Table 3).

After damage, decay can spread around the wound. As decay generally spread in the majority of the damaged stem [START_REF] Metzler | Comparing Norway spruce and silver fir regarding impact of bark wounds[END_REF]Gill, 1992a), we made the simplistic assumption that decay developed in 100% of the damaged stems. Additionally, the rate of decay spread was modeled using the model of [START_REF] Löffler | Zur ausbreitung der wundfäule in der fichte[END_REF]. This model was developed for Norway spruce in Germany. It predicts the length of the decay column (dl) as a function of the tree Kraft status (Ks), damage width (w in cm), damage length (l in cm), time elapsed since the damage (p in years) and mean ring width (rw in mm). The tree Kraft status is a discrete variable used to measure tree social status [START_REF] Kraft | Beiträge zur lehre von den durchforstungen, schlagstellungen und lichtungshieben[END_REF]. It ranges between 1 (dominant tree) and 5 (suppressed tree).

Here also, very different models have been proposed ( Čermák et al., 2004a,b;[START_REF] Čermák | Stem decay by Stereum sanguinolentum after red deer damage in the Českomoravská vrchovina Highlands[END_REF]. We verified that the length of the simulated decay column was consistent with local field observations at the final cut, i.e., with column decay being generally shorter than 3 m and very rarely reaching 4 m (Gheysen T., pers. comm.).

Simulation scenarios

We simulated the even-aged sylviculture of Norway spruce stands following the local guidelines previously adjusted using the same forest stand dynamics model and through a long iteration process with field experts [START_REF] Perin | De nouvelles normes sylvicoles pour les futaies pures équiennes d'épicéa et de douglas en appui à la gestion de la forêt publique en Wallonie[END_REF].

The simulated plantation density was 2000 trees/ha, of which 90% were expected to survive until the first thinning. The first thinning was simulated 17, 19, 22, 25 or 29 years after the plantation in sites whose site index (SI) is 33, 30, 27, 24 or 21 m, respectively. The type of all simulated thinnings (T in Eq. 13) is 0.15. The first thinning reduced stand density to 1280 trees per hectare with a stochastic coefficient fixed at S = 0.692. Before the first thinning, the trees are pruned to a height of 2 m as it is usually done in the study area to ease tree marking. The second, third and subsequent thinnings reduced stand rdi to 0.49, 0.53 and 0.55 with S = 0.3. The thinnings were simulated every 6 years which is the usual cutting cycle length for spruce stands in the study area.

Four protection treatments were tested : without any protection ("the unprotected treatment"), protecting all trees individually at the first thinning ("bark-scraping all"), selecting and protecting 400 crop trees/ha at the first thinning ("bark-scraping 400"), and protecting all trees with a fence (erected at plantation, "the fence treatment"). Different meth- The last solution, bark-scraping, was the cheapest and the most frequently implemented in the study area at the time of writing this manuscript, and so this solution was selected. It consists in scraping the bark generally with a Gerstner plane. The scraping must be light so as not to damage the cambium while stimulating a resin flow (Fig. 3). This solution might, however, not be as effective as more costly ones [START_REF] Ueckermann | Modifizierung der mechanisch-biologischen Maßnahmen zum Schälschutz der jüngsten Fichtenaltersstufe[END_REF]. Simulations were carried out for stands using five values of site index (21, 24, 27, 30, 33 m), 11 values of bark-stripping rates (BSR) ranging between 0 and 10% and four protection treatment levels. As the model included stochastic terms, for each combination, 5 repetitions were simulated leading to 1,100 simulations.

Financial assessment

The financial losses caused by bark-stripping damage were assessed by quantifying the costs (C) and revenues (R) of timber production during a full rotation. The revenues and costs of hunting activities were not considered.

The revenues were computed by multiplying the merchantable volume of trunks (up to a diameter of 7 cm) by unit prices. The volume of the decayed trunk sections was computed using local allometric relationships (Section 2.3.4) and the height reached by the column decay (Section 2.4). The timber unit prices were computed using data of the public sales of timber in 2021 in Southern Belgium. Of all sales, we selected 499 public sales (426,986 m 3 in total) whose volume comprised at least 85% of Norway spruce. We excluded sales originating from sanitary cuts. The unit price was then estimated by girth class (Table 4) [START_REF] Sanchez | Analyse des prix de ventes des bois en forêts publiques[END_REF]. The unit price for the decayed trunk sections could not be estimated with public sale records and it was set at 5e/m 3 , the estimated market price (FNEF, 2022).

The costs considered were those usually incurred following plantation of Norway spruce in Southern Belgium. The costs of the plantation and subsequent operations are given in Table 5. The cost of fencing a 1 ha plantation was set at 6000 e. The cost of bark-scraping a tree to protect it against bark-stripping was set at 1.34 e/tree. [50 -60[ 14.9 [60 -70[ 25.4 [70 -80[ 35.0 [80 -90[ 43.6 [90 -100[ 51.2 [100 -110[ 57.8 [110 -120[ 63.4 [120 -130[ 68.0 [130 -140[ 71.6 [140 -150[ 74.2 [150 -160[ 75.8 [160 -170[ 76.4 [170 -180[ 76.0 [180 -190[ 74.6 [190 -200[ 72.2 > 200 68.8 For each simulation, the cash flows were summarized by computing the net present value (NPV ∞ , Eq. 19) which, for a given discount rate (r), is the difference between the discounted revenues and costs over an infinite time horizon. If the NPV ∞ is greater than zero, then the discounted revenues are greater than the discounted costs and that the investment generates profits. Different rotation lengths were tested and the selected optimum rotation length was the one maximizing NPV ∞ . The tested rotation lengths corresponded to the different thinning dates (Section 2.5), i.e., every 6 years from the first thinning to 100 years.

NPV = ∑ n i=0 R (i) -C (i) (1 + r) i NPV ∞ =NPV. (1 + r) n (1 + r) n -1 (19)
The opportunity cost of bark-stripping damage, i.e. the discounted loss due to bark-stripping damage, was computed with the unprotected simulations as NPV ∞,BSR=0 -NPV ∞ , where NPV ∞,BSR=0 was the average of NPV ∞ for the simulations with the same site index and BSR = 0%. To discuss the magnitude of the opportunity costs, we computed the corresponding annuities by multiplying the opportunity costs by the discount rate. Such annuities can be interpreted as the amount that need to be paid annually to offset the cost of the bark-stripping damage.

These computations critically depend on the discount rate, which is the key variable to estimate the discounted (or present) value of future cash flows. For such projects, there is still scant agreement on what the appropriate rate might be [START_REF] Gollier | The long-run discount rate controversy[END_REF]. We therefore avoided a single value and we repeated the computations with different rates (1%, 2%, 3% and 4%).

In addition, we verified that computing the opportunity cost of bark-stripping damage over a single rotation (i.e. computing the difference between NPV rather than between NPV ∞ ) led to the same conclusions whatever the discount rate, and barely affected our estimates of the opportunity cost (Suppl. Table S1).

Statistical analyses

Simulation results were mostly compared visually because the differences were generally well marked. Even so, some models were fitted to synthesize the results.

The relationships between the opportunity cost of barkstripping damage in response to stand site index, discount rate, and bark-stripping rate was modeled with a linear model whose intercept was zero and whose slope depended on the discount rate and site index (Eq. 20).

Opportunity cost

= (β 1 + β 2 • r + β 3 • SI) • BSR (20)
Similarly, the NPV ∞ was modeled with a linear model whose intercept and slope depended on the site index, discount rate, bark-stripping rate and protection treatment denoted T in Eq. 21.

NPV ∞ = (α 1,T + α 2 • SI + α 3 • r) + (β 1,T + β 2 • r) • BSR (21)

Results

Stand characteristics over time

Over the 100 simulated years, the main characteristics of the simulated stands depended primarily on site index (SI) and little on treatment or bark-stripping rate (Fig. 4). For a given site index, the timber volume productivity and, the intensity and type of the thinnings were thus very similar across the treatments and bark-stripping rates. Trees grew faster in the most productive site with a high site index. At stand initialization (i.e. at 8 year old), the mean tree diameter ranged between 3.2 cm (SI = 21 m) and 6.0 cm (SI = 33 m). Mean productivity ranged between 11.7 m 3 /ha/year (SI = 21 m) and 17.9 m 3 /ha/year (SI = 33 m).

The number of trees with bark-stripping damage increased with bark-stripping rate (BSR) depending also on the treatment and plantation age (Fig. 5). The number of bark-stripped trees generally peaked in 25-year-old stands and then slowly declined as the damaged trees were thinned. For example, in unprotected stands with SI = 27 m, the proportion of damaged trees was about 31% and 81% with BSR = 2% or 8%, respectively.

At clear-cut, the number of trees containing decayed wood increased with the bark-stripping rate and depended on the treatment (Fig. 6). When the bark-stripping rate was low (e.g. BSR = 2%), most of the harvested trees were healthy (Fig. 5) but the proportion of harvested trees with damage increased sharply with BSR. In unprotected plantations, at least half of the harvested trees were damaged if BSR ≥ 5% and more than 80% of trees were damaged if BSR ≥ 10%.

At clear-cut, the height of the decay column in the damaged trees ranged between 0.6 m and 5.6 m with an average of 2.1 m (Fig. 6C). This variability resulted mostly from the stochastic terms of the models used to predict the damage dimensions (Section 2.4).

Losses in timber volume and value

For all treatments except fencing, the total volume of harvested undamaged trunk sections (during thinnings and clearcut) decreased with bark-stripping rate (Fig. 7). The reduction was obviously stronger in unprotected plantations than in plantations with individual protections. The volume of undamaged trunk section was very similar in the two scenarios with individual protections (bark-scraping all and bark-scraping 400).

Individual protections were installed at the first thinning. As the first thinnig occured later in the stands of low productivity (low site index), individual protection were installed later in the low productivity stand than in the high productivity stand. In the low productivity stands, the proportion of damaged trees could then already be high when individual protections were applied (bark-scraping all and bark-scraping 400). The effect of these protections thus became limited (Suppl. Fig. S2).

Net present value

Optimum rotation time

The optimum rotation value varied substantially with the site index (SI), discount rate (r), and bark-stripping rate (BSR, Suppl. Fig. S1). Optimum rotation was shorter for the most productive sites (high site index) and/or when using a high discount rate. For example, when r = 2% and BSR = 0%, the mean optimum rotation was 65.6 and 96.3 years in stands with SI = 33 m and SI = 21 m, respectively.

The optimum rotation also depended, though to a lesser extent, on the interaction between the treatment and the barkstripping rate. For example, when r = 2%, BSR = 0% and SI = 27 m, the optimum rotation was 76 years in unprotected stands and 82 years in fenced stands. In the same conditions but with BSR = 10%, the optimum rotation did not differ between treatments (Suppl. Fig. S1, 82 years for all treatment). Generally, the optimum rotation tended to be longer when the bark-stripping rate was high. When the bark-stripping rate was low, the optimum rotation of fenced stands was longer than that in the other scenarios. The opposite was nevertheless observed when bark-stripping rate was high and discount rate low.

Moreover, combining a low site index (SI = 21) m and a low discount rate (r = 1%) gave an optimum rotation that was longer than the tested ones (100 years, Suppl. Fig. S1). An optimum rotation of 100 years was nevertheless considered for these cases.

Opportunity cost of bark-stripping damage

The opportunity cost of bark-stripping damage ranged between 0 and 19, 445 e/ha. For example, with r = 2%, in unprotected stands of average fertility (SI = 27 m), the net present value was about 16,492 e/ha if no damage occurred (BSR=0%). It fell to 11,048 e/ha with the highest simulated level of damage (BSR=10%) as the opportunity cost was 5,443 e/ha (i.e., 37% of NPV ∞ without damage).

The opportunity cost of bark-stripping damage depended on the bark-stripping rate, discount rate, and site index (Fig. 8). More precisely, the opportunity cost increased linearly with bark-stripping rate, and the slope of this relationship increased with the site index and decreased with the discount rate (Table 6). The opportunity cost expressed in absolute value increased with site index while its the relative value decreased with site index (Fig. 8, 9 and Suppl. Fig. S5).

Net present value across treatments

The net present value increased with site index, decreased with discount rate, and varied across treatments. When no damage occurred (BSR = 0%), the net present value of the fenced treatment was lowest (α 1, fα 1,u = -7.651 • 10 3 , p < 0.001, Table 7) and the net present values of the scenarios with individual protections were only slightly, and not significantly BSR ∩ SI 3.10 • 10 3 1.88 • 10 2 < 0.001 lower than that of the unprotected treatment (e.g. α 1,s4α 1,u = -1.373, p = 0.860, Fig. 10).

The net present value decreased linearly with the barkstripping rate for most treatments (β 1,u = -1.308 • 10 5 , p < 0.001) except for the fenced treatment, which was obviously unaffected by the bark-stripping rate (Table 7, Fig. 10 and Suppl. Fig. S3). For all scenarios except fencing, the slopes of this relationship varied little across treatments when the discount rate was low (r = 0.01%), whereas the slopes of the unprotected scenario were markedly steeper than that of the other scenarios when the discount rate was high (r ≥ 0.03%, Fig. 10,Table 7).

The net present value of the fence scenario was higher than that of the other scenarios only when r = 1, SI < 27 and BSR > 5%. The scenarios with individual protections had a higher net present value than the unprotected scenario when SI ≥ 27m. The difference between these scenarios increased with r (Fig. 10, Suppl. Fig. S3).

Discussion

Model assumptions

Modeling such a complex system inevitably entails making simplistic assumptions, but, it allows hypotheses to be tested while controlling many factors. In particular, our simulations did not allow for any changes in the growth or mortality of Norway spruce trees due to climate change or any other disturbances (Cukor et al., 2019a;[START_REF] Šn ¸epsts | Cervid bark-stripping is an explicit amplifier of storm legacy effects in Norway spruce (Picea abies (L.) Karst.) stands[END_REF]. The models of tree growth and mortality were calibrated with empirical data collected between 1971 and 2010 in healthy pure softwood plantations. However, such plantations have been increasingly affected by climate change, warmer temperatures, bark beetle outbreaks [START_REF] Jönsson | Spatio-temporal impact of climate change on the activity and voltinism of the spruce bark beetle, Ips typographus[END_REF] and droughts [START_REF] Lévesque | Drought response of five conifer species under contrasting water availability suggests high vulnerability of Norway spruce and European larch[END_REF]. Nevertheless, the extent of such damage could be reduced if spruce is progressively replaced by other tree species and if winters become milder [START_REF] Ligot | Modeling recent bark stripping by red deer (Cervus elaphus) in south Belgium coniferous stands[END_REF][START_REF] Candaele | Mitigation of bark stripping on spruce: the need for red deer population control[END_REF][START_REF] Konôpka | Influence of tree species and size on bark browsing by large wild herbivores[END_REF]. The combined effects of climate changes and changes in forest structure and composition need to be addressed and clarified in future studies.

Bark-stripping is one only type of damage caused by ungulate populations. Other damage can be expected when ungulates are abundant. Browsing of juvenile trees and antler rubbing on tree trunks can further reduce timber yield [START_REF] Ward | Modelling the cost of roe deer browsing damage to forestry[END_REF]. In this work, the cost of this additional damage was considered constant, although we could expect it to Additionally, we considered that bark-stripping damage and the associated decay would not affect tree growth or survival. While the vast majority of damaged trees survive, very rare exceptions can occur when the bark is removed over the entire stem circumference on small trees and in highly impacted areas [START_REF] Gill | The economic implications of deer damage: a review of current evidence[END_REF]. Furthermore, contrasting effects of bark-stripping damage on tree growth have been reported in the literature. These effects depend on wound size and position (Gill, 1992a;[START_REF] Mäkinen | Increment and decay in Norway spruce and Scots pine after artificial logging damage[END_REF]. It seems likely that a growth reduction can be observed in a damaged tree, but such growth reduction may be limited (except when the wound is very large) to a period of about 10 years after the wounding [START_REF] Gill | The economic implications of deer damage: a review of current evidence[END_REF]Cukor et al., 2019b) and can ultimately be negligible [START_REF] Gill | The economic implications of deer damage: a review of current evidence[END_REF]. However, damaged stems can have reduced mechanical resistance ( Čermák et al., 2004a; Šn ¸epsts et al., 2022). Bark damage impairs tree hydraulics and may induce physiological drought and increase tree susceptibility to drought (Cukor et al., 2019b,a;[START_REF] Vacek | Bark stripping, the crucial factor affecting stem rot development and timber production of Norway spruce forests in Central Europe[END_REF], bark-beetle attacks, and wind damage [START_REF] Snepsts | Stem damage modifies the impact of wind on Norway spruces[END_REF][START_REF] Šn ¸epsts | Cervid bark-stripping is an explicit amplifier of storm legacy effects in Norway spruce (Picea abies (L.) Karst.) stands[END_REF][START_REF] Krisans | Effect of barkstripping on mechanical stability of Norway spruce[END_REF].

We assumed that very few bark-stripping damage occurred in plantations younger than 8 years old (the simulations started in stands 8 years old) in line with several studies [START_REF] Vospernik | Probability of bark stripping damage by red deer (Cervus elaphus) in Austria[END_REF][START_REF] Jerina | Red deer (Cervus elaphus) bark stripping on spruce with regard to spatial distribution of supplemental feeding places[END_REF][START_REF] Gill | The economic implications of deer damage: a review of current evidence[END_REF][START_REF] Konôpka | Influence of tree species and size on bark browsing by large wild herbivores[END_REF]. At 8 years old, the mean tree diameter ranged between 3.2 cm (SI = 21m) and 6.0 cm (SI = 33m), and so it seems reasonable to assume that bark-stripping damage seldom occurs on such small trees because their stem will not be stiff or accessible enough [START_REF] Gill | The economic implications of deer damage: a review of current evidence[END_REF]. Furthermore, in the study area, all trees in spruce plantations are generally pruned to a height of 2 m before the first thinning, and this operation, which is carried out earlier in the most productive stands, may hypothetically affect the behavior of ungulates. Nevertheless, we simulated damage occurrence without taking pruning into account because its effect is still poorly documented.

We also assumed that rot developed in all bark-stripping wounds, although it could develop in a smaller proportion of the wounded stems. This may have resulted in an overestimation of the losses induced by bark-stripping damage. The proportion of the wounds in which rot develops ranges widely across studies and environmental conditions. In a literature review, Gill (1992a) estimated this range at between 73.0% and 99.7%. Similar values were recently found in Germany (Metzler et al., 2012, 93%) but lower values were found in Latvia (Burn ¸eviča et al., 2016, 13-50%) and in the Czech Republic ( Čermák and Strejček, 2007, 68%). The rate of infection probably depends on wound size and environmental conditions [START_REF] Vasaitis | Wound occlusion and decay in Picea Abies stems[END_REF][START_REF] Vasiliauskas | Damage to trees due to forestry operations and its pathological significance in temperate forests: a literature review[END_REF]Cukor et al., 2019a). Although, infection was assumed to occur in every wound, the decay spread was modeled taking into account wound size and environmental conditions [START_REF] Löffler | Zur ausbreitung der wundfäule in der fichte[END_REF]. This model was deemed valid in the study area, but according to other studies carried out in different environmental conditions, the decay may spread slightly faster ( Čermák et al., 2004a;[START_REF] Čermák | Stem decay by Stereum sanguinolentum after red deer damage in the Českomoravská vrchovina Highlands[END_REF].

How badly does bark-stripping harm timber production and profitability?

The cost of bark-stripping damage can be substantial. It ranged between 0 and over 100% depending mostly on the rate of bark-stripping damage, discount rate, and, to a lesser extent, site index (Fig. 8 and Suppl. Fig. S3). With high barkstripping rate (BSR = 10%) and without protection, around 85% of the trees harvested at the final cut contained decayed wood and the volume of this decayed wood accounted for 15% of stand volume. According to [START_REF] Vacek | Bark stripping, the crucial factor affecting stem rot development and timber production of Norway spruce forests in Central Europe[END_REF] and [START_REF] Heyninck | Pertes en volume et en valeur d'épicéas adultes touchés par des dégâts de cervidés[END_REF], for one damaged tree, the decayed wood volume could instead account for 30-40 % of stem volume. These studies using other methodological approaches thus suggest that our predictions of the decayed wood volume at the stand level could be underestimated. However, the study of [START_REF] Vacek | Bark stripping, the crucial factor affecting stem rot development and timber production of Norway spruce forests in Central Europe[END_REF] was conducted on younger spruce stands in the Czech Republic and the study of [START_REF] Heyninck | Pertes en volume et en valeur d'épicéas adultes touchés par des dégâts de cervidés[END_REF] analyzed a very small sample of Norway spruce stems (n = 17) in a single site in Belgium.

Our results can be used to estimate the average cost of bark-stripping in the study area, making assumptions about the average bark-stripping rate, site index, and discount rate. Within the study area,the mean bark-stripping rate is about 4% [START_REF] Gheysen | A regional inventory and monitoring setup to evaluate bark peeling damage by red deer (Cervus elaphus) in coniferous plantations in Southern Belgium[END_REF]. This rate also corresponds to the damage rate tolerated by the forest administration. In Southern Belgium, Norway spruce yields 14.7 m 3 /ha on average, which corresponds to a site index of 27 m. If we choose a discount rate of r = 2%, which seems reasonable for a project maturing over 50 years [START_REF] Gollier | The long-run discount rate controversy[END_REF], the average opportunity cost of bark-stripping damage will be 2647 e/ha which corresponds to an annuity of 53 e/ha/year. These figure are high, corresponding to a loss of net revenue of 19%. This cost is also of the same order of magnitude as the plantation costs or hunting leases.

Cost also depends on site fertility. Like [START_REF] Gill | The economic implications of deer damage: a review of current evidence[END_REF] for Sitka spruce plantations, we found that the relative cost of bark-stripping damage was greater in the least productive stands (Suppl. Fig. S5). In the most productive stands, the relative cost was lower, but the absolute cost was nevertheless higher (Fig. 8). For example, with BSR=4% and r =2%, the opportunity cost was about 2391 e/ha (35 % of NPV ∞ ) in the least productive stands and 3880 e/ha (17 %) in the most productive stands (with SI=33 and 21 m).

With our virtual experiment, we made a comprehensive assessment of the cost of bark-stripping damage over a full cycle of timber production. In particular, using Eq. 20 and the fitted parameters (Table 6), we can predict the opportunity cost of bark-stripping damage for any given values of the discount rate, site index and bark-stripping rate. Such information can be particularly helpful in making the management of Norway spruce plantations and deer populations more effective [START_REF] Ward | Modelling the cost of roe deer browsing damage to forestry[END_REF].

How should forest management be adapted?

Rotation should be kept unchanged or only slightly lengthened

Forest managers may be tempted to shorten the rotation in stands where bark-stripping damage is observed because they are eager to start a new rotation. However, we found that the optimum rotation should either remain the same or be only slightly lengthened (e.g., 6 years more, Suppl. Fig. S1) irrespective of the discount rate used. In old stands, new bark-stripping damage becomes less likely [START_REF] Vospernik | Probability of bark stripping damage by red deer (Cervus elaphus) in Austria[END_REF][START_REF] Jerina | Red deer (Cervus elaphus) bark stripping on spruce with regard to spatial distribution of supplemental feeding places[END_REF][START_REF] Gill | The economic implications of deer damage: a review of current evidence[END_REF] and the spread of the decay is very low ( Čermák et al., 2004a;[START_REF] Čermák | Stem decay by Stereum sanguinolentum after red deer damage in the Českomoravská vrchovina Highlands[END_REF][START_REF] Löffler | Zur ausbreitung der wundfäule in der fichte[END_REF]. Considerations about the rotation length should therefore first take into account the effects of the discount rate, site index and certainly also, though not examined here, the risk of future disturbance to stands [START_REF] Zimová | Reducing rotation age to address increasing disturbances in Central Europe: potential and limitations[END_REF] and other forest ecosystem services [START_REF] Sing | A review of the effects of forest management intensity on ecosystem services for Northern European temperate forests with a focus on the UK[END_REF].

Cheap protections can still be helpful

Fencing to protect trees against bark-stripping damage is unlikely to be cost-effective. In most cases, its cost was appreciably greater than the opportunity cost of the damage (Fig. 8, 10 and Suppl. Fig. S3). Moreover, we probably underestimated its cost (6000 e/ha) because the fence maintenance costs were not counted. Fences in forests usually last about 15 years, much less than the time trees need protection against bark-stripping damage (about 40 years) [START_REF] Gill | The economic implications of deer damage: a review of current evidence[END_REF]. Fences were cost-effective only when the chosen discount rate was low (r ≤1%), the bark-stripping rate high and/or likley when saplings need also to be protected from browsing [START_REF] Jensen | Shrubs protect oak seedlings against ungulate browsing in temperate broadleaved forests of conservation interest: A field experiment[END_REF]. Cheaper protections might be cost-effective particularly in the most productive stands. In this study, we assumed that such individual protections were placed at the time of the first thinning (at 17-29 years of age) [START_REF] Perin | De nouvelles normes sylvicoles pour les futaies pures équiennes d'épicéa et de douglas en appui à la gestion de la forêt publique en Wallonie[END_REF]. At that time, the trees are big enough to be pruned and their bark scraped to protect them against bark-stripping. Nevertheless, a substantial proportion of the trees may already be damaged at that age (Gill, 1992b;[START_REF] Candaele | Mitigation of bark stripping on spruce: the need for red deer population control[END_REF][START_REF] Vospernik | Probability of bark stripping damage by red deer (Cervus elaphus) in Austria[END_REF]. In the simulations with high bark-stripping rates, individual protections could be fitted on already damaged trees, so having no effect or a very small effect on future timber value. This occurred mostly in the less productive stands and when the bark-stripping rate was very high (BSR > 7%, Suppl. Fig. S4). In the most productive stand, the first thinning occurs earlier [START_REF] Perin | De nouvelles normes sylvicoles pour les futaies pures équiennes d'épicéa et de douglas en appui à la gestion de la forêt publique en Wallonie[END_REF], and enough healthy crop trees can probably be found making individual protection a valid solution. In less fertile stands, individual protection should be fitted several years before the first thinnings.

Conclusion

Levels of wild ungulate populations have usually been adjusted to the damage levels, with limited regard to the actual cost of such damage. The model we propose in this study can be used to assess the cost of bark-stripping damage balancing long-term revenues against short-term costs of protection measures and long-term costs of bark-stripping damage. Knowing the true cost of bark-stripping damage is essential to improving the economic efficiency of deer and forest management. Protecting plantations against bark-stripping damage with fences was found unlikely to be worthwhile. By contrast, individual protections placed on crop trees could be helpful, particularly in the most fertile stands. Loss of revenue depended greatly on the factors tested: we estimated of the average damage cost could be about 53e/ha/year, reducing timber yield by 19%. Bark-stripping is one of the several factors affecting the proffitablity of Norway spruce plantations. To properly guide forest management, the other factors (e.g. browsing, changing environmental conditions, pest outbreak) that sometimes interact with bark-stripping should also be considered. The barkstripping damage should also be further investigated for unevenaged and mixed stands.

Figure 1 .

 1 Figure 1. Map of the study area as illustrated by Candaele et al. (2021).

Figure 2 .

 2 Figure 2. Design of the sampling method as illustrated by Candaele et al. (2021). The plots are installed on the nodes of a square grids that fall into Norway spruce plantations of 8-36 years old (sensitive stands, A). At each selected nodes, six trees were measured in three circular plots (B).

Figure 3 .

 3 Figure 3. Within the study area, spruce trees are often protected from bark damage by ungulates by lightly scraping (e.g. with a Gerstner plane) the bark to stimulate a resin flow. The scars remain visible during several years.

Figure 4 .Figure 5 .

 45 Figure 4. Changes in the main characteristics of the simulated stands: density, mean quadratic diameter, basal area, standing volume and relative density index (rdi). For clarity, only the simulations with averaged productivity stands are shown (SI = 27 m).

Figure 7 .Figure 8 .Figure 9 .

 789 Figure 7. Undamaged volume and value of the trees thinned and harvested across treatments and bark-stripping rates. For clarity, only the simulations with averaged productivity stands are shown (SI = 27m).

Table 1 .

 1 Equations of the stand dynamics and bark-stripping damage. The symbols and abbreviations are described in table 2.

	Process

Table 2 .

 2 Symbols and abbreviations.

	Symbols	Description
	age	Stand age
	Hdom y	Stand top height at y years of age (m)
	Hdom 50 = SI Site index, stand top height at y years of
		age (m)
	N	stand density (trees/ha)
	gbh i	Girth of tree i at 1.3 m (cm)
	dbh	Tree diameter at 1.3 m (cm)
	∆gbh i	Annual girth increment of tree i
		(cm/year)
	cg	Quadratic mean of tree girth (cm)
	rdi	Stand density index
	N max	Maximum tree density (trees/ha)
	rdi max	Maximum stand density index
	Gdom	Mean girth of the 100 largest trees/ha
		(cm)
	BA rtithm
	cut damaged	Density of damaged thinned trees
		(trees/ha)
	cut healthy	Density of healthy thinned trees
		(trees/ha)
	N healthy	Stand density of healthy trees, before
		thinning (trees/ha)
	N damaged	Stand density of damaged trees, before
		thinning (trees/ha)
	BSR	Bark-stripping rate
	f summer	Summer damage probability function
	f winter	Winter damage probability function
	τ summer,s	Summer bark-stripping rate in stand s
	τ winter,s	Winter bark-stripping rate in stand s
	ε i	Model residual
	dl	Length of the decay column (cm)
	Ks	Tree Kraft status, social position
	w	Damage width (cm)
	l	Damage length (cm)
	p	Time elapsed since the damage (years)
	rw	Mean ring width (mm)
	µ	Distribution mean
	θ	Distribution standard deviation
	NPV	Net present value over a finite (n years)
		time horizon (e/ha)
	NPV ∞	Net present value over an infinite time
		horizon (e/ha)
	r	Discount rate

y Stand basal area at year y (m 2 /ha) g10 i Tree girth at 10% of tree heigth (cm) g h Tree girth at h m (cm) a h , b h , c h Fitted parameters of the taper equation T Thinning type S Thinning stochastic coefficient s i

Score of tree i used in the thinning algo-

Table 3 .

 3 Distribution parameters of bark-stripping damage height, width and length.

	Variables	Season	Distribution	Parameter ± std. error
	Height (cm)	All	Weibull(k, δ )	4.692 ± 0.034 112.421 ± 0.232
	Width (cm)	All	Log-N (µ, θ 2 ) 1.582 ± 0.006	0.634 ± 0.004
	Length (cm) Summer Log-N (µ, θ 2 ) 2.834 ± 0.022	0.922 ± 0.015
	Length (cm)	Winter	Log-N (µ, θ 2 ) 2.236 ± 0.007	0.704 ± 0.005
	ods have been proposed to protect trees individually against
	bark-stripping damage : lattice sleeve, bark protection net and
	bark scraping. All these protections are usually implemented
	at the time of the first thinning and generally after pruning.

Table 4 .

 4 Timber price by girth class estimated from public timber sales in 2021, in Southern Belgium.

	Girth class (cm) Price (e/m 3 )
	< 40	0.0
	[40 -50[	3.3

Table 5 .

 5 Costs common to all simulation scenarios.

	Cost	Year Price (e/ha)
	Plantation 0	2400
	Weeding	1	640
	Weeding	3	640
	Weeding	5	640
	Pruning	17	1219

Table 6 .

 6 Parameter estimates, standard errors and p-value of the fitted linear model of the opportunity cost of bark-stripping damage (Eq. 20) in response to the bark-stripping rate (BSR), site index (SI) and discount rate (r).

	Param.	Effect	Estimate Std. error	p
	β 1	BSR	1.05 • 10 5	5.43 • 10 3 < 0.001
	β 2	BSR ∩ r -4.75 • 10 6	7.13 • 10 4 < 0.001
	β 3

Table 7 .

 7 Parameter estimates, standard errors, and p-values of the fitted linear model of the net present value (NPV ∞ , Eq. 21) in response to the bark-stripping rate (BSR), site index (SI), discount rate (r), and protection treatment. The levels of the protection treatment are denoted u for the unprotected treatment, s for bark-scraping of all trees, s4 for bark-scraping of 400 crop trees/ha, and f for fencing. The symbol ∩ indicates an interaction between two explanatory variables.

	Parameter	Effect	Estimate Std. error	p
	α 1,u	Treatment	2.468 • 10 4 1.245 • 10 3 < 0.001
	α 1,s4 -α 1,u	Treatment	-1.373 • 10 2 7.762 • 10 2	0.860
	α 1,s -α 1,u	Treatment	-1.216 • 10 3 7.762 • 10 2	0.117
	α 1, f -α 1,u	Treatment	-7.651 • 10 3 7.762 • 10 2 < 0.001
	α 2	SI	1.463 • 10 3 3.457 • 10 1 < 0.001
	α 3	r	-1.764 • 10 6 2.454 • 10 4 < 0.001
	β 1,u	BSR	-1.308 • 10 5 1.392 • 10 4 < 0.001
	β 1,s4 -β 1,u BSR ∩ Treatment	2.041 • 10 4 1.312 • 10 4	0.198
	β 1,s -β 1,u	BSR ∩ Treatment	2.540 • 10 4 1.312 • 10 4	0.053
	β 1, f -β 1,u BSR ∩ Treatment	6.569 • 10 4 1.312 • 10 4 < 0.001
	β 2	r ∩ BSR	2.596 • 10 6 4.149 • 10 5 < 0.001

Availability of data and material

An archive ZIP file containing the data generated by the simulations, an R scripts, and a read-me text file are available in the supplementary materials.

Code availability

The model is open-source and can be downloaded (hdl.handle.net/ 2268/260468).

Unprotected

Bark-scraping 400 Bark-scraping all Fence

Declarations Funding

The study was funded by the public service of Wallonia (SPW) under the funding: "Accord-cadre de recherches et vulgarisation forestières" (ACRVF 2009(ACRVF -2014(ACRVF and 2014(ACRVF -2019)).


P Develter

J Bosse

O Rabaste

P Forster

J.-P Ovarlez

FALSE ALARM REGULATION FOR OFF-GRID TARGET DETECTION WITH THE MATCHED FILTER

Keywords: Radar detection, Off-Grid, GLRT, Matched Filter, False Alarm Regulation

In the state-of-the-art, the Probability of False Alarm (P F A )threshold relationship for the popular Matched Filter (MF) is often derived assuming that unknown non-linear parameters lie on a grid. However, these parameters vary continuously in practice. This is known as the off-grid case. In this article, an asymptotic P F A -threshold relationship for the popular Matched Filter is derived in the off-grid case under complex white Gaussian noise hypothesis using expected Euler characteristics. This asymptotic relationship fits very well with corresponding Monte-Carlo trials in the moderate to low P F A regime.

INTRODUCTION

Detection of signals embedded in noise is a classical problem encountered in many fields [START_REF] Kay | Fundamentals of statistical processing[END_REF]. Signals that depend on unknown parameters are usually dealt with using a Generalized Likelihood Ratio Test (GLRT) where the unknown parameters are replaced with their Maximum Likelihood Estimates (MLE) in the Likelihood Ratio Test. However, analytical MLE solutions are not always available. In this scenario, those parameters are usually supposed to lie on a discrete set, called the grid. Generally, this grid is chosen so as to have statistically independent tests. However, the parameters of interest have no reason to fall precisely on the grid. Thus, the signals of interest will always be mismatched with the signals under test. This off-grid mismatch may induce a loss in performance. This is the case when the GLRT coincides with the well-known Matched Filter, which can present up to a 3dB loss per parameter dimension.

Detection under mismatched signal models has been addressed for mismatch lying in a cone [START_REF] Besson | Detection of a signal in linear subspace with bounded mismatch[END_REF], quadratically constrained mismatch [START_REF] De Maio | Fractional QCQP with applications in ML steering direction estimation for radar detection[END_REF], among other types of mismatch [START_REF] Bandiera | Advanced Radar Detection Schemes Under Mismatched Signal Models[END_REF][START_REF] Liu | Robust detection in mimo radar with steering vector mismatches[END_REF][START_REF] Besson | Adaptive detection of a signal known only to lie on a line in a known subspace, when primary and secondary data are partially homogeneous[END_REF].

Whereas off-grid mismatch has been explored extensively in sparse reconstruction [START_REF] Tang | Compressed sensing off the grid[END_REF][START_REF] Lasserre | Bayesian sparse Fourier representation of off-grid targets with application to experimental radar data[END_REF] and direction of arrival estimation [START_REF] Dai | Root sparse Bayesian learning for off-grid DOA estimation[END_REF], it has received less attention in detection (see, for example [START_REF] Bosse | Adaptive subspace detectors for off-grid mismatched targets[END_REF][START_REF] Develter | On the false alarm probability of the normalized matched filter for off-grid target detection[END_REF]).

The most natural and efficient way to address the off-grid issue is to consider the true GLRT where the unknown parameters lie on a continuous space instead of an unrealistic discrete grid. There are two main issues with this procedure: it can be computationally heavy and its theoretical test statistics depend on the maximum of a continuum of non-independent variables, which may be difficult to evaluate. In [START_REF] Develter | On the false alarm probability of the normalized matched filter for off-grid target detection[END_REF], we derived a P F A -threshold relationship for the off-grid Normalized Matched Filter. In this paper, our contribution focuses on the case of the off-grid radar Matched Filter with respectively one (Doppler or angle) and two (both Doppler and angle) unknown parameters. Based on results on random Gaussian fields [START_REF] Adler | On excursion sets, tube formulas and maxima of random fields[END_REF], we propose a new asymptotic P F A -threshold relationship in the case of complex circular white Gaussian noise.

Section 2 presents the signal model and the off-grid GLRT detector analyzed in this paper. In Section 3, we derive the corresponding theoretical P F A -threshold relationship and its validation with Monte-Carlo simulations in Section 4.

Notations: Matrices are in bold and capital, vectors in bold. For any matrix A or vector, A T is the transpose of A and A H is the Hermitian transpose of A. I is the identity matrix and CN (µ, Γ) is the circular complex Normal distribution of mean µ and covariance matrix Γ. ⊗ is the Kronecker product. a r and a i are the complex and imaginary parts of the complex vector a.

PROBLEM FORMULATION

Our problem of interest is to detect a potential known signal s(ξ) ∈ C N of unknown amplitude, depending on unknown parameter ξ and embedded in noise. Formally, this reduces to solving the following hypothesis testing problem:

H 0 : r = n , H 1 : r = a s(ξ) + n , (1)
where n ∈ C N is the noise vector distributed according to a complex circular Gaussian distribution n ∼ CN (0, Γ) with known covariance matrix Γ and a ∈ C is the complex amplitude of the target. In the rest of the paper, the number of unknown parameters will be fixed to one or two. In the case of a single parameter (e.g. Doppler shift ξ = θ and N pulses), the signal s(θ) can generally be expressed as:

s(θ) = 1 √ N 1, e 2iπθ , . . . , e 2iπ(N -1)θ T (2)
In the case of two unknown parameters (ξ = (θ, µ) T ) as, for example, for Space-Time Adaptive Processing (STAP), the signal s(θ, µ) ∈ C N P can be represented by the joint Doppler θ (N pulses) and the direction of arrival µ (P sensors) parameters as:

s(ξ) = s(θ) ⊗ s(µ), (3)
where s(θ) ∈ C N and s(µ) ∈ C P are defined as in [START_REF] Besson | Detection of a signal in linear subspace with bounded mismatch[END_REF]. When ξ is known and the complex amplitude a is unknown, the GLRT of this problem denoted as oracle GLRT reduces to the well-known MF test [START_REF] Scharf | Signal detection in Gaussian noise of unknown level: an invariance application[END_REF][START_REF] Scharf | Matched subspace detectors[END_REF]:

s(ξ) H Γ -1 r 2 s(ξ) H Γ -1 s(ξ) H1 ≷ H0 w 2 . (4)
In that case, the corresponding general P F A -threshold relationship is given by the following equation:

P F A = exp (-w 2 ). (5)
Since the true target parameters are unknown, tests are generally evaluated for parameter ξ values lying on a discrete grid. However, in reality, target parameters will rarely fall precisely on the grid. For typical Fourier resolution grid spacing, loss of the matched filter can reach up to 3dB per dimension (and thus up to 6dB in STAP), degrading detection performance, as illustrated in Figure 1.

The true GLRT solution, named off-grid GLRT or off-grid MF, is to evaluate (4) for all ξ varying continuously in the search domain D, referred as a cell (interval [k/N, (k+1)/N ] in 1D) in the sequel:

max ξ∈D s(ξ) H Γ -1 r 2 s(ξ) H Γ -1 s(ξ) H1 ≷ H0 w 2 . ( 6
)
The resulting empirical detector performance is shown in Figure 1. It is close to the oracle GLRT (which perfectly knows the position of targets) while the Matched Filter exhibits a detection loss due to the off-grid targets. To guarantee a given value of P F A , the off-grid GLRT requires setting the threshold w 2 in (6). For that purpose, we need to know the statistics of the maximum of a non-independent random variable continuum. To our knowledge, there is no known solution to this problem in our context. The aim of the next section is to fill this gap in the white noise scenario (Γ = I).

AN ASYMPTOTIC P F A -THRESHOLD RELATIONSHIP FOR THE OFF-GRID MATCHED FILTER

In this section, we derive a new asymptotic P F A -threshold relationship for the off-grid MF under white noise, using expected Euler characteristics: studies [START_REF] Adler | On excursion sets, tube formulas and maxima of random fields[END_REF][START_REF] Siegmund | Testing for a signal with unknown location and scale in a stationary Gaussian random field[END_REF] provide a general theoretical framework for solving this kind of problem.

Estimating the P F A for a single parameter steering vector

In the case of range or Doppler detection, the MF response can be seen as the following 1D random field Y (θ):

Y (θ) = s(θ) H r 2 . (7)
In order to use the results of [START_REF] Adler | On excursion sets, tube formulas and maxima of random fields[END_REF], we need to introduce the following Gaussian random field:

X(α, θ) = (γ 1 (θ) cos α + γ 2 (θ) sin α ) T r, (8)
where

γ 1 (θ) = s r (θ) s i (θ) , γ 2 (θ) = -s i (θ) s r (θ)
, r = r r r i is a 2N -real-valued noise vector following a centered Gaussian distribution of covariance I 2N /2 and α is the phase of the product s(θ) H r in [0, 2π]. It can be shown [START_REF] Develter | On the false alarm probability of the normalized matched filter for off-grid target detection[END_REF] that max α∈[0,2π] X(α, θ) = Y (θ), so that the P F A to characterize in this paper can be written as:

P F A = P max θ∈D Y (θ) > w 2 = P max α,θ∈[0,2π]×D X(α, θ) > w .
(9) X is easily shown to be stationary and we denote by Λ the covariance matrix of its gradient ∂X ∂α , ∂X ∂θ .

The excursion set A w (X) associated with X for a threshold w is defined as the set of parameters such that X(α, θ) exceeds w [START_REF] Adler | On excursion sets, tube formulas and maxima of random fields[END_REF]:

A w (X) = {(α, θ), X(α, θ) > w} . (10)
In [START_REF] Adler | Random fields and geometry[END_REF], it is shown that E(φ(A w (X))), the expected Euler characteristic of A w (X) is a precise estimation of the probability of X(α, θ) exceeding w and so the desired PFA: for some c > 1. In this paper, we will not go into details about the Euler characteristic φ(.) of a set: it is a topological characterization of the set, and for 2D sets, it counts the number of connected components minus the number of holes. For more details, readers can refer to the corresponding chapters of [START_REF] Adler | Random fields and geometry[END_REF]. Our first result is based on the following fundamental theorem for random fields on R 2 :

|P F A -E(φ(A w (X)))| < O exp -c w 2 , (11)
Theorem 3.1 [START_REF] Adler | On excursion sets, tube formulas and maxima of random fields[END_REF] Let X be a zero-mean, isotropic stationary random field on the parameter space T ⊂ R 2 with gradient covariance matrix Λ = λ I and assume that ∂T , defined as the boundary of T , is continuously differentiable except at an at most finite number of points. Then

E(φ(A w (X))) = |T | ρ 2 (w) + |∂T | 2 ρ 1 (w) + φ(T )Ψ w σ , ( 12
)
where σ is the variance of X (constant since X is stationary),

ρ k (w) = exp(-w 2 /2σ 2 )λ k/2 (2π) (k+1)/2 σ k H k-1 w σ with H k the k-th
Hermite polynomial and Ψ denotes the tail of the standard Gaussian distribution functions. || denotes the volume.

Applying this general result to our specific problem, it can be shown that the following corollary holds:

Corollary 3.1.1 The P F A -threshold relationship of the off- grid MF tested on a cell D = [θ 1 , θ 2 ]
with steering vectors of size N is:

P F A = π(N 2 -1) 3 (θ 2 -θ 1 ) w + 1 exp (-w 2 ) + O(exp (-c w 2 )) for some c > 1. ( 13
)
When D = [0, 1], only the first term in the parenthesis should remain: this special case was given in [START_REF] Hayward | Cfar detection of targets with unknown doppler shifts[END_REF] for Doppler detection.

As in the case of the off-grid NMF [START_REF] Develter | On the false alarm probability of the normalized matched filter for off-grid target detection[END_REF], the P F A of the off-grid MF is expressed as the sum of a constant term, equal to the P F A of the on-grid MF given in [START_REF] Liu | Robust detection in mimo radar with steering vector mismatches[END_REF] and corresponding to the behavior of the random field at the boundaries of the search domain plus another term proportional to the length of the cell.

Let us give a sketch of the proof of the previous corollary. Our result is achieved when derivating the mean Euler characteristic of the excursion set of X defined in [START_REF] Bosse | Adaptive subspace detectors for off-grid mismatched targets[END_REF] and apply [START_REF] Develter | On the false alarm probability of the normalized matched filter for off-grid target detection[END_REF] in order to get the P F A . First, note that in our case σ 2 = 1/2. The search domain for the target parameter θ is an interval D = [θ 1 , θ 2 ]. Since for all θ ∈ D, X(0, θ) = X(2π, θ), T is shaped akin to an hollow cylinder in R 3 as represented on Figure 2a: it has no boundaries along the α-axis. It results |T | = 2π(θ 2 -θ 1 ), |∂T | = 4π and φ(T ) = 0 since the Euler characteristic of a cylinder is zero. As the attentive reader may have noticed, our process X is not isotropic as required in the hypothesis of the theorem, but what matters in this hypothesis is that X should be stationary on ∂T and its gradient covariance matrix Λ should be proportional to identity matrix i.e. Λ = λ I. It is possible to have such a matrix Λ rewriting s in a centrosymmetric form and using an adequate change of variables ν = π N 2 -1 3 θ.

Note that we could derive this relationship only because X is stationary under white noise. Unfortunately, this method cannot be applied to colored noise.

Estimating the P F A for STAP detection

In the case of STAP detection, the MF test quantity can also be rewritten as a Gaussian random field, now depending on three unknown parameters (angle, Doppler shift, and phase):

X(α, ξ) = (γ 1 (ξ) cos α + γ 2 (ξ) sin α ) T r, (14)
where the vector s used in γ 1 and γ 2 is now the steering vector for STAP detection defined in (3), and r is a 2N P realvalued noise vector. This calls for the use of Theorem 3.3.2 in [START_REF] Adler | On excursion sets, tube formulas and maxima of random fields[END_REF], in which the mean Euler characteristic of sufficiently "nice" random fields defined on a subset of R 3 is given:

Theorem 3.2 [START_REF] Adler | On excursion sets, tube formulas and maxima of random fields[END_REF] Let X meet the same conditions as in theorem 3.1, except it is now defined on T ⊂ R 3 and ∂T is twice Fig. 3: Theoretical relationship (13) (resp. ( 16)) (red) along an empirical (blue) obtained under white noise using 5 10 7 Monte Carlo draws with N = 10 (resp. N = 8 and P = 4). The on-grid relationship (5) (yellow) is shown for comparison purposes.

continuously differentiable except at edges or creases of finite length and vertices where the edges meet. Then:

E(φ(A w (X))) = |T | ρ 3 (w) + |∂T | 2 ρ 2 (w) + H(∂T ) π ρ 1 (w) + φ(T )Ψ w σ , (15)
where H(∂T ) is a curvature integral that equals the mean curvature if ∂T is smooth and involves all of the mean curvature, the lengths of the edges, and the angles between them otherwise: details are given in Lemma 8 of [START_REF] Worsley | Estimating the number of peaks in a random field using the Hadwiger characteristic of excursion sets, with applications to medical images[END_REF]. The other terms are defined as in Theorem 3.1.

Applying the previous theorem to X(α, ξ) leads to:

Corollary 3.2.1 The P F A -threshold relationship of the offgrid MF for STAP detection tested on a cell D

= [θ 1 , θ 2 ] × [µ 1 , µ 2 ]
with steering vectors of size N P is:

P F A = π 6 ∆θ∆µ (N 2 -1)(P 2 -1) 2w 2 -1 + π 3 ∆θ N 2 -1 + ∆µ P 2 -1 w + 1 e -w 2 +O(e -c w 2 ) for some c > 1, (16)
with ∆θ = θ 2 -θ 1 and ∆µ = µ 2 -µ 1 .

Once again, when D = [0, 1] × [0, 1], only the first term in the bracket remains. Three terms appear in the developed formula: the first term accounts for the acceptance zone inside the surface of the 2D cell, while the second and third terms are the same as in [START_REF] Scharf | Signal detection in Gaussian noise of unknown level: an invariance application[END_REF]. The second term translates what happens at the edge of the cells and is similar to what happens on a 1D cell. The corners of the cell explain the last term that is, once again, identical to the P F A for the on-grid MF given in [START_REF] Liu | Robust detection in mimo radar with steering vector mismatches[END_REF].

To give an idea of the proof, the shape of the parameter space is a bit more complicated in this case. The search domain for ξ is a

2D cell D = [θ 1 , θ 2 ] × [µ 1 , µ 2 ].
As X is still periodic in α, T can be viewed as a rectangular cuboid of dimensions 2π, ∆θ, ∆µ where its faces along the phase axis are stuck together, as drawn in Figure 2b. This leads to |T | = 2π∆θ∆µ and |∂T | = 4π(∆θ + ∆µ). The quantity H(∂T ) π is equal to the mean length of the four edges, i.e., 2π, and finally, we have again φ(T ) = 0. The previous remark on the isotropy of X still holds here.

NUMERICAL RESULTS

In order to verify (13) (resp. ( 16)), we computed an empirical P F A -threshold relationship with 5 × 10 7 Monte Carlo draws in the Fourier resolution cell [0, 1/N ] (resp. [0, 1/N ] × [0, 1/P ]), with steering vectors of fixed size N = 10 (resp. N = 8 and P = 4). The results are shown on Figure 3a (resp. 3b). In both cases, the relationship seems to be very accurate, except for P F A close to one where there is a noticeable mismatch.

CONCLUSION

In this work, we analyzed the statistics of the off-grid GLRT, which is more robust to the presence of off-grid targets than the classical Matched Filter. We have derived its P F Athreshold relationship under white noise for Doppler or angle and joint Doppler-angle STAP radar detection, but the relation is valid for all applications dealing with Fourier steering vectors as can be widely encountered in spectral analysis.

The next natural step would be to derive a similar relationship for colored noise distribution. This, however, is more involved due to the non-stationarity of the corresponding random fields.

Fig. 1 :

 1 Fig.1: STAP detection performance for off-grid targets drawn at random in a cell: (blue) oracle detection, (yellow) GLRT (6) and (red) Matched Filter (4). P F A = 10 -6 , N = 8, P = 4.

  One unknown parameter: in this case, T is a hollow cylinder. Two unknown parameters: in this case, T can be viewed as a rectangular cuboid with two opposite faces stuck together.

Fig. 2 :

 2 Fig. 2: Representations of the search domains T .