Influence of climatic conditions on Normalized Difference Vegetation Index variability in forest in Poland (2002–2021)

The influence of climate change on forest condition is noticeable. Forest ecosystem stress caused by climate change has already been manifested in several parts of Europe, including Poland. Thus, the main objective of this paper is to investigate for the entire area of Poland a long‐term trend and variability of forest greenness expressed as the Normalized Difference Vegetation Index (NDVI), derived from two decades (2002–2021) of remote sensing Moderate Resolution Imaging Spectroradiometer (MODIS) data. In the next step, selected meteorological elements – temperature (T), precipitation (P) and evapotranspiration (ETo), derived from ERA5‐Land reanalysis—were used to determine the influence of climatic conditions on the variability of NDVI in forests. The study documents the general greening of forests in Poland in 2002–2021. The greening is mostly visible in central‐eastern Poland, where the annual mean NDVI increased by 0.030 in 20 years, while it is weaker in the Baltic coast and in the southern edges of Poland (increase by 0.009 in 20 years). Overall, the positive, statistically significant trends in annual NDVI prevail over the negative, statistically significant trends and account for 32.5% of forest area, whereas the negative trends account for 3.9%. The study indicates an overall moderate impact of meteorological elements on variability of NDVI in forests in Poland. The most important factors affecting forest condition are P and ETo. The strongest correlations between NDVI and P and ETo reach 0.55 and are located in central Poland, in the form of a belt from western to eastern borders.


| INTRODUCTION
Forests are the predominant terrestrial ecosystem on Earth, covering approximately 31% of the world's land area (FAO and UNEP, 2020) and accounting for 3/4 of the gross primary production of the Earth's biosphere (Pan et al., 2013).
They contribute to regulating climate at many different levels, from local to global scale.At the local scale, forest's presence in a particular place can, for example, induce rainfall, while deforestation leads to increased occurrence of droughts (Wright et al., 2017).At the same time, forests play an important role in carbon sequestration and evaporative cooling processes on the globe.Hence, the forests' impact on mitigation of the climate change is indisputable.
On the other hand, the influence of climate change on forest condition is also noticeable.The global warming observed over the past decades has an impact on increased evapotranspiration and water shortage that consequently increase the severity and duration of droughts worldwide (Dai, 2013;Forzieri et al., 2017).Numerous studies (Forzieri et al., 2016;Haarsma et al., 2013) reported that climate change may increase the frequency of extreme weather events, not only droughts, but also devastating wind storms and wildfires in many regions, and hence may cause severe, long-term damage to forest ecosystems (Bryn & Potthoff, 2018;Hofgaard et al., 2012;Karlsen et al., 2017;Morin et al., 2018).The prolonged and repeatable severe weather conditions-like droughts-gradually impact the forest structure, enhancing the forest vulnerability to insects, pests or fungus pathogen attacks.Several studies highlighted that multiple stress factors affect forest resistance to negative phenomena (Camarero & Gutierrez, 2004;Ciais et al., 2005;Desprez-Loustau et al., 2006;Fuhrer et al., 2006).In Europe, forests are in large proportion monocultural and thus represent the simplified ecosystem, which is more vulnerable to drought stress and insect outbreaks (Klapwijk & Björkman, 2018;Schuler et al., 2017).At the same time, it has to be stressed that forests in Europe will experience an increase of extreme drought events over the next decades (IPCC, 2021).
Satellite remote sensing permits to assess the vegetation vigour based on the spectral properties of the plant.The vegetation reflectance in the visible wavelengths is controlled by plant pigments, while in infrared, it is driven by the internal cell structure of plants.Healthy plants heavily absorb the red spectrum of electromagnetic radiation and reflect near-infrared radiation.As plants undergo stress, less red radiation is absorbed and less near-infrared radiation is reflected.The most commonly used tool for assessing vegetation (including forests) is the normalized difference vegetation index (NDVI) (Adole et al., 2016;Buras et al., 2020;Huang et al., 2021;Rouse et al., 1974;Soubry et al., 2021).NDVI is a normalized transform of the near-infrared to red reflectance ratio, designed to standardize vegetation index values to between À1 and +1 (Didan & Munoz, 2019).When derived from thoroughly calibrated satellite-borne sensors, it remains a reliable ecological indicator (Huang et al., 2021).NDVI has a long, 50-year history in the research of vegetation vigour, but other vegetation products-enhanced vegetation index (EVI), gross primary production (GPP) and net primary production (NPP) among others-also gain popularity in recent years.EVI formula, apart from red and near-infrared radiation, uses also the blue radiation, in order to stabilize the index value against variations in aerosol concentration levels (Didan & Munoz, 2019).According to the research results, NDVI is more chlorophyll-sensitive, while EVI is more sensitive to canopy structural variations (Huete et al., 2002).Therefore, NDVI is more likely to reflect changes in leaf colouration as, e.g., in the course of premature leaf senescence under drought, whereas EVI may better reflect early leaf shedding (Buras et al., 2020).Although NDVI tends to saturate at high biomass levels (Gao et al., 2000), research shows that it reflects natural vegetation coverage better than EVI (Li et al., 2010).
Daily coverage of the data provided by the Moderate Resolution Imaging Spectroradiometer (MODIS), proved to be an extremely useful dataset for mapping and monitoring the land cover changes at the national and global scale since 2000.Several spectral indices-NDVI among othersderived from a time series of MODIS data coupled with climatic variables were examined for meteorological and agricultural drought detection and monitoring (Keshavarz et al., 2014;Park et al., 2019;Rhee et al., 2010) and for forest drought assessment (Buras et al., 2020;Schuldt et al., 2020) and forest mapping (Barka et al., 2019).The applied coupling methods were often based on single and multiple linear regressions between vegetation indices and climate elements.In Poland, there are number of studies on forest status and condition derived predominantly at local scale (Bartold, 2016;Bartold & Bochenek, 2019;Bochenek et al., 2015;Ochtyra et al., 2016;Piekarski & Zwolinski, 2014;Siwicki & Ufnalski, 1998), but none of them are regarding the area of whole country.
Thus, the main objectives of this research are (1) to investigate long-term trend and variability of forest greenness expressed as the NDVI values, derived from two decades of MODIS data, and (2) to couple the NDVI with climatic variables in order to determine the influence of selected meteorological elements on the variability of NDVI in forests.All this was made for the entire area of Poland and for the time period 2002-2021.In order to assess which forest type-coniferous or broadleaved-is more vulnerable to the variability of meteorological conditions, the analyses were carried out for all forests, as well as for coniferous and broadleaved forests, at the country scale and at the pixel level.

| STUDY AREA
Poland is the ninth-largest country in Europe, and the largest country in central Europe-it covers the area of 312,696 km 2 and extends from 49 to 54 50 0 northern latitude and from 14 07 0 to 24 09 0 eastern longitude.The climate of Poland is determined primarily by its location in central Europe, in mid-latitude zone.This affects both the amount of incoming solar radiation and main features of atmospheric circulation.
Due to the mainly latitudinal arrangement of the land relief in Europe (lack of orographic barriers), the Polish climate is mostly influenced by the western transfer of air masses, and thus indirectly by the Atlantic Ocean.The country's warm temperate climate is characterized mainly by thermally differentiated winters (mean annual temperature from À3.5 C in the north-eastern regions and in the sub-mountain and foothill regions to 1.5 C in the western regions), warm summers (mean annual temperature from 14.5 C in foothill regions to 19.5 C in the central part of the country) and the mean annual precipitation sums from 450 mm in the central regions to 1200 mm in the mountains (1991-2020) (Tomczyk & Bednorz, 2022).
Forests cover almost 9.5 million ha, that is, 29.6% of the area of Poland (Raport o stanie, 2022).They are mostly located in the northern and north-western parts of the country, as well as in the mountains in the southern edges (Figure 1b).Coniferous forests slightly prevail, as they account for 50.4% of the forest area, while the broadleaved habitats account for 49.6% (Raport o stanie, 2022).However, apart from mountain areas, where the share of spruce, fir and beech in stands species composition is larger, in most parts of the country, the predominant species is pine-it covers 58% of the forests area (Raport o stanie, 2022).The reason for that is that pine monocultures were mostly located on the areas with the poorest, often former agricultural soils.After the pine, the following species dominate in the species composition: birch, oak, spruce and beech, but their area share ranges from 6 to 8% (Raport o stanie, 2022).
Polish forests span across various climatic, geological and hydrological conditions.They cover areas with different land relief and different natural vegetation types.Thus, the country area is divided into 8 main nature-forest lands (Figure 1a), characterized by similar physical-geographical features.The natural conditions of forest production in each land are formed mostly by diverse climatic conditions (i.e., mean annual temperature, annual precipitation sum and length of growing season) and also by diverse topography and hydrological and geological conditions (Zielony & Kliczkowska, 2012).From the forest management point of view, the country area is also divided into 429 Forest Inspectorates (Figure 1b).

| Meteorological elements
To visualize the climatic conditions, we used gridded data originating from ERA5-Land reanalysis (Muñoz-Sabater, 2019;Muñoz-Sabater et al., 2021).The hourly data representing meteorological elements, which are generally known to have a major influence on forest productivity dynamics (Chu et al., 2019;Liu et al., 2015;Yang et al., 2019), that is, 2-m temperature (T, in C), precipitation (P, in mm) and evapotranspiration (ETo, in mm) was downloaded for the period 1971-2021.The P and ETo were additionally used to compute the climatic water balance (CWB = P À ETo) (Thornthwaite, 1948).Spatial extent of the data was 49 N to 55 N latitude and 14 E to 25 E longitudecorresponding to the spatial extent of Poland.The resolution of reanalysis data is 0.1 Â 0.1 .In the next step, the mean annual values of T, P, Eto and CWB were prepared for each grid cell of the ERA5-Land reanalysis.

| MODIS-NDVI vegetation index
We used a vegetation index (VI)-NDVI-derived from MODIS onboard Terra and Aqua satellites-products MOD13Q1 and MYD13Q1 (Didan, 2021a(Didan, , 2021b)).The MOD13Q1 and MYD13Q1 products were downloaded for the growing season (April-October) of 2002-2021, because Aqua started on 4th July 2002, while Terra has many data gaps in 2001.Using both Terra and Aqua data streams together provides a quasi-8-day product time series, because the data are processed 8 days out of phase at a 16-day interval (Didan & Munoz, 2019).MOD13Q1 and MYD13Q1 products are provided at the spatial resolution of 250 m.Each granule consists of 4800 Â 4800 pixels, so three granules were needed to cover the area between 49 N and 55 N latitude and 14 E and 25 E longitude-corresponding to the spatial extent of Poland and resulting in 1545 granules needed in total to cover the time period 2002-2021.
Together with the NDVI product, the corresponding pixel reliability and day of year products were used.
Based on the pixel reliability information, only the pixels indicated as good or marginal quality were retained.In the next step, each of the retained pixels was allocated to the respective month on the basis of the day of year information.Then a monthly maximum NDVI was calculated for each of the retained pixels.This approach is based on the logic that low-value observations are either erroneous or have less vegetation vigour for the period under consideration (Holben, 1986).
Having monthly maximum NDVI values, the annual NDVI values were prepared.There are many different approaches to calculating annual NDVI, from which the most common is calculating the mean out of maximum monthly values (Ghebrezgabher et al., 2020;Mao et al., 2012).However, in this paper, the annual NDVI was prepared as the median of the maximum monthly NDVI values of the growing season for the year in question.Similar approach was used in, for example, Romania (Pr av alie et al., 2022) and Germany (West et al., 2022).The rationale behind using median is that it represents the situation of the vegetation over the whole growing season better than, for example, sole maximum value.At the same time, median is relatively robust to outliers, which can occur, for example, due to sanitary cuts throughout the year.This is especially important when comparing annual NDVI values to mean annual values of meteorological elements.

| Forest masks
Forest mask was prepared on the basis of Corine Land Cover (CLC) 2006, 2012and 2018databases (CLMS, 2021).The CLC database provides information on land cover in 44 classes over European countries.The CLC uses a Minimum Mapping Unit (MMU) of 25 ha for areal phenomena and a minimum width of 100 m for linear phenomena (CLMS, 2021).As forest we considered three land cover classes: broadleaved forests (311), coniferous forests (312) and mixed forests (313).The forest layers for 2006, 2012 and 2018 were intersected and the areas that remained forest for the three periods composed the forest mask.For each MODIS pixel (grid cell 250 m Â 250 m), the percentage of forest coverage was calculated.The pixels with ≥50% of forest cover were taken for further analysis.In total, the final forest mask contained 1,849,881 pixels, representing the area of 115,616 km 2 .
Additionally, the masks for broadleaved forest and coniferous forest were used.They were prepared following the same steps, as used for forest mask, except that only the pixels with ≥80% of particular forest type were considered.The threshold of 80% coverage of broadleaved or coniferous forest was used to assure the homogeneity of forest pixels.Following these selection criteria, 174,243 pixels were retained as the broadleaved forest mask and 798,777 pixels were retained as the coniferous forest mask, representing the area of 10,890 and 49,924 km 2 , respectively.

| Statistical analyses
In order to determine the multi-annual variability of NDVI and meteorological elements, as well as the influence of the latter on NDVI, several research methods were used.In the first step, the spatial distribution of mean annual NDVI was presented.Then, the spatial distribution of the slopes of the trend in changes in annual NDVI was presented, and the significance of these slopes was assessed (α = 0.05).Next, the pixels with the significant slopes were categorized into four classes: (1) a strong positive slope, which is bigger than the median of all positive statistically significant slopes, (2) a moderate positive slope, which is smaller than the median of all positive statistically significant slopes, (3) a strong negative slope, which is smaller than the median of all negative statistically significant slopes and (4) a moderate negative slope, which is bigger than the median of all negative statistically significant slopes.All this was carried out at the pixel level.
In the next step, spatially averaged values of NDVI in forests were analysed.The spatially averaged values of NDVI were calculated as area averages of all NDVI values in the MODIS grid cells (i.e., all pixels) within the respective forest masks in Poland.The equations of a linear trend for all investigated forest types were generated and the statistical significance of these trends was assessed using the F test (α = 0.05).Also, the mean annual, spatially averaged NDVI was calculated for three time periods: whole period 2002-2021 and two decades-2002-2011 and 2012-2021.The mean annual NDVI prepared in this way was calculated for all forest types and the nature-forest lands.
The analysis of meteorological elements was prepared in the similar manner as for NDVI.The spatial distribution of mean annual T, P, ETo and CWB was presented.Also, the spatial distribution of the slopes of the trend in changes in annual T, P, ETo and CWB was presented, and the significance of these slopes was assessed (α = 0.05).All this was carried out at the pixel level.Next, spatially averaged values of meteorological elements were analysed and the courses of mean annual, spatially averaged T, P, ETo and CWB were prepared.The equations of a linear trend for these meteorological elements were generated and the statistical significance of these trends was assessed using the F test (α = 0.05).It was made for the same time frame as for NDVI-2002-2021, and additionally for the longer period 1971-2021, in order to strengthen the information coming from the climate trend analysis.
The influence of the variability of meteorological elements on NDVI was assessed with the methods of spatial correlations.To this end, the mean annual T, P, ETo and CWB were resampled to fit the MODIS grid cells (using nearest neighbour interpolation method) and masked for all investigated forest types and the nature-forest lands.The strength of the correlation between the mean annual NDVI and meteorological elements in forests was assessed using the Pearson correlation coefficient, expressed by the formula r ¼ cov xy= S x S y , where cov xy is the covariance in the bivariate distribution of the variables x (mean annual values of a respective meteorological element, averaged spatially over the entire reference area) and y (mean annual NDVI, averaged spatially over the entire reference area), S x and S y are the standard deviations in the marginal distributions of the variables x and y respectively.All this was made for all forest types and nature-forest lands reference areas.The significance of linear correlations calculated in this way was assessed at two significance levels of α = 0.05 and α = 0.1.Spatial distribution of the correlation coefficient in forests was also prepared-this time the correlation coefficient was calculated at the pixel level.In order to better visualize the spatial distribution of the correlation coefficients between NDVI and meteorological elements, the bivariate choropleth maps were used, with Forest Inspectorates' areas (Figure 1b) as primary fields.In this method, each primary field is coloured according to the level of two features.For each Forest Inspectorate, the number of pixels with statistically significant positive and negative correlation coefficient was calculated.The colour on the final map indicates twofold information: the first feature is the percentage of forest area covered with pixels with significant positive correlations, and the second feature is the percentage of the forest area covered with pixels with significant negative correlations.

| Variability of NDVI
Multi-annual mean NDVI in Polish forests equals 0.804, but it differs across the country and between forest types.The biggest values occur in the sub-mountain and foothill regions in the south (0.834 for Land VII and 0.855 for Land VIII), and also in the northern regions (Figure 2 and Table 1).Central regions have lower values of annual mean NDVI (0.779 for Lands III and IV).Broadleaved forests have slightly bigger values of NDVI (0.841) than coniferous forests (0.791) (Table 1).
In the period 2002-2021, a general greening of forests is noticeable.The mean annual NDVI increased by 0.017 in total, during the period 2002-2021.The biggest increase was in the central Poland and reached 0.030 in total, while the lowest was in the Baltic coast and in the Sudeten Mountains (0.009 in total).Even though the NDVI increase is rather small in absolute terms, the linear trend in changes in mean annual NDVI is statistically significant for all investigated forest types (Figure 3).Similarly, all forest types and all nature-forest lands have bigger values of mean NDVI in the second decade (2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020)(2021) than in the first (2002-2011) (Table 1).In fact, the strong and moderate positive trends account for 32.5% of forest area and appear mostly in the central-eastern and southern parts of Poland (Figure 4).
In these regions, the mean annual NDVI increased by 0.020 to 0.030 over 20 years.On the contrary, there are also places with statistically significant negative trends, which account for 3.9% of forest area.Negative trends appear mostly in the north-western Poland, for example, in Szczecin Coastland (Figure 4.1) and large pinedominated forest complex in the Note c river basin (Figure 4.3), as well as in the large part of the Pisz Forest in northern Poland, dominated by the pine forest (Figure 4.2) and in many single, separate locations.Negative trends presented in Figure 4.4 are not related to forest dieback.These negative trends can be explained by a heavy tornado that occurred on the 11th of August 2017 (Ho sciło & Lewandowska, 2018).The pixels with negative trend are located on the edges of the affected  area.It has to be noticed that heavily damaged forest was excluded from our analysis, because in the CLC 2018, it was classified as transitional woodland and shrub (class 324).

| Variability of meteorological elements
Considering meteorological elements-T, P, ETo and CWB-the important climatic trends are noticeable.The spatial distribution of mean annual T reflects the direction of increase of continentality of climate-the biggest temperatures (>10 C) are in western Poland, while the lowest (>7 C), apart from the mountains, are in the north-eastern region.The spatial distribution of P shows relatively large mean annual sums in the Baltic coast (>750 mm) and in the sub-mountain and foothill regions in the south.Lower sums of P are in the central Poland (600-650 mm).The mean annual ETo has rather homogenous distribution over Poland, with mean values around 600 mm, while the CWB reflects mostly the distribution of P. The positive values of CWB are observed over most of the country, meaning that more water is delivered to the surface than evaporates.However, there are some small, separate locations in the central Poland, where the CWB is negative (Figure 5).The slope of the trend in changes in mean annual T is positive over the entire country, which means that the mean annual T increased by 1 to 1.6 C during the period 2002-2021 (Figure 5b).The largest, statistically significant average increase in T was observed in the southern and eastern parts of the country.In general, the increase in T is confirmed by the multi-annual trend in the longer period of 1971-2021 (Figure 6a), although the trend for the last 20 years is steeper, with a slope of 0.066 CÁyear À1 for 2002-2021 compared with 0.039 CÁyear À1 for 1971-2021.The statistically significant increase in ETo is reported mostly in central and eastern regions (Figure 5f), while the slopes of the trend in changes in P and CWB appear insignificant in the whole country (Figure 5d, h).Interestingly, the direction of the trends of P and CWB is positive (though insignificant) in last 20 years, but negative (and significant for CWB) in the longer period 1971-2021 (Figure 6b, d).

| Relationship between NDVI and meteorological elements
Statistical analysis of correlation between mean annual NDVI and meteorological elements reveals some interesting relationships.At the country scale, the correlation coefficient between NDVI and T is weak and statistically insignificant for each forest type and each nature-forest land (Table 2).However, a close-up look at the significant correlations at the pixel level shows that statistically significant negative correlations occur mostly in western Poland (Figure 7a).Additional cluster of Forest Inspectorates with more than 10% of the forest area covered with statistically significant negative correlations is in the eastern edge of Poland, the Białowieża Forest region (Figure 7a).
Correlations between NDVI and the other meteorological elements show an entirely different picture (Table 2).At the country scale, strong and statistically significant positive correlations between NDVI and P are in central and western regions (0.45 in Land III, 0.49 in Land IV and 0.53 in Land V), as well as in the south-western region, in the Sudeten Mountains (0.55 in Land VII).In Sudetes, there is a very strong correlation between NDVI and P for the broadleaved forests (0.65 in Land VII).Coniferous forests, which occupy more than a half of the forested areas in Poland show a moderate correlation between NDVI and P (0.38, α = 0.1) with regard to the whole country.The lowest values of correlation coefficient are in the Carpathians (south-eastern region).Interestingly, there is a cluster of Forest Inspectorates with more than 10% of statistically significant negative correlations (Figure 7b).
Similarly to P, ETo also correlates well with NDVI.Here, the strongest correlations of 0.55 (Land III) and 0.51 (Land IV) are observed in the central regions.Similarly, as in the case of P, coniferous forests show the moderate correlation between NDVI and ETo (0.39, α = 0.1) for the whole Poland.Sudetes is quite an interesting region, because the correlation between ETo and NDVI is negative (and insignificant).A closeup look at the significant correlations at the pixel level shows that in Sudetes, there is a cluster of Forest Inspectorates with more than 10% of the forest area covered with statistically significant negative correlations (Figure 7c).
The relationship between NDVI and CWB, an element derived from P and ETo, shows similar correlation coefficients to previously described between NDVI and P.However, here the correlations are mostly weaker than the ones obtained with P.
In turn, a close-up look at the significant correlations between NDVI and P, ETo and CWB at the pixel level shows a similar pattern in all three cases.Forest Inspectorates with many significant positive correlations are clustered in the central belt, from western to eastern borders (Figure 7b-d).The northern regions as well as the sub-mountain and foothill regions in the south have mostly insignificant correlations.
This study examines the long-term trend in forest greenness expressed as NDVI and in meteorological elements-T, P, ETo and CWB-and illustrates the correlation between variability of NDVI and meteorological elements.In Poland, values of the slopes of the trends in changes in meteorological elements, especially T and P, are similar to the ones reported elsewhere.Numerous studies report similar, significant temperature rise in Poland (Degirmendži c et al., 2004;Marsz & Styszy nska, 2022) and in Europe (Krauskopf & Huth, 2020;Pr av alie et al., 2022;van der Schrier et al., 2013).In turn, the long-term trend in changes in annual sums of P reported in this study is insignificant for the long period of 1971-2021.Similar, insignificant (positive or negative) trends were also reported in other long time periods, for example, in 1951-2000, the positive insignificant 2.9 mmÁyear À1 trend in P over Poland (Degirmendži c et al., 2004), while in 2001-2018, the negative insignificant À1.3 mmÁyear À1 trend in P over Poland (Ziernicka-Wojtaszek & Kopci nska, 2020).Despite the lack of clear trends in annual sums of P, many researchers highlight the increasing variability of annual sums, which results in an increase in the frequency of both drought and floods (Ziernicka-Wojtaszek & Kopci nska, 2020).An important feature regarding spatial distribution of P over Poland is a characteristic decrease of the mean annual precipitation sums that is observed in central Poland-this phenomenon is known as a rain shadow (Tomczyk & Bednorz, 2022).Combined with relatively homogenous values of Eto, it results in negative values of CWB in some separate locations in central Poland, which makes these regions especially vulnerable to vegetation disturbances.
This study documents the general greening of forests in Poland in the period 2002-2021.It is in line with the results of many research studies worldwide (Ding et al., 2020;Guo et al., 2018;Liu et al., 2015;Yang et al., 2019) and in Europe (Kempf, 2023;Pr av alie et al., 2022).Kempf (2023) reported a strong, statistically significant greening in central Europe in 2001-2021, with the strongest positive trends of NDVI occurring in winter.Guo et al. (2018) also reported positive trends of NDVI in European countries in 1982-2015, but these trends were not strong.Liu et al. (2015) observed the increase in NDVI by 0.46 Â 10 À3 per year from 1982 to 2012 globally with decadal variations.For many regions of the world, including central Europe, they observed a greeningbrowning-greening trend over the periods 1982-2004, 1995-2004 and 2005-2012, respectively.According to their study (Liu et al., 2015)   For instance, in Romania, a country with similar climate and similar features of the natural environment, the percentage of positive, statistically significant trends equals only 4% in the time period 1987-2018(Pr av alie et al., 2022)).As the time period of our study is much shorter (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020)(2021), it may suggest that the greening process accelerates in the 21st century, especially after the year 2005 or 2006 (Liu et al., 2015;Piao et al., 2006;Piao et al., 2011).This study indicates an overall impact of meteorological elements on variability of forest greenness in Poland.P and ETo are the most important factors affecting the forest condition, while the correlation between NDVI and T is rather weak at the country scale.Similar results can be found in many other studies (Liu et al., 2015;Piao et al., 2014;Yang et al., 2019), for example, Yang et al. (2019) showed that over 60% of vegetation evolution can be explained by rainfall, while temperature explains only 15%.Getting into more details, the growth response of, for example, Scots pine-the most popular coniferous tree species in Poland-showed a significant correlation with summer P, while decreasing association with T (Misi et al., 2019).The strongest correlations between NDVI and P (and NDVI and ETo) occur in central Poland.Many significant positive correlations cluster in the central belt, from western to eastern borders, while northern and sub-mountain and foothill regions in the south have mostly insignificant correlations.This pattern reflects the locations of the decreased amounts of the soil moisture (Somorowska, 2022;Urban et al., 2022), induced by decreased amounts of P. All this makes the central Poland's forests especially vulnerable to evaporative stress and precipitation-dependent to a very high extent.What can rise doubts are the mostly positive correlations between NDVI and ETo, meaning that the more water evaporates, the better the forest condition is.However, in Poland ETo is strictly related to the amount of the P.
In general, at the country scale, the correlation between NDVI and T is weak and insignificant.However, this result can be biased due to the spatial variability of the Polish forest.The weakened spatial relationship between trends in NDVI and T is supported also by Piao et al. (2014) and Liu et al. (2015).They showed that the relationship between the NDVI and interannual temperature variability has weakened substantially in the northern mid and high latitudes since 1980, possibly due to extreme climate events.Indeed, temperature rise increases the frequency of extreme weather events, droughts, wind storms and wildfires, and hence may cause severe, long-term damage to forest ecosystems, by enhancing the forest vulnerability to insects, pests or fungus pathogen attacks (Bryn & Potthoff, 2018;Hofgaard et al., 2012;Karlsen et al., 2017;Morin et al., 2018).Unlike for agricultural lands, where the spectral response to the influencing factor is quick, in the case of forests this process might be very extended in time, so a significant time lag in correlation between NDVI and meteorological element can occur.The lags in the forests' response to the negative changes in meteorological elements and hence the lags in the maximum values of correlation coefficients are also within the scope of researchers' interests (Mao et al., 2012;Moreira et al., 2019).
Another source of uncertainty is how the meteorological means were calculated.Among many ways to calculate climatology of meteorological elements (e.g., means out of summer or growing seasons months), this paper uses annual means.It should be remembered that the climatology calculated on the basis of different months might have different and even contrasting effect on the vegetation condition (Chen et al., 2018;Kern et al., 2022;Wang et al., 2003).
One also should not forget that meteorological elements are important, but not the only factors driving the forest greening/browning trends.Of course, temperature rise can serve as an indirect factor of vegetation growth, by, for example, increasing the atmospheric CO 2 concentration, which in turn is likely to cause an increase in vegetation productivity (Liu et al., 2015;Piao et al., 2006).However, there are other factors, such as forest site condition, forest habitats, the age structure of forests, reforestation, forest regrowth or forest management practices, that may drive greening trends.On the other hand, the reasons for browning trends can be non-climatic or anthropogenic as well, for example, insect outbreak or pathogen attacks.

| CONCLUSION
The results presented in this paper show in detail the variability of annual NDVI in Polish forests in the period 2002-2021, as well as the correlations between NDVI and selected meteorological elements, which are known to have a major influence on forest productivity dynamics: T, P, ETo and CWB.
Among forests in Poland, the biggest values of mean annual NDVI are observed in the sub-mountain and foothill regions in the south and also in the northern regions, while central regions have lower NDVI.In general, broadleaved forests have slightly bigger values of NDVI than coniferous forests.However, a division into two classes (broadleaved/coniferous) is rather coarse, and a more detailed classification into species should be investigated in the future.In the period 2002-2021, a general greening of the forests is noticeable, particularly visible in central-eastern Poland.
This study indicates an overall moderate impact of meteorological elements on variability of NDVI in forests in Poland.The most important factors affecting forest condition are P and ETo.The strongest correlations between NDVI and P and ETo are in central Poland, in the form of a belt from western to eastern borders.It is the region where forests are especially vulnerable to evaporative stress.Surprisingly, the correlation between NDVI and T is weak and insignificant in the whole country, although there are some significant negative correlations in western Poland.Still, one needs to remember that there might be significant time lags in maximum values of correlation coefficient between NDVI and meteorological elements, as forests' response to the negative changes in meteorological conditions is slow.That is why further research is needed, based on increasingly longer measurement data series.
The research presented in this study fills the knowledge gap on if and how strong is the response of Polish forest ecosystem to climate variability.However, the obtained results do not limit to Poland only.Identifying changes and variability in the condition of forests-the predominant terrestrial ecosystem on Earth-is particularly important in the context of the recent climate change and forests' impact on mitigation of the climate change.

F
I G U R E 3 The course of mean annual NDVI in Polish forests (all forests-grey line, broadleaved forests-red line, coniferous forests-green line) in the period 2002-2021.The graph also shows the linear trends (dashed lines), the trend line equations, the coefficient of determination R 2 and the statistical significance level p. F I G U R E 4 Slopes of the trend in changes in mean annual NDVI (NDVIÁyear À1 ) in Polish forests in the period 2002-2021.The values of the trend's slopes are divided into 5 classes, according to their statistical significance at the significance level of α = 0.05 and slope's steepness (1 insignificant and 4 significant classes).Selected areas are enlarged below the main picture.F I G U R E 5 Mean annual T (a), P (c), ETo (e), CWB (g) and the slopes of the trend in changes in these elements (b, d, f, h, respectively) over Poland during the period 2002-2021.The shading indicates the statistically insignificant slopes at the significance level of α = 0.05.F I G U R E 6 The course of mean annual T (a), and mean annual sums of P (b), ETo (c) and CWB (d) over Poland in the period 1971-2021.The graphs also show the linear trends in the period 1971-2021 (dashed lines), and the linear trends in the period 2002-2021 (solid lines), the trend line equations, the coefficient of determination R 2 and the statistical significance level p.
in central-eastern Poland, where the annual mean NDVI increased by 0.030 in 20 years, while it is weaker in northern, western and southern Poland.At the same time, the biggest values of mean annual NDVI are in the sub-mountain and foothill regions in the south, and also in the northern regions, suggesting that the spatial differentiation of the forest condition over Poland is gradually decreasing.Overall, the positive, statistically significant trends in annual NDVI prevail over the negative, statistically significant trends, and account for 32.5% of forest area, whereas the negative trends account for 3.9%.It is much bigger proportion than in other studies.

F
I G U R E 7 Percentage of forest area covered with pixels with statistically significant (α = 0.05) positive and negative correlations between NDVI and T (a), P (b), ETo (c), CWB (d) in individual Forest Inspectorates in the period 2002-2021.
T A B L E 1 , the 2005-2012 trend of NDVI for central Europe was 1.352 Â 10 À3 per year, which is comparable with 0.088 Â 10 À3 per year obtained in this study for 2002-2021.The greening is mostly visible T A B L E 2 Correlation coefficients between NDVI and individual meteorological elements in Poland and in each nature-forest land in the period 2002-2021.Statistically significant values at the significance level of α = 0.05 (α = 0.1) are in bold (italic). Note: