Conceptual models of dissolved carbon fluxes in a two-layer stratified lake: interannual typhoon responses under extreme climates

Extreme climates affect the seasonal and interannual patterns of carbon (C) distribution in lentic ecosystems due to the regimes of river inflow and thermal stratification. Typhoons rapidly load substantial amounts of terrestrial C into smaller subtropical lakes (i.e., Yuan-Yang Lake in Taiwan, hereafter referred to as YYL), renewing and mixing the water column. We developed a conceptual dissolved C model and hypothesized that allochthonous C loading and river inflow intrusion may affect the dissolved inorganic C (DIC) and dissolved organic C (DOC) distributions in a small subtropical lake under these extreme climates. A two-layer conceptual C model was developed to explore how the DIC and DOC fluxes respond to typhoon disturbances on seasonal and interannual timescales in YYL while simultaneously considering autochthonous processes such as algal photosynthesis, remineralization, and vertical transformation. To compare the temporal patterns of fluxes between typhoon years (2015–2016) and non-typhoon years (2017–2018), monthly field samples were obtained and their DIC, DOC, and chlorophyll a (Chl a) concentrations measured. The results demonstrated that net ecosystem production was 3.14 times higher in typhoon years than in non-typhoon years. These results suggested that a loading of allochthonous C was the most crucial driver of the temporal variation in C fluxes in typhoon years because of changes in physical and biochemical processes, such as photosynthesis, mineralization, and vertical transportation. However, the lowered vertical transportation rate shaped the seasonal C in non-typhoon years due to thermal stratification within this small subtropical lake.

1 Introduction

The Sixth Assessment Report from the Intergovernmental Panel for Environmental Change (Masson-Delmotte, 2021) suggests that not only is air temperature going to increase by at least about 1.5 ^∘C by 2050 but high-intensity storms and drought events will also become more frequent as a result of global warming and climate change. With respect to freshwater ecosystems, extreme climates may change the mixing regimes of freshwater columns (Kraemer et al., 2021; Maberly et al., 2020; Woolway et al., 2020) and induce heat wave events (Woolway et al., 2021a, b), droughts (Marcé et al., 2019), and floods (Woolway et al., 2018). Freshwater ecosystems store around 0.32–1.8 Pg C yr⁻¹, which is approximately equivalent to shallow coastal areas, thereby providing carbon (C) to food webs that support human resources by, for example, acting as processing hot spots in regional C cycles (Aufdenkampe et al., 2011; Cole et al., 2007; Engel et al., 2018; Lauerwald et al., 2015; Raymond et al., 2013).

Carbon fluxes in small lakes (lake area < 1 km²) are sensitive to climate change due to the ease with which C mixes within water columns (Doubek et al., 2021; MacIntyre et al., 2021; Winslow et al., 2015). Moreover, storms induce dramatic changes in thermal stratification and water inflows (Lin et al., 2022; Olsson et al., 2022b; Vachon and Del Giorgio, 2014; Woolway et al., 2018). River inflows and wind turbulence have been shown to release allochthonous C from sediments into the water column after storm events in small stratified lakes (Bartosiewicz et al., 2015; Czikowsky et al., 2018; Vachon and Del Giorgio, 2014). These small lakes account for 25 % to 35 % of the total area of the Earth's surface lakes (Cole et al., 2007; Downing et al., 2006; Raymond et al., 2013). Compared with the case in larger lakes, our understanding of C fluxes in small lakes remains uncertain, as small lakes have usually been ignored in calculations of C fluxes on a global scale (Cole et al., 2007; Raymond et al., 2013). Thus, elucidation of the C fluxes in small lakes under extreme weather conditions is the key to optimizing estimations of global C fluxes in extreme climates.

Understanding the influences of physical, hydrological, and biogeochemical processes on the fate of C fluxes in smaller lake ecosystems is challenging work (Aufdenkampe et al., 2011; Cole et al., 2007; Raymond et al., 2013; Tranvik et al., 2009; Vachon et al., 2021; Woolway et al., 2018). The physical and biogeochemical regimes under climate change remain uncertain, such as biological compositions, mixing regimes, morphometric characteristics, and air–water energy fluxes (evaporation and transpiration) (Woolway et al., 2020). The dissolved inorganic carbon (DIC) concentration is an important factor in estimating CO₂ fluxes within lake ecosystems (Smith, 1985). Among C fluxes in a freshwater body, the partial pressure of CO₂ (pCO₂), defined as CO₂ emission across the air–water interface, is affected by DIC, water temperature, wind speed, and pH (Jähne et al., 1987; Smith, 1985). River inflows, sediment C burial, and heterotrophic respiration in the water column contribute to DIC dynamics in lakes (Hope et al., 2004; Vachon et al., 2021); simultaneously, autotrophic organisms, such as plankton and submerged vegetation, capture DIC via photosynthesis (Amaral et al., 2022; Nakayama et al., 2020, 2022). Moreover, calcification and mineralization may consume dissolved oxygen within water, inducing uncertainty in pCO₂ estimation (Hanson et al., 2015; Lin et al., 2022; Nakayama et al., 2022). Dissolved organic carbon (DOC) might contribute to CO₂ emission from lake water to the atmosphere through mineralization and remineralization within lake ecosystems (Hanson et al., 2015; Sobek et al., 2005). In subtropical freshwater ecosystems, the DOC concentration is a vital factor in describing variances in mineralization and remineralization rates for dissolved C (Lin et al., 2022; Shih et al., 2019).

Typhoons might significantly impact C distributions within the water columns in subtropical regions (Chiu et al., 2020; Lin et al., 2022). Kossin et al. (2013) investigated global storm events with an accumulated rainfall of about 50 mm, which accounts for approximately 10 %–40 % of precipitation in a subtropical typhoon event. Some studies found that extreme rainstorms not only impact the dissolved carbon in large lakes and catchments due to weathering (Sun et al., 2021; Zhou et al., 2023) but also that typhoon disturbances quickly mix, renew, or dilute the water in small subtropical lakes (Kimura et al., 2012, 2017; Lin et al., 2022). However, the complex interactions between biogeochemical and physical regimes for autochthonous and allochthonous C introduce uncertainty in elucidating the complete patterns between typhoons and dissolved C concentrations in small subtropical lakes. This uncertainty hinders our understanding of the seasonal and interannual variations in DIC and DOC concentrations (Lin et al., 2022). Thus, to understand the seasonal regimes and to estimate C fluxes in subtropical lakes, we investigated the variations in DIC and DOC due to typhoon disturbances.

Typhoons' effects on C fluxes have previously been studied in a small, two-layer stratified subtropical lake, Yuan-Yang Lake (YYL) in Taiwan (Chiu et al., 2020; Jones et al., 2009; Lin et al., 2021, 2022). Jones et al. (2009) used a conceptual hydrology model and sensor data to estimate CO₂ emission in YYL during typhoon disturbances that occurred in October 2004: 2.2–2.7 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ of CO₂ was released into the atmosphere. CO₂ emissions into the atmosphere were recorded at around 3.0–3.7 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ in summer and autumn 2015 because of substantial loads of terrestrial C via river inflows after strong typhoons in YYL (Chiu et al., 2020). In particular, vertical mixing, thermal stratification, and river retention regimes were found to be essential physical processes in the C fluxes in YYL (Lin et al., 2021, 2022). The results of these studies suggest that river intrusion and thermal stratification are key factors shaping the seasonal and interannual patterns of C fluxes during typhoon disturbances. River intrusion controls not only the C fluxes, algal biomass, and nutrient loading but also influences the length of stratification and hydraulic retention times (Lin et al., 2021, 2022; Maranger et al., 2018; Nakayama et al., 2020; Olsson et al., 2022a, b; Zwart et al., 2017; Vachon and Del Giorgio, 2014). Therefore, we hypothesized that allochthonous C loading and river inflow intrusion might affect DIC and DOC distributions. Further, autochthonous processes in small subtropical lakes, such as algal photosynthesis, remineralization, and vertical transportation, must also be considered. Here, we tested our hypothesis by developing two-layer conceptual C models to assess C flux responses to typhoon disturbances in small subtropical lakes.

Figure 1 Sampling locations and bathymetry maps of Yuan-Yang Lake (YYL). The dark green rectangle shows the buoy station, which is located at the deepest site of the lake. The red arrows and white cross show the river mouths of the inflows and outflows, respectively. The gray triangle shows the location of the water depth sensor.

2 Materials and methods

2.1 Study site

YYL is a shallow (mean water depth: 4.3 m), oligotrophic (total phosphorous: 10–20 µg P L⁻¹; total nitrogen: 20–60 µg N L⁻¹) subtropical mountain lake (Chou et al., 2000; Tsai et al., 2008; Wu et al., 2001) on Chilan Mountain at around 1640 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ in northern central Taiwan (24.58^∘ N, 121.40^∘ E) (Fig. 1). Its water is brown due to its humic acid content (colored dissolved organic matter: 20–50 ppb quinine sulfate equivalent; specific ultraviolet absorbance at 254 nm assessed by a portable fluorometer – model C3, Turner Designs, Sunnyvale, CA, USA; mean pH: 5.4). YYL is surrounded by old-growth trees such as Chamaecyparis formosensis, Chamaecyparis obtusa var. formosana, and Rhododendron formosanum Heiml (Chou et al., 2000). Precipitation is over 3000 mm yr⁻¹, and typhoon precipitation contributes up to half of the total precipitation input to YYL annually (Chang et al., 2007; Lai et al., 2006). Due to the rapid renewal of the waterbody, the water retention time (or residence time) was around 4.4 d during Typhoon Megi from 27 September to 1 October 2016 (Lin et al., 2022). The water surface temperature ranges from 15 to 25 ^∘C during March to August, and the water column overturns in September (Kimura et al., 2012, 2017; Lin et al., 2021). The concentrations of DIC, DOC (Lin et al., 2021), total nitrogen, and total phosphate (Chiu et al., 2020; Tsai et al., 2008) as well as bacterial compositions (Shade et al., 2011) increase within YYL from autumn to winter. YYL has been registered as a long-term ecological study site by the Ministry of Science and Technology (MOST) of Taiwan since 1992 and became part of the Global Lake Ecological Observatory Network (GLEON) in 2004.

2.2 Water sampling and chemical analysis

We collected water quality samples (DOC, DIC, and Chl a) at water depths of 0.04, 0.50, 1.00, 2.00, and 3.50 m at the buoy site (Fig. 1). From January 2015 to December 2018, we measured the water surfaces of six river inflows and one outflow each month using a horizontal van Dorn bottle (2.20 L, acrylic) (Fig. 1). These liquid samples were collected using a portable hand pump and glass microfiber filter papers (47 mm GF/F, nominal pore size 0.70 µm; Whatman, Maidstone, Kent, UK) to obtain filtrate samples. Water samples were stored at around 4 ^∘C in a refrigerator until analysis. Samples were analyzed using an infrared gas detector to detect DIC and DOC concentrations with persulfate digestion (model 1088 Rotary TOC autosampler; OI Analytical, College Station, TX, USA). The filter papers were kept refrigerated in opaque bottles at around −25 ^∘C in a refrigerator until the samples were analyzed. In the laboratory, the filter papers were extracted with methanol to obtain the Chl a concentration using a portable fluorometer (model 10-AU-005-CE; Turner Designs, Sunnyvale, CA, USA), with specific wavelengths of 430 nm (blue) and 662 nm (red). All analysis was completed within 72 h of exposure to light in order to reduce degradation.

Table 1 Comparison of Yuan-Yang Lake's rainfall and hydrological records in typhoon and non-typhoon years.

SD denotes standard deviation.

2.3 Data analysis and numerical modeling

Three water quality variables (DIC, DOC, and Chl a) were compared between different layers (upper and lower layers), years (typhoon and non-typhoon years), and seasons (spring, summer, autumn, and winter). First, we separated our investigation data into typhoon years and non-typhoon years, as described in Sect. 2.3.1. Next, to simulate the DIC and DOC concentrations under extreme weather scenarios in YYL, we developed a conceptual equations model to generate continuous DIC and DOC data for the upper and lower layers, as shown in Sect. 2.3.2. This also helped us understand the transportation, photosynthesis, and remineralization rates between seasons and between typhoon and non-typhoon years.

2.3.1 Typhoon and non-typhoon years

We collected meteorological data from a tower located about 1.0 km from YYL (Lin et al., 2021, 2022). Data on rainfall (model N-68; Nippon Electric Instrument, Tokyo, Japan) and wind speed (model 03001; R. M. Young, Traverse City, MI, USA) were stored in a datalogger (model CR1000; Campbell Scientific, Logan, UT, USA) every 10 min. River discharge (Q_in, m³ d⁻¹) was estimated every 10 min using the rainfall data and a water depth meter (model HOBO U20; Onset Computer, Bourne, MA, USA) at the end of a river inflow (Fig. 1) using the Manning formula. Transparency was estimated using Secchi disk data measured at a certain interval in the time frame from 10:00 to 14:00 LT (local time, GMT + 08:00).

Four strong typhoons were recorded using wind speed and rainfall meteorological parameters. As Table 1 shows, these events contributed a total of 2254 mm of precipitation in all 24 months of 2015 and 2016; this accounted for 35.6 % of the typical rainfall across 2 years (> 3000 mm yr⁻¹). However, no typhoon rainfall was recorded at YYL in 2017 or 2018; the total precipitation was around 1398 mm yr⁻¹, below half of average years. The annual average wind speed from 2015 to 2016 was 0.09 m s⁻¹ higher than in 2017 and 2018 (Table 1). Despite no significant difference in the average water depth between 2017 and 2018, the average discharge in 2015 and 2016 was 774 m³ d⁻¹ higher than in 2017 and 2018 (Table 1). Thus, we considered 2015 and 2016 to be typhoon years and 2017 and 2018 to be non-typhoon years.

2.3.2 Conceptual two-layer DIC and DOC model

Nakayama et al. (2010) successfully developed a conceptual two-layer dissolved oxygen model based on strong wind turbulence at Tokyo Bay. Lin et al. (2021) pointed out that thermal stratification that inhibits vertical C flux between the upper and lower layers in shallow stratified lakes makes it possible to develop conceptual two-layer C models (Lin et al., 2022; Nakayama et al., 2022). The phytoplankton and remineralization effects on DIC and DOC fluxes ($\mathrm{d}\mathrm{DIC}/\mathrm{d}t$ and $\mathrm{d}\mathrm{DOC}/\mathrm{d}t$, respectively, $\mathrm{mg}\mathrm{C}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$) were considered in a conceptual two-layer equation model as shown in Eqs. (1)–(4). The fluxes in the upper layer (from the water surface to 2.5 m water depth) were calculated as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{rl}{V}_{\mathrm{U}}{\displaystyle \frac{\mathrm{d}{\mathrm{DIC}}_{\mathrm{U}}}{\mathrm{d}t}}=& O_{\mathrm{U}}{\mathrm{DIC}}_{\mathrm{R}}-Q_{\text{out}}{\mathrm{DIC}}_{\mathrm{U}}-V_{\mathrm{U}}{\mathit{\alpha }}_{\text{PU}}{\mathrm{Chl}}_{\mathrm{U}}\\ & +V_{\mathrm{U}}{\mathit{\alpha }}_{\text{MU}}{\mathrm{DOC}}_{\mathrm{U}}+A_{\mathrm{I}}{w}_{\mathrm{I}}({\mathrm{DIC}}_{\mathrm{L}}-{\mathrm{DIC}}_{\mathrm{U}})\\ & +Q_{\mathrm{L}}{\mathrm{DIC}}_{\mathrm{L}}-{\displaystyle \frac{A_{\mathrm{s}}{F}_{{\mathrm{CO}}_{\mathrm{2}}}}{C_{\mathrm{U}}}}+{\mathrm{Pa}}_{\mathrm{U}},\end{array}\text{(2)}& {\displaystyle }\begin{array}{rl}{V}_{\mathrm{U}}{\displaystyle \frac{\mathrm{d}{\mathrm{DOC}}_{\mathrm{U}}}{\mathrm{d}t}}=& O_{\mathrm{U}}{\mathrm{DOC}}_{\mathrm{R}}-Q_{\text{out}}{\mathrm{DOC}}_{\mathrm{U}}\\ & -V_{\mathrm{U}}{\mathit{\alpha }}_{\text{MU}}{\mathrm{DOC}}_{\mathrm{U}}\\ & +A_{\mathrm{I}}{w}_{\mathrm{I}}({\mathrm{DOC}}_{\mathrm{L}}-{\mathrm{DOC}}_{\mathrm{U}})\\ & +Q_{\mathrm{L}}{\mathrm{DOC}}_{\mathrm{L}}+{\mathrm{Pb}}_{\mathrm{U}}.\end{array}\end{array}$$

Table 2 Parameters of the two-layer conceptual model in Yuan-Yang Lake.

Those in the lower layer (from 2.5 to 4.0 m water depth) were calculated as follows:

$$\begin{array}{}\text{(3)}& {\displaystyle }\begin{array}{rl}{V}_{\mathrm{L}}{\displaystyle \frac{\mathrm{d}{\mathrm{DIC}}_{\mathrm{L}}}{\mathrm{d}t}}=& Q_{\mathrm{L}}{\mathrm{DIC}}_{\mathrm{R}}-V_{\mathrm{L}}{\mathit{\alpha }}_{\text{PL}}{\mathrm{Chl}}_{\mathrm{L}}+V_{\mathrm{L}}{\mathit{\alpha }}_{\text{ML}}{\mathrm{DOC}}_{\mathrm{L}}\\ & +A_{\mathrm{I}}{w}_{\mathrm{I}}({\mathrm{DIC}}_{\mathrm{U}}-{\mathrm{DIC}}_{\mathrm{L}})-Q_{\mathrm{L}}{\mathrm{DIC}}_{\mathrm{L}}\\ & +{\displaystyle \frac{A_{\mathrm{B}}{\mathrm{BF}}_{\mathrm{DIC}}}{C_{\mathrm{U}}}}+{\mathrm{Pa}}_{\mathrm{L}},\end{array}\text{(4)}& {\displaystyle }\begin{array}{rl}{V}_{\mathrm{L}}{\displaystyle \frac{\mathrm{d}{\mathrm{DOC}}_{\mathrm{L}}}{\mathrm{d}t}}=& Q_{\mathrm{L}}{\mathrm{DOC}}_{\mathrm{R}}-V_{\mathrm{L}}{\mathit{\alpha }}_{\text{ML}}{\mathrm{DOC}}_{\mathrm{L}}\\ & +A_{\mathrm{I}}{w}_{\mathrm{I}}({\mathrm{DOC}}_{\mathrm{U}}-{\mathrm{DOC}}_{\mathrm{L}})-Q_{\mathrm{L}}{\mathrm{DOC}}_{\mathrm{L}}\\ & +{\mathrm{Pb}}_{\mathrm{L}}.\end{array}\text{(5)}& {\displaystyle }{V}_{\text{total}}=V_{\mathrm{U}}+V_{\mathrm{L}}\text{(6)}& {\displaystyle }{Q}_{\text{in}}=O_{\mathrm{U}}+Q_{\mathrm{L}}\end{array}$$

In the above expressions, as shown in Table 2, the total lake volume (V_total, 53 544 m³) comprises the upper layer (V_U, 45 456 m³) and the lower layer (V_L, 8808 m³) (Eq. 5), the lake surface area (A_s) is 36 000 m², and the bottom of the lake area (A_B) is 3520 m². The interface is 2.5 m vertically, and the interface area (A_I) is 7264 m² in YYL. The water depth varied from 4.56 to 4.66 m during the typhoon period (Chiu et al., 2020; Lin et al., 2022). Therefore, we can assume that changes in the lake volume and area were negligible. The coefficient C_U, with a value of 1000, was used to establish a standard unit for $F_{{\mathrm{CO}}_{\mathrm{2}}}$ ($\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), considering the air–water CO₂ exchange using Fick's law, as follows:

$$\begin{array}{}\text{(7)}& F_{{\mathrm{CO}}_{\mathrm{2}}}=k_{{\mathrm{CO}}_{\mathrm{2}}}\cdot K_{H}p{{\mathrm{CO}}_{\mathrm{2}}}_{\mathrm{water}}-p{{\mathrm{CO}}_{\mathrm{2}}}_{\mathrm{air}}.\end{array}$$

Here, $k_{{\mathrm{CO}}_{\mathrm{2}}}$ is the gas transfer velocity from empirical wind speed equations (Cole and Caraco, 1998; Jähne et al., 1987; Smith, 1985; Wanninkhof, 1992); K_H is Henry's coefficient calculated using water temperature empirical equations (Plummer and Busenberg, 1982); $p{{\mathrm{CO}}_{\mathrm{2}}}_{\mathrm{air}}$ (µatm) is the CO₂ partial pressure in the atmosphere using air pressure data (Lin et al., 2021, 2022), where the atmospheric CO₂ concentration is assumed to be 400 ppm; and $p{{\mathrm{CO}}_{\mathrm{2}}}_{\mathrm{water}}$ (µatm) is the CO₂ partial pressure at the water surface at around 0.04 m water depth from water quality data (temperature, pH, and DIC concentrations). The empirical equation (Cai and Wang, 1998) was also used by Lin et al. (2021). $F_{{\mathrm{CO}}_{\mathrm{2}}}$ contributed approximately half of the net ecosystem production (NEP) across the water surface to the atmosphere for YYL (Lin et al., 2021). Further, because sediment carbon may be an important flux into shallow subtropical lakes, the sediment C flux (BF_DIC or BF_DOC, mg C L⁻¹) in the lower layer was considered (Lin et al., 2022).

We assumed that the river discharge and outflow discharge (Q_out, m³ d⁻¹) were in a quasi-steady state (Q_in=Q_out), divided into upper-layer discharge (O_U, m³ d⁻¹) and lower-layer discharge (Q_L, m³ d⁻¹) (Eq. 6). Lin et al. (2021) showed that the buoyancy frequencies in YYL were 0.011 ± 0.004 s⁻¹, 0.013 ± 0.004 s⁻¹, 0.006 ± 0.003 s⁻¹, and 0.007 ± 0.004 s⁻¹ for spring (March–May), summer (June–August), autumn (September–November), and winter (December–February), respectively, inhibiting the vertical profile of DIC mixing due to stratification. We estimated the percentages of O_U and Q_L based on the buoyancy frequency following Lin et al. (2021, 2022). O_U values were 75 %, 80 %, 45 %, and 50 % of Q_in for spring, summer, autumn, and winter, respectively, and Q_L values were 25 %, 20 %, 55 %, and 50 % of Q_in for the corresponding respective seasons (Table 2).

Extreme weather events might induce stronger seasonal thermal stratification from spring to summer and longer overturns from autumn to winter, thereby changing C distribution and transportation within waterbodies (Kraemer et al., 2021; Olsson et al., 2022a; Woolway et al., 2020). Thus, we attempted to simulate extreme climate scenarios; we shifted the ratio of Q_in for each season and tested the river intrusion hypothesis. We established two extreme conditions: “Level 1” and “Level 2”. Level 2 is the more extreme condition: O_U is 80 % (spring), 85 % (summer), 50 % (autumn), and 50 % (winter) of Q_in; Q_L is 20 % (spring), 15 % (summer), 50 % (autumn), and 50 % (winter) of Q_in. Level 1 is the condition between the present and the Level 2 condition: O_U is 77 % (spring), 82 % (summer), 47 % (autumn), and 50 % (winter) of Q_in; Q_L is 23 % (spring), 18 % (summer), 53 % (autumn), and 50 % (winter) of Q_in (Table 2).

The contributions of photosynthesis production depended on the chlorophyll a (Chl a) concentration (Chl_U or Chl_L, mg L⁻¹) and on the absorption coefficients in the upper layer (α_PU, d⁻¹) and lower layer (α_PL, d⁻¹). The coefficients of DOC remineralization rates in the upper layer (α_MU, d⁻¹) and lower layer (α_ML, d⁻¹) were also considered in the conceptual models. Pa_U, Pa_L, Pb_U, and Pb_L are constants in the conceptual models. To obtain unknown values (α_PU, α_MU, α_PL, α_MU, α_ML, w_I, BF_DIC, BF_DOC, Pa_U, Pa_L, Pb_U, and Pb_L), we applied multiple linear regression analysis. Further, these unknown values were tested by trial and error to obtain the parameters of the “Best-fit” condition (Nakayama et al., 2022). The same parameters of the Best-fit condition were used to obtain the extreme conditions for Level 1 and Level 2. We used the coefficient of determination (R²) and the Nash–Sutcliffe efficiency (NSE) coefficient (Nash and Sutcliffe, 1970) to quantify the performance of the equation model with DIC and DOC sampling data (observation data) for each simulation as follows:

$$\begin{array}{}\text{(8)}& \mathrm{NSE}=\mathrm{1}-{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n({\text{Obs}}_i-{\text{Sim}}_i{)}^{\mathrm{2}}}{{\sum }_{i=\mathrm{1}}^n({\text{Obs}}_i-\stackrel{\mathrm{‾}}{\text{Obs}}{)}^{\mathrm{2}}}},\end{array}$$

where Obs represents observation data of DIC and DOC concentrations, and Sim represents Best-fit scenario data for the conceptual model.

2.3.3 DIC and DOC fluxes

Net ecosystem production (NEP) was defined as the difference between gross primary production (GPP) and ecological respiration (ER), NEP = GPP − ER, due to photosynthesis and respiration via biota (Dodds and Whiles, 2020). Given that we assumed that the C fluxes were dependent on the river inflows in YYL (Fig. 1), we estimated the NEP via end-member analysis using the C concentration of the river inflow and outflow (Lin et al., 2021; Nakayama et al., 2020), following Eqs. 9–12. The upper-layer NEP of the DIC flux (mg C d⁻¹) was obtained from Eq. (1) as follows:

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& \mathrm{Upper}{\mathrm{flux}}_{\mathrm{DIC}}=C_{\mathrm{U}}{\mathit{\alpha }}_{\text{PU}}{\mathrm{Chl}}_{\mathrm{U}}-C_{\mathrm{U}}{\mathit{\alpha }}_{\text{MU}}{\mathrm{DOC}}_{\mathrm{U}}\\ & -C_{\mathrm{U}}{\displaystyle \frac{A_{\mathrm{I}}{w}_{\mathrm{I}}({\mathrm{DIC}}_{\mathrm{L}}-{\mathrm{DIC}}_{\mathrm{U}})}{V_{\mathrm{U}}}}-C_{\mathrm{U}}{\displaystyle \frac{{\mathrm{Pa}}_{\mathrm{U}}}{V_{\mathrm{U}}}}\\ =& C_{\mathrm{U}}{\displaystyle \frac{O_{\mathrm{U}}{\mathrm{DIC}}_{\mathrm{R}}+Q_{\mathrm{L}}{\mathrm{DIC}}_{\mathrm{L}}-Q_{\text{out}}{\mathrm{DIC}}_{\mathrm{U}}}{V_{\mathrm{U}}}}\\ & -{\displaystyle \frac{A_{\mathrm{S}}}{V_{\mathrm{U}}}}{F}_{{\mathrm{CO}}_{\mathrm{2}}}\\ =& C_{\mathrm{U}}{\displaystyle \frac{\mathrm{1}}{t_{\text{rU}}}}\left({\displaystyle \frac{O_{\mathrm{U}}}{Q_{\text{in}}}}{\mathrm{DIC}}_{\mathrm{R}}+{\displaystyle \frac{Q_{\mathrm{L}}}{Q_{\text{in}}}}{\mathrm{DIC}}_{\mathrm{L}}-{\mathrm{DIC}}_{\mathrm{U}}\right)-F_{\mathrm{C}}\\ & t_{\text{rU}}={\displaystyle \frac{V_{\mathrm{U}}}{Q_{\text{in}}}}.\end{array}\end{array}$$

The upper-layer flux of the DOC flux ($\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$) was estimated from Eq. (2):

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& \mathrm{Upper}{\mathrm{flux}}_{\mathrm{DOC}}=C_{\mathrm{U}}{\mathit{\alpha }}_{\text{MU}}{\mathrm{DOC}}_{\mathrm{U}}\\ & -C_{\mathrm{U}}{\displaystyle \frac{A_{\mathrm{I}}{w}_{\mathrm{I}}({\mathrm{DOC}}_{\mathrm{L}}-{\mathrm{DOC}}_{\mathrm{U}})}{V_{\mathrm{U}}}}-C_{\mathrm{U}}{\displaystyle \frac{{\mathrm{Pb}}_{\mathrm{U}}}{V_{\mathrm{U}}}}\\ & =C_{\mathrm{U}}{\displaystyle \frac{O_{\mathrm{U}}{\mathrm{DOC}}_{\mathrm{R}}+Q_{\mathrm{L}}{\mathrm{DOC}}_{\mathrm{L}}-Q_{\text{out}}{\mathrm{DOC}}_{\mathrm{U}}}{V_{\mathrm{U}}}}\\ & =C_{\mathrm{U}}{\displaystyle \frac{\mathrm{1}}{t_{\text{rU}}}}\left({\displaystyle \frac{O_{\mathrm{U}}}{Q_{\text{in}}}}{\mathrm{DOC}}_{\mathrm{R}}+{\displaystyle \frac{Q_{\mathrm{L}}}{Q_{\text{in}}}}{\mathrm{DOC}}_{\mathrm{L}}-{\mathrm{DOC}}_{\mathrm{U}}\right).\end{array}\end{array}$$

The lower-layer flux of DIC flux ($\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$) was estimated from Eq. (3):

$$\begin{array}{}\text{(11)}& \begin{array}{rl}\mathrm{Lower}{\mathrm{flux}}_{\mathrm{DIC}}=& C_{\mathrm{U}}{\mathit{\alpha }}_{\text{PL}}{\mathrm{Chl}}_{\mathrm{L}}-C_{\mathrm{U}}{\mathit{\alpha }}_{\text{ML}}{\mathrm{DOC}}_{\mathrm{L}}\\ & -C_{\mathrm{U}}{\displaystyle \frac{A_{\mathrm{I}}{w}_{\mathrm{I}}({\mathrm{DIC}}_{\mathrm{U}}-{\mathrm{DIC}}_{\mathrm{L}})}{V_{\mathrm{L}}}}\\ & -{\displaystyle \frac{A_{\mathrm{B}}{\mathrm{BF}}_{\mathrm{DIC}}}{V_{\mathrm{L}}}}-C_{\mathrm{U}}{\displaystyle \frac{{\mathrm{Pa}}_{\mathrm{L}}}{V_{\mathrm{L}}}}\\ =& C_{\mathrm{U}}{\displaystyle \frac{Q_{\mathrm{L}}({\mathrm{DIC}}_{\mathrm{R}}-{\mathrm{DIC}}_{\mathrm{L}})}{V_{\mathrm{L}}}}\\ =& C_{\mathrm{U}}{\displaystyle \frac{\mathrm{1}}{t_{\text{rL}}}}{\displaystyle \frac{Q_{\mathrm{L}}}{Q_{\text{in}}}}({\mathrm{DIC}}_{\mathrm{R}}-{\mathrm{DIC}}_{\mathrm{L}})\\ & t_{\text{rL}}={\displaystyle \frac{V_{\mathrm{L}}}{Q_{\text{in}}}}.\end{array}\end{array}$$

The lower-layer flux of DOC flux ($\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$) was estimated from Eq. (4):

$$\begin{array}{}\text{(12)}& \begin{array}{rl}\mathrm{Lower}{\mathrm{flux}}_{\mathrm{DOC}}=& C_{\mathrm{U}}{\mathit{\alpha }}_{\text{ML}}{\mathrm{DOC}}_{\mathrm{L}}\\ & -C_{\mathrm{U}}{\displaystyle \frac{A_{\mathrm{I}}{w}_{\mathrm{I}}({\mathrm{DOC}}_{\mathrm{U}}-{\mathrm{DOC}}_{\mathrm{L}})}{V_{\mathrm{L}}}}\\ & -{\displaystyle \frac{A_{\mathrm{B}}{\mathrm{BF}}_{\mathrm{DOC}}}{V_{\mathrm{L}}}}-C_{\mathrm{U}}{\displaystyle \frac{{\mathrm{Pb}}_{\mathrm{L}}}{V_{\mathrm{L}}}}\\ =& C_{\mathrm{U}}{\displaystyle \frac{Q_{\mathrm{L}}({\mathrm{DOC}}_{\mathrm{R}}-{\mathrm{DOC}}_{\mathrm{L}})}{V_{\mathrm{L}}}}\\ =& C_{\mathrm{U}}{\displaystyle \frac{\mathrm{1}}{t_{\text{rL}}}}{\displaystyle \frac{Q_{\mathrm{L}}}{Q_{\text{in}}}}({\mathrm{DOC}}_{\mathrm{R}}-{\mathrm{DOC}}_{\mathrm{L}}).\end{array}\end{array}$$

Thus, the respective total fluxes of DIC and of DOC are as follows:

$$\begin{array}{}\text{(13)}& {\displaystyle }{\mathrm{Flux}}_{\mathrm{DIC}}={\displaystyle \frac{V_{\mathrm{U}}\mathrm{Upper}{\mathrm{flux}}_{\mathrm{DIC}}+V_{\mathrm{L}}\mathrm{Lower}{\mathrm{flux}}_{\mathrm{DIC}}}{V_{\text{total}}}},\text{(14)}& {\displaystyle }{\mathrm{Flux}}_{\mathrm{DOC}}={\displaystyle \frac{V_{\mathrm{U}}\mathrm{Upper}{\mathrm{flux}}_{\mathrm{DOC}}+V_{\mathrm{L}}\mathrm{Lower}{\mathrm{flux}}_{\mathrm{DOC}}}{V_{\text{total}}}}.\end{array}$$

Here, F_C is $\frac{A_{\mathrm{S}}}{V_{\mathrm{U}}}{F}_{{\mathrm{CO}}_{\mathrm{2}}}$, and t_rU and t_rL are residence times (d)T in the upper and lower layers, respectively. These parameters were used for the Best-fit condition, as shown in Table 2.

3 Results

3.1 Measurement data (monthly DIC, DOC, and Chl a concentrations) in typhoon and non-typhoon years

The comparisons between the two typhoon years (2015 and 2016) revealed no significant differences in DIC, DOC, or Chl a concentrations between the upper and lower layers; however, all of these parameters differed significantly between the layers in the non-typhoon years of 2017 and 2018 (Fig. 2). This was due to typhoon-induced mixing and lower thermal stratification between the upper and lower layer (Lin et al., 2021, 2022). Overall, the average DIC_U and DIC_L were 2.06 and 3.66 mg C L⁻¹, respectively; the average DOC_U and DOC_L were 5.87 and 8.02 mg C L⁻¹, respectively; and Chl_U and Chl_L were 2.13 and 18.5 µg C L⁻¹, respectively. In typhoon years, the average DIC_U and DIC_L were 2.34 and 4.07 mg C L⁻¹, respectively; the average DOC_U and DOC_L were 6.10 and 8.38 mg C L⁻¹, respectively; and the Chl_U and Chl_L were 2.38 and 12.2 µg C L⁻¹, respectively (Fig. 2). In non-typhoon years, the average DIC_U and DIC_L were 1.81 and 3.28 mg C L⁻¹, respectively; the average DOC_U and DOC_L were 5.66 and 7.67 mg C L⁻¹, respectively; and Chl_U and Chl_L were 1.89 and 24.4 µg C L⁻¹, respectively (Fig. 2).

Figure 2 Comparisons of (a) DIC, (b) DOC, and (c) Chl a between upper (DIC_U, DOC_U, and Chl_U) and lower (DIC_L, DOC_L, and Chl_L) layers, grouped by typhoon and non-typhoon years. Bullet points show the water sampling data. We used a t test to obtain the p values (blue text). Total sampling numbers are 41 (n=41) for each measurement from January 2015 to December 2018; n=20 in typhoon years and n=21 in non-typhoon years.

Figure 3 Seasonal variations in (a) upper-layer DIC (DIC_U), (b) lower-layer DIC (DIC_L), (c) upper-layer DOC (DOC_U), (d) lower-layer DOC (DOC_L), (e) upper-layer Chl a (Chl_U), and (f) lower-layer Chl a (Chl_L), grouped by typhoon and non-typhoon years. The bullet points show the water sampling data. To determine seasonality, we used a one-way ANOVA to obtain the p values. The significance of the p values is defined as follows: “ns” (not significant) represents p values ≥ 0.05, “^∗” shows p values from 0.05 to 0.01, and “${}^{\ast \ast }$” represents p values from 0.01 to 0.001.

ANOVA results indicated no significant differences in DIC concentrations among seasons during typhoon years (p values ≥ 0.05), suggesting a lack of statistically significant variation in DIC data across seasons (Fig. 3a, b). However, the DOC concentration showed significant differences between seasons in typhoon years (Fig. 3c, d). No significant differences between Chl_U and Chl_L were observed among the seasons (Fig. 3e, f), whereas the standard deviation (SD) values of DIC and DOC were higher in summer and autumn (Fig. 3) due to terrestrial C loading (Chiu et al., 2020). In summer, the SD values of DIC_U and DOC_U were 3.51 mg C L⁻¹ and 3.69 mg C L⁻¹, respectively (Fig. 3a, c, e). In autumn, DIC_L and DOC_L had the highest SD values (4.06 and 4.17 mg C L⁻¹, respectively) (Fig. 3b, d). Notably, the maximums of DIC_U and DOC_U were 7.06 and 15.6 mg C L⁻¹, respectively, and those of DIC_L and DOC_L were 10.9 and 19.8 mg C L⁻¹, respectively, in typhoon years (Fig. 3a, b, c, d).

Figure 4 Pearson correlation coefficients of the DIC, DOC, and Chl a concentrations for upper-layer and lower-layer DIC (DIC_U and DIC_L, respectively), upper-layer and lower-layer DOC (DOC_U and DOC_L, respectively), and upper-layer and lower-layer Chl a (Chl_U and Chl_L, respectively) during (a) typhoon years and (b) non-typhoon years. Black-crosses show insignificant values (p values are > 0.05).

Positive Pearson correlations of 0.45 to 0.80 were observed between the DOC and DIC in typhoon years (Fig. 4a). In the non-typhoon years, the upper-layer DIC was the only variable that was negatively correlated with DOC in the upper and lower layers (Fig. 4b). DIC in the lower layer was positively correlated with the Chl_L (Fig. 4) due to the abundant respiration in the lower layer (Lin et al., 2021; Tsai et al., 2008).

Table 3 Best-fit parameters of a two-layer conceptual model of DIC and DOC for Yuan-Yang Lake from 2015 to 2018.

Figure 5 Continuous daily DIC and DOC data for the (a, c) upper layer (DIC_U and DOC_U, respectively) and (b, d) lower layer (DIC_L and DOC_L, respectively) using the conceptual equation model under extreme climates from 2015 to 2018. Blue lines are the original Best-fit scenario data, for which the parameters of the DIC model in non-typhoon years are as shown in Table 3, and gray dots show water sampling (observation) data for each month (n=41) from January 2015 to December 2018. Orange regions show the Level 1 scenario and pink regions show the Level 2 scenario.

3.2 Performance of simulation data in the conceptual two-layer DIC and DOC models

The results for typhoon years demonstrated that DIC_U was around 1.5–5.0 mg C L⁻¹ (Fig. 5a, b) and DIC_L was around 5.0 mg C L⁻¹ (Fig. 5d). In the non-typhoon years (2017–2018), the NSE values of DIC_U and DIC_L were 0.61 and 0.70, respectively. On the other hand, the DOC values fit our observation data (R² values are 0.91 and 0.46; the NSE coefficients from Eq. 8 are 0.95 and 0.73) (Fig. 5c, d; Table 3). The parameters for the conceptual two-layer DIC and DOC models showed different regimes between typhoon and non-typhoon years (Table 3). In typhoon years, the photosynthesis absorption rate coefficients (α_PU and α_PL) were negative (photosynthesis < respiration) for each season. YYL was a C source due to a large allochthonous C loading during typhoons; the respiration was elevated by around 30- to 150-fold from summer to autumn. However, the values of the transportation coefficients (w_I) were higher in autumn than in the other seasons (Table 3). Further, the higher remineralization rates during typhoon disturbances from summer to autumn resulted in positive α_MU and α_ML. In non-typhoon years, the remineralization rates were negative (Table 3).

Figure 6 Interannual (a) Upper flux_DIC, (b) Lower flux_DIC, (c) Flux_DIC, (d) Upper flux_DOC, (e) Lower flux_DOC, and (f) Flux_DOC values ($\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$) grouped by typhoon and non-typhoon years. The Obs condition (black boxes) shows the observation data, as in Fig. 5. The Fit condition (blue boxes) shows the Best-fit scenario data using the conceptual two-layer carbon model, as in Fig.5. Level 1 (yellow boxes) and Level 2 (red boxes) show the extreme scenarios, as in Fig. 6. For the definitions of fluxes, the reader is referred to Sect. 2.3.3. Positive values refer to a carbon sink, whereas negative values show carbon released within YYL.

We simulated the responses of DIC and DOC to typhoons using conceptual two-layer C models. The results showed that the DIC was more sensitive to typhoon disturbances than DOC under the Level 1 and Level 2 scenarios (Fig. 5). Overall, the C level declined in the upper layers but increased in the lower layers (Fig. 5). DIC and DOC in the upper layer tended to decline from 1.0 mg C L⁻¹ (Level 1) to 2.0 mg C L⁻¹ (Level 2) (Fig. 5a, c); however, they increased to 10.0 and 20.0 mg C L⁻¹ in the lower layer under the Level 1 and Level 2 scenarios, respectively (Fig. 5b, d). Under extreme weather conditions, Level 2 usually shifted to different typhoon responses for each season (Figs. S1 and S2 in the Supplement) due to extreme river intrusions and the strength of the thermal stratification (Lin et al., 2021). DIC changed more than DOC under the Level 1 and Level 2 scenarios (Fig. 5) because photosynthesis, transportation, and remineralization rates may crucially affect the seasonal and interannual patterns of DOC as well.

Table 4 Seasonal carbon flux averages ($\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$) for each season in Yuan-Yang Lake. Positive values represent a carbon sink, whereas negative values show that carbon was released. F_C denotes carbon emission across water to air using the empirical equations method (Lin et al., 2021). The bold font represents average values from spring to winter.

3.3 Interannual responses of the NEP to typhoons under extreme weather scenarios

We used field observation data (Fig. 5) to estimate the C fluxes in Table 4. The typhoon disturbances in summer and autumn played an important role in promoting C release by YYL (Table 4). Overall, YYL released 245 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$ of DIC and 415 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$ of DOC during the typhoon years; during non-typhoon years, it released 51.7 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$ of DIC and 22.8 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$ of DOC (Table 4). The average F_C was 219 and 133 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$ released from YYL into the atmosphere in typhoon and non-typhoon years, respectively, which was 1–2 times larger than Flux_DIC (Table 4). In summer, the upper layer exhibited declines in both the DIC and DOC concentrations, with DIC being approximately 3.7 times higher in typhoon years than in non-typhoon years (Table 4). In autumn, 216 mg C d⁻¹ of upper-layer DIC was released in typhoon years, but 46.1 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$ of upper-layer DOC was produced. The upper-layer Flux_DIC was negative in autumn in typhoon years, as 268 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$ more F_C was released than in non-typhoon years. In addition, the lower layer exhibited the largest release of C into the outflow in typhoon years; however, the flux in the lower layer was more than twice as high in the summer compared with the autumn of those years (Table 4). The average total Flux_DIC was a release of approximately 3.14 times more C in typhoon years than in non-typhoon years. The average total Flux_DOC showed an increase of 62.3 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$ of DOC between typhoon and non-typhoon years due to the over 10-fold higher flux in the upper layer (Table 4).

We compared the fluxes with different model conditions (Best-fit, Level 1, and Level 2), as shown in Fig. 6, demonstrating that the responses of Flux_DIC to typhoons differed dramatically between the Level 1 and Level 2 scenarios (Fig. 6a, b, c); in particular, the Upper Flux_DIC released more C in typhoon years and absorbed more C in non-typhoon years than Obs (Fig. 6a). Not only were the absolute values of Flux_DIC over 3 times higher in typhoon years than in non-typhoon years (Table 4), but the SD was higher in typhoon years as well (Fig. 6). However, DOC fluxes changed less under the Level 1 and Level 2 scenarios (Fig. 6d, e, f), a finding that is consistent with our continuous DOC data (Fig. 5c, d).

Figure 7 Contributions of typhoon disturbances to (a) DIC and (b) DOC fluxes ($\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$). The contributions of typhoons to C fluxes (ΔFlux_DIC and ΔFlux_DOC) represented the intervals of Flux_DIC (or Flux_DOC) in typhoon years and non-typhoon years (e.g., ΔFlux_DIC = Flux_DIC in typhoon years − Flux_DIC in non-typhoon years).

We not only attempted to establish the contributions of thermal stratification and river intrusion to DIC and DOC fluxes using the Level 1 and Level 2 scenarios in this conceptual model but we also examined the contributions of typhoon disturbances to DIC and DOC fluxes (Fig. 7). We found that DIC and DOC fluxes were both released (ΔFlux_DIC and ΔFlux_DOC, respectively) under the Level 1 and Level 2 scenarios (Fig. 7). Typhoon disturbances contributed −102.2 and −62.3 mg C m³ d⁻¹ with respect to Flux_DIC and Flux_DOC in measurement (observation) data (Obs), respectively (Fig. 7). The averages of ΔFlux_DIC under Best-fit conditions decreased by 7.0 % (Level 1) and 30.0 % (Level 2); ΔFlux_DOC decreased by 56.9 % (Level 1) and 118 % (Level 2), respectively (Fig. 7).

4 Discussion

4.1 Biochemical and physical differences in DIC and DOC fluxes between typhoon and non-typhoon years in YYL

The total precipitation was 35.6 % higher in typhoon years than in non-typhoon years (Table 1). Water retention and typhoon-induced upwelling controlled the dynamics of DIC and DOC during the summer and autumn (Chiu et al., 2020; Jones et al., 2009; Tsai et al., 2008, 2011). The absence of typhoon-induced upwelling affected water quality data differences between the upper and lower layers (Chiu et al., 2020; Lin et al., 2022; Tsai et al., 2008, 2011). DIC, DOC, and Chl a concentrations differed significantly between the upper and lower layers in non-typhoon years (Fig. 2). Further, the abundance of microorganisms can lead to intensive respiration in the lower layers during the non-typhoon period in YYL; for example, an anoxic condition at the hypolimnion may decrease the efficiency of C mineralization and remineralization in non-typhoon years (Carey et al., 2022; Chiu et al., 2020; Lin et al., 2022; Shade et al., 2010; Shade et al., 2011). Thus, the thermal stratification and anoxic condition may have been controlled by the seasonal and interannual patterns of DIC and DOC fluxes in non-typhoon years (Fig. 5; Tables 3, 4).

Figure 8 Interactions of measurement data in typhoon years and non-typhoon years. The black circles represent typhoon years, gray triangles are non-typhoon years, and the solid lines represent the linear regression line in typhoon years (black lines) and non-typhoon years (gray lines).

We found significantly positive correlations between DOC and DIC concentrations in typhoon years (Fig. 4); this was caused by substantial amounts of C loading into YYL during the strong typhoon period in 2015 (Fig. 8a, b) due to the occurrence of the last drought year in 2014 (Chiu et al., 2020). Not only did prolonged (or hysteresis) effects lead to dramatic C emissions in 2015 due to the drought year in 2014 (Chiu et al., 2020) but the rapid water retention and C loading during strong typhoon periods (Lin et al., 2022) also induced vigorous algal biomass production in the lower layer from May to September 2017 (Fig. 8d). Conversely, without the typhoon-induced mixing and refreshing of the water column, DIC and DOC concentrations in non-typhoon years might be mineralization dominated (Chiu et al., 2020; Hanson et al., 2015; Lin et al., 2022; Vachon et al., 2021), which could result in the positive linear relationship between DIC and DOC concentrations (Figs. 4b, 8a) as well as negative remineralization rates in non-typhoon years (Table 3). Therefore, these hydraulic, biogeochemical processes and the hysteresis effects might describe different patterns in the measurement data between typhoon years and non-typhoon years.

Thermal stratification and allochthonous C loading may drive the responses of fluxes to typhoons in YYL. The absolute values of fluxes were higher in typhoon years than in non-typhoon years (Table 4). We found that precipitation from typhoons loaded large amounts of allochthonous C into YYL during summer and autumn, which might explain the higher fluxes in autumn compared with the other seasons (Table 4). Typhoons dramatically changed the seasonal and interannual patterns of DIC fluxes due to river intrusion (Figs. 5a, b, S1), which proves to our hypothesis and concurs with the results of previous studies (Chiu et al., 2020; Lin et al., 2021, 2022). In summer, the DOC and DIC concentrations were spatially different between layers, observed as a two-layer system within the water column, because the upper and lower layers were inhibited due to strong thermal stratification (Lin et al., 2021, 2022), thereby resulting in positive upper DIC and DOC fluxes and negative lower DIC and DOC fluxes (Table 4).

Because of the absence of typhoon-induced mixing and allochthonous C loading, the total fluxes were lower in non-typhoon years than those in typhoon years (Table 4). In typhoon years, our results showed that typhoon-induced upwelling and loading increased by 102.2 $\mathrm{mg}\mathrm{DIC}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$ and 62.3 $\mathrm{mg}\mathrm{DOC}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$ in YYL (Table 4). Additionally, CO₂ emission (F_C) was 43 % higher (∼ 83 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$) in typhoon years than in non-typhoon years (Table 4). Therefore, we found that typhoon disturbances control DIC loading and C emissions in YYL, consistent with our previous studies (Chiu et al., 2020; Lin et al., 2021, 2022). Simultaneously, bio-photochemical mineralization and degradation may play a key role in shaping C fluxes, as colored DOC reduces ultraviolet radiation (UVR) and active photosynthetic radiation (PAR) (Allesson et al., 2021; Chiu et al., 2020; Schindler et al., 1996; Scully et al., 1996; Williamson et al., 1999), resulting in higher light intensity and water temperature in summer and consuming 3.7 times more DIC and DOC than in the other seasons (Table 4). These results suggested that the allochthonous C loading and light duration might be the most crucial factor for DIC and DOC fluxes in typhoon years. Conversely, the transportation rate shaped the seasonal C concentrations due to thermal stratification in non-typhoon years.

4.2 Model limitations under the extreme weather scenarios

Water temperature might be a crucial driver in controlling C fluxes in YYL (Chiu et al., 2020; Lin et al., 2021, 2022). We found that the fluxes and $F_{{\mathrm{CO}}_{\mathrm{2}}}$ in summer were usually higher than in winter (Tables 3, 4) due to the higher levels of photosynthesis, remineralization, and thermal stratification strength (Lin et al., 2021, 2022). With the conceptual two-layer C models (Table 3), photosynthesis absorption (α_PU and α_PL), remineralization (α_MU and α_ML), and transportation (w_I) represented the seasonal variations in DIC and DOC data well. These parameters of the conceptual two-layer C models appeared in reasonable patterns (Table 3). The higher remineralization and photosynthesis rates resulted in higher absolute values of fluxes in the autumn of typhoon years (Tables 3, 4). In non-typhoon years, the photosynthesis rates contributed to the total fluxes (Tables 3, 4). Thus, the conceptual two-layer C models suitably characterizes the seasonal and interannual responses of DIC and DOC fluxes to typhoons in YYL.

Under the extreme weather event scenarios (Level 1 and Level 2), the release of DIC and DOC was 30 %–118 % higher within YYL, considering thermal stratification and river intrusion. (Fig. 7). In non-typhoon years, Flux_DIC was higher than Flux_DOC (Fig. 6c, f) under the Level 1 and Level 2 scenarios or events for each season (Figs. S1, S2). However, the results showed that DOC fluxes were more sensitive than DIC fluxes to the impact of typhoon disturbances (Fig. 7). For autumn, our results show that the transportation rate parameter in typhoon years was over 10-fold higher than in non-typhoon years (Table 3) and that the DOC concentration was dominated by mineralization in typhoon years (Fig. 8), as biogeochemical processes within lakes, such as respiration, mineralization, and sediment burial, may impact DOC fluxes (Bartosiewicz et al., 2015; Hanson et al., 2015; Maranger et al., 2018). Simultaneously, the physical processes, such as typhoon-induced mixing and autumn overturns, upwelled the sediment and lower-layer DOC into the surface water in YYL (Kimura et al., 2017; Lin et al., 2022). Therefore, we suggest that autumn overturns and mineralization after typhoon disturbances might be vital with respect to influencing the vertical distribution of DOC concentrations in typhoon years.

Ejarque et al. (2021) successfully developed a conceptual one-layer model of DOC and DIC, considering bacterial respiration, photomineralization, and degradation in a temperate mountain lake. In addition, Nagatomo et al. (2023) revealed the significance of “freshwater carbon” capture due to photosynthesis and suggested that DIC might be underestimated if submerged vegetation is ignored. We suggest that photobiochemical processes (such as photomineralization) and submerged vegetation should be considered in the upper layer to clarify and validate the responses of the total C fluxes under extreme climates in a two-layer stratified lake from the aspect of a freshwater carbon ecosystem.

5 Conclusions

We successfully developed two-layer conceptual C models to obtain continuous DIC and DOC data in YYL and to simulate extreme conditions. Our conceptual two-layer C models revealed that allochthonous and autochthonous processes both accounted for C flux responses to typhoon disturbances. Without typhoons, thermal stratification was the primary driver of seasonal and interannual patterns of DIC and DOC. In typhoon years, the changes in seasonal river intrusion regimes in YYL resulted in a 3-fold higher total Flux_DIC than in non-typhoon years. However, our model should be improved for application to extreme climate scenarios by considering other processes within a lake, such as sediment burial, photodegradation processes, and anoxic conditions. The present results suggest that physical processes (river intrusion and vertical transportation) and biogeochemical processes (mineralization, photosynthesis, and respiration) in a small subtropical lake account for the C flux responses to typhoons on seasonal and interannual timescales.