Wind Damage and Temperature Effect on Tree Mortality Caused by Ips typographus L.: Phase Transition Model

: The aim of this study was to develop methods for constructing a simple model describing tree mortality caused by Ips typographus L. using a minimum number of variables. We developed a model for areas spanning natural mountain forests in the Tatra National Park (Slovakia) and the Šumava National Park (Czech Republic), and in managed Czech forests located in four areas varying in environmental conditions. The model describes the time series of tree mortality dynamics caused by I. typographus using two submodels: a long-term dynamics submodel, and a short-term dynamics autoregressive distributed lag(ADL) model incorporating a two year delay and temperature variable averaged over the April-May period. The quality of ﬁt for our models ( R 2 value) ranged from 0.87 to 0.91. The model was formulated to capture the average monthly temperature effect, a key weather factor. We found that for high-elevation stands located at least 1000 ma.s.l., forest damage was predominantly inﬂuenced by May temperatures. For lower-elevation managed forests with warmer climates, the weather effect was insigniﬁcant.


Introduction
Outbreaks of xylophagous insects, one of the key economically important pests affecting coniferous stands, are typical for forests of Scandinavia, North of the European part of Russia, and mountain forests of Central Europe [1][2][3][4][5][6][7][8][9][10]. Over the past decades, extensive bark beetle outbreaks occurred across large areas in Central Europe, Western United States, and Canada causing extensive tree mortality rates and affecting forest ecosystem functioning. Bark beetle outbreaks often follow various natural disturbances such as    Tatranská Javorina is a protected district of the Tatra NP located on the north-oriented slopes. The deep valleys of study area are dominated by autochthonous mountain Norway spruce.
Populations of I. typographus in the study area are largely univoltine, although favourable weather conditions may allow beetles to complete a sister generation [47].
The strategy of non-intervention has not been consistently followed since bark beetle outbreaks erupted at the beginning of the 1990s. Tatra NP pest management approaches were modified several times in the study area, changing from a non-intervention regime prior to 1994, to intensive forest protection management, including application of trap trees, insecticides, and sanitation cutting in the period of 1995-1996. During 1997-2003, implementation of pest control measures again deviated from non-management approach to sanitary cuttings and intensive application of pheromone trap barriers [48]. From 2004 onward, 1650 ha (55%) of studied forest area gained highest protection level and has not been subjected to any human interventions. Two bark beetle outbreaks occurred in area over the last decades: first was in 1990s, and second after 2004. Regarding case of this study area, data and data acquisition procedure are described in more detail in [4].

Šumava National Park
The second non-management study area is in Šumava NP (Czech Republic). The forests, mainly spruce and mixed forests, cover more than 85% (54,439 ha) of the area. Zonal spruce forests (i.e., bog woodlands and mire spruce woods) grow naturally in elevations above 1050 ma.s.l., and occur in mire areas and inversion locations. Large areas of spruce forests have been subjected to significant natural disturbances in last few decades [49,50]. Šumava NP is characterized by a diverse mosaic of old-growth forests, windthrow areas, forests impacted by bark beetles, and areas influenced by past traditional forestry.
Local climate is transitory, with oceanic and continental climate influences throughout year, minute temperature variations, and relatively high rainfall.
Long-lasting debates on zoning and optimal strategies to protect mountain spruce natural stands in this area [51,52] escalated after hurricanes Kyrill and Emma occurred in 2007 and 2008, respectively. Subsequent bark beetle outbreak during 2008-2012 substantially affected NP forests. In January 2007,~20% of forest area was granted a non-intervention status. Only environmentally sound forest protection activities, such as debarking of lying or standing trees, with all dead wood retained at a site, were allowed on 10% of forest area. In the rest of the Šumava NP territory, to a certain extent, common forest management activities of varying intensity could be implemented. However, chaotic political interventions and management instability eroded this concept [53].

Karlovy Vary Division
The Karlovy Vary division is in western Bohemia and includes most of the Doupovské Mountain. The average annual air temperature is 6 • C, and the average annual total precipitation is 600-800 mm [54]. The forest stands are composed of 60% conifers (mainly Norway spruce). Spruce bark beetle has been in the endemic state a long time here. More important is damage by wind, the most serious of which was in 2007 by hurricane Kyrill.

Horní Planá Division
Horní Planá division is in the South Bohemian region. Spruce stands dominate in the study area. The stands are managed by the Military Forests and Farms of the Czech Republic. The average annual air temperature is 4-5 • C, and the average annual total precipitation is 800-1000 mm [54].
Tree mortality caused by bark beetles has existed at a low level for a long time. However, situation changed after the hurricane Kyril windstorm in January 2017, which destroyed ten-fold of annual harvest volume overnight, providing an ample breeding base for bark beetle development.

Lipník nad Bečvou Division
Lipník nad Bečvou division, another of the Czech army's training polygon since 1946, is in northern Moravia. The average annual air temperature does not exceed 5-6 • C. The average annual total precipitation ranges from 700 to 800 mm [54], but only in the last decade. Wooded area occupies 85% of the division, with spruce-monodominant stands covering~23,000 ha. A low degree of static stability of stands is responsible for frequent wind disturbances occurring in this area. Furthermore, climate change and wind-induced severe bark beetle outbreak led to a collapse of spruce forests in this area. Hořovice division, is a region in central Bohemia, with an average annual air temperature of 5-6 • C, and average annual total precipitation of 600-800 mm [54]. A prevailingly forested Brdy area occupies 95% of the division territory. Forest stands are dominated by Norway spruce covering~28,000 ha, or about 70% of the area. The bark beetle population has existed in an endemic condition for a long period of time. However, the consequences of hurricane Kyril were similar to those observed in the Horní Planá division. Despite a lower amount of damaged wood (one-fold of total annual harvest volume), the windfall created abundant breeding material for bark beetle proliferation.

Tree Mortality Data
The volume of wind-damaged wood has been recognized as an important driver of spruce ecosystem dynamics because acutely stressed, uprooted, or windfallen trees provide food resources for bark beetle development and population amplification [6,[22][23][24][25][26]. As direct assessment of bark beetle population density is complicated, we used a time series of the volume V (t) of bark beetles-induced wood damage (tree mortality) in year t as an indirect indicator of I. typographus population dynamics.
Wind and beetle-caused wood damage (tree mortality) data were sourced from field surveys conducted by organizations that manage forests in our study areas in Tatra NP and Šumava NP, as well as from the Military Forest and Farms state enterprise.

Meteorological Data
To ensure the accuracy of weather condition estimates, we used the data obtained by meteorological stations located inside, or in close proximity to, our study area (Table 1). For each study year, daily records of temperature were averaged over monthly periods. To understand the potential effects of climatic conditions on tree mortality, we used mean monthly air temperature.

Modelling Design
We considered outbreaks of forest insects by analogy with first order phase transitions (for example, boiling) in physical systems. For physical systems, first order phase transitions are characterized by a potential function-a function of free energy, which has two local minima F 1 and F 2 at some states of the object, separated by a local maximum F 12 -a potential barrier. The system can be in one of the states when it has a minimum energy value. During a phase transition under the influence of external factors, for example, changes in the temperature of the environment, the physical system passes from one stable state X 1 to another stable state. To define potential function, an approach based on the statistical study of the residence times of the system in different states can be used [55].
To design an outbreak model using this approach, we introduced system state function F(X), which we defined as an inverse value 1 f (X) of distribution density function f (X) of population density X values over a long time period, T. This function was sufficient and could be realized for all possible states of a system. If probability of a certain state X m realization for a sufficiently long observation time T is small, then value of function F(X m ) will be large. If state X m is realized frequently during period T, then value of F(X m ) is small.
Function F(X) potentially depends on a large number of various factors and it is difficult to write down dependences of F(X) on these factors. However, F(X) can be classified according to quantity, values, and positions of local minima and maxima. Function F(X) with one global minimum at the value X = X 1 will characterize a population with one stable state. Existence of a system near a stable state X 1 is associated with implementation of negative feedbacks in the system when its state deviates from value of X 1 . A state function with two local minima F(X 1 ) and F(X 2 ) (where X 1 << X 2 ) and one local maximum F(X 12 ) between local minima will be characterized as a population with possible transitions from state to state. If the value of F(X 12 ) is large, this state is observed very Forests 2022, 13, 180 6 of 21 rarely and system goes from state X 1 to X 2 and back, skipping over state F(X 12 ). Using the approach describing a first-order phase transition model proposed by [56], a state function can be written as an expansion in a Taylor series in some generalized indicator of the system state-order parameter q: where Z is an external factor influencing system state; a, b, c are some constants of expansion. A general view of the state function F(X) − F 0 for a bistable system with two stable states X 1 and X 2 is shown in Figure 2. ized as a population with possible transitions from state to state. If the value of F is large, this state is observed very rarely and system goes from state 1 X to X back, skipping over state ) ( 12 X F . Using the approach describing a first-order transition model proposed by [56], a state function can be written as an expansi Taylor series in some generalized indicator of the system state-order parameter q where Z is an external factor influencing system state; a, b, c are some constants of sion. A general view of the state function 0 ) ( F X F − for a bistable system with t ble states 1 X and 2 X is shown in Figure 2. Local minima of a state function can be determined from condition dF dq = 0. After differentiation of Equation (1), we obtained a critical value of q, and when this was reached a phase transition occurs: The action of such environmental factors as weather conditions, parasites and predators, amount of available food, within the framework of the first order phase transition model, can be considered as an action of an external field h: If external field h > 0 increases for a bistable state function (3), then depth of left local maximum will decrease proportionally to the external field (see Figure 2). In this case, at a certain value of the external field, a local minimum at X 1 will disappear, system will transfer to a state with a density X 2 , and an outbreak will occur. A change in sign of the external field and/or an increase in the influence of regulatory factors (for example, an increase in herbivore pressure, and a decrease in volume and quality of available food) will lead to a reverse jump in the state of the population, cessation of the outbreak, and system's return to a stable state.
When analysing fluctuations in state of insect populations, two types of fluctuations can be distinguished: high-frequency and low-frequency. In the case of high-frequency fluctuations, when state of population changes little over time, population will always be near one of stable states X 1 or X 2 . Low-frequency density fluctuations characterize situations when transitions from one stable state to another occur, and characteristic time of such fluctuations-so-called Kramer's time [57]-corresponds to the average time between two adjacent outbreaks.
It is rather difficult to collect data on changes in population density around X 1 , although probability of a population jump from X 1 to X 2 state depends on the amplitudes of these fluctuations. Oscillations around X 2 value can be measured fairly well due to their large absolute values. In this case, so-called ADL (autoregressive distributed lag) models can be used for modelling [15,46]: component of Equation (4), characterizing influence of populations state in k previous years (i.e., we can talk about a lag in time series of population), W(j) are external factors (weather conditions, food supply) in j-th year; k, u, a j , and b m are constants, k is order of autoregression and characterizes period of past seasons influence on current population state, n is period during which external factors affect current state of population. We used following data fitting quality indicators with model (6) [15]:

1.
A determination coefficient R 2 , characterizing proportion of variance explained by model. The adjusted determination coefficient (adj.R 2 ) is used to account for influence of variables number. In this case, variables number for different models is the same, so R 2 is shifted in relation to adj.R 2 by the same value.

2.
t-test for estimation difference of model coefficients from zero.

3.
F-test to test hypothesis that all coefficients of linear model are equal to zero.
Coefficients a j = ∂q(t) ∂q(t−j) characterize sensitivity of current population state to changes in previous seasons; coefficients b j = ∂q(t) ∂W(t−j) characterize sensitivity of current population state to the influence of external factors.
In Equation (4), following variable values are known: on the left side of equation, time series {X(i)} of system state dynamics, and on the right side-series {X(i − j} (in fact, the same time series of "forest-insects" state taken with some delay j), and a number of indicators W(i) selected from minimum discrepancies condition of model series to the time series {X(i)} of population count data. In fact, this means that Equation (4) should be considered as a linear regression equation, for which estimation of orders of k and u on the right side are needed; and then, using standard procedure for finding coefficients of linear regression equations, calculate coefficients a j and b m .
It is possible to use methods of autocorrelation analysis [58][59][60][61] and calculate partial autocorrelation function (PACF) to estimate order k of the AR-component of model. However, such a calculation is performed correctly only if the studied time series is stationary and standard deviations do not change over time [58,62]. If these conditions are not met, and there is a trend in population density values and (or) a time-varying range of density fluctuations, then it is necessary to transform the observed time series without losing the properties of the series of interest. Hence, it could be transformed into a Linear Time Invariant (LTI) series. A LTI series should have a time-constant mean, a time-constant mean variance, and a time-constant oscillation frequency of variable. To assess the quality of the model, four criteria were used: value of determination coefficient R 2 should be close to 1, values of regression coefficients in Equation (4) should be significant according to t-test, transformed series of accounting data, and model series should be synchronous (synchronicity is assessed by characteristics of cross-correlation functions between these series). We can say that it adequately described the dynamics of the studied population if the model meets all these criteria.

Data Processing
To construct ADL models to reduce variance of population state values, we transformed count data to a logarithmic scale of densities. Low-frequency filtering allowed us to identify long-term trends in densities, and high-frequency filtering was used to reduce possible counting errors. To isolate the high-frequency component, i.e., to select oscillations with higher frequencies, we applied a Hann filter [63] in the following functional form: where τ = 1 is the time between two adjacent counts of population state, represented by a formula with weighting coefficients 0.24, 0.52, and 0.24 for values ln (X − 1), ln X(t), and ln X(t + 1), respectively: Thus, to model dynamics of population using data from long-term records of tree mortality (wood damage), researcher can use, on the one hand, first-order phase transition models to explain outbreak phenomenon itself and, on the other hand, ADL-models to describe changes in the population state during an outbreak. Parameters of ADL model can be calculated using statistical package Statistica 10.

Model and Model Coefficients
As it is rather difficult to estimate population density for xylophagous populations, an indirect indicator is usually used-area S or volume V of bark beetle attacked wood in one year. In this case, to construct the state function F(ln V) instead of F(X), where F(ln V) = 1 f (ln V) and f (ln V)-distribution density function of values ln V. For a long time period T 1 , outbreaks are not observed, volume of wood attacked by bark beetle is small, probability f (V 1 ) ∝ T 1 T (were T-general time of stand observation) that stand is in an inter-outbreak state is high, and then value F(ln V 1 ) = 1 f (ln V 1 ) is small. Because the transition from a state with a low V 1 value to a state with a high tree mortality caused by bark beetle (volume of wood attacked by bark beetle), level V 2 occurs quickly, and F(ln V 12 ) value for this transition state V 12 is large.
Using data on the volume of wood attacked by bark beetles, we constructed a state function of "forest-insects" system as an inverse function of density distribution of tree mortality volume by caused bark beetle f (ln V) over a long observation period T (Figure 3).
cause the transition from a state with a low V1 value to a state with a high tree mortality caused by bark beetle (volume of wood attacked by bark beetle), level V2occurs quickly, and ) (ln 12 V F value for this transition state V12 is large. Using data on the volume of wood attacked by bark beetles, we constructed a state function of "forest-insects" system as an inverse function of density distribution of tree mortality volume by caused bark beetle ) (ln V f over a long observation period T (Figure 3).  As seen in Figure 3, outbreak phase can be characterized by value of ln V 2 , a potential barrier is characterized by value of ln V 12 , and phase of rarefied state can be characterized by value of ln V 1 .
The partial autocorrelation function (PACF) was used as an indicator for connection between ln V(t) and ln V(t − k). The PACF calculation is correct only for stationary time series, therefore, to start calculations at Karlovy Vary division study area, we selected linear trend ln VTR(t) = A − Bt of series {ln V (t)} and then will calculate PACF for time series {{∆ ln V(t)} = {ln V(t) − ln VTR(t)} (Figure 4).

EVIEW
10 of 23    As follows from Figure 5, the order of autoregression k = 2. Current fluctuations in the system state are not random and depend on two previous seasons. To assess the in- As follows from Figure 5, the order of autoregression k = 2. Current fluctuations in the system state are not random and depend on two previous seasons. To assess the influence of external factors (weather or volume of windfallen trees) on the system state, it is difficult to use a model of phase transitions and autoregression. Therefore, to assess effect of weather, we will choose values with maximum coefficient of determination R 2 for ADL model. Under these conditions, ADL model for first differences ∆ ln V(t) will have the following form: In this case, the maximum value of R 2 is achieved at the following values of coefficients of ADL model for investigated plot ( Table 2). Dynamics of tree mortality by I. typographus at the Karlovy Vary division study area and model calculations are shown in Figure 6.    Figure 6.  The synchronism of the empirical and modelled time series of wood damage caused by bark beetle is estimated using a cross-correlation function (CCF) as shown for Karlovy Vary division study area (Figure 7). Similar calculations and figures for the stationary component of time series, the partial cross-correlation function, the empirical and modelled time series of tree mortality caused by bark beetle, the cross-correlation function of empirical and model time series of damaged forest, were carried out for other study areas. Those data are available both on the intensity of wind damage and monthly weather conditions. The quantified values of ADL model parameters, which considered the effect of windthrows and air temperature in certain months for other study areas, are shown in Table 3. Similar calculations and figures for the stationary component of time series, the partial cross-correlation function, the empirical and modelled time series of tree mortality caused by bark beetle, the cross-correlation function of empirical and model time series of damaged forest, were carried out for other study areas. Those data are available both on the intensity of wind damage and monthly weather conditions. The quantified values of ADL model parameters, which considered the effect of windthrows and air temperature in certain months for other study areas, are shown in Table 3. Table 3. Model parameters to estimate the volume of wood damaged by bark beetle (tree mortality) during an outbreak in different study areas (for abbreviations refer to Table 2).

Variable Coefficient Std. Error t-Test p-Value
Tatra National Park  Parameter V in t − 2 proved to be significant in nearly every study area, except the Karlovy Vary division study area. Damage in the previous year was significant in every study area.
A series of field data on the dynamics of tree mortality caused by bark beetle is shown for the SA Tatra NP, Šumava NP, Horní Planá, Lipník nad Bečvou, and Hořovice divisions ( Figure 8).
Two outbreaks were observed in the studied period in Tatra NP. The first outbreak occurred in the mid-1990s, and the second outbreak was observed after 2004. The tree mortality remained relatively high for a long period of time. Such outbreaks in a period of less than 10 years are unusual for higher elevation forests. In Šumava NP, tree mortality caused by bark beetles remained relatively unchanged during the study period, although it was maintained at high levels.
The data for Horní Planá resembles a typical epidemic curve, with low fluctuation of tree mortality.
Tree mortality in the Lipník nad Bečvou study area resembles a typical permanent high level of tree mortality. The mature spruce stands were destroyed by the year 2018.
Tree mortality in SA Hořovice was much more stable.
Karlovy Vary division study area. Damage in the previous year was significant in eve study area. A series of field data on the dynamics of tree mortality caused by bark beetle is show for the SA Tatra NP, Šumava NP, Horní Planá, Lipník nad Bečvou, and Hořovice divisio (Figure 8).

Model Performance
According to calculation of PACF, order of autoregression for model of forest attacked by bark beetle populations is equal to 2 (Table 3). That is, current level of tree mortality caused by bark beetle is determined by levels of tree mortality caused by bark beetle in previous two years, and a positive feedback is observed between tree mortality caused by bark beetle in past year and current year: higher level of tree mortality caused by bark beetle in past year causes a higher level in current year. On the contrary, there is a negative feedback between tree mortality caused by bark beetle in year before last and current yearhigher level of tree mortality caused by bark beetle in year (t − 2), lower level is in year t. These relationships can be explained by increasing population density at a high level of tree mortality caused by bark beetle per year (t − 1) that enhances food resources available for next beetle generation, intensifying mass attacks on host trees. On the contrary, if food supply decreases during two seasons, trees damaged in year (t − 2) become unsuitable for colonization in year t. Model (6) indicates that an important driver of beetle-induced tree mortality dynamics is the volume of wind-fallen trees. Windfalls commonly occur in winter and mostly depend on tree growing conditions. Trees felled in the current year become a food resource for bark beetles in the following year. Weather conditions affect the dynamics of damage in different ways. For example, in the Šumava NP trial plot, it was found that current level of tree mortality caused by bark beetle significantly depends on average May temperature of year (t − 2). Sign of damage level from weather changes sensitivity coefficient is positive. It indicates that warmer weather in May is associated with a higher level of tree mortality caused by bark beetle two years later. The lag between temperature rise and level of tree mortality caused by bark beetle can be related to the physiological response of trees, manifested in decline of its resistance capacities after wind damage. In Karlovy Vary study area and Tatra NP study area, average monthly temperature in May in year t was significant. In remaining managed study areas at lower elevations, monthly average temperatures were not statistically significant.
The question arises: under what conditions will outbreak stop? From perspective of proposed approach, to stop an outbreak, it is necessary to reduce effects of external field to an extent that a minimum of local fluctuations in the current level of tree mortality caused by bark beetle is reached. In other words, outbreak ceases under impact of decline in bark beetle population density. Since model (6) is a second-order autoregression model (Yule model), spectral density p(f ) of plant damage time series according to model (6) will be described by Equation (7) [13,64,65]: Empirical and theoretical spectral density of tree mortality time series caused by bark beetle calculated according to Equation (7) are shown in Figure 9.
Frequencies f max = 0.142, which correspond to the maximum spectral density of time series of 3.64 for tree mortality caused by bark beetle, are close to the empirical and theoretical spectra and wave length L = 1 f max = 1 0.142 ≈ 7 years as indicated in Figure 9. As can be seen from Equation (7), the cyclicity of outbreaks depends on the values of parameters of ADL model. If the value of a 1 (susceptibility of ∂ ∆ ln V(t) ∂ ∆ ln V(t−1) volume of wood attacked by bark beetle in year t to the volume of wood attacked by bark beetle in year (t − 1)) increases, then the frequency of the spectrum maximum will shift towards lower frequencies, and the maximum of the spectral power will increase ( Figure 10, curve 2). If, on the contrary, the value of a 1 decreases, then the frequency of the spectrum maximum will increase, and the value of the maximum spectral power will decrease (Figure 10  As can be seen from Equation (7), the cyclicity of outbreaks depends on the values of parameters of ADL model. If the value of a1(susceptibility of wood attacked by bark beetle in year t to the volume of wood attacked by bark beetle in year (t − 1)) increases, then the frequency of the spectrum maximum will shift towards lower frequencies, and the maximum of the spectral power will increase ( Figure 10, curve 2). If, on the contrary, the value of a1 decreases, then the frequency of the spectrum maximum will increase, and the value of the maximum spectral power will decrease ( Figure  10, curve 3).

Estimation of Wind Damage and Temperature Effect on Tree Mortality Caused by Bark Beetle
How strongly do such factors as the volume W(t) of trees felled by the windfall and the temperature T(t) of the environment influence the dynamics of tree mortality caused by bark beetle? If we assume that there are no wind impacts on the stands and W(t) = 0, and changes in air temperature do not affect fluctuations in the dynamics of bark beetle population and level of connected tree mortality, then the stability of the tree mortality process caused by bark beetle will depend on the a 1 and a 2 components of model (6) [64]: Using data in Tables 2 and 3, we estimated realization of inequalities (Equation (8)) ( Table 4). Table 4. Values of inequalities (Equation (8)) and stability margin for series of studied tree mortality caused by bark beetle {∆ln V(t)} at the different survey areas.

Parameters and Conditions
In accordance with Table 4, for all study areas, except for Šumava, conditions (8) are satisfied. That is, a 1 and a 2 components in Equation (6) can be considered stable. Table 4 shows values of stability margin for different study areas. Also from Table 4, values of stability margin are small, and in general it can be concluded that in absence of a windthrow and elevated air temperatures, stability of tree mortality volumes caused by bark beetle and, therefore, dynamics of bark beetle population is rather small, and a loss of stability for bark beetle population dynamics is possible.
In addition to conditions (8), for AR models of forest insects population dynamics, there are important characteristics that make it possible to assess changes in population state under various external environment transformations, and proper population characteristics, such as stability margin η [66]. Stability margin characterizes proximity of this point to the boundaries of stability zone. For discrete systems, Mikhailov criterion and Mikhailov hodograph are used [67] to estimate the stability margin.
To assess stability of a 1 and a 2 components in Equation (6), it is necessary to write down the characteristic equation D(z) = z 2 − a 1 z + a 2 , replace variable z: z = 1+ν 1−ν and build Mikhailov hodograph D(jν) = D(z)| z = exp(jν) . (where j = √ −1). To construct the Mikhailov hodograph, the real and imaginary parts of polynomial are separated D(jν). The stability margin η is the circle radius centred at the point z = 0,which can be written inside Mikhailov hodograph [68].
By definition, value of margin stability η ≥ 0 and smaller value of η, provide a greater likelihood of a "breakdown" and a loss of system stability under external influences. Despite seemingly complicated calculation procedure, stability margin is estimated using a simple program in MATLAB package [68], and only values of coefficients of AR components of Equation (6) are required for calculations.

Model
Tree mortality by bark beetle can be considered as transitions between two stable (or metastable) states in the bistable system of the pest population: a state with a low-density level X 1 (and therefore a low level of damage V by the pest), and a state with a high level of population density X 2 (and therefore with a high level of tree mortality). In this case, the level of forest damage by the windfall plays the role of the external field [46]. The proposed model makes it possible to identify key factors in the dynamics of insect infestation in forests-the current values of damage after a windfall, average temperature of April or May, tree mortality caused by bark beetle in the previous two years. The proposed model of tree mortality does not require knowledge of bark beetle populations and tree state characteristics. Forest damage by wind, tree mortality by bark beetles and weather data are sufficient to build a model.
Further research is needed to model parameters of local tree mortality. A calculation showed that tree mortality by a bark beetle is not random and depends on both the history of past tree mortality and the current weather. However, using the proposed model for a short-term forecast of tree mortality dynamics is difficult and, in this case, it is necessary to predict future weather in a crucial seasonal time period for Ips typographus [34]. However, according to the data of the past tree mortality and the weather conditions in April-May, the risk of damage can be estimated by the order of autoregression k depending on the values of derivatives and maximum potential of a bistable system [46,58]. It should be noted that the proposed ADL models do not describe a long-term state with a low population density due to the lack of field data for these periods. To construct a general model of bark beetle population dynamics (or dynamics of forest tree mortality associated with it), data are needed both in the state of outbreak and in a stable state with a low level of density. In this case, it is possible to construct a more general model of population dynamics (and tree mortality) as analogues of a first-order phase transition model in physical systems with two stable states and oscillations around each of the states.

Model Coefficients: Environmental and Internal Factors
The coefficients of the model differ in the key weather factor: average monthly temperature. These differences may be due to differences in landscape characteristics and elevation of the study areas. For forest stands at high elevation (not lower than 1000 ma.s.l.) (Tatra National Park study area, Šumava National Park study area, and Horní Planá division), damage is influenced by the weather in May. As shown at Table 3 the t-test value and the p-level coefficients at variable T are significant for p ≤ 0.047 (Tatra National Park), p = 0.065 (Šumava National Park), p = 0.138 (Horní Planá division). For lower-lying forests in the Lipník nad Bečvou division study area (500-650 m a.s.l.) and in the Hořovice division study area (600-800 m a.s.l.), with warmer climates, the model coefficients for weather are insignificant (see t-test values and p-level with variable T in Table 3 for Lipník nad Bečvou and Hořovice divisions), that is, for these territories the damage intensity does not depend on the weather conditions. It should be noted that average temperatures of spring months for the test areas of Tatra National Park, Šumava NP, and Horní Planá division study area, are from 4.8 to 5.1 • C. For Lipník nad Bečvou division study area and Hořovice division study area at lower elevations, average temperatures for spring months were approximately 0.5 • C higher and ranged from 5.4 to 5.5 • C.
Comparing coefficients for variable ln V(i − 1), which characterize susceptibility ∂ ln V(i−1) of current tree mortality caused by bark beetle ln V(i) to the previous year's damage ln V(i − 1), it can be said that in forest stands at high elevations (approximately 1000 ma.s.l.), values of these coefficients were maximum, whereas for Lipník nad Bečvou study area (elevation approximately 500 ma.s.l.) values of these coefficients were minimal. Differences in population dynamics of I. typographus related to the different elevation of study areas were described by [69] for areas of Slovak Republic and by [70] for areas of Czech Republic. Lipník nad Bečvou study area is located at an especially (relatively)  [17,71]. Under acute or chronic drought stress the defence abilities of spruce are minimal [72]. The coefficients of the variable ln W for all plots, except for Tatra National Park, are approximately equal. The small value of this coefficient for the Tatra National Park study area indicates that the volume of trees which fell under windthrow did not affect tree mortality caused by bark beetles. This finding could possibly be explained by the dominance of high elevations and steep northern-oriented slopes in the study area. Such conditions are not favourable for colonization of downed trees by Ips typographus [73].

Conclusions
Models that account for autoregressive properties of population densities are well known [74][75][76][77]. However, the combination of phase transition models and ADL-models make it possible to describe different stages of forest insect's outbreaks. According to the first-order phase transition model, outbreak occurs because of system transition from stable state with low density X 1 to a state of high density X 2 . An outbreak occurs either with fluctuations around X 1 (this is possible if height of the potential barrier between states X 1 and X 2 is small), or under the influence of an external field, when state X 1 disappears. For stands damaged by Ips typographus, wood volume felled by wind acts as an external "field". In the absence of wind, the damage level fluctuates around X 1 value. After exposure to an external "field" and a jump into the outbreak state with a volume of damaged wood near X 2 , level of damaged wood fluctuates around X 2 (with effect of weather and windthrow) according to the proposed ADL-model.
Phase transition models of the first kind and an ADL-model were previously used to model population dynamics and outbreaks of phyllophagous insects [46]. For these models, pest population densities were used as a variable. Variables such as the density of parasites or predators (which are extremely difficult to account in the field) were not used for modelling and were replaced by negative feedback from previous density of simulated population.
The modelled system of a stand damaged by bark beetle differs from such classical systems as "predator-prey" or "resource-consumer" (where a population of bark beetle is considered as a consumer, and trees as an available resource). The system is described by only one variable-volume of damaged wood, and there is no data about bark beetle density. The concept of an outbreak as a first-order phase transition and an ADL-model of simulated variable oscillations (dead wood volumes) near the metastable state X 2 allows us to solve the problem of a lack of field measurement data.
Calculations performed within the framework of the proposed model have shown that the use of a fairly small amount of data from past seasons and the meteorological parameters make it possible to predict current forest damage caused by insects. As follows from the values of determination coefficients R 2 , the ADL model accountedfor at least 87% variance in the dynamics of forest stand damage. The proposed model has the potential to improve forecasting bark beetle population systems.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.