p CO 2 Dynamics of Stratiﬁed Reservoir in Temperate Zone and CO 2 Pulse Emissions During Turnover Events

: This study explores the dynamic changes in the partial pressure of CO 2 ( p CO 2 ) with depth, and the temporal variations of CO 2 net atmospheric ﬂux (NAF) in a stratiﬁed reservoir. A total of 16 ﬁeld campaigns were conducted from the summer stratiﬁcation to fall turnover period in 2017. A random forest (RF) model was developed to estimate the p CO 2 using concurrently measured water quality variables. The results showed that the vertical distribution of p CO 2 and associated temporal variations of the NAF are closely related to the stratiﬁcation strength of the reservoir. The reservoir surface p CO 2 was supersaturated (1542 µ atm) in summer (July 11), but this decreased to undersaturation as algae grew. Meanwhile, dissolved CO 2 continuously accumulated below the reservoir mixed-layer due to the thermal stratiﬁcation barrier and organic-rich ﬂoodwater intrusion. Vertical mixing began instantly as the stratiﬁcation strength began to weaken in mid-October, and the surface p CO 2 increased sharply up to 1934 µ atm. Consequently, the NAF drastically increased to 3235 mg − CO 2 m − 2 · day − 1 , which implies that the NAF changes seasonally and large CO 2 pulsing occurs during the turnover events. The results provide valuable information about p CO 2 variability and physical mixing processes, as well as carbon budget estimation in stratiﬁed reservoirs, and offer an improved understanding of these phenomena.


Introduction
Since 1750, concentrations of atmospheric greenhouse gases (GHGs) such as CO 2 , CH 4 , and N 2 O have increased continuously due to human activities. As of 2010, total GHG emissions were 49 Gt CO 2 -eq yr −1 ; in 2011, the atmospheric CO 2 concentration was 391 ppm, which constitutes a 40% increase over pre-industrial levels [1]. Conversion of GHGs emissions into global warming potential (GWP 100 ) produces the following component contributions: 76% of total GWP 100 from CO 2 , 16% from CH 4 , and 6.2% from N 2 O [2]. CO 2 not only accounts for the highest proportion of GHG emissions, but also exerts the greatest radiative forcing, which affects warming. Thus, CO 2 is the most important species in global carbon cycle analyses.
Most global GHG inventory studies have divided the terrestrial and oceanic biospheres to estimate carbon absorption and emission. Currently, sources of GHGs from the terrestrial biosphere are categorized according to IPCC guidelines [3] as relating to energy; industrial processes and product use (IPPU); agriculture, forestry, and other land use (AFOLU); or waste. Inland waters such as rivers, lakes, and reservoirs are not considered separately [4]. However, the CO 2 net atmospheric flux (NAF) through the interface between inland waters and the atmosphere is an important component of the global carbon cycle [5,6]. Previous work has suggested that about 50% of the carbon absorbed by the terrestrial biosphere is emitted from freshwater, which covers about 2% of the Earth's surface [4,7]. Cole et al. [8] collected 4665 samples from 1835 lakes to examine the dissolved CO 2 partial pressure (pCO 2 ) in the lake surface; 87% of them were supersaturated instead of in equilibrium with the atmosphere. Algesten et al. [9] calculated the carbon loss by mineralization and sedimentation and net CO 2 emission to the atmosphere for 79,536 lakes and 21 rivers in Scandinavian watersheds. They reported that 30-80% of the organic carbon flowing into the freshwater ecosystem was lost from the lakes, with mineralization in the lakes and subsequent emission to the atmosphere (mainly by diffusion) being the most important carbon loss process. Cole et al. [5] showed that the carbon NAF in freshwater was greater than that in surrounding areas in terms of emissions per unit area. Roehm et al. [10] reported that the pCO 2 in the Eastmain River, Quebec, Canada was 400-1800 µatm, i.e., continuously at supersaturation.
The large CO 2 emissions in inland waters arise from the fact that organic carbon (OC) precipitates about 50 times faster in inland waters than in oceans and accumulates in sediments, becoming a key factor in CO 2 supersaturation in the lake [5,[11][12][13][14]. OC reserves per unit area in reservoirs formed by artificial dams are much higher than those in natural lakes [15]. More than one million dams have been constructed and operated worldwide to supply water, control flooding, and generate hydroelectric power [16]. Dams interrupt the flow of natural rivers and increase the area of the reservoir, which accelerates the emission of carbon accumulated in reservoir sediments to the atmosphere through diverse pathways including diffusion, ebullition, and degassing processes [17,18]. According to Deemer et al. [19], as well as the sedimentation and decomposition of organic matter from outside the reservoir, artificial reservoirs have GHG emission sources which do not exist in natural lakes, including emissions from deposits exposed due to large fluctuations in the water level, and spillway emissions and emissions from hydroelectric power generation. Bevelhimer et al. [20] quantified GHG emissions from six hydropower reservoirs in the southwestern United States via selected emission mechanisms, namely diffusive emissions, ebullition, and tailwater degassing, and found diffusive emissions to be dominant. Under increased rainfall, flooding may transport large amounts of organic matter into the artificial reservoir from the upper watershed; this organic matter is accumulated in the lake and decomposed to emit final products such as gaseous CO 2 and CH 4 , which maintain the CO 2 supersaturation in the hypolimnion [21,22].
However, most of the reported carbon emission from reservoirs have been estimated using limited field data [5,23], for which the methodological and analytical uncertainties are quite high [6,24]. In particular, because carbon emissions from reservoir surface water vary widely over time and space [19,25], data measured at a specific time and place cannot represent overall CO 2 fluxes, which are very temporally and spatially dynamic [26,27]. In deep lakes and reservoirs, mixing occurs between the upper and lower layers when thermal stratification is destroyed, and the accumulated supersaturated CO 2 is spread to the surface layer. However, the CO 2 pulse emission mechanism that occurs during the turnover period (in which thermal stratification is destroyed) has not been sufficiently studied [6,19]. This CO 2 pulse emission process is closely related to thermal stratification and physical mixing processes in stratified reservoirs. However, data on CO 2 concentration vertical distributions in stratified reservoirs is limited; thus, large-scale carbon emissions after the turnover period have not been accurately accounted for in reservoir carbon emissions estimations [28][29][30].
The purposes of this study were to investigate variations in pCO 2 with water depth, and temporal changes of the CO 2 NAF in a stratified reservoir located in a mid-latitude monsoon climate zone, to elucidate the dynamic changes in pCO 2 and water quality variables with depth stemming from changes in the thermal stratification structure before and after the turnover period. A water thermometer chain was installed with probes at given depths to monitor changes in thermal stratification strength in the reservoir, and continuous measurements were made. Variations in pCO 2 with depth were measured directly using a nondispersive infrared (NDIR) method; water quality variables that can affect pCO 2 were measured simultaneously. In addition, a stepwise multiple linear regression (MLR) model and a random forest (RF) model for estimating pCO 2 were developed using measured water quality variables as explanatory variables, and compared with pCO 2 values calculated via thermodynamic equilibrium theory. These estimated pCO 2 values were analyzed against the measured values, according to the variations of stratification strength and water quality, and used for estimating CO 2 NAF based on the thin boundary layer method. The results of this study are expected to enhance our understanding of pCO 2 variability dynamics in stratified reservoirs and the mechanisms behind CO 2 pulse emissions during turnover periods.

Site Description
Daecheong Reservoir (36.5 • N, 127.5 • E), which constitutes the study area (Figure 1), was formed when the Daecheong Dam was built in the middle of the Geum River in Korea; the reservoir has a long, narrow shape. The dam fulfills multiple functions, providing water supply, flood control, and hydroelectric power generation. The total drainage area is 3204 km 2 and the surface area of the reservoir is 72.8 km 2 . The total reservoir capacity is 1490 million m 3 , and the effective storage capacity is 790 million m 3 . The reservoir has a mid-latitude monsoonal climate and exhibits a typical warm monomictic mixing regime [7], with full vertical mixing in winter and a steady stratification for the rest of the year. The average depth of the Daecheong Reservoir is about 20 m. During the study period, the maximum depth at the dam front was 42 m. The annual precipitation during the study year (2017) was 1138.0 mm, which is slightly greater than the average annual precipitation of 1110.4 mm over the last 10 years.
Water 2018, 10, x FOR PEER REVIEW 3 of 18 calculated via thermodynamic equilibrium theory. These estimated pCO2 values were analyzed against the measured values, according to the variations of stratification strength and water quality, and used for estimating CO2 NAF based on the thin boundary layer method. The results of this study are expected to enhance our understanding of pCO2 variability dynamics in stratified reservoirs and the mechanisms behind CO2 pulse emissions during turnover periods.

Site Description
Daecheong Reservoir (36.5° N, 127.5° E), which constitutes the study area (Figure 1), was formed when the Daecheong Dam was built in the middle of the Geum River in Korea; the reservoir has a long, narrow shape. The dam fulfills multiple functions, providing water supply, flood control, and hydroelectric power generation. The total drainage area is 3204 km 2 and the surface area of the reservoir is 72.8 km 2 . The total reservoir capacity is 1490 million m 3 , and the effective storage capacity is 790 million m 3 . The reservoir has a mid-latitude monsoonal climate and exhibits a typical warm monomictic mixing regime [7], with full vertical mixing in winter and a steady stratification for the rest of the year. The average depth of the Daecheong Reservoir is about 20 m. During the study period, the maximum depth at the dam front was 42 m. The annual precipitation during the study year (2017) was 1138.0 mm, which is slightly greater than the average annual precipitation of 1110.4 mm over the last 10 years.

Sampling and Analysis
Water quality is measured regularly (monthly) at four sites (R1-R4) in the Daecheong Reservoir; the site locations are shown in Figure 1. R4 was selected as the study site for its depth and strong thermal stratification during summer (June to August). The study period spanned July 2017 to November 2017 (five months) and included the summer season, when thermal stratification is strong, and fall, when turnover occurs due to vertical mixing. A buoy was installed at site R4 (36°28′28″ N, 127°29′16″ E), and the field study was performed at this fixed site. The reservoir depth varied between 37 and 42 m during the study period.

Sampling and Analysis
Water quality is measured regularly (monthly) at four sites (R1-R4) in the Daecheong Reservoir; the site locations are shown in Figure 1. R4 was selected as the study site for its depth and strong thermal stratification during summer (June to August). The study period spanned July 2017 to November 2017 (five months) and included the summer season, when thermal stratification is strong, and fall, when turnover occurs due to vertical mixing. A buoy was installed at site R4 (36 • 28 28" N, 127 • 29 16" E), and the field study was performed at this fixed site. The reservoir depth varied between 37 and 42 m during the study period. To investigate the thermal stratification structure in the reservoir, a thermometer chain was installed at R4, and continuous water temperature monitoring was performed at given water depths. Water thermometers (HOBO Water Temp Pro, Onset Co., Bourne, MA, USA) were installed at depths of 0, 2, 4, 6,8,10,12,14,16,18,21,24,27,30, and 33 m, and water temperature was measured at 10 min intervals. As both meteorological and hydrological data were needed to analyze the characteristics of the hydrodynamic mixing in the reservoir, data were collected from the Korea Meteorological Administration (KMA). Temperature, precipitation, wind direction, and wind speed data were obtained from an automated weather station located 1.2 km south of the monitoring point, and average local pressure (hPa) and solar radiation (MJ m −2 ) data were obtained from the Daejeon meteorological station located 15 km southwest. Hydrological and dam operational data (inflow, outflow, and reservoir water level) were collected from the K-water database [31].
A total of 16 field campaigns were conducted during the study period, especially during the turnover period, when the pCO 2 temporal variability in the reservoir increased as the water body became unstable. In these field campaigns, pH, salinity (ppt), electronic conductivity (Cond, µS·cm −1 ), water temperature (Tw, • C), and dissolved oxygen (DO, mg·L −1 ) were measured using multi water quality sensors YSI-6600 and YSI PRO (YSI, Yellow Springs, OH, USA). Sampling was performed according to the field sampling manual [32] at 2 to 5 m intervals at up to a 30 m depth. The water samples were collected so as to completely fill a 1 L sterilized bottle with no air bubbles and were immediately transferred to the laboratory for analysis. The pH was measured with YSI pH1200 (YSI, USA), which has a measurement error of ± 0.02 at pH 7 and ± 0.01 at pH 4-10. Water quality variables, including chemical oxygen demand (COD, mg·L −1 ), total organic carbon (TOC, mg·L −1 ), dissolved organic carbon (DOC, mg·L −1 ), chlorophyll a (Chl-a, mg·m −3 ), and turbidity (Turb, NTU), were analyzed according to the Standard Method [33]. The pH measured at the site was used as a reference, and the pH value measured in the laboratory was used for the analysis. pCO 2 (µatm) was measured with water depth using an in situ NDIR detector (Pro-Mini CO 2 , Pro-Oceanus Systems Inc., Bridgewater, NS, Canada), which measures the wet mole fraction of gaseous CO 2 in equilibrium with the surrounding water and uses an infrared (IR)-based detector to ensure accurate CO 2 measurements at low and high concentrations. The IR detector measures the gas mole fraction of CO 2 in the measuring cell, and the output is corrected for the total dissolved gas pressure (TDGP) change in the cell. The pressure sensor in the detector cell outputs a gas pressure, which is used internally to calibrate the CO 2 measurement values. The measurement range of the NDIR detector was set to 0-6000 µatm, and the error of the measurements was ± 2% (120 µatm) of the maximum range.

Indirect pCO 2 Estimation and Development of Statistical Models
To investigate the dynamic changes in pCO 2 with water depth, pCO 2 was directly measured in the field using the NDIR method. However, when only water quality data were available, pCO 2 was calculated indirectly using two different methods: first, the CO2SYS program, into which the measured pH and alkalinity (Alk) values were entered [34]; second, statistical models (MLR and RF), with water quality variables and pCO 2 as the explanatory and dependent variables, respectively. pCO 2 was estimated via a best-fit statistical model due to the existence of previous studies indicating that the pH/Alk calculation is quite sensitive to the accuracy of the measured pH [35][36][37].
The CO2SYS program analyzes carbonate systems based on thermodynamic equilibrium theory and requires input data such as pH, Alk, DIC, water temperature, and salinity. In this study, the water temperature measured in the field and the pH and Alk values measured in the laboratory were used. Carbonate analysis was performed with salinity = 0 (freshwater) samples, which reflected the freshwater characteristics of the reservoir under investigation.
The MLR and RF models used the pCO 2 measured 12 times at 0, 5, 10, 15, and 20 m as the dependent variable and the station pressure and water quality variables (pH, Tw, DO, Cond, Turb, Alk, DOC, COD, TOC, and Chl-a) as explanatory variables. The MLR model was developed using the R olsrr package [38,39], and the RF model was developed using the RF package [40]. The MLR model assumes that there is a linear relationship between one dependent variable and several independent variables and estimates the regression coefficient based on the given learning data. The MLR model was developed using a stepwise method that evaluates the performance of the given prediction model by excluding insignificant explanatory variables one by one. The RF model is an ensemble version of classification tree analysis and is used to model complex nonlinear relationships. The RF model randomly constructs trees with slightly different characteristics and determines the final predicted results from several decision trees using the average or majority voting method. The number of trees (n tree ; a major variable) was 500, and the number of randomly extracted variables was set to p/3, where p is the number of variables.
The performance of the MLR and RF models was evaluated via the root mean squared error (RMSE), adjusted coefficient of determination (Adj-R 2 ), and Akaike information criteria (AIC) values. Adj-R 2 values closer to 1 and smaller AIC and RMSE values indicate better model performance.

Analysis of Reservoir Thermal Stratification and Stability
Physical properties such as the thermal stratification structure and water body stability were analyzed during the study period using the rLakeAnalyzer [41]. The input data for the rLakeAnalyzer include reservoir topography, reservoir water level, wind speed, salinity, and water temperature with water depth. The thermocline depth (m), Schmidt stability (S t ), and Lake number (LN) were calculated to schematize the water temperature with depth in a time series and identify the water body thermal stratification structure and strength. S t is an indicator of water body stratification strength, and LN, an indicator of the ratio of internal vertical mixing strength, is derived from stratification strength and wind force [42,43]. Larger values of both indicators signify a more stable water body. The stabilization index of the water body was linked with the environmental variables profile data collected in the field campaigns, and the vertical distribution of pCO 2 was analyzed according to stratification strength.
The results observed at specific water depth (Figure 2b-f) also showed a strong negative correlation between measured pCO 2 and DO; similar results were found between DO and pCO 2 (R 2 = 0.74, p < 0.001) in Janauca Lake in the Amazon River flood plain [44]. However, in that case, the DO concentration in water increased and the concentration of dissolved CO 2 decreased due to active algal photosynthesis, which differs from the pattern observed in this study. The concentration of Chl-a (median = 1.1 mg·m −3 ) was very low in every layer except the surface layer in August and September, and the DO concentration in the reservoir mid-layer was decreased. Typically, DO concentration is influenced by gas exchange with the atmosphere, consumption by organic matter aerobic decomposition, and algal photosynthesis and respiration. Considering that measurements were taken primarily during the daytime during the survey period, DO reduction may have resulted largely from the decomposition of organic matter. Therefore, the strong inverse correlation between pCO 2 and DO and increase in pCO 2 in the water body may have been due to the aerobic decomposition of organic matter arising predominantly from the watershed.
Water 2018, 10, x FOR PEER REVIEW 6 of 18 0.74, p < 0.001) in Janauca Lake in the Amazon River flood plain [44]. However, in that case, the DO concentration in water increased and the concentration of dissolved CO2 decreased due to active algal photosynthesis, which differs from the pattern observed in this study. The concentration of Chl-a (median = 1.1 mg·m −3 ) was very low in every layer except the surface layer in August and September, and the DO concentration in the reservoir mid-layer was decreased. Typically, DO concentration is influenced by gas exchange with the atmosphere, consumption by organic matter aerobic decomposition, and algal photosynthesis and respiration. Considering that measurements were taken primarily during the daytime during the survey period, DO reduction may have resulted largely from the decomposition of organic matter. Therefore, the strong inverse correlation between pCO2 and DO and increase in pCO2 in the water body may have been due to the aerobic decomposition of organic matter arising predominantly from the watershed. The results from this study are consistent with those from previous studies [45][46][47][48], which demonstrate that DOC (allochthonous DOC) is an important source of CO2 supersaturation in terrestrial freshwater. Jonsson et al. [49] reported that pCO2 was increased due to aerobic decomposition of allochthonous DOC from the watersheds of 16 lakes in northern Sweden. However, the positive correlation shown between pCO2 and DOC in Jonsson et al. [49] was only observed herein at a water depth of 15 m (Figure 2e). This discrepancy may have arisen because most of the watershed DOC, which was transported via floodwaters, was decomposed before reaching the R4 site in front of the dam. The positive correlation between pCO2 and TOC, which is only evident in the mid-layer of the reservoir (15 and 20 m) where floodwaters containing organic matter were introduced, supports this inference.

Variations in Stratification Structure and Thermal Stability
Temporal changes in atmospheric temperature and the LN and St values analyzed by the rLakeAnalyzer are shown with temporal-depth changes in water temperature, as measured by the thermometer chain, in Figure 3a-d. Figure 3d also shows the dynamic changes in the metalimnion thickness and thermocline depth. The St and LN values were estimated using water temperaturewith-depth data, which were measured at 10 min intervals from July 11 to December 31, 2017. During the study period, the atmospheric temperature ranged from a minimum of -9.1 °C to a maximum of 31 °C. LN ranged from 0.1 to 123.9, and St ranged from 9.6 to 2823.8 J·m −2 .
(e) (f) The results from this study are consistent with those from previous studies [45][46][47][48], which demonstrate that DOC (allochthonous DOC) is an important source of CO 2 supersaturation in terrestrial freshwater. Jonsson et al. [49] reported that pCO 2 was increased due to aerobic decomposition of allochthonous DOC from the watersheds of 16 lakes in northern Sweden. However, the positive correlation shown between pCO 2 and DOC in Jonsson et al. [49] was only observed herein at a water depth of 15 m (Figure 2e). This discrepancy may have arisen because most of the watershed DOC, which was transported via floodwaters, was decomposed before reaching the R4 site in front of the dam. The positive correlation between pCO 2 and TOC, which is only evident in the mid-layer of the reservoir (15 and 20 m) where floodwaters containing organic matter were introduced, supports this inference.

Variations in Stratification Structure and Thermal Stability
Temporal changes in atmospheric temperature and the LN and S t values analyzed by the rLakeAnalyzer are shown with temporal-depth changes in water temperature, as measured by the thermometer chain, in Figure 3a-d. Figure 3d also shows the dynamic changes in the metalimnion thickness and thermocline depth. The S t and LN values were estimated using water temperature-with-depth data, which were measured at 10 min intervals from July 11 to December 31, 2017. During the study period, the atmospheric temperature ranged from a minimum of -9.1 • C to a maximum of 31 • C. LN ranged from 0.1 to 123.9, and S t ranged from 9.6 to 2823.8 J·m −2 .
At the beginning of the experiment on July 11, the air temperature was 26.7 • C and the S t was 2383 J·m −2 , displaying stable thermal stratification. The maximum stratification strength occurred on August 12, with a S t of 2674 J·m −2 . Then, a large amount of water entered the reservoir after five consecutive days of rainfall (125.5 mm), resulting in a temporary decrease in stratification strength. At this time, the LN value decreased sharply from 105.2 to 16.8, then increased to 120.6. The floodwaters propagated into the reservoir along the upper layer of the thermocline (the equal density layer) formed before the rainfall, increasing the depth of the existing thermocline from 15 m to 25 m and inducing an abrupt change in the thickness of the metalimnion layer.
As the air temperature dropped below 10 • C from October 30 onwards, the thickness of the metalimnion layer decreased sharply from 8.8 to 4.4 m, and the S t weakened to 790.2 J·m −2 . On November 27, the S t had decreased to~300 J·m −2 , and the water layer was completely mixed up to a depth of 30 m, thereby removing most of the thermal stratification. As the atmospheric temperature decreased, the stratification became weaker, and the thermocline depth also tended to At the beginning of the experiment on July 11, the air temperature was 26.7 °C and the St was 2383 J·m −2 , displaying stable thermal stratification. The maximum stratification strength occurred on August 12, with a St of 2674 J·m −2 . Then, a large amount of water entered the reservoir after five consecutive days of rainfall (125.5 mm), resulting in a temporary decrease in stratification strength. At this time, the LN value decreased sharply from 105.2 to 16.8, then increased to 120.6. The floodwaters propagated into the reservoir along the upper layer of the thermocline (the equal density layer) formed before the rainfall, increasing the depth of the existing thermocline from 15 m to 25 m and inducing an abrupt change in the thickness of the metalimnion layer.
As the air temperature dropped below 10 °C from October 30 onwards, the thickness of the metalimnion layer decreased sharply from 8.8 to 4.4 m, and the St weakened to 790.2 J·m −2 . On November 27, the St had decreased to ~300 J·m −2 , and the water layer was completely mixed up to a depth of 30 m, thereby removing most of the thermal stratification. As the atmospheric temperature decreased, the stratification became weaker, and the thermocline depth also tended to deepen gradually. In addition, the correlation between atmospheric temperature and stratification strength (St) was high, with a coefficient of 0.95.

Comparisons of Measured and Estimated pCO 2 Values
A total of 58 pCO 2 values (µatm) measured with water depth were compared with the pCO 2 values estimated by CO2SYS (Figure 4). The pCO 2 values estimated by CO2SYS were 49.4% lower than the average measured pCO 2 . The RMSE between the estimated and measured pCO 2 values was 1531 µatm for the entire dataset, 541 µatm for the surface layer (water depth = 0 m), and 2038 µatm for the bottom layer (water depth = 20 m). Thus, the deviation increased as the water depth increased (Figure 4). Moreover, the coefficient of determination (R 2 ) was 0.43, indicating low significance. These results are consistent with previous studies performed at similar pH and DIC concentration ranges [32,35]. As mentioned in previous studies, this error between the measured and calculated pCO 2 values can be attributed primarily to inaccurate pH measurements, and the error is greater at lower pH [35]. Furthermore, the dynamic variation in dissolved CO 2 concentration-due to the decomposition of organic matter contained in floodwater introduced into the mid-layer of the reservoir-differs from the thermodynamic equilibrium assumption used by CO2SYS [6,47].
for the bottom layer (water depth = 20 m). Thus, the deviation increased as the water depth increased (Figure 4). Moreover, the coefficient of determination (R 2 ) was 0.43, indicating low significance. These results are consistent with previous studies performed at similar pH and DIC concentration ranges [32,35]. As mentioned in previous studies, this error between the measured and calculated pCO2 values can be attributed primarily to inaccurate pH measurements, and the error is greater at lower pH [35]. Furthermore, the dynamic variation in dissolved CO2 concentration-due to the decomposition of organic matter contained in floodwater introduced into the mid-layer of the reservoir-differs from the thermodynamic equilibrium assumption used by CO2SYS [6,47]. Thus, MLR and RF models were developed from major environmental variables to better predict pCO2 in water. The MLR model used a stepwise method for all explanatory variables, which included air pressure (Pa), Tw, Cond, DO, pH, COD, TOC, DOC, and Alk. Similar methods to those used in this study have been applied previously [35], however the previous study developed MLR and RF models using water quality data from only the surface layer of the water bodies (n = 10), and hence Chl-a and Turb-which herein showed little variation except in the surface layer-were determined to be major factors. In the current statistical evaluations, Chl-a and Turb contributed little to the prediction performance of the model; thus, they were excluded from the explanatory variables. The selected optimal MLR model included explanatory variables Tw, DO, pH, and DOC (Equation 1). The Adj-R 2 and AIC values were 0.866 and 577.5, respectively. (1) The significance of the explanatory variables in the RF prediction model can be placed in the following order: DO > pH > Tw > Pa > Cond > Alk > DOC > COD > TOC. To prevent overfitting, cross validation was used to select the most parsimonious model; DO, pH, Tw, Pa, Cond, and DOC were selected as predictors, and Adj-R 2 was 0.976, which indicates very high prediction performance ( Figure 5). Thus, MLR and RF models were developed from major environmental variables to better predict pCO 2 in water. The MLR model used a stepwise method for all explanatory variables, which included air pressure (P a ), Tw, Cond, DO, pH, COD, TOC, DOC, and Alk. Similar methods to those used in this study have been applied previously [35], however the previous study developed MLR and RF models using water quality data from only the surface layer of the water bodies (n = 10), and hence Chl-a and Turb-which herein showed little variation except in the surface layer-were determined to be major factors. In the current statistical evaluations, Chl-a and Turb contributed little to the prediction performance of the model; thus, they were excluded from the explanatory variables. The selected optimal MLR model included explanatory variables Tw, DO, pH, and DOC (Equation (1)). The Adj-R 2 and AIC values were 0.866 and 577.5, respectively. 11, 435.74 (1) The significance of the explanatory variables in the RF prediction model can be placed in the following order: DO > pH > Tw > P a > Cond > Alk > DOC > COD > TOC. To prevent overfitting, cross validation was used to select the most parsimonious model; DO, pH, Tw, P a , Cond, and DOC were selected as predictors, and Adj-R 2 was 0.976, which indicates very high prediction performance ( Figure 5).
The pCO 2 measured by the NDIR method ranged from 215 to 6000 (mean = 2008) µatm, and the predicted values from the MLR and RF models ranged from −290 to 4679 (mean: 1797) µatm and 425 to 5525 (mean: 1895) µatm, respectively. The errors between the values estimated by the MLR and RF models and the measured values were 18.5% and 13.0%, respectively, which were much improved over the error for CO2SYS estimated values (49.4%) based on thermodynamic equilibrium theory ( Figure 5). Thus, the unmeasured pCO 2 values were estimated using the RF model, which featured the lowest error and superior prediction performance.
to 5525 (mean: 1895) µatm, respectively. The errors between the values estimated by the MLR and RF models and the measured values were 18.5% and 13.0%, respectively, which were much improved over the error for CO2SYS estimated values (49.4%) based on thermodynamic equilibrium theory ( Figure 5). Thus, the unmeasured pCO2 values were estimated using the RF model, which featured the lowest error and superior prediction performance.

Dynamic Vertical Variations of pCO2 and Water Quality Variables
Temporal changes in daily precipitation (mm), daily mean wind speed (m·s −1 ), and Chl-a concentration in the surface layer (mg·m −3 ) are shown in Figure 6; the 12 campaign dates are denoted (a-l), respectively. Chl-a concentration data collected by K-water are also presented. Vertical profiles of pH, DOC, DO, and Tw for each survey date are shown in Figure 7. Variations of the pCO2 with depth are visualized along with the pCO2 predicted by the RF model and measured pCO2. The pCO2 measurement instrument failed after the July 27 measurement; the month of August was not included in the analysis.
On July 11 (Figure 7a), at the beginning of the study, thermal stratification was present. The pCO2 in the surface layer was 1542 µatm, and thus supersaturated rather than in equilibrium with the atmosphere. The pH was high nearer the surface, with values of 8.5 or more to a depth of 10 m, and relatively low (7.5) at 20 m. The vertical distribution of pCO2 was inversely correlated with pH, and the correlation coefficient was greater at a depth of 20 m than at the surface layer. DO sag phenomena occurred on July 20 ( Figure 7b) and July 27 ( Figure 7c); DO drastically decreased at depths of 5-15 m (the upper part of the thermocline was at 15 m). The pCO2 values were in the range of 3127-3155 µatm, and the DOC concentration was 2.448-2.738 mg·L −1 , higher than the values at other water depths. The increase in pCO2 in the upper part of the thermocline was closely related to the rainfall event. Due to the 28.5 mm of rainfall experienced from July 14 to 15, a mid-layer density current containing organic matter flowed down to a depth of 10 m (i.e., in the upper part of the thermocline, which was equal in density current); this was related to the reduction in DO and pH and increase in pCO2 from organic matter decomposition. The maximum pCO2 values (3638-3895 µatm) were observed at depths of 10-12 m (above the thermocline) on July 27 (Figure 7c) after 25 mm of rainfall and again on July 24. Generally, heavy rainfall transports large amounts of organic matter from the watershed; pH decreases and pCO2 increases as organic matter decomposes in the water [50][51][52]. These results are consistent with a study by Jonsson et al. [47] showing high correlations between DOC concentration and pCO2 in 16 lakes in northern Sweden. In this study, aqueous pCO2 was

Dynamic Vertical Variations of pCO 2 and Water Quality Variables
Temporal changes in daily precipitation (mm), daily mean wind speed (m·s −1 ), and Chl-a concentration in the surface layer (mg·m −3 ) are shown in Figure 6; the 12 campaign dates are denoted (a-l), respectively. Chl-a concentration data collected by K-water are also presented. Vertical profiles of pH, DOC, DO, and Tw for each survey date are shown in Figure 7. Variations of the pCO 2 with depth are visualized along with the pCO 2 predicted by the RF model and measured pCO 2 . The pCO 2 measurement instrument failed after the July 27 measurement; the month of August was not included in the analysis.
On July 11 (Figure 7a), at the beginning of the study, thermal stratification was present. The pCO 2 in the surface layer was 1542 µatm, and thus supersaturated rather than in equilibrium with the atmosphere. The pH was high nearer the surface, with values of 8.5 or more to a depth of 10 m, and relatively low (7.5) at 20 m. The vertical distribution of pCO 2 was inversely correlated with pH, and the correlation coefficient was greater at a depth of 20 m than at the surface layer. DO sag phenomena occurred on July 20 ( Figure 7b) and July 27 ( Figure 7c); DO drastically decreased at depths of 5-15 m (the upper part of the thermocline was at 15 m). The pCO 2 values were in the range of 3127-3155 µatm, and the DOC concentration was 2.448-2.738 mg·L −1 , higher than the values at other water depths. The increase in pCO 2 in the upper part of the thermocline was closely related to the rainfall event. Due to the 28.5 mm of rainfall experienced from July 14 to 15, a mid-layer density current containing organic matter flowed down to a depth of 10 m (i.e., in the upper part of the thermocline, which was equal in density current); this was related to the reduction in DO and pH and increase in pCO 2 from organic matter decomposition. The maximum pCO 2 values (3638-3895 µatm) were observed at depths of 10-12 m (above the thermocline) on July 27 (Figure 7c) after 25 mm of rainfall and again on July 24. Generally, heavy rainfall transports large amounts of organic matter from the watershed; pH decreases and pCO 2 increases as organic matter decomposes in the water [50][51][52]. These results are consistent with a study by Jonsson et al. [47] showing high correlations between DOC concentration and pCO 2 in 16 lakes in northern Sweden. In this study, aqueous pCO 2 was supersaturated throughout the various DOC ranges (0.4-9.9 mg·L −1 ), which may have been caused by allochthonous DOC from the watershed and heterotrophic respiration in the reservoir.
On September 6 (Figure 7d), the thermocline was at a depth of 25 m, which was lower by about 10 m than that on July 27. The depths at which DO concentration was depleted and the pCO 2 value was 2000 µatm or higher were extended to 10-30 m. The vertical distributions of pCO 2 , DO, and pH changed significantly as the floodwater, which included large amounts of organic matter, flowed into the upper part of the thermocline (25 m depth) after 57 mm of precipitation fell on July 31. The pCO 2 of the surface layer was 346 µatm, which was undersaturated with regard to the atmospheric equilibrium state, but tended to significantly increase with increasing water depth. In particular, at a depth of 10 m, pCO 2 had increased considerably to 4022 µatm and maintained similar concentrations up to 30 m in depth. DO was minimized at 10 m with a value of 1.59 mg/L, and the DO concentrations were low, with values of 3.87 mg/L or less down to 30 m. At this time, the maximum wind speed was 4.5 m·s −1 , and the pCO 2 in the surface layer seemed to be in equilibrium with the atmospheric concentration as a strong wind blew. In contrast, the drastic decrease in pCO 2 in the surface layer in September was closely related to an increase in algal biomass after August (Chl-a concentration on August 2 = 1.0 mg·m −3 → September 7 = 7.8 mg·m −3 ) (Figure 6). Dense algae in the surface layer consume CO 2 during the day and increase the DO [53]. This is consistent with previous studies suggesting that plankton metabolism significantly influences pCO 2 , and that higher DO concentrations are associated with lower CO 2 emissions from the reservoir surface layer [20,44].  On September 6 (Figure 7d), the thermocline was at a depth of 25 m, which was lower by about 10 m than that on July 27. The depths at which DO concentration was depleted and the pCO2 value was 2000 µatm or higher were extended to 10-30 m. The vertical distributions of pCO2, DO, and pH changed significantly as the floodwater, which included large amounts of organic matter, flowed into the upper part of the thermocline (25 m depth) after 57 mm of precipitation fell on July 31. The pCO2 of the surface layer was 346 µatm, which was undersaturated with regard to the atmospheric equilibrium state, but tended to significantly increase with increasing water depth. In particular, at a depth of 10 m, pCO2 had increased considerably to 4022 µatm and maintained similar concentrations up to 30 m in depth. DO was minimized at 10 m with a value of 1.59 mg/L, and the DO concentrations were low, with values of 3.87 mg/L or less down to 30 m. At this time, the maximum wind speed was 4.5 m·s −1 , and the pCO2 in the surface layer seemed to be in equilibrium with the atmospheric concentration as a strong wind blew. In contrast, the drastic decrease in pCO2 in the surface layer in September was closely related to an increase in algal biomass after August (Chl-a concentration on August 2 = 1.0 mg·m −3 → September 7 = 7.8 mg·m −3 ) ( Figure 6). Dense algae in the surface layer consume CO2 during the day and increase the DO [53]. This is consistent with previous studies suggesting that plankton metabolism significantly influences pCO2, and that higher DO concentrations are associated with lower CO2 emissions from the reservoir surface layer [20,44].
On September 12 (Figure 7e), the pCO2 in the surface layer was 410 µatm, similar to the atmospheric pCO2 level. Meanwhile, the pCO2 values at depths of 14-25 m were high, ranging from 4460 to 6000 µatm due to the decomposition of organic matter in floodwater. On October 10, the pCO2 in the surface layer was at its minimum value, of 215 µatm. However, the pCO2 in the mid-layer (depths of 14-25 m) maintained high values of 3000 µatm or more, possibly due to the heavy rainfall event (86.5 mm·day −1 ) on September 11. Drastic variations in DO, pH, and pCO2 were observed below the mixed-layer of the reservoir as large amounts of organic matter were introduced with the rainwater inflow and attendant organic matter decomposition.
As the thermal stratification strength weakened, the thermocline became deeper. On October 17 (Figure 7g), the pCO2 in the surface layer was 410 µatm, similar to the atmospheric concentration. However, the pCO2 was 2000 µatm or higher at depths of 10 m or less, gradually increasing with depth. On October 23 (Figure 7h), the pCO2 in the surface layer increased drastically to 1934 µatm, and similar pCO2 values (1751-2388 µatm) were observed from the surface layer to a depth of 25 m as the atmospheric temperature dropped and the thermocline weakened after October 10. At 30 m, On September 12 (Figure 7e), the pCO 2 in the surface layer was 410 µatm, similar to the atmospheric pCO 2 level. Meanwhile, the pCO 2 values at depths of 14-25 m were high, ranging from 4460 to 6000 µatm due to the decomposition of organic matter in floodwater. On October 10, the pCO 2 in the surface layer was at its minimum value, of 215 µatm. However, the pCO 2 in the mid-layer (depths of 14-25 m) maintained high values of 3000 µatm or more, possibly due to the heavy rainfall event (86.5 mm·day −1 ) on September 11. Drastic variations in DO, pH, and pCO 2 were observed below the mixed-layer of the reservoir as large amounts of organic matter were introduced with the rainwater inflow and attendant organic matter decomposition.
As the thermal stratification strength weakened, the thermocline became deeper. On October 17 (Figure 7g), the pCO 2 in the surface layer was 410 µatm, similar to the atmospheric concentration. However, the pCO 2 was 2000 µatm or higher at depths of 10 m or less, gradually increasing with depth. On October 23 (Figure 7h), the pCO 2 in the surface layer increased drastically to 1934 µatm, and similar pCO 2 values (1751-2388 µatm) were observed from the surface layer to a depth of 25 m as the atmospheric temperature dropped and the thermocline weakened after October 10. At 30 m, the DO concentration was as low as 3.1 mg·L −1 , and the pCO 2 was 3640 µatm, the maximum value. These results indicate that the CO 2 accumulated in the mid-and lower layers of the reservoir was transferred to the surface layer by vertical mixing as the thermal stratification in the reservoir weakened (S t = 887 J·m −2 , LN = 9.3). October 23 (Figure 7h

Temporal Variations of the CO 2 NAF
The CO 2 NAF between the reservoir surface and the atmosphere was estimated using the thin boundary layer method [54,55]. The difference between the measured pCO 2 values in the water and the atmosphere was used as the concentration gradient. The gas exchange coefficient was derived as a function of wind speed and water temperature using the Wanninkhof and Knox [56] equation. The NAF shows (Figure 8) large temporal variations. There was an overall negative NAF (i.e., absorption of CO 2 from the atmosphere) from July to mid-October, when the primary production of the reservoir increased with water temperature. This dramatically changed to a positive NAF (i.e., emission of CO 2 from the reservoir) after the turnover event. On November 21, 2017, the NAF increased up to 3235 mg-CO 2 m −2 ·day −1 , and further to 16,519 mg-CO 2 m -2 ·day −1 in January 2018 (not shown in Figure 8), which is much higher than the average NAF of 1400 mg-CO 2 m −2 ·day −1 from other temperate reservoirs reported in Louis et al. [57]. These results suggest that the CO 2 NAF in stratified reservoirs changes seasonally. Moreover, a large CO 2 pulse emission could initiate at the reservoir surface during the turnover period.
Water 2018, 10, x FOR PEER REVIEW 13 of 18

Temporal Variations of the CO2 NAF
The CO2 NAF between the reservoir surface and the atmosphere was estimated using the thin boundary layer method [54,55]. The difference between the measured pCO2 values in the water and the atmosphere was used as the concentration gradient. The gas exchange coefficient was derived as a function of wind speed and water temperature using the Wanninkhof and Knox [56] equation. The NAF shows (Figure 8) large temporal variations. There was an overall negative NAF (i.e., absorption of CO2 from the atmosphere) from July to mid-October, when the primary production of the reservoir increased with water temperature. This dramatically changed to a positive NAF (i.e., emission of CO2 from the reservoir) after the turnover event. On November 21, 2017, the NAF increased up to 3235 mg-CO2 m −2 ·day −1 , and further to 16,519 mg-CO2 m -2 ·day −1 in January 2018 (not shown in Figure 8), which is much higher than the average NAF of 1400 mg-CO2 m −2 ·day −1 from other temperate reservoirs reported in Louis et al. [57]. These results suggest that the CO2 NAF in stratified reservoirs changes seasonally. Moreover, a large CO2 pulse emission could initiate at the reservoir surface during the turnover period.

Conclusions
Reservoirs have been reported as major emitters of GHGs, but the temporal variability of pCO2 in stratified reservoirs and CO2 pulse emissions mechanisms during fall turnover events have not been sufficiently studied. In this study, variations in pCO2 vertical distribution with stratification strength were examined in a monomictic reservoir located in a temperate region. pCO2 vertical mixing and supersaturation in the surface layer were explored, along with the changing thermal stratification structure before and after a turnover period in which thermal stratification was destroyed.
Calculations of pCO2 using pH and Alk values based on the thermodynamic equilibrium theory (in CO2SYS) showed increased error relative with increasing water depth; CO2SYS underestimated the measured values by 49.4%. This error can be attributed to inaccurate pH measurements and dynamic variation of the dissolved CO2 concentration due to the decomposition of organic matter, which differs from the thermodynamic equilibrium assumption used by CO2SYS. The MLR and RF models that were developed to estimate pCO2 in the reservoir showed greatly improved errors of 18.5% and 13.0%, respectively, compared to the mean measured pCO2 value.
The pCO2 distributions with water depth during the study period were closely related to rainfall events and the strength of the thermal stratification. Until the end of July, pCO2 was supersaturated in the surface layer and not in equilibrium with the atmosphere. However, after August, pCO2 rapidly

Conclusions
Reservoirs have been reported as major emitters of GHGs, but the temporal variability of pCO 2 in stratified reservoirs and CO 2 pulse emissions mechanisms during fall turnover events have not been sufficiently studied. In this study, variations in pCO 2 vertical distribution with stratification strength were examined in a monomictic reservoir located in a temperate region. pCO 2 vertical mixing and supersaturation in the surface layer were explored, along with the changing thermal stratification structure before and after a turnover period in which thermal stratification was destroyed.
Calculations of pCO 2 using pH and Alk values based on the thermodynamic equilibrium theory (in CO2SYS) showed increased error relative with increasing water depth; CO2SYS underestimated the measured values by 49.4%. This error can be attributed to inaccurate pH measurements and dynamic variation of the dissolved CO 2 concentration due to the decomposition of organic matter, which differs from the thermodynamic equilibrium assumption used by CO2SYS. The MLR and RF models that were developed to estimate pCO 2 in the reservoir showed greatly improved errors of 18.5% and 13.0%, respectively, compared to the mean measured pCO 2 value.
The pCO 2 distributions with water depth during the study period were closely related to rainfall events and the strength of the thermal stratification. Until the end of July, pCO 2 was supersaturated in the surface layer and not in equilibrium with the atmosphere. However, after August, pCO 2 rapidly decreased to below saturation as algae grew. In contrast, in the middle layer of the reservoir, floodwater carrying organic matter flowed into the upper part of the thermocline (the equal density layer after the rainfall event), and the middle layer DO concentration decreased due to the aerobic decomposition of the organic matter. Subsequent increases of the pCO 2 were observed at the same layer. During periods of strong thermal stratification in the reservoir, high pCO 2 values were maintained in the middle layer without spreading to the surface layer. However, after mid-October, when the stratification began to weaken, water instantly mixed between the upper and lower layers, producing relatively constant pCO 2 , pH, DO, and DOC from the surface to a depth of 28 m. As a result, the pCO 2 in the surface layer increased sharply, and the pCO 2 in the deep layer decreased, reaching equilibrium. These results suggest that large-scale CO 2 pulse emissions were initiated at the reservoir surface during the turnover period.
The results of this study are important for understanding the carbon cycle in stratified reservoirs located in temperate zones, and particularly the "pulsed emissions" of carbon during turnover periods. However, as the study was conducted near the front of the dam, the results may not be representative of the whole domain of the Daecheong Reservoir. To analyze spatial variability, this study will be expanded to further locations in the Daecheong Reservoir in the future. In addition, the temporal and spatial changes of GHGs in the reservoir will be examined by applying numerical hydrodynamic and water quality modeling techniques, in order to comprehensively evaluate the processes of carbon emission from the stratified reservoir.

Conflicts of Interest:
The authors declare no conflict of interest.