Introduction

Eddy covariance (EC) technology is a kind of micro-meteorological observation technology for measuring, long-term, fluxes of material and energy between the ground’s surface and the atmosphere. Spatial and temporal evaluations of distributions across a contribution area and the reliability of the magnitude of the gathered data depend primarily on footprint function or source proportion (Nicolini et al. 2017). In particular, material flux variations are affected by non-uniform features of the underlying surface, such as land use, vegetation and biomass, soil type, and leaf surface index. The footprint function model can be used to analyze the effects of surface structure conditions such as roughness, measured height, and micro-meteorological conditions (wind speed, turbulence intensity, and atmospheric stability) on the effective source area of the material flux (Arriga et al. 2017).

With the development of EC, many achievements have been reached in the study of flux contribution areas, resulting in the proposals of various models. Commonly used models include analytical models, Lagrangian random diffusion models, large eddy simulation models, and closed models (Pandey et al. 2011). Comparative analyses of the advantages and disadvantages of each model have been performed (e.g., Guo et al. 2018; Wang et al. 2014a; Arriga et al. 2017). Thus, use of the footprint theory to interpret and analyze measurements of carbon exchanges between the atmosphere and ecosystems in complex terrain areas has progressed (Neftel et al. 2008; Sogachev and Dellwik 2017; Fry et al. 2018). In 2001, Kormann and Meixner proposed an analytical footprint model (the KM model) based on the advection–diffusion model, which is also used for basic research in the dynamic characteristics of CO2 fluxes over days, months, seasons, or years and to analyze the effects of land use changes on flux traces. In these studies, positive values indicate that CO2 is released to the atmosphere, and negative values indicate an atmospheric CO2 sink, the atmospheric CO2 drawdown during hydrochemical erosion of soluble substances or photosynthesis of plants. The Agroscope Reckenholz Tanikon (ART) footprint tool (2007-03-13; Version 1.0) software developed by the Swiss Federal Agricultural Research Center (Neftel et al. 2008) is widely used where the terrain is complex (Arriga et al. 2017; Guo et al. 2018; Zhang et al. 2010; Gong et al. 2017; Meyer et al. 2015). Its primary function is to analyze the effects and contributions of various underlays on the carbon flux, and its purpose is to provide a reliable basis for examining the utilization and governance of carbon sources and sinks.

Current studies and assessments of carbon fluxes are being conducted in human-affected areas such as farmland (Wang et al. 2014a; Wu et al. 2017), cities, and surrounding areas (Guo et al. 2018). Various wetlands (Zhang et al. 2008; Yang et al. 2017), grasslands (Arriga et al. 2017), forests (Dou et al. 2018; Ma et al. 2015), oceans (Tortell et al. 2012; Ito et al. 2017; Wrobel 2017), deserts (Meyer et al. 2015), and snowpacks (Berryman et al. 2018; Zhao et al. 2015), where human influence is relatively small, have also been studied. Achievements are numerous.

In recent years, EC has been applied to glaciers, which cover about 11% of the global land area. Most of these studies focus on energy balance research (Guo et al. 2011; Litt et al. 2015; Lund et al. 2017; Yao et al. 2014), examining variations in carbon fluxes in the glacier area using hydrochemistry analysis (Krawczyk and Bartoszewski 2008; Donnini et al. 2016; Galeczka et al. 2016; Feng et al. 2012; Wang et al. 2016). The complex terrain and changes in the underlying surface affect instrumentation stability, thus making long-term monitoring of CO2 flux difficult in glacier areas.

Mountain glaciers, important components of the cryosphere, are not only recorders and early warning devices of climate change (Wang et al. 2008), but also important drivers of climate change. As a whole, global glaciers are currently in a state of retreat. According to Farinotti et al. (2015), 97.52% of the glaciers of the Tianshan Mountains in China are retreating. On a background of glacier ablation, the key factor in exploring hydrochemical erosion in cold regions is the change in atmospheric CO2 flux at the gas/liquid interface. This not only has significance for the study of hydrochemical erosion but also provides references for global carbon cycle research. Therefore, EC was used to measure and monitor the atmospheric CO2 flux in an ablation area of Koxkar Glacier in western Mt. Tianshan, China. The distribution characteristics of the CO2 flux contribution area were analyzed using the ART footprint tool during various typical periods to evaluate the reliability of EC in monitoring CO2 flux when the underlying surface moraine areas are complex. The footprint tool was also used to provide reference data for estimating and analyzing CO2 carbon sources and sinks in a mountain glacier area.

Methods

Site description

Koxkar Glacier is located on the southern slope of the Tomur peak in the Western Tianshan Mountains of China. The upper limit is Koxkar peak, elevation 6342 m. The effective positive difference from glaciation is more than 3000 m. The material balance line is at an altitude of 4550–4700 m (Zhang et al. 2006), and the area is 82.89 km2 (the area of ablation is about 30.6 km2, and its length is 19.0 km), 70.2% of the runoff control area above the hydrologic station (118.12 km2). The ice volume is 15.79 km3 (Han et al. 2010). On the ice surface is a wide range of foul layers and moraine, the latter accounting for about 83% of the total ablation area. However, the distribution of moraine thickness is quite different; it is discontinuous above 3700 m but exceeds 2 m at the glacier snout.

Paleogene mudstones, siltstones, and gravel were in the ridges and snout of the glaciers, and quaternary moraine was on its surface. A large number of marine terrigenous clastic rocks and carbonatites were distributed below 3400 m (Wang et al. 2010, 2016), providing sufficient source material for glacial hydrochemical erosion. Multi-year monitoring at the end of Koxkar Glacier indicates that the average annual runoff is about 0.9–1.2 × 108 m3 a−1, the majority (~ 94.5% of the annual total) occurring in the warm season (May–October)(Han et al. 2010; Wang et al. 2010).

Instrumentation

From July 22, 2015 to July 18, 2017, a relatively flat and open ablation area (N 41°43′, E 080°09′, 3212 m) covered by moraine was selected for CO2 flux monitoring (Fig. 1) for two reasons: (1) 83% of the ablation area was covered by glacial debris and (2) the moraine, 0.72 m in thickness, could effectively mitigate glacial ice melting and provide a stable observation site. Selecting a strong ablation area for monitoring on bare ice or where the moraine area is less predominant on the glacier surface could yield better representation of the data, but the dynamic surface underneath could affect instrument stability or cause damage to the instrument if it were to fall. An open-path EC system, which uses a three-dimensional ultrasonic wind and temperature meter [CSAT3; Campbell Scientific Inc. (CSI), USA], was used to measure the three-dimensional wind speed, and an open-path CO2/H2O infrared gas analyzer (Li-7500; LiCor Inc, USA) was used to measure CO2/H2O concentration pulsation. Additionally, an HMP45C temperature-humidity probe was used to measure temperature and relative humidity (Sogachev and Dellwik 2017). Utilizing a data acquisition device (CR1000 with an additional 16 Gb compact flash memory card; CSI), the 10 Hz original data were automatically collected and stored. After online virtual temperature correction and air density fluctuation corrections, averages were calculated and stored every 30 min.

Fig. 1
figure 1

Observation location and surrounding terrain (grid 50 × 50 m) in the debris area of Koxkar Glacier where CO2 flux was monitored

Full size image

To eliminate or reduce the influence of topography on the EC system, sensor height was reduced (i.e., a smaller effective observation area was used); therefore, the sensor for the infrared gas analyzer was mounted at a height of 2.0 m. Simultaneously, precipitation observations were made, environmental gradient observation systems were utilized to monitor wind speed and direction, temperature, and humidity at heights of 2 m, 4 m, and 10 m, and a CNR1 four-component net radiometer was used at a height of 3.5 m in the camp area.

Data quality and processing

According to the theory of EC, CO2 flux after correcting for sensible and latent heat fluxes is:

$${F_{{\text{c-WPL}}}}=\overline {{w\prime \rho \prime }} +\mu \frac{{\overline {{{\rho _{\text{c}}}}} }}{{\overline {{{\rho _{\text{d}}}}} }}\overline {{w\prime \rho \prime }} +\left( {1+\mu \frac{{\overline {{{\rho _{\text{v}}}}} }}{{\overline {{{\rho _{\text{d}}}}} }}} \right)\frac{{\overline {{{\rho _{\text{c}}}}} }}{{\overline {T} }}\overline {{w\prime {T_{\text{a}}}}} ,$$
(1)

where \({F_{{\text{c-WPL}}}}\) is CO2 flux modified by the WPL correction (mmol m−2 s−1), \(T~\)and \({T_{\text{a}}}\) are the virtual and air temperatures (°C), respectively, \(w\) is the vertical wind speed (m s−1), \(\mu\) is the ratio between the molecular masses of dry air and water vapor, and \({\rho _{\text{c}}}\), \({\rho _{\text{d}}}\), and \({\rho _{\text{v}}}\) are the densities of CO2 gas, dry air, and water vapor (mmol m−3), respectively (Wang et al. 2007).

When the atmosphere is relatively stable or the turbulent mixing effect is weak, CO2 dissolved or released by the surface water cannot reach sensor height; it accumulates below the sensor. The storage term, Fs, is:

$${F_{\text{s}}}=\mathop \smallint \limits_{0}^{z} \frac{{\overline {{\partial c}} }}{{\partial t}}{\text{d}}z,$$
(2)

where \(z\) represents the measurement and integration height, and \(\frac{{\overline {{\partial c}} }}{{\partial t}}\) represents the time rate of change of the CO2 mixing ratio.

Mountain glaciers are distributed primarily in high-altitude mountainous regions because conditions for glacier development are limited. Based on the material conservation equation, vertical and horizontal advection items customarily account for less than 5% of the net exchange (Nicolini et al. 2018; Wohlfahrt and Galvagno 2017; Wharton et al. 2017; Wang et al. 2014a). For the purposes of this study, and because these observations are relatively complicated, vertical and horizontal advection items are ignored, yielding a simplified equation for net glacial exchange (NGE) of CO2 in Koxkar Glacier:

$${\text{NGE}}={F_{{\text{c-WPL}}}}+{F_{\text{s}}},$$
(3)

All flux data and micro-meteorological data were subject to data quality control. In the former, pre-processing before analysis, including outlier removal, three-dimensional coordinate rotation, and a WPL correction, were performed (Wharton et al. 2017).

Field flux data are easily influenced by weather conditions such as precipitation, dew, and fog, as well as the instrument itself; irrational data must be removed (Wharton et al. 2017; Gu et al. 2003). After data quality control, 92.3%, 86.8%, and 81.2% of the total data remained for periods of snow accumulation (October to February), snow melting (March to May), and glacial melting (June to September), respectively. When missing data spanned an interval of 2 h or less, an internationally-accepted linear interpolation method was applied (Zhu et al. 2006; Anita et al. 2010). When those intervals exceeded 2 h, a linear relationship between carbon flux (\({F_{{\text{c-WPL}}}}\)) and the environmental factor net radiation (\({\text{NR}}\)) was applied for interpolation, Table 1. The equations passed the test of significance, where p < 0.01.

Table 1 The relationship between carbon flux (\({F_{{\text{c-WPL}}}}\), µmol m−2 s−1) and net radiation (\({\text{NR}}\),W m−2) in various typical phases
Full size table

ATR-footprint model

To examine the reliability of EC when investigating CO2 fluxes in a debris-covered area with complex topography, the ART footprint tool was used based on the KM analytical model (Kormann and Meixner 2001) to analyze the contribution area during three typical periods. The ART footprint tool was used to determine the two-dimensional footprint density function \(\varphi (x,y)\), calculating the five parameters (A, B, C, D, and E) that define the footprint function. With these five parameters, the footprint function was defined as:

$$\varphi (x,y)={D_{xy}} \cdot C \cdot \exp \left( { - \frac{B}{x}} \right) \cdot {x^{ - A}},$$
(7)

where \({D_{xy}}\) describes the Gaussian crosswind distribution:

$${D_{xy}}=\frac{1}{{\sqrt {2\pi } \cdot D \cdot {x^E}}} \cdot {\text{exp}}\left( { - \frac{{{y^2}}}{{2 \cdot {{(D \cdot {x^E})}^2}}}} \right).$$
(8)

Parameters A, B, and C thus describe the crosswind integrated distribution, whereas parameters D and E specify the crosswind distribution. The specific calculation methods for parameters A through D are described by Neftel et al. (2008). The footprint contributions of the user-defined fields were computed. For each grid point lying inside any quadrangle (50 × 50 m; see Fig. 1), its contribution to the footprint is calculated, and the sum yields the footprint contribution of the field. The footprint contributions of the quadrangles are given as percentages of the integral over the considered domain. During data processing, 100 cells of 100 × 100 m, 70 × 70 m, 50 × 50 m, and 25 × 25 m were compared. At these scales, 91.06%, 90.33%, 90.17%, and 87.54%, respectively, of the CO2 fluxes of the glacial ablation season were explained using the ART footprint tool. Ice cliffs are common glacier landforms in glacier areas that are covered by moraine. They are distributed in the NNW direction about 130 m from the EC system, necessitating the selection of 50 × 50 m quadrangles. Output parameters determined by the ART footprint tool included zero plane displacement, footprint coverage, and regional footprint contribution rate (Neftel et al. 2008; Guo et al. 2018).

Results

Weather and hydrology

Weather

During the observation period, the temperature of the meteorological field near the hydrological section ranged from − 21.2 to 21.7 °C, averaging 0.7 °C. During the cold seasons (October–February) and controlled by high pressures from Mongolia (Yao et al. 2014), it was cold and dry, averaging − 6.7 °C and 54.7% relative humidity; precipitation was supplied primarily by moist air from the Atlantic and Arctic Oceans (Shen et al. 2003), accounting for only 6.07% of the total annual precipitation. Development of regional low-pressure air masses brought the ablation period of snow cover and averages of 1.1 °C and 59.3% relative humidity. Precipitation occurred primarily during the melting period, accounting for 70.2% of the annual rainfall; the average temperature was 9.7 °C, peaking at 21.7 °C, strongly accelerating melting. The relative humidity was 61.0%.

With the effects of the alpine topography, surface wind directions were predominantly NW and NNW (Fig. 2), consistent with the findings of Wang et al. (2014b) and, at 3700 m above sea level, of Han et al. (2010). On the other hand, significant differences were found in the frequency of wind and wind speed during periods of snow accumulation, snow ablation, and glacial melting (Fig. 2). While snow accumulated, winds were affected by high pressures in Mongolia and Siberia, and they were predominantly NW and NNW (Han et al. 2010), accounting for 45.04% of the measurements. Average wind speed was 2.56 ms−1; wind speeds did not exceed 4 m s−1 87.38% of the time and exceeded 8 m s−1 0.73% of the time. Winds that were NNW remained mostly constant, but NW winds were significantly less during ablation compared to snow accumulation periods, accounting for 11.43% and 9.86% of the winds during snow-melting and glacial ice-melting periods, respectively. Conversely, SSW winds were more significantly more frequent. Remarkably, average wind speeds were significantly greater, 2.94 and 3.21 m s−1, during snow and glacial ice melting periods, respectively. Strong winds, those exceeding 8 m s−1, accounted for 2.48% of the winds during glacial ice melting periods.

Fig. 2
figure 2

Windrose figures of Koxkar Glacier on the southern slope of the Tianshan Mountains. a Snow accumulation season, b snow melting season, and c glacial ice ablation season

Full size image

Hydrology

Water levels were monitored using a Campbell SR50A sensor and corrected by applying manual readings. Total runoffs were found to be 7.51 × 107 m3 and 9.08 × 107 m3 from July 21 to November 10, 2015 and from March 31 to September 30, 2016. Moreover, from June to August 2016, total runoff was 6.41 × 107 m3, 18.34% less than that found by Han et al. (2010), and the temperature was 0.21° lower than that during the summer of 2008. The decrease in runoff conformed to the degree-day model (Li et al. 2012). Unfortunately, a large error in runoff data was caused by the frozen river surface from November to the following March. Discontinuous data was a result of the damaged hydrologic section after October 2016.

Fluxes in CO2 (F c_wpl)

At Koxkar Glacier, CO2 flux ranged from − 408.95 to 81.58 mmol m−2 day−1, averaging − 58.68 mmol m−2 day−1. This value is not only greater than that found for the sea ice area in the Amundsen Sea in Antarctica (Tortell et al. 2012), but also much greater than the annual CO2 flux budget, which averages between − 1.0 and − 2.4 mmol m−2 day−1 at the air–sea interface (Ito et al. 2017). It is similar to the annual CO2 sequestration in wet areas (− 560 to − 980 g m−2 year−1, equivalent to − 34.87 to − 61.02 mmol m−2 day−1; Table 2) (Fortuniak et al. 2017). Thus, the Koxkar Glacier ablation area was mostly a sink for CO2. Average CO2 fluxes were − 87.05 mmol m−2 day−1 and − 38.41 mmol m−2 day−1 in the warm (May–September) and cold seasons, respectively. A similar trend was found in the seasonal variation of CO2 flux on the underlying surface with vegetation, particularly in farmland (Fleischer et al. 2016).

Table 2 Variations in CO2 fluxes during various typical periods, in mmol·m−2 day−1
Full size table

Also, although the average temperature was − 6.5 °C during snow accumulation, and snow and ice melting nearly ceased, atmospheric CO2 flux was − 28.15 mmol m−2 day−1, and − 67.19 mmol·m−2 day−1 during the day, indicating that the study area remained a carbon sink. This was because sedimentation of suspended matter in the atmosphere created heavy pollution in some of the snow layer. During daytime, strong solar radiation caused a small amount of melting in the surface snow, and hydrochemical reactions took up atmospheric CO2 by elution of soluble matter (Niu et al. 2017; Mouri 2016) (e.g., Eqs. 9, 10).

$${\text{CaC}}{{\text{O}}_3}+{\text{C}}{{\text{O}}_2}+{{\text{H}}_2}{\text{O}} \to {\text{C}}{{\text{a}}^{2+}}+2{\text{HCO}}_{3}^{ - },$$
(9)
$${\text{CaA}}{{\text{l}}_2}{\text{S}}{{\text{i}}_2}{{\text{O}}_{8{\text{(s)Anorthite}}}}+2{\text{C}}{{\text{O}}_{2{\text{(aq)}}}}+2{{\text{H}}_2}{\text{O}} \to {\text{C}}{{\text{a}}^+}+2{\text{HCO}}_{3}^{ - }+{{\text{H}}_2}{\text{A}}{{\text{l}}_2}{\text{S}}{{\text{i}}_2}{{\text{O}}_{{\text{8(s)weathered}}}},$$
(10)

The CO2 fluxes observed during various typical periods remarkably differed. During snow melting, it was − 74.49 mmol m−2 day−1, much greater than that during snow accumulation. The specific reason might be that, in addition to elution in the layers of snow, the chemical reaction of soluble substances absorbing CO2 could have resulted in snowmelt water runoff. During the snow melting, CO2 flux was slightly lower than during glacial ice melting. This might be because the average temperature during snow melting was 0.5 °C, and the surface ice was almost frozen, inhibiting infiltration of the melting water, thus restricting chemical reactions, such as carbonate hydrolysis, inside debris (e.g., Eq. 9).

Because of rising temperatures and increased solar radiation, snow and ice melting increased, increasing hydrochemical erosion and leading to atmospheric CO2 sinking (negative CO2 flux) in the daytime. Inversely, mean nighttime CO2 fluxes were positive during snow accumulation (Table 2). Snowmelt water under the influence of solar radiation would refreeze at night (average temperature, − 7.0 °C), and some dissolved CO2 in the liquid meltwater would be released, producing positive flux.

CO2 storage (F s)

Daily atmospheric CO2 storage in the debris area of the Koxkar Glacier was between − 6.21 and 6.50 mmol m−2 day−1, averaging 2.66 × 10− 3 mmol m−2 day−1, or less than 1% of NGE, well below the ratio of CO2 storage to net ecosystem flux in grasslands, forests, and wetland ecosystems (Fortuniak et al. 2017; Rains et al. 2016). Reasons for this were primarily: (1) the instrument was set up only 2.0 m above the debris, and the smaller space limits CO2 storage, (2) porosity produced by glacial movement possibly reduced CO2 storage in the debris, and (3) compared with forests (Zhang et al. 2010), grasslands, and wetlands (Fortuniak et al. 2017), vegetation or microbial respiration is likely to be lower.

Discussion

Effective contribution area of CO2

Glacial topography is complex, given large differences in elevation, a large number of moraine hills, ice cliffs, and super-glacial rivers and lakes. Therefore, selecting the location of flux observation points and determining the effective range of instrument monitoring are the foundations of precise analyses and studies of regional carbon fluxes. In this study, the ART footprint tool was used to analyze CO2 flux data for three typical phases during a study period beginning in 2015 and ending in 2017. The farthest point distribution of the CO2 flux footprint is shown in Fig. 3, where significant differences are evident between the farthest points and distribution directions in 30-min data. The farthest points were located primarily in the major wind direction near the glacier middle streamline. The farthest distances varied with the period: glacier melting period (1499.13 m) > snow melting period (760.65 m) > snow accumulation period (573.14 m), and this was similar to wind speed findings (Fig. 2). Our results are, therefore, consistent with the synergistic effect of winds on CO2 flux found in urban and wet areas (Guo et al. 2018; Jia et al. 2010).

Fig. 3
figure 3

Distribution of points farthest from the eddy covariance station. a Snow accumulation period, b snow melting period, and c glacier ablation season

Full size image

Influenced by the terrain and wind speed and direction, the farthest distribution points during snow accumulation were concentrated between 110° and 345° (Fig. 3). This differed from the major wind direction, which was NW (Fig. 2). This might have resulted from relatively high temperatures at lower altitudes, where less wind was found (Han et al. 2008), the melting time of snow was shortened under solar radiation, and atmospheric CO2 was drawn down during elution of ions in the snow layer. The farthest distribution points were concentrated between 180° and 315° during snow melting, corresponding to wind speed and direction, thus indicating that the distribution direction of the farthest point tended to be gradually concentrated by increases in wind speed. These findings suggest that influences from grassland CO2 fluxes at the ridges and at the glacier end could be ignored, and CO2 fluxes observed at the eddy station were affected primarily by hydrochemical reactions of soluble substances at the ice (moraine)–gas interface.

Also, in 30-min data, flux contribution rates exceeded 75% for 95.74%, 92.45%, and 84.71% of all data during snow accumulation, snow melting, and glacial melting periods, respectively. Mean values for total contribution rates accounted for 93.30%, 91.39%, and 90.17% of the data in various periods, respectively. Flux footprint contribution areas were predominantly within 150 m of the EC station (cells were numbered 1, 6, 26, 51, 52, 56, 76, 77, 81, 82, 86, 87, 91, and 92; Fig. 4). The center point of the farthest cell (no. 92) was 190.39 m from the EC station, and the average maximum contribution rate was only 0.71% across the three typical periods. This also shows that those contribution areas with significant influences on CO2 flux were concentrated, confirming that grassland CO2 flux around the glaciers had little effect on observations at the EC station.

Fig. 4
figure 4

Spatial distribution of the CO2 footprint contribution in 100, 50 × 50 m cells around the EC station. a Snow accumulation period, b snow melting period, and c glacier ablation season

Full size image

Strong ice ablation and hydrochemical erosion was observed in cells 18, 19, 20, 22, 23, and 86, containing exposed ice cliffs, but these were far from the EC observation station, and their average contribution rate was less than 1.65%. Thus, their influence could also be ignored.

Furthermore, the atmosphere in the boundary layer was divided using the Monin–Obukhov length (L in Eq. 11) into steady state (L > 0) and non-steady state (L < 0) conditions. Results are shown in Table 3. The proportion of daytime spent in the steady state was 35.71% ± 8.92%, compared to 49.99% ± 20.94% during the night. This proportion was significantly lower than the value recorded at night during glacier ablation (69.90%), showing that the valley wind formed by solar radiation differences was affected by the terrain and the underlying surface, causing significant atmospheric turbulence in the boundary layer during the day. This finding is consistent with those of Wang et al. (2014b) who analyzed the microclimate of the Koxkar Glacier. In other words, during the day, the atmosphere was in unstable stratification, and the turbulence was strong. At night, the atmosphere at the near-surface appeared to be inverted, and turbulence was relatively weak. Meanwhile, under atmospheric steady-state conditions, the average percentages of the integral over the considered domain was 79.09% ± 1.84% during the day, slightly greater than found at night (78.27% ± 1.48%), but significantly lower than found during non-steady state conditions (89.95% ± 0.16% during the day, 89.70% ± 0.33% at night).

Table 3 Comparisons of the major parameters of the CO2 flux footprint under stable and non-stable atmospheres
Full size table
$$L= - \frac{{u_{*}^{3}{T_{\text{v}}}}}{{kg{Q_{{\text{vo}}}}}},$$
(11)

where \(k\) is the von Karman constant, \(g\) is gravitational acceleration, \({Q_{{\text{vo}}}}\) is the temperature of the near-surface layer, \({u_{\text{*}}}\) is the frictional velocity, and \({T_{\text{v}}}\) is the covariance of the vertical turbulent velocity and the temperature fluctuation, which is proportional to the turbulent heat flux (Sogachev and Dellwik 2017).

A detailed analysis showed that the distance of the farthest point in the CO2 footprint for each typical phase under steady-state atmospheric conditions averaged 202.61 ± 69.33 m, significantly greater than that found under non-steady state conditions (68.55 ± 10.34 m). Similarly, the farthest point (1491.13 m) under steady-state conditions was 3.95-fold the maximum distance found under non-steady state conditions. This is consistent with the findings of Wu et al. (2017), who found that the flux contribution area under steady-state conditions was larger than the range under non-steady state conditions in the winter wheat farmland ecosystem of the North China Plains.

Influences of temperature and precipitation

During our study period, a significant (p < 0.01) negative correlation was found between NGE and air temperature in the moraine area of Koxkar, yielding an R value of − 0.64, similar to results reported for forest vegetation (Arriga et al. 2017; Fortuniak et al. 2017; Sogachev and Dellwik 2017). The former was caused by the confiscation of additional atmospheric CO2 by hydrochemical reactions increased by the melting of ice and snow as temperatures increased; the latter was caused by higher temperatures, which promotes enzyme activities in plants, thus enhancing net productivity and increasing carbon absorption capacity (Li and Zhang 2015). In fact, the linear relationship for the glaciated area was good only on days free of precipitation (Fig. 5a). When precipitation was observed, the correlation coefficient between NGE and air temperature was only 0.36 (n = 275; Fig. 5b). Continuous heavy precipitation events caused glacial temperatures to fall, resulting in decreases in ice and snow melting, thus restricting the hydrolysis of soluble ions and affecting the confiscation of atmospheric CO2. However, a relationship between NGE and precipitation was not obvious (Fig. 5c).

Fig. 5
figure 5

Variations in NGE (mmol m−2 day−1) influenced by a daily temperature (°C) on days free of precipitation, b daily temperature (°C) on days with recorded precipitation, and c precipitation (mm)

Full size image

When precipitation was less than 8.8 mm on a given day, the variation of NGE with precipitation was not significant, possibly the result of the effective precipitation in Koxkar being concentrated over a 16- to 19-h period, whereby 68.2% of the precipitation fell within 2 h. Based on this, a small amount of precipitation over a short period of time cooled the study area to less than 0.16 °C, thereby rendering glacial melting intensity and hydrochemical erosion to relatively weak influences. On the other hand, after eluting aerosols during precipitation, dissolved atmospheric CO2 could participate in surface hydrochemical reactions. Although this could inhibit or affect changes in NGE resulting from less precipitation or shorter periods of precipitation, this effect was non-significant.

During the 34 days where daily precipitation exceeded 8.8 mm, the exponential fitting equation for NGE and precipitation yielded a correlation coefficient of 0.42 (Fig. 5c). Using linear regression, the coefficient drops to 0.27. On the whole, NGE increased with increasing precipitation, attributable not only to heavy precipitation, but also to long durations of precipitation. In one case, 24 h of continuous precipitation decreased the regional temperature by more than 0.92 °C compared with the previous day. Additionally, precipitation on August 5, 2015 was 13.2 mm, and the temperature fell 1.5 °C; precipitation on September 13, 2016 was 39.0 mm, and the temperature fell 1.4 °C; the precipitation on June 1, 2017 was 19.2 mm, and the temperature fell 1.9 °C. As a result, snow and ice melting slowed, hydrochemical erosion weakened, and CO2 drawdown at the gas–liquid interface was reduced. On the other hand, the precipitation carried large amounts of dissolved CO2 to the surface to participate in hydrochemical reactions, possibly inhibiting the CO2 drawdown. Perhaps CO2 was released into the atmosphere as a result of evaporation or freezing.

Influence of runoff volume

No significant relationship between total runoff and NGE was found on Koxkar Glacier, not only because NGE changes were affected by many atmospheric factors, but also because the complexity of the runoff water sources could cause changes in physical and chemical reactions in regional bodies of water. This was reflected primarily in: (1) snow and ice ablation and NGE undergoing the combined effects of temperature, precipitation, radiation, and other factors. In addition, NGE was affected by the amount of glacier ablation and the type, intensity, and composition of soluble matter (Gibson et al. 2017), (2) atmospheric CO2 as an open system, whereby changes in the water–gas interface flux are not only controlled by hydrochemical reactions and erosion media (the amount of glacial meltwater), but also by the solubility of CO2, the type and intensity of hydrochemistry, and the carbon cycle. The intensity of hydrochemical erosion in the glacier area was relatively complex and included the hydrolysis of carbonate and Na/K feldspar. The H+ generated by regional sulfide oxidation, as shown in Eq. (12), could also suppress the CO2 sink (Wang et al. 2010; Krawczyk and Bartoszewski 2008), (3) the englacial and subglacial runoff processes were complicated, and CO2 exchanged with the outside was limited, possibly inhibiting chemical erosion and CO2 drawdown (Wang et al. 2016; Krawczyk and Bartoszewski 2008), (4) the debris, which was affected by glacial melting, porosity, moisture content, and chemical composition, could also result in changes in the types and intensities of hydrochemical reactions (Krawczyk and Bartoszewski 2008; Brown 2002) that affected NGE.

$$4{\text{Fe}}{{\text{S}}_2}+15{{\text{O}}_2}+8{{\text{H}}_2}{\text{O}} \to 2{\text{F}}{{\text{e}}_2}{{\text{O}}_3}+8{\text{SO}}_{4}^{{2 - }}+16{{\text{H}}^+}$$
(12)

The relationship between daily runoff and NGE was further analyzed in the absence of precipitation, yielding a significant negative correlation between the two on a monthly statistical scale (Fig. 6). The daily runoff during snow ablation (from April to May) was the smallest, averaging 1.8 ± 0.24 m3·s−1. On the other hand, NGE was greatest compared to daily runoff, yielding a graphical slope of − 237.20 (Table 4). The daily runoff in the early stage of snow accumulation (from October to November) was 3.10 ± 0.56 m3 s−1, yielding a graphical slope of − 69.82. The daily runoff in the early stage of glacier ablation (June) was 3.76 ± 1.89 m3 s−1 (slope − 29.92) but during the late stage (August to September), it was 7.08 ± 4.33 m3 s−1 (slope − 27.55). At peak glacier ablation (July) daily runoff also peaked (9.24 ± 1.73 m3 s−1; slope − 15.02); only 6.33% of this was within the snow melting period. Therefore, the graphical slope of the relationship between daily runoff and NGE for the Koxkar Glacier area followed a trend, whereby it was greatest during the snow melting period followed by the snow accumulation period, the early glacial melting period, the late glacial melting period, and finally, the peak glacial melting period. Because of the damping effect of the snow body itself and the undeveloped hydrological channels under and in the ice during snow melting, runoff formation was slower than during ice melting, providing sufficient time for the full reactions of soluble substances and enhancement in CO2 sinking. Moreover, this was verified by the relatively small slope between NGE and runoff in the absence of snow, and a certain percentage of snow melting during early (June) and late (August to September) glacier ablation periods compared to the peak period (July).

Fig. 6
figure 6

Variations in NGE (mmol·m−2 day−1) influenced by daily runoff (m3·s−1) in various glacier melting periods in the Koxkar Glacier basin

Full size image
Table 4 Regression equations determined between net glacier-system CO2 exchange (mmol·m−2day−1) and runoff (m3·s−1) during various snow and ice melting periods
Full size table

Furthermore, when analyzing the same runoff in the Koxkar Glacier area, significant differences might be present in NGE at different stages. For example, when runoff was 8.9 m3 s−1, NGE was approximately − 61.59 mmol m−2 day−1 in the late melting period (August to September) and − 121.59 mmol m−2 day−1 at the peak of glacial melting (July), as shown in Fig. 6. This was primarily the result of the complex water sources and runoff formation processes: (1) the runoff sources were not limited to precipitation, groundwater, and glacial ice meltwater, but also included snow meltwater (Gibson et al. 2017; Jones et al. 2017). Differences in the soluble ionic compositions of various water sources and changes in ice and snow ablation areas, which can occur monthly, particularly in the latter, led to differences in the effective soluble ionic compositions of runoff waters. For example, a large number of terrigenous clastic rocks of the palaeo-ocean and carbonatites distributed around 3400 m above sea level, and the chemical reactions are shown in Eqs. (9) and (10) could occur, temporarily consuming atmospheric CO2. However, pyrite appeared in the moraine around 3700 m above sea level, as shown in Eq. (12), and the H+ generated could inhibit atmospheric CO2 sink, leading to large differences in the CO2 sink from hydrochemical erosion in runoff formation. (2) The effects of variability in hydrogeometric channels within ice over time, as well as the amount of dissolved or stored CO2 in the channel, could control the type and intensity of chemical reactions (Kociuba 2017), thus affecting NGE.

Conclusions

  1. 1.

    (1) Atmospheric CO2 flux averages − 58.68 mmol m−2 day−1 in the Koxkar Glacier area of western Mt. Tianshan, China. Even in winter, when ice and snow melt stagnates, the water chemistry reaction of elution in the snow layer caused by short-term solar radiation during the daytime can lead to small amounts of CO2 drawdown, indicating that the Koxkar Glacier area is a significant carbon sink.

  2. 2.

    (2) The CO2 flux footprint contribution areas were primarily within 150 m of the EC station and contribute, on average, 93.30%, 91.39%, and 90.17% of the total during snow accumulation, snow melting, and glacial melting, respectively. This shows that the contribution areas with significant influences on CO2 flux were concentrated, demonstrating that grassland CO2 flux around the glaciers have little effect at the EC stations.

  3. 3.

    (3) Under steady-state atmospheric conditions, the average proportions of CO2 flux in the major grid contribution areas are slightly higher in the day than at night, but significantly lower than in non-steady state conditions. The distance to the farthest point of the CO2 footprint under steady-state atmospheric conditions averages more than the average distance under non-steady state conditions (68.55 ± 10.34 m). Our data show that CO2 flux is affected primarily by hydrochemical erosion reactions in the glacier area.

  4. 4.

    (4) In areas of glacial melting, a good negative correlation between NGE and air temperature exists on precipitation-free days. Precipitation events can cause decreases in temperature, and when precipitation exceeds 8.8 mm in a given day, it can affect the confiscation of atmospheric CO2.