Meteorological Tower Observed CO₂ Flux and Footprint in the Forest of Xiaoxing’an Mountains, Northeast China

Since the time of the industrial revolution in the 1860s, with rapid development of the world economy and continuous improvement of human living standards, the global climate system has faced increasingly severe challenges, such as the problem of accelerated global warming (Nisbet and Myers, 2007). In the 21st century, the trend of global warming has intensified and the hottest weather ever recorded in the past hundreds of years has occurred over much of the world. The frequent occurrence of the El Niño phenomenon has also exacerbated the occurrence of climate-related disasters causing enormous economic losses for many countries. The Sixth Assessment Report of the United Nations Intergovernmental Panel on Climate Change highlighted that the sharp increase in atmospheric greenhouse gas concentration is the primary cause of global warming, and that CO₂ is one of the greenhouse gases that have greatest impact on global warming (Azar and Johansson, 2022). Forest ecosystems represent an important part of the terrestrial carbon cycle. Forest plants consume CO₂ and release O₂ through the process of photosynthesis. Therefore, forests have a strong carbon sink function throughout the year (Curtis et al., 2002), represent the largest carbon pool of all terrestrial ecosystems, and play an extremely important and irreplaceable role in maintaining the regional ecological environment and global carbon balance (Pan et al., 2011).

The eddy covariance method (Mauder et al., 2013; Charuchittipan et al., 2014), which has become recognized as a technical approach suitable for direct observation of CO₂ emission and absorption by a forest ecosystem (Falge et al., 2002; Ma et al., 2016; Baldocchi et al., 2018), has been adopted in scientific research all over the world (Richardson et al., 2010; Pita et al., 2013). However, owing to the influence of observation height, topography, wind direction, and wind speed, the value of the flux measured by the sensors on a flux tower represents only a certain limited surrounding area in practical application (Chen et al., 2009). Therefore, when performing flux research, it is necessary to determine the spatial representation of the flux measured by the sensors on a flux tower (Zhang and Wen, 2015).

Many models for calculation of the flux contribution area have been developed, e.g., the KM (Kormann and Meixner, 2001), FSAM (Flux Source Area Model; Schmid and Lloyd, 1999), Hsieh (Hsieh et al., 2000), and Horst–Weil models (Horst and Weil, 1992). The range of the flux contribution area increases with increase of the observation height under the same atmospheric stability conditions (Liu et al., 2014), and the flux source area has obvious seasonal variation (Gong et al., 2015). The source area distribution in the growing season has been shown to be greater than that in the nongrowing season in earlier studies of forest ecosystems (Wei et al., 2012). However, other studies have shown that the source area in the nongrowing season is greater than that in the growing season (Zhao et al., 2005), particularly when the atmosphere is unstable (Wu et al., 2012), and when the length of the flux source area under stable conditions is substantially longer than that under unstable conditions (Chen et al., 2016). Kljun et al. (2015) proposed a three-dimensional flux footprint model based on a scale analysis algorithm, which has fast operation speed and is suitable for calculation of long-term continuous time series of flux source area with satisfactory application effect (Ji et al., 2020). In recent years, among the various models available for calculating the contribution percentage of different types of underlying surface flux in a source area, the Agroscope Reckenholz–Tänikon (ART) footprint tool, developed by Neftel et al. (2008) and based on the KM model, has been widely used (Arriga et al., 2017; Wang and Xu, 2018).

The Xiaoxing’an Mountains are located in the temperate continental monsoon climate zone of Northeast China and cover an area of 12.06 million ha. The region includes more than 5 million ha of coniferous and broad-leaved mixed forest (CBF). The forest volume is approximately 450 million m³ and the forest coverage rate is 80.6%. The forest represents the largest and most complete virgin Korean pine forest (KPF) in Asia, and it plays an important role in China’s forest carbon sink. Thus, the region of the Xiaoxing’an Mountains has great potential for carbon sequestration.

Owing to the sparsity of the flux observation network in China, there is only one flux observation tower in the Xiaoxing’an Mountains. Currently, research on the CO₂ flux, spatial representation, and flux source area of this underlying surface in China is lacking. The findings of this study will fill this knowledge gap and provide a basis for assessment of the carbon source/sink potential of CBF in the Xiaoxing’an Mountains.

This study has selected eddy observation data from the Wuying National Climate Observatory acquired from January 2015 to November 2017 to explore the characteristics of the CO₂ flux of the forest ecosystem on different timescales. The characteristics of CO₂ spectra and cospectra, and distributions of the CO₂ flux source area using the Kljun model and ART footprint tool, are investigated. The study intends to improve understanding of the land–air CO₂ exchange activities in the Xiaoxing’an Mountains of Northeast China and provide references for further study of CO₂ flux in other ecosystems.

Wuying is located on the southern slopes of the Xiaoxing’an Mountains, and it has the largest and most intact virgin KPF belt in Asia. It is one of the most important areas of productive forest in China, with a total forest reserve of 11.23 million m³ and a forest coverage rate of 93.2%. In the region, annual average temperature is 0.6°C, annual precipitation is 610.7 mm, the annual number of sunshine hours is 2196.0 h, the length of the frost-free period is 117 days, and the accumulated temperature ≥ 10°C is 2141.8°C. The coniferous species are mainly Pinus koraiensis, and the broadleaved species are mainly Tilia amurensis, Quercus mongolica, Fraxinus mandshurica, Ulmus japonica, Ulmus laciniata, Populus ussuriensis, Acer mono Maxim, and Acer triflorum (Sun et al., 2021).

The observation site is located at the Wuying National Climate Observatory of the China Meteorological Administration. The average height of nearby trees is approximately 23–25 m, the terrain is slightly undulating (high to the north and low to the south), and the KPF is distributed mainly to the east, while the CBF is distributed mainly to the west (Fig. 1). A 70-m meteorological gradient observation tower has been built at the observation site, and an eddy measurement system is installed on the tower at the height of 50 m. This measurement system comprises a data collector (Li-7550, LI-COR Biosciences, USA), three-dimensional ultrasonic anemometer (WindMaster, Gill Instruments Ltd., USA), and a CO₂/H₂O analyzer (Li-7500, LI-COR Biosciences, USA). It is used for monitoring CO₂, H₂O, momentum, sensible heat, latent heat, three-dimensional wind speed, and other elements between the ground and atmosphere, and for the study of the carbon cycle of forest ecosystems.

Figure 1.  Location and surrounding land use types of the study site in the Xiaoxing’an Mountains (CBF: coniferous and broadleaved mixed forest, KPF: Korean pine forest, CR: construction land and roads).

This study used EddyPro 7.0 to process the raw data (frequency: 10 Hz) from the eddy measurement instrument into 30-min averages. Outliers with unreasonable values and/or obvious errors were eliminated. Then, linear detrending, coordinate rotation (secondary rotation), time lag correction, frequency response correction, ultrasonic virtual temperature correction, and density effect correction were performed. Finally, a series of turbulence statistical operations were conducted to determine the average, pulsation, variance, covariance, and flux (Wang et al., 2016). The power spectra of each component were obtained after the original data were processed by corresponding low-frequency filtering, pulsation calculation of each component, and fast Fourier transform. The high-frequency region of the obtained power spectra was processed by using low-pass filtering, moving average, and turbulence spectra normalization to obtain turbulence spectra and cospectra. The full calculation process is detailed in Sun et al. (2021).

In this study, the EddyPro software was used to screen the 30-min flux value output using the following procedure. (1) Eliminate data associated with instrumental error. (2) Exclude data acquired of 1 h before and 1 h after precipitation for CO₂ flux. (3) The threshold of CO₂ flux in this study was set as −3 to 3 mg m⁻² s⁻¹ (Li et al., 2005; Ji et al., 2019), and flux data exceeding this threshold were eliminated. (4) According to the stability test and development turbulence condition test, the quality of the flux data was evaluated, and the data were divided into three levels: 0, 1, and 2, corresponding to best, good, and poor, respectively. In this study, level 2 data were eliminated, which left an average of 68% (67% in 2015, 69% in 2016, and 69% in 2017) of the original data for further research following the quality control process.

Following quality control, the flux data were interpolated by using the average diurnal variation method (Falge et al., 2001; Baldocchi, 2003). The specific steps of the method are as follows. (1) Interpolate the average value for 14 days before and 14 days after the missing time of daytime flux. (2) Interpolate the average values of 7 days before and 7 days after the missing time of nighttime flux.

According to the actual climate conditions in the study area, the seasons are divided into spring (March, April, and May), summer (June, July, and August), autumn (September, October, and November), and winter (December, January, and February).

In this study, only the annual change of CO₂ flux was studied using the 3-yr data record; all other elements were studied using the 1-yr data record from December 2016 to November 2017.

The formula for calculation of CO₂ flux can be expressed as follows:

where F_c is the CO₂ flux (mg m⁻² s⁻¹); $ {w}' $ and $ {\rho }_{c}' $ are the fluctuation values of vertical wind speed and CO₂ density, respectively, and $ \overline{{w}'{\rho }_{c}'} $ represents the covariance of the two.

To accurately reflect the CO₂ flux characteristics of the forest underlying surface, it is necessary to consider the CO₂ flux stored under the canopy and introduce the net ecosystem exchange (NEE, i.e., carbon flux):

where NEE represents the change of ecosystem carbon storage caused by photosynthesis, biological and abiotic respiration, or consumption of CO₂ in the atmosphere by the plants in the ecosystem (a negative value of NEE indicates that the forest absorbs CO₂ in the atmosphere, while a positive value indicates that CO₂ is released into the atmosphere) (Fortuniak et al., 2017); F_s is the CO₂ storage flux, which indicates the accumulation and consumption of CO₂ by the underlying surface; z_h is the height of the eddy measuring instrument (50 m in this study); $ \Delta c/\Delta t $ is the rate of increase of CO₂ concentration with time; and z is the height above the ground surface (m).

The Kljun footprint model developed by Kljun et al. (2015) establishes a coordinate system with the observation point defined as the origin (0, 0), and the positive direction of the x axis represents the upwind distance, such that the following expression holds:

where F_c(0, 0, z_m) is flux density (per unit area); z_m is the effective dynamic height of the observation point [z_m = z_h − d₀, where z_h is the height of the eddy measuring instrument and d₀ is the zero plane displacement; here, d₀ adopts the calculated value (Table 1) in the literature (Sun et al., 2022)]; R denotes the integration domain, which is generally 10%–90%; Q_c is the source or sink integrated over the unit area; and f(x,y) is the footprint function, which is the conversion function of the carbon source or sink. It represents the contribution rate density of a point (x, y) on the surface to the observed value at z_m. The footprint function f(x,y) can be expressed as follows:

where $\bar{{f}^{y}}\left(x\right)$ is the crosswind-integrated footprint function and D_y is the crosswind dispersion function. Assuming that the crosswind dispersion function has Gaussian characteristics, f(x,y) can be written as follows:

where y is the crosswind distance from the centerline (i.e., the x axis), and σ_y is the standard deviation of the crosswind distance. The parameters related to ${\bar{f}^{y}}$ include the effective dynamic height z_m, atmospheric boundary layer height h, average wind speed at height z_m $\bar{u}\left({z}_{{\rm{m}}}\right)$, friction wind speed u_*, and roughness length z₀. Among the above parameters, h is related to L (the Monin–Obukhov length), and the detailed calculation method refers to the literature (Kljun et al., 2015); z₀ adopts the calculated value (Table 1) in the literature (Sun et al., 2022); while all other parameters are obtained through EddyPro. Through dimensional analysis (П Theorem), the above variables can be composed into the following dimensionless parameter group:

where k is the von Karman constant (k = 0.4), and Ψ_M is the stability correction function, which can be expressed as follows:

with $\chi ={\left(1-19{z}_{{\rm{m}}}/L\right)}^{1/4}$, where z_m/L > 0 is the atmospheric stability condition and z_m/L < 0 is the atmospheric instability condition.

The nondimensional form of the crosswind-integrated footprint, F^y*, can be written as a yet unknown function φ of the nondimensional upwind distance X*. Thus, F^y*= φ$ \left({X}^{*}\right) $, with ${X}^{*}={\varPi }_{2}{\varPi }_{3}{\varPi }_{4}^{-1}$ and ${F}^{*}={\varPi }_{1}{\varPi }_{3}^{-1}{\varPi }_{4}$, and the fitting function $ {\widehat{F}}^{y*}\left({\widehat{X}}^{*}\right) $ of $ {F}^{y*}\left({X}^{*}\right) $ can be obtained through parametric fitting:

where a, b, c, and d are the fitting parameters related to the roughness (z₀), with values of 1.4524, −1.9914, 1.4622, and 0.1359, respectively (Kljun et al., 2015), which are calculated using the backward Lagrangian stochastic particle dispersion model. This study used the online tool FFPonline provided by Kljun et al. (2015) (http://footprint.kljun.net/index.php) to calculate the flux footprint, and then research the flux footprint characteristics of the Xiaoxing’an Mountains forest ecosystem.

Routine observational data recorded at the Wuying National Meteorological Station were used for the analysis described in Section 3.1. Other parts of the analysis used the eddy observation data of the gradient observation tower at the Wuying National Climate Observatory, which were processed and output by EddyPro. The normalized and exponentially binned (co)-spectra of the 30-min output were used for the turbulence spectra. The data used for the flux footprint calculation included z_h, d₀, z₀, L, σ_y, u_*, and wind direction.

Figure 2 shows the changes of daily average temperature, sunshine, and precipitation in the study area in the Xiaoxing’an Mountains from January 2015 to November 2017. The trends of change of temperature, sunshine, and precipitation are broadly similar with obvious seasonal changes, which are more (higher) in the growing season (May–September) and less (lower) in the nongrowing season (October to April of the following year). The precipitation in the growing season accounts for 80% of the total annual precipitation, and the radiation in the growing season accounts for 47% of the total annual radiation.

Figure 2.  Annual changes of daily average temperature, sunshine hours, and precipitation at the study area in the Xiaoxing’an Mountains.

Over an entire year, the forest ecosystem of the Xiaoxing’an Mountains acts as a carbon sink (Fig. 3). The annual change of NEE (2015–2017) shows a “V” shape in the growing season (May–September) and a “U” shape in the nongrowing season (October to April of the following year). The maximum NEE in each year was −9.26 gC m⁻² day⁻¹ (28 June 2015), −9.11 gC m⁻² day⁻¹ (12 June 2016), and −8.71 gC m⁻² day⁻¹ (31 July 2017). The curve of cumulative NEE of the 3 years shows a downward trend with rapid decline in the growing season and a slow decline in the nongrowing season. From 2015 to 2017, the annual total NEE was −756.84, −834.73, and −629.37 gC m⁻², respectively (average: −740.31 gC m⁻²). Owing to lack of data from December 2017, the total annual NEE in 2017 was calculated from December 2016 to November 2017. However, the cumulative NEE from January to November 2015 (−712.75 gC m⁻²), January to November 2016 (−812.06 gC m⁻²), and January to November 2017 (−606.70 gC m⁻²) reveals that the order of carbon sequestration capacity is 2016 > 2015 > 2017.

Figure 3.  Annual changes of net ecosystem exchange (NEE) capacity at the study area in the Xiaoxing’an Mountains.

The difference in carbon sequestration capacity between the 3 years might be related to variations in temperature, sunshine, and precipitation. The average temperature, number of sunshine hours, and amount of precipitation during May–August were 17.15°C, 895 h, and 444.4 mm, respectively, in 2015, 16.99°C, 805.5 h, and 557.4 mm, respectively, in 2016, and 17.48°C, 926.9 h, and 396.8 mm, respectively, in 2017. In the May–August period, the amount of precipitation was greatest, average temperature was lowest, and number of sunshine hours was fewest in 2016, while the amount of precipitation was smallest, average temperature was highest, and number of sunshine hours was greatest in 2017. Corresponding to the carbon sequestration capacity, it can be found that enhanced precipitation is conducive to photosynthesis and absorption of CO₂ by plants. When the air temperature is too high, to prevent excessive transpiration of water, the stomata of plant leaves will be partially closed and photosynthesis will be weakened, which is not conducive to the absorption of CO₂ and results in low NEE value.

The daily change of NEE in the study area in the Xiaoxing’an Mountains is shown in Fig. 4. It can be seen that the greatest concentration of blue and red areas is in the growing season, indicating negative values during daytime [0600–1800 LST (local standard time)] and positive values at night (1800–0600 LST the following day). In the nongrowing season, dark blue areas appear at noon in winter, while light blue areas are evident at all other times. It shows that CO₂ exchange between the forest ecosystem and the atmosphere is strongest in the wet season, whereby the system is characterized as a strong carbon sink in the daytime and a strong carbon source at night. In winter, the system generally acts as a weak carbon sink in the daytime and a weak carbon source at night.

Figure 4.  Daily change of NEE at the study area in the Xiaoxing’an Mountains.

Owing to the different sunrise/sunset times, substantial changes in solar radiation and the growth period of trees, and obvious changes in turbulence intensity in different seasons, the daily change of NEE varies seasonally (Fig. 5). After sunrise, with the increase of solar radiation, trees absorb CO₂ through the process of photosynthesis, whereby CO₂ is transmitted downward from the atmosphere above the canopy through turbulent exchange, NEE rapidly decreases below zero, and the net CO₂ absorption rate of the forest underlying surface increases. After sunset, owing to the respiration of trees, the forest emits CO₂ to the atmosphere, and NEE approaches zero and gradually becomes positive.

Figure 5.  Average daily change of NEE at the study area in the Xiaoxing’an Mountains.

In summer, convective activity is strong during the day, the solar altitude angle is large, radiation is intense, and the photosynthesis process is most active, resulting in the lowest value of NEE (−0.235 mgC m⁻² s⁻¹ at 0830 LST). Conversely, the respiration of trees is strongest at night, resulting in the highest value of NEE (0.076 mgC m⁻² s⁻¹ at 0430 LST). The diurnal changes of NEE in spring, autumn, and winter are broadly similar. The lowest and highest values of NEE in spring appear at 0900 and 0430 LST, respectively, ranging from −0.083 to 0.022 mgC m⁻² s⁻¹; in autumn, the lowest (−0.074 mgC m⁻² s⁻¹) and highest (0.034 mgC m⁻² s⁻¹) values appear at 1000 and 0100 LST, respectively; and in winter, the lowest (−0.074 mgC m⁻² s⁻¹) and highest (0.017 mgC m⁻² s⁻¹) values appeared at 1100 and 0800 LST, respectively. The appearance of the peaks in the morning rather than in the noon might reflect the inhibition of leaf photosynthesis by high vapor pressure deficit (Huang et al., 2011; Jia et al., 2015).

The monthly change of NEE increased gradually during June–October (Fig. 6), and reached its highest value (−0.0025 mgC m⁻² s⁻¹) in October, indicating that the carbon sequestration capacity of the forest ecosystem of the Xiaoxing’an Mountains was weakest in autumn. The change of NEE was minimal during October–April, but it gradually decreased from April to June. In June, plants grow vigorously in response to the increased solar radiation energy, and the strengthened photosynthesis leads to increased absorption of atmospheric CO₂ and the greatest carbon fixation capacity (−0.0612 mgC m⁻² s⁻¹). The average value of NEE in October was close to zero, which might reflect the decrease in solar radiation in autumn, weakening of the photosynthesis process of trees during daytime, and strong respiration of trees and soil microorganisms at night.

Figure 6.  Average monthly change of NEE at the study area in the Xiaoxing’an Mountains (IQR: interquartile range).

Eddies with different sizes superimposed on each other form turbulence, and the relative strength of eddies with different scales is defined as turbulence spectra. The turbulence spectra show the contributions of various types of turbulent vortices to the turbulent energy, and the covariance spectra between atmospheric parameters are often used to analyze the contributions of turbulent vortices of different scales to the corresponding turbulent flux (Kaimal and Finnigan, 1994). According to Kolmogorov’s theory, in the inertial subrange of the near-surface layer, turbulence spectra follow the −2/3 law and the covariance follows the −4/3 law (Kaimal et al., 1972). The spectral peak frequency is the frequency corresponding to the maximum of the turbulent spectral energy. The spectral peak wavelength corresponds to the scale of the eddies that contribute most to the turbulent energy; the larger the spectral peak wavelength, the larger the corresponding scale of the eddies (Sun et al., 2021).

We analyzed the CO₂ spectral characteristics of the forest canopy (Fig. 7). It was found that the normalized spectra of CO₂ follow the −2/3 law at the higher frequencies, but that the spectral line dispersion relative velocity spectra under different stability conditions is large (Sun et al., 2021), indicating that atmospheric stability has great impact on CO₂ concentration in the study area in the Xiaoxing’an Mountains. The spectral peak wave-length under different atmospheric stability conditions is approximately 172–398 m. The normalized cospectra of the vertical wind component (w) and CO₂ (cw) do not completely follow the −4/3 law in the higher frequencies, and the slope of the fitting line is closer to −1 (Fig. 8). The peak wavelength of the cw cospectra is approximately 172–301 m under different atmospheric stability conditions. The CO₂ spectra and cw cospectra curves gra-dually move toward the higher frequencies with increasing atmospheric stability, and the frequency of the spectral peaks also increases accordingly, indicating that the scale of the eddies that contribute most to the turbulent energy decreases with increase of atmospheric stability.

Figure 7.  Normalized spectra of CO₂ (c) change with normalized frequency under different atmospheric stability conditions (x axis nz/u is normalized frequency).

Figure 8.  Normalized cospectra of cw change with normalized frequency under different atmospheric stability conditions (x axis nz/u is normalized frequency).

Figure 9 shows the distributions of wind direction and wind speed during daytime and nighttime in 2017 near the observation station. During the observation period, the study area experienced mostly southwesterly winds throughout. Wind speed during daytime was generally 2–6 m s⁻¹ at most, and winds with speed > 8 m s⁻¹ occurred rarely; at night, wind speed was generally 0–4 m s⁻¹, indicating that the observation period was mainly dominated by weak winds. The distributions of wind direction and wind speed in each season are similar to those of the entire year (Table 2). In spring, the southwesterly wind was dominant (55.4%), while a northeasterly wind occurred least (12.34%). In summer, the southwesterly wind was dominant (43.72%), while the northwesterly wind occurred least (8.3%). In autumn, the southwesterly wind was dominant (60.65%), while the northwesterly wind occurred least (7.24%). In winter, the southwesterly wind was dominant (76.2%), while the northeasterly wind occurred least (6.95%).

Figure 9.  Distributions of annual wind direction and wind speed (WS; m s⁻¹) at the study area in the Xiaoxing’an Mountains: (a) daytime and (b) nighttime.

In this study, we used the Kljun model to calculate the flux footprint of the 10%–90% contribution rate at 30-min intervals. Through analysis of the annual source area distribution (Fig. 10), it was found that the extent of the source area when the atmosphere had unstable stratification was substantially smaller than that when the atmosphere had stable stratification.

Figure 10.  Annual CO₂ flux footprint at the study area in the Xiaoxing’an Mountains in 2017: (a) unstable atmospheric conditions and (b) stable atmospheric conditions.

The size of the source area with a flux contribution rate of 90% was 2.68 km², and the length of the source area was 565.29–1339.17 m in the case of unstable stratification. It is because turbulent movement in the vertical direction is intense, vertical transmission of material is rapid, and the flux measured by the observation instrument is derived from a close distance in the upwind direction, resulting in the small source area. When the atmosphere has stable stratification, the source area with a flux contribution rate of 90% was 20.5 km², and the length of the source area was 1432.24–4317.06 m. In this case, turbulent movement of the air is weak, and the vertical diffusion speed of material is slow; therefore, the observation instrument can be affected by eddy currents from more distant places, resulting in the large source area (Ji et al., 2020). Moreover, it can be found from Fig. 10 that the main flux was derived from the southwest direction under either stable or unstable conditions. It was found that the main factors affecting the flux contribution area are the wind direction and the wind speed (Fig. 9), consistent with the research of Ji et al. (2020).

Analysis of the distribution of flux source areas with a contribution rate of 10%–90% in different seasons revealed large seasonal differences in flux contribution area (Fig. 11). The maximum length of the source area in the different seasons was toward the southwest (180°–270°), consistent with the distribution of the source area throughout the year, and the shape of the source area was consistent with the wind direction characteristics in the different seasons (Table 1). The length of the source area did not exceed 2000 m under atmospheric unstable conditions, and the order of the range of the source area was as follows: winter > spring > autumn > summer. The length of the source area did not exceed 6000 m under stable conditions, and the order of the range of the source area was as follows: winter > spring > autumn > summer, consistent with that under unstable conditions. It indicates that the distribution of the source area in the study area under different atmospheric conditions is greater in the nongrowing season than in the growing season, and the underlying reason is that turbulence activity is strong in the growing season and weak in the nongrowing season.

Figure 11.  CO₂ flux footprint during different seasons at the study area in the Xiaoxing’an Mountains: (a), (c), (e), and (g) unstable conditions in spring, summer, autumn, and winter, respectively, and (b), (d), (f), and (h) stable conditions in spring, summer, autumn, and winter, respectively.

As can be seen from the distribution of the crosswind integral function (Fig. 12), the peak value of the flux contribution rate in different seasons was located near the sensor. The peak distance was 30.79 m under atmospheric unstable conditions and 153.95–184.74 m under stable conditions. The crosswind footprint function value decreased markedly with the change from unstable conditions to stable conditions. Under unstable conditions, strong turbulence activity causes faster material transport in the vertical direction and reduces the extent of transport in the horizontal direction (Ji et al., 2020). The distance of lateral transport by wind under stable conditions is markedly greater than that under unstable conditions. Under stable conditions, the scale of the eddies is small, and the sensor can detect eddies derived from distances further from the tower, which is consistent with the finding in relation to the CO₂ spectral characteristics that the scale of the eddies that contribute most to the turbulent energy decreases with increase of atmospheric stability.

Figure 12.  Diagrams of crosswind integral function at the study area in the Xiaoxing’an Mountains: (a) spring, (b) summer, (c) autumn, and (d) winter.

The area to the east of the gradient tower at the study area in the Xiaoxing’an Mountains is dominated by KPF and the area to the west is dominated by CBF. Different types of trees make different contributions to the total flux; therefore, it is necessary to quantify the flux contribution rate of different tree species. We selected a 2000-m × 2000-m study area (Fig. 1), divided the area into 24 quadrilateral segments according to different land types, calculated the contribution rate of the different segments to the total flux using the ART footprint tool (Kljun et al., 2002, 2003; Neftel et al., 2008), and overlaid vegetation types to obtain the flux contribution rate of the different vegetation types throughout an entire year (Fig. 13). The calculation results show that the annual total contribution rate of the selected areas is 97%, and that the order of the contributions of the source areas from large to small is as follows: CBF > KPF > CR > water. Because the west side of the observation station is predominantly CBF and the east side is mostly KPF (Fig. 1), and because the flux contribution area mentioned above is mainly located to the southwest of the tower, the proportion of the contribution of CBF is relatively large, i.e., up to 3.4 times that of KPF. The construction land in the study area, which comprises mostly buildings, roads, and footpaths associated with scenic spots in the Wuying National Forest Park, makes a contribution of only 3.6% to the CO₂ flux; thus, it can be neglected.

Figure 13.  Proportion of source area contribution to total flux by land use type.

The Xiaoxing’an Mountains have the largest and most complete virgin KPF in Asia, which offers great potential for carbon sequestration. To understand the carbon source/sink characteristics of the forest ecosystem of the Xiaoxing’an Mountains, this study used eddy observation data from the Wuying National Climate Observatory acquired from January 2015 to November 2017 to analyze the CO₂ flux and source area distribution characteristics. It was found that the annual change of NEE in the forest ecosystem of the Xiaoxing’an Mountains presents a “V” shape in the growing season (May–September) and a “U” shape in the nongrowing season (October to April of the following year). The diurnal change of NEE exhibits single-peak characteristics with decline during daytime and increase at night. In the humid and warm seasons (May–September), CO₂ exchange between the ground and the air is strongest, and the NEE in each month is less than zero. The carbon sequestration capacity is strongest in June and weakest in October, which is different from the finding reported in other literature indicating that winter is the season with weakest carbon sequestration (Huang et al., 2011; Jia et al., 2015). This discrepancy might reflect the effects of decrease in solar radiation in autumn, weakening of the photosynthesis process of trees during daytime, and the strong respiration of trees at night.

The annual average NEE is −740.31 gC m⁻² yr⁻¹, lower than that reported for the mangrove forest of the Mai Po Nature Reserve (China) but substantially higher than that of other forest ecosystems (Table 3). It is also different from the finding of Huang et al. (2019) who reported that NEE decreases with increasing latitude. This indicates that the forest ecosystem of the Xiaoxing’an Mountains has considerable carbon sequestration capacity. The NEE differences between different forest ecosystems might be related to subjective factors such as different data quality control methods and data interpolation methods, while the objective factors may be related to differences in regional climatic characteristics, landforms, and type, density, and growth stage of vegetation.

By analyzing the CO₂ spectra and cw cospectra, it was found that the CO₂ spectra meet the −2/3 power rate in the higher frequencies, while the cw cospectra do not completely accord with the −4/3 power rate, and the fitting slope is closer to −1. The peak wavelengths of the CO₂ spectra and cw cospectra are 172–398 and 172–301 m, respectively, which are similar to those of temperature spectra and humidity spectra (Sun et al., 2021).

Research on the distributions of wind direction, wind speed, and flux source area near the observation tower revealed that southwesterly winds prevail throughout most of the year. The main flux measurements are derived from the southwest, and the shape of the source area is consistent with the characteristics of wind direction and wind speed irrespective of the condition of atmospheric stability. It is clear that wind direction and wind speed are the main factors affecting the flux contribution area, which is consistent with the findings of Ji et al. (2020). The source area in the nongrowing season is larger than that in the growing season, and the order of the source area is as follows: winter > spring > autumn > summer, which is consistent with the research of Zhao et al. (2005). The source area at the time of stable atmospheric stratification is substantially larger than that under unstable atmospheric conditions. Under unstable atmospheric stratification, turbulence activity is intense, material exchange in the vertical direction is rapid, and diffusion is fast, resulting in a small range of the flux source area. Under stable atmospheric conditions, the observation instrument can be affected by eddy currents from more distant areas and thus the range of the source area is correspondingly larger, consistent with previous conclusions (Wu et al., 2012; Chen et al., 2016; Ji et al., 2020).

The maximum distance of the 90% flux source area is not more than 6000 m under stable atmospheric conditions and 2000 m under unstable atmospheric conditions. In previous studies, the maximum distance of the source area obtained by Ji et al. (2020) was 7000 m, that found by Wei et al. (2012) was 4500 m, and that reported by Gong et al. (2015) was not more than 2500 m. The furthest extent of the source area shows notable differences between different forest ecosystems, which might be related to the height of the observation instrument and the canopy height. In this study, the observation height was 50 m and the canopy height was 23 m, the observation height and canopy height of Ji et al. (2020) were 40 and 15 m, respectively, those of Wei et al. (2012) were 25 and 16 m, respectively, and those in Gong et al. (2015) were 38 and 11 m, respectively. Some studies have shown that the higher the measurement height, the larger the range of the high-throughput contribution area and the further the maximum distance extends from the source area (Liu et al., 2014). However, the maximum detection distance at the height of 50 m in this study was smaller than the maximum distance of 40 m in Ji et al. (2020), which might be related to the difference in height between the measuring instrument and the canopy. In this study, the height difference between the two was 27 m, i.e., the observation height was approximately twice the canopy height, while the height difference in Ji et al. (2020) was 25 m, i.e., the observation height was approximately 2.5 times the canopy height. Although the actual height difference is similar in both cases, the ratio between the two is different, which might account for the difference in the maximum distance of the source area. Although the underlying surface in the other two studies was a forest ecosystem, the different topography and forest types might also have had an impact on the flux source area, which is a subject that should be discussed in depth in future research.

It was found that the contribution of the source area around the observation station to the total flux varied from large to small in the following order: CBF > KPF > CR > water. The contribution rate of CBF is 3.4 times that of KPF. This is because the flux contribution area is mainly located to the southwest of the observation tower, and the contribution area to the east direction is very small; consequently, the source area of CBF in the west is much larger than that of KPF in the east. Additionally, although the contribution rate of CR is small (only 3.6%), which could be neglected in theory, the CO₂ emission generated by human activities will cause underestimation in the assessment of the carbon sink capacity of the forest ecosystem of the Xiaoxing’an Mountains.

In this study, we examined the distribution characteristics of CO₂ flux and source area in the Xiaoxing’an Mountains. Owing to improper instrument maintenance, a large number of data were missing, which meant that only data acquired during 2015–2017 were relatively complete for research purposes. To improve understanding of the overall variation law of matter and energy exchange between the ground and the atmosphere in the forest ecosystem of the Xiaoxing’an Mountains, it will be necessary to study multiyear observational data. Certain differences were found between the area of focus of this study and other forest ecosystems, and although plausible explanations were proposed, their verification was prevented owing to lack of detailed data. Consequently, further research will be needed in the future. Additionally, owing to lack of detailed information on land use types, the forest types around the observation station were crudely classified as either CBF or KPF. When data availability allows, further investigation will be required to determine whether there might be some deviation in the contribution rate of the source area.

The authors thank James Buxton, MSc, for editing the English text of a draft of this manuscript.