THREE-YEAR VARIABILITY OF ENERGY AND CARBON DIOXIDE FLUXES AT CLEAR-CUT FOREST SITE IN THE EUROPEAN SOUTHERN TAIGA

. Forest clearing strongly influences the energy, water and greenhouse gas exchange at the land surface - atmosphere interface. To estimate effects of clear cutting on sensible ( H ), latent heat ( LE ) and CO 2 fluxes the continuous eddy covariance measurements were provided at the recently clear-cut area situated in the western part of Russia from spring 2016 to the end of 2018. The possible effects of surrounding forest on the air flow disturbances and on the spatial pattern of horizontal advection terms within the selected clear-cut area were investigated using a process-based 3D momentum, energy and CO 2 exchange model. The modeling results showed a very low contribution of horizontal advection term into total turbulent momentum fluxes at flux tower location in case of the southern wind direction. The results of field flux measurements indicated a strong inter-and intra-annual variability of energy and CO 2 fluxes. The energy budget is characterized by higher daily and monthly LE fluxes throughout the entire period of measurements excepting the first two months after timber harvest. The mean Bowen ratio ( β=H/LE ) was 0.52 in 2016, 0.30 - in 2017 and 0.35 - in 2018. Analysis of CO 2 fluxes during the first year following harvest showed that the monthly CO 2 release at the clear-cut area consistently exceeded the CO2 uptake rates. The mean net ecosystem exchange ( NEE ) in the period was 3.3±1.3 gC∙m -2 ∙d -1 . During the second and the third years of the flux measurements the clear-cut was also a prevailed sink of CO 2 for the atmosphere excepting short periods in June and in the first part of July when daily CO 2 uptake was higher than CO 2 release rates. The mean NEE rates averaged for the entire warm period of corresponding years were 1.2±2.3 gС∙m -2 ∙d -1 in 2017 and 2.8±2.5 gC∙m -2 ∙d -1 in 2018, respectively. The mean ratio between gross primary production ( GPP ) and ecosystem respiration ( TER ) was 0.58 in 2016, 0.84 - in 2017 and 0.74 - in 2018.


INTRODUCTION
Clear-cutting is one of the most widespread logging practice, which can substantially transform energy, water vapor and CO 2 exchange between forest ecosystems and the atmosphere and affect the climate system at multiple scales. Effects of different logging practices as well as natural forest disturbances on vegetation-atmosphere interaction is a key topic of various experimental and modeling studies provided over the recent decades (Amiro et al. 2010;Masek and Collatz 2006;Olchev et al. 2009;Radler et al. 2010;Coursolle et al. 2012;Ma et al. 2013;Molchanov et al. 2017). Clear-cutting influences surface albedo, net radiation, surface roughness and consequently the energy flux partitioning into sensible (H), latent (LE) and soil heat fluxes (McCaughey 1985;Amiro 2001;Rannik et al. 2002;Kowalski et al. 2003). Amiro et al. (2006) and Williams et al. (2013) reported on significant decrease of evapotranspiration rate at clear-cut sites that can be observed during several years after logging. At the same time some studies showed that the Bowen ratio (β=H/LE) at the clear-cut sites can be quickly recovered (within one growing season) to pre-disturbance values (Matthews et al. 2017). Most of the recent studies focused on effects of clear-cutting on ecosystem-atmosphere interaction were dedicated to analysis of CO 2 balance changes in forest ecosystems induced by the clear-cut harvesting (Amiro et al. 2010;Aguilos et al. 2014;Grant et al. 2010;Paul-Limoges et al. 2015;Rodrigues et al. 2010;etc.). Numerous studies reported that clear-cutting lead to forest ecosystem change from CO 2 sink to CO 2 source in annual or growing season balances. Aggregated analysis of FLUxNET data sets for territory of the USA and Canada provided by Amiro et al. (2010) showed that forest ecosystems after clear-cutting require usually about 20 years to restore their ability to be a CO 2 sink in the annual balance. The time that is necessary for ecosystem after logging to reach a compensation point between carbon uptake and release rates is varied depending on local climate and geographical conditions. Also, most of available studies use a hypothesis that clear-cutting lead to the substan-tial decrease of gross primary production (GPP). The ecosystem respiration (TER) rate in general doesn't change substantially after the forest disturbance mainly due to rapid compensation of the decreased autotrophic respiration contribution into TER by increased heterotrophic one (Amiro et al. 2010;Paul-Limoges et al. 2015;Williams et al. 2014).
Most of previous experimental studies were performed at disturbed forest ecosystems in North America and Europe but effects of clear-cutting on ecosystem-atmosphere exchange in Russian boreal forests are still very poorly investigated and represented in a couple of experimental and modeling studies only (e.g. Machimura et al. 2008;Mamkin et al. 2016Mamkin et al. , 2019Molchanov et al. 2017). Boreal forests cover large areas in Russia and they are very sensitive to different natural and anthropogenic disturbances (zamolodchikov et al. 2017). It makes very necessary to investigate the consequences of the forest disturbances for biogeophysical and biogeochemical climate regulation functions of the forest ecosystems, first of all in the context of prediction the carbon balance of Russian forests and their influence on climate system. To derive the long-term variability of CO 2 fluxes in forest ecosystems an eddy covariance measuring technique is mainly used (Aubinet et al. 2012;Burba et al. 2013). The eddy covariance flux measurements are usually provided in undisturbed forest ecosystems over uniform vegetation canopy. The forest damaging can obviously lead to strong air flow disturbances that make very difficult the application of eddy covariance techniques for the flux measurements in such areas. Vegetation heterogeneity can lead to disruption of the basic assumptions used in the method and hence, to large uncertainties in measured fluxes. To improve the accuracy of the energy, water vapor and greenhouse gas (GHG) flux estimates the effects of non-homogeneous canopy should be taken into account in the flux analysis (Belcher et al. 2011). During the last time the influence of the clear-cuts on the air flows has been investigated using mathematical models in several studies (e.g. Frank and Ruck 2008; Sogachev et al. 2005; 199 GES 02|2019 Mamkin et al. 2016;Levashova et al. 2017). Within the framework of our study the temporal variability of the energy, water vapor and CO 2 fluxes during the three growing seasons following harvest using continuous eddy covariance flux measurements was analyzed and the possible influence of forest edges on the air flows within the clear-cut using 3D hydrodynamic turbulent exchange model was assessed.

Site description
The flux measurements were performed at the recently clear-cut forest ecosystem in the sustainable management zone of the Central Forest Biosphere Reserve (CFBR) in Tver region in the western part of Russia (56.44° N,33.05° E, 250 m a.s.l.) (Fig. 1). CFBR is located in the south-western part of Valdai Hills and its territory belongs to the humid continental climate zone (Dfb type according to the Köppen-Geiger classification scheme) (Peel et al. 2007;Kuricheva et al. 2017). The Climate Moisture Index (CMI) (Willmott and Feddema 1992), calculated as the ratio between annual precipitation and potential evapotranspiration, ranged between 0.3-0.4 that corresponded to moderately wet surface moisture conditions (Novenko et al. 2018). The forest vegetation is consisted of typical species of southern taiga e.g. Norway spruce (Picea abies), European white birch (Betula pubescens) and Eurasian aspen (Populus tremula) (Knohl et al. 2002).
The area of experimental clear-cut site is about 4.5 ha and it is situated at a flat plain with well-drained sod-pale podzolics soils. Organic carbon content in the upper soil horizons varied between 2.73 and 5.79%. It is surrounded by mixed forest stand with Norway spruce (Picea abies), European white birch (Betula pendula) and Eurasian aspen (Populus tremula). The site was clear-felled in March-April 2016. After logging the large amount of harvest residuals, stumps and litter were remained on the ground (Mamkin et al. 2016(Mamkin et al. , 2019. During the first months after logging the site was free of any vegetation. Active vegetation regeneration started after the ground defrosting in the second half of April. In June-August the vegetation was mainly represented by herbaceous vegetation with dominated starwort (Stellaria graminea), sowtit (Fragaria vesca) and wood-sour (Oxalis acetosella), as well as by a small number of juvenile aspen (Populus tremula). In the first half of August the leaf area index (LAI) of vegetation increased to 2.5 m 2 •m -2 and the grassy and woody vegetation reached 70-90 cm height while the height of the surrounding forest

Meteorological and eddy-covariance flux measurements
The 3 m tower for flux measurements at the clear-cut site was installed immediately after the timber harvest in April 2016. The tower location was chosen taking into account the dominating southern wind direction in spring and summer (Mamkin et al. 2016). The eddy covariance equipment includes openpath CO 2 /H 2 O gas analyzer LI-7500A (LI-COR Inc., USA) and 3-D ultrasonic anemometer WindMaster Pro (Gill Instruments, UK). The instruments were mounted on the tower at the height of 2.4 m above the ground.
The air temperature, relative humidity and precipitation rates were measured by the weather transmitter WxT 520 (Vaisala Inc., Finland) that was installed at the height of 2 m above the ground. Global and reflected short-wave solar radiation as well as longwave incoming and outgoing radiation was measured using 4-component radiometer NR01 (Hukseflux Thermal Sensors, The Netherlands) at the height of 1.9 m. To obtain the temperature and volumetric water content (SWC) of the upper soil layer four reflectom-eters CS655 (Campbell Sci. Inc., USA) were installed around the tower in the soil at the 10 cm depth. Soil heat flux was measured using three heat flux sensors HFP01SC (Hukseflux Thermal Sensors, The Netherlands) installed in the soil at the 5 cm depth.
The eddy covariance data was sampled at the frequency of 10 Hz using Analyzer interface unit LI-7550 (LI-COR Inc., USA). Meteorological parameters were collected with data logger CR3000 (Campbell Sci. Inc., USA) at the frequency of 0.

Data post-processing
All steps of data post-processing were performed according to guidelines for data analysis (Aubinet et al. 2012;Burba et al. 2013). Flux calculation from the raw data was carried out for 30-min time intervals using EddyPro data processing software (LI-COR Inc., USA), with implementing of all necessary corrections and statistical tests. The footprint characteristics were estimated using the model suggested by Kljun et al. (2004). Quality check included 0-2 flag policy (Mauder and Foken 2006). The fluxes were determined with corresponding storage terms that was estimated according to Migliavacca et al. (2009). After data post-processing all fluxes with flag 2, as well as the fluxes with flag 0 and 1, containing the spikes and measured under e.g. rainfall and dew events, weak turbulence and low wind, were removed from the final data sets. The u*-filtering procedure was implemented for estimation of net ecosystem exchange (NEE) of CO 2 . The mean threshold values of u* were 0.086 m·s -1 for 2016, 0.197 m·s -1 -for 2017 and 0.064 m·s -1 -for 2018, respectively. The u*-filtering, gap filling of the flux data and partitioning of NEE into TER and GPP was per-
A three-dimensional hydrodynamic model for momentum, water vapor and carbon dioxide exchange between a non-uniform land surface and the atmosphere A three-dimensional (3D) hydrodynamic model uses a 1.5-order closure scheme and it is based on a system of averaged Navier-Stokes and continuity equations for the mean wind speed components, as well as on the reaction-advection-diffusion equation for H 2 O and CO 2 transfer within the atmospheric surface layer (Garratt 1992;Wyngaard 2010;Mukhartova et al. 2015). The 3D wind speed components, and concentration of any green house gases (GHG), e.g. CO 2 , c = c (x, y, z,t) , are considered as functions of horizontal coordinates x, y, and vertical coordinate z, where x is the coordinate along the prevailed wind direction, and z is equal to the height above the ground surface. The model uses Reynolds's decomposition and expresses the wind speed and concentrations of any considered GHG as sums of their mean values and deviations: and c, -their deviations. In case of neutral atmospheric stratification the averaged Navier-Stokes and continuity equations can be written as: where ρ 0 is the density of dry air, δP is the mean pressure deviation from the hydrostatic distribution, ! F cor and ! F d are the Coriolis and drag force acting within the vegetation cover. The Coriolis force components can be written as: where Ω is angular velocity of the Earth rotation (7.29·10 -5 rad s -1 ) and is the unit vector along the axis of rotation. Taking into account the significance of the term η Z = sinψ ( ψ is the geographic latitude), the components of ! F cor can be written as: The drag force can be expressed as follows: where c d is the dimensionless drag coefficient (we assume in our study that c d =0.4), LAD is the leaf area density (m 2 •m -3 ), that describes the total area of vegetation elements (leaves, branches, tree trunks) per unit volume.
The 1.5-closure scheme assumes that the turbulent fluxes can be expressed using the turbulent kinetic energy E and turbulent exchange coefficient K as follows: where C µ is a dimensionless model parameter, and ε is the dissipation rate for turbulent kinetic energy.
The values of E and ε are found from the system of differential equations (Sogachev and Panferov 2006;Mukhartova et al. 2017;Olchev et al. 2017): where φ (φ=ε•E -1 ) is the supplemented function characterizing the scale of turbulence. The dimensionless constants σ E =σ φ =2 introduce the Prandtl number for turbulent  kinetic energy and turbulent Schmidt number for the function φ, respectively. The function P E describes the shear generation of turbulent kinetic energy and can be expressed as: The coefficients C φ1 =0.52 and C φ2 =0.8 are model constants. The term Δφ describes the increase of energy dissipation caused by air flow interaction with vegetation elements. In the first approximation it can be expressed as follows:

Meteorological conditions
The three year measuring period from spring 2016 to the end of 2018 is characterized by various weather conditions (Fig. 2)

Modelling of the wind field and turbulent patterns within the clear-cut area
To describe the spatial pattern of wind speeds and atmospheric fluxes, as well as to estimate the possible influence of forest edges on the air flows at the flux tower location the 3D hydrodynamic turbulent exchange model was applied. The main attention in our modeling experiments was paid to descriptions of the 3D wind and flux distributions within and Analysis of the spatial variability of horizontal and vertical wind speed components showed significant heterogeneity of the air flows within and around the clear cut area. It is also enhanced by complex shape of the clear-cut boundary (Fig. 1). Results of numerical experiments showed that the largest effect of the clear-cut on the air flows is mainly appeared along the windward forest edge whereas at the leeward forest edge this effect is less pronounced (Fig. 3). The leeward part of the clear-cut is characterized by prevailed downward air flows (with negative w) whereas the upward air flows (w is positive) are appeared at its windward part (Fig. 3). It is very important to point out that the vertical wind velocity around the tower location as well as its horizontal gradient at 3-4 m height above the ground surface was relatively small (close to 0 m•s -1 ) (Fig. 3c) that can be used as one of criteria indicating sufficient reestablishing of the air flow at tower location after its disturbance at the leeward forest edge.
To estimate effects of the air flow disturbances at forest edges to wind and turbulent conditions within the clear-cut area we also analyzed the spatial patterns of the horizontal and vertical momentum fluxes. Negligible horizontal turbulent flux is a requirement for representative eddy covariance measurements. Fig. 4 shows an example of vertical distribution of momentum fluxes along profile crossing the clear-cut are from south to north and passing through the tower location. Maximum horizontal turbulent fluxes are detected at leeward part of the clear-cut area and at windward forest edge. Moreover the high anomalies of horizontal turbulent fluxes are also observed above the clear-cut area in the atmospheric layers higher than 20-30 meters above the ground. Area of minimum (close to zero) horizontal turbulent fluxes is situated in the northern (windward) part of the clear-cut and its boundaries overlap the flux tower location (Fig. 4a). Analysis of vertical turbulent fluxes along selected profile at the height of flux tower is also showed their very low variability at the tower location. All these factors confirm obviously the representativeness of provided eddy covariance flux measurements under southern wind directions. In case of the northern winds the influence of the forest edge on wind and turbulence patterns at flux tower location is much higher that can obviously lead to significant uncertainties in flux estimations using the eddy covariance technique.

Temporal variability of energy fluxes
Analysis of the eddy covariance flux data showed that the monthly LE fluxes exceeded H fluxes during the entire three-year period of flux measurements (Fig. 5). The total LE fluxes integrated over the period from 06 May to 18 October are also higher than total sums of H fluxes (Table 1) (Table 1).
It is important to point out that the lowest values of LE were measured during the first year after the timber harvest. The previously conducted comparisons of the energy fluxes in undisturbed forest and at the clear-cut area (Mamkin et al. 2019) showed that the clear-cutting led to decrease of LE by 30% in the first growing season following harvest. Moreover, the H also decreased due to the harvesting by 22%. The decreasing trend in turbulent energy fluxes due to clear-cutting can be explained by increased albedo and consequently reduced surface available energy (available energy is a difference between net radiation and sum of soil heat flux and canopy energy storage term). It should be also taken into account that the low LAI values influence the transpiration and evapotranspiration rates. Moreover, the comparison showed that in spite of lower H and LE fluxes measured at the clear-cut site the mean seasonal β for the clear-cut and undisturbed forest was almost the same (Mamkin et al. 2019).
In year 2017 and 2018 LE grew by ~20% following active regeneration of grassy and woody vegetation at the clear-cut area, and β decreased from 0.52 to 0.30-0.35. The main differences between LE fluxes were observed between the first and the second years after the harvest although the weather conditions in the second and the third years were more contrasting. Fast recovery of β value to pre-disturbance values have been discussed in several studies (e.g. Amiro et al. 2006;Matthews et al. 2017;Williams et al. 2013). Particularly, Williams et al. (2013) also reported about decreased summer and autumn H and increased LE for period from the first to the third years following harvest at the clear-cut of Norway spruce (Picea abies) forest in Massachussets (USA).

Temporal variability of CO 2 fluxes
Assessment of the temporal variability of NEE over the entire period of flux measurements at the clear-cut area showed that the CO 2 fluxes are characterized by a large variability mainly governed by weather conditions (first of all, incoming solar radiation, air temperature) and amount of regenerated photosynthesizing vegetation (Fig. 6) Grant et al. 2010). It is mainly governed by GPP/TER ratio that is varied in our study between 0.1-0.9. The GPP/TER ratio is usually grown with vegetation recovery but this process can be interrupted by the short-term periods of GPP/TER decreasing caused by changes of key meteorological parameters (e.g. the high air temperature in 2018 resulted in strong increase of TER rate and decrease of NEE). The same effects were observed by Paul-Limoges et al. (2015) between the second and the third years following harvest at the Douglas-fir clear-cut in the Vancouver island in Canada and by Pypker and Fredeen (2002) between 5 th and 6 th years at the spruce-fir clear-cut   . 2015) reported that the clear-cut of poplar (Populus dettoides) in subtropical China became a CO 2 sink at the end of the first growing season following harvest. It is obvious that to explain adequately the temporal and spatial variation of NEE rate, effect of all possible factors that control GPP and TER rates at the clear-cut (including analysis of environmental condition, vegetation and soil properties, biomass of debris, etc.) should be investigated within aggregated experimental and modeling studies.

CONCLUSIONS
The results of our experimental and modeling studies showed a strong influence of clear-cutting on the energy and water vapor fluxes between forest ecosystem and the atmosphere. It is manifested in disturbance of the spatial air flow patterns, in change of microclimatological conditions and the energy and CO 2 fluxes. Decreased net radiation and higher albedo in summer period are resulted in lower LE and H fluxes. Sufficient soil moisture and regenerated vegetation promote higher LE fluxes comparing with H ones. Obtained NEE dynamics are consistent with the hypothesis that clear-cutting turns forest ecosystems from CO 2 sink to CO 2 source for the atmosphere for several years after logging. NEE was reached maximal values in the first year after the harvest and minimal -in the second one. GPP increased from the first to the third year while the TER decreased in the second and increased in the third year, respectively. Variability of GPP/TER ratio is well corresponded to results obtained in other post-clear felled boreal forest ecosystems of the young successional stages. Representativeness of our eddy covariance flux measurements were controlled by results of numerical experiments using a 3D hydrodynamic model. It was shown that the area at tower location under prevailed southern wind direction is characterized by low horizontal and vertical momentum fluxes as well as by small horizontal gradient of the vertical wind speed component. Disturbing effect of the forest edges on our flux measurements can be therefore neglected. The results obtained in the present study can be applicable for predicting the influence of deforestation on the regional weather conditions in boreal ecozone and for estimating its consequences for the climate system.

ACkNOwLEDGEMENTS
The study was funded by the RFBR and Russian Geographical Society according to the research project № 17-05-41127. It was also partially supported by the Presidium of the Russian Academy of Sciences (programs № 51 «Climate change: causes, risks, consequences, problems of adaptation and regulation» and № 41 "Biodiversity of natural systems and biological resources of Russia".