Simulating CO₂ dynamics in Brackish Lake Obuchi, Japan: Low pCO₂ maintenance through diverse nutrient inputs

Coastal waters play a crucial role as a pathway for the transport of terrestrial chemical constituents to marine environments. Carbon cycling in coastal areas is essential for understanding climate change. Here, we used a three-dimensional hydrodynamic-ecosystem model to investigate CO₂ dynamics and its controlling factors in Lake Obuchi, a brackish lake in Aomori Prefecture, Japan. The model successfully reproduced the seasonal variations in water temperature, salinity, nutrients, chlorophyll-a, and carbonate system components in the lake. The simulation results show that the partial pressure of CO₂ (pCO₂) in Lake Obuchi remained below atmospheric levels throughout the year, suggesting that the lake functions as a CO₂ sink. The pCO₂ variability is dominated by spatially uniform seasonal patterns (90% of variance) driven by temperature changes, while summer periods exhibit distinct spatial patterns (6.0% of variance) that vary with river discharge magnitude. In Lake Obuchi, pCO₂ is maintained at low levels because of the carbonate system components of water masses formed by mixing river water and seawater. Additionally, primary production within these water masses further reduces pCO₂. This primary production is sustained by nitrogen-rich river water, combined with seasonal phosphorus supply from sediments during summer and autumn and from the sea during winter and spring. Our findings reveal that the seasonal combination of these nutrient sources maintains year-round primary production, which in turn contributes to the reduction in pCO₂ in Lake Obuchi.

In the global carbon cycle, coastal waters play a crucial role in connecting terrestrial and marine environments. Coastal waters are areas where the interactions of terrestrial matter, oceanic influences, and human impacts create complex dynamics. Although coastal waters cover only about 7% of the Earth’s surface, their importance in the carbon cycle exceeds that proportional to their area (Cai 2011).

Coastal waters can be categorized into several types based on their characteristics, including estuaries, brackish waters, tidal wetlands, and continental shelves (Borges and Abril 2012). Each of these categories has unique characteristics and affects the carbon cycle differently. For example, continental shelf waters and tidal wetlands are considered to be net CO₂ sinks (Cai 2011; Bauer et al. 2013), whereas estuaries are generally considered to be CO₂ sources, primarily because of the decomposition of organic matter supplied by rivers and the physical mixing of CO₂-rich river water with seawater (Frankignoulle et al. 1998).

However, studies have shown that CO₂ dynamics in estuaries and brackish waters are not straightforward. For instance, Jiang et al. (2008) reported that some eutrophic estuaries serve as CO₂ sinks because of enhanced photosynthetic activity of phytoplankton in nutrient-rich environments, which promotes CO₂ uptake. Kubo et al. (2022) demonstrated that in Lake Hamana, Japan, the influx of low-carbon, high-nitrogen water from sewage treatment facilities has led the lake to become a CO₂ sink. Herrmann et al. (2015) proposed that the ratio of dissolved inorganic nitrogen to total organic carbon in inflowing water can serve as an indicator of the metabolic balance in these water bodies. Water residence time is also crucial for understanding CO₂ dynamics in estuaries and brackish waters, as longer residence times tend to make the metabolic balance of the water body more apparent (Li et al. 2023). Indeed, Maher and Eyre (2012) investigated three estuaries in Australia with relatively long residence times (> 1 month) and demonstrated that they were autotrophic and serve as CO₂ sinks. Thus, carbon cycling in estuaries and brackish lakes depends on factors such as the chemical properties of inflowing river water and the residence time in these systems, indicating that further research in various coastal environments can deepen our understanding.

Lake Obuchi (Fig. 1; 40.9553N, 141.3619E) represents an ideal study site for investigating CO₂ dynamics in brackish waters because its hydrodynamics, water quality, and biogeochemical processes have been well documented over the past two decades (Ueda et al. 2000, 2004, 2006, 2007, 2011, 2017; Tsubono et al. 2021), and recent CO₂ measurements suggest that it acts as a CO₂ sink (Higashi et al. in press). Located in Rokkasho Village, Aomori Prefecture, Japan, Lake Obuchi receives freshwater input from the Futamata River and discharges water into the Pacific Ocean through the Obuchi River. The lake has a surface area of 3.7 km², average depth of 2.5 m, and maximum depth of 4.5 m (Ueda et al. 2006). Tides drive seawater into Lake Obuchi. During summer, a halocline develops with surface salinity of 10–15 and bottom water salinity of 20–30 (Ueda et al. 2000). When winds exceeding 10 m/s continue for approximately 3 h or more, the water column becomes fully mixed, temporarily eliminating the density stratification (Ueda et al. 2000). The average water residence time in Lake Obuchi is estimated to be 25 to 30 days (Ueda et al. 2011). In addition, a hydrodynamic transport model estimated the half-life of substances released into the lake as 11.5 to 15.9 days (Tsubono et al. 2021).

Lake Obuchi is known to have moderate nutrient concentrations. According to surveys conducted from 2004 to 2015, the annual mean and standard deviation of nitrate concentration was 6.7 ± 6.6 μmol/L, whereas the phosphate concentration was 0.26 ± 0.29 μmol/L (Ueda et al. 2017). Due to the inflow of nitrogen-rich river water from the Futamata River, primary production in Lake Obuchi is considered to be mainly limited by phosphorus (Ueda et al. 2017). In summer and autumn, density stratification leads to decreased dissolved oxygen concentrations in the bottom water, resulting in the release of phosphate and ammonium from sediments (Ueda et al. 2000) and the formation of a nutricline.

Net community production in Lake Obuchi has been measured as 7–27 μmolC/L/day in the surface water and 32–63 mmolC/m²/day in the water column (Higashi et al., in press). The annual mean and standard deviation of chlorophyll-a (Chl-a) concentration is 6 ± 5 μg/L (Ueda et al. 2017), and diatoms dominate throughout the year (Ueda et al. 2004). Higashi et al. (in press) measured partial pressure of CO₂ (pCO₂) in Lake Obuchi for approximately 30 days in total in 2018 and 2019, demonstrating that the lake was undersaturated relative to the atmosphere. Observational data are temporally limited, however, and the year-round CO₂ dynamics and their controlling factors in Lake Obuchi have not been fully studied.

The objective of this study is to develop a three-dimensional numerical model that simulates hydrodynamics, ecosystem processes, and cycling of chemical constituents in Lake Obuchi to clarify the year-round CO₂ dynamics and their controlling factors. Lake Obuchi is particularly suitable for investigating mechanisms through numerical simulation because of the availability of past observational data on hydrodynamics and water quality, including nutrients, as well as recent direct CO₂ measurements. In addition to the observational data reported by Higashi et al. (in press), this study presents new observational data. Section 2 describes the numerical model and observational data used in this study. Section 3 presents the comparison between numerical simulation results and observational data. Section 4 discusses the controlling factors of pCO₂ dynamics in Lake Obuchi. The final section provides a summary of the conclusions.

2.1 Numerical model

To simulate the hydrodynamics in Lake Obuchi, we employed the Regional Ocean Modeling System (ROMS; Shchepetkin and McWilliams 2005). This model is identical to the one previously used to estimate the time required for substance dilution in Lake Obuchi (Tsubono et al. 2021). ROMS solves the Reynolds-averaged Navier–Stokes equations under the Boussinesq and hydrostatic approximations and is a free-surface model that explicitly represents sea surface height variations. The computational grid was set to 40 meshes in both east–west and north–south directions (Fig. 1c). The model’s vertical structure was configured with five layers. ROMS employs terrain-following coordinates (Haidvogel et al. 2000), allowing calculations across five vertical layers regardless of total water depth. Bathymetric surveys were conducted from 2016 to 2017, and the results were linearly interpolated spatially to establish the model’s bathymetry (Fig. 1c). To accurately reproduce the average salinity within the lake, the bathymetric depth was adjusted once during the model setup phase as a static correction using Green’s function method (Tsubono et al. 2021) before the simulation runs and then remained constant throughout the simulations. The western, northern, and southern boundaries of the domain were set as closed boundaries, whereas the eastern boundary was set as open. Although Lake Obuchi partially freezes during winter, this model does not account for ice effects. The ROMS time step was set to 5 s, with the barotropic component being solved with 15 time-splits. Harmonic horizontal viscosity and diffusion were employed, with the horizontal viscosity coefficient set to 0.5 m²/s and the horizontal diffusion coefficient to 0.05 m²/s. The background vertical viscosity coefficient was set to 1 ×ばつ 10^–6 m²/s and the background vertical diffusion coefficient to 1 ×ばつ 10^–5 m²/s. The K-profile parameterization (Large et al. 1994) was used as the surface mixing closure scheme.

To calculate the carbon cycle in Lake Obuchi, we conducted simulations by coupling the marine ecosystem model Biogeochemical Elemental Cycling (BEC; Moore et al. 2013) with ROMS (Misumi et al. 2021). BEC is a nutrient-phytoplankton-zooplankton-detritus type model that simulates five nutrients (phosphate, nitrate, ammonium, silicate, iron), three types of phytoplankton (small phytoplankton, diatoms, and nitrogen fixers), and one type of zooplankton (Fig. S1). Production by calcifiers is treated implicitly as a portion of small phytoplankton production.

The parameters of the BEC model used in this study are shown in supplemental materials. Briefly, based on observational data showing diatom dominance in Lake Obuchi (Ueda et al. 2004), the model parameters were adjusted to enhance the competitive advantage of diatoms. Specifically, we reduced the maximum growth rate and photosynthetic efficiency of small phytoplankton, while increasing the nutrient uptake efficiency of diatoms by lowering their half-saturation constants. Additionally, we increased grazing pressure on small phytoplankton. Parameters were set to suppress calcifier blooms.

The BEC model employs the formulation of Dunne et al. (2007) for organic matter burial. This formulation expresses the burial flux (F_{POC_burial}) from organic carbon flux reaching the seafloor (F_{POC_bottom}) based on observational data of sediment accumulation rates and organic carbon preservation efficiency as: F_{POC_burial} = F_{POC_bottom} ×ばつ (0.013 + 0.53 ×ばつ F²_{POC_bottom})/(7.0 + F²_{POC_bottom}), where both F_{POC_bottom} and F_{POC_burial} have units of mmolC/m²/d. This relationship represents increasing burial efficiency with increasing organic carbon supply to the bottom waters. The BEC model introduces a scaling parameter (parm_POMbury) to the burial flux calculated by this formulation to adjust the model according to resolution and applied environments. In this study, we set this value to 0.5.

The BEC model varies calcium carbonate burial rates using the lysocline as a threshold. In this simulation, since no depths exceed the lysocline, all calcium carbonate reaching the bottom is buried. The BEC model formulates opal burial rates to vary from 4 to 20% depending on the opal flux reaching the bottom, with 20% in high-productivity regions where flux exceeds 2.0 mmol SiO₂/m²/d (Ragueneau et al. 2000). In this simulation, diatoms constantly dominate and productivity is consistently high compared to the open ocean throughout the year, so approximately 20% of opal reaching the bottom is buried.

Surface boundary conditions were based on meteorological data (temperature, specific humidity, precipitation, atmospheric pressure, and wind speed) observed in 2017 at the Institute for Environmental Sciences, adjacent to the northern side of Lake Obuchi (Fig. S2). For shortwave radiation, where observations were not available, we used data from the nearest grid point of the Normal Year Forcing dataset (Large and Yeager 2004). Net longwave radiation was calculated using a parameterization included in ROMS based on Berliand and Berliand (1952).

River discharge was based on daily averages of hourly measurements from the Futamata River during 2018–2021 (Fig. S3). Small streams exist on the northern and southern sides of Lake Obuchi (Fig. 1c). We distributed the observed discharge among the Futamata River and the northern and southern streams with weights of 0.6, 0.2, and 0.2, respectively. Hereafter, inflow from the Futamata River and these two streams is collectively referred to as "Futamata River inflow." Water temperature for the Futamata River was set to the adiabatic condition, and salinity was set to zero. For nutrients and carbonate system components in the Futamata River, we used constant concentrations based on median values calculated from in-situ measurements taken in various months from 2017 to 2021 (Table 1).

During summer and autumn, the sediments in Lake Obuchi become anoxic, leading to the release of phosphate and ammonium from the sediments into the bottom water (Ueda et al. 2000). To account for nutrient supply from the sediments, nudging with a time constant of 3 days was applied to the observed concentrations of phosphate and ammonium (Table 2) for model grid cells with total water depth greater than 3 m. Since this phenomenon is observed only during a limited period, the nudging was applied only from July to October. Given that iron is not considered a limiting factor because of its abundance in the brackish lake, the dissolved iron concentration was nudged to 10 nmol/L for all bottom grid cells.

For lateral boundary conditions at the ocean side, radiation condition with nudging was applied. The nudging was performed toward monthly climatological values of seawater temperature, salinity, and nutrient concentrations from the nearest grid point of the World Ocean Atlas 2013 (Locarnini et al. 2013; Zweng et al. 2013; Garcia et al. 2014a, b) (Fig. S4). Similarly, total dissolved inorganic carbon (DIC) and total alkalinity (TA) were nudged toward annual mean values from the nearest grid point of GLODAP v2 (Lauvset et al. ×ばつ 1° GLODAP version 2. Earth Syst Sci Data 8:325–340" href="/index.cgi/larger-text/https://progearthplanetsci.springeropen.com/articles/10.1186/s40645-025-00751-1#ref-CR19" id="ref-link-section-d71058211e1351">2016), 2022 μmol/L and 2251 μmol/L, respectively. Tidal forcing was based on interpolated results from eight tidal constituents (M2, S2, N2, K2, K1, O1, P1, Q1) of TPXO8 (Egbert and Erofeeva 2002), with amplitudes adjusted to match tidal amplitudes at the nearby Mutsu-Ogawara Port. Additionally, annual and semi-annual tidal constituents were incorporated based on harmonic analysis of hourly tidal height observations from 2010 to 2020 at the Mutsu-Ogawara tide-gauge station, located approximately 4 km south of the Obuchi River mouth.

Initial conditions were set to a state of rest, with water level at 0 m and temperature, salinity, and other chemical constituents set to the same values as the lateral boundary conditions. Biological parameters, such as phytoplankton carbon, were initialized to very low values (e.g., 0.1 μmol/L for phytoplankton carbon for all the functional types). The model was run for 8 years using repeated 1-year cycles of the surface boundary conditions, river discharge, and lateral boundary conditions, and the simulated data for the final 4 years were used for the analysis.

2.2 Observational data

We compared observational data gathered from the central part of Lake Obuchi (Station C in Fig. 1c) in 2016–2022 with the simulated results (Fig. 2). Those observational data from 2018 and 2019 have already been presented by Higashi et al. (in press). We deployed a mooring system at the surface of Station C for continuous CO₂ monitoring. Additionally, during the deployment and recovery of the mooring system, we collected water samples at depths of 0, 1, 2, and 3.5 m and analyzed nutrient concentrations, carbonate system components, and other parameters. The measurements of samples from 0- and 1-m depths were treated as surface layer data, whereas those from 2- and 3.5-m depths were treated as bottom layer data. For simulation results, values from the uppermost and lowermost layers were treated as surface and bottom layer data.

The CO₂ concentration (xCO₂) was measured using a portable CO₂ analyzer (CO2-09, Kimoto Electronic Co.) connected to a battery (12 V, 115Ah; M31MF, AC Delco) and solar panel (160 W; AT-MA160E, JA Solar) via a charge controller (SR-ML2430, LVYUAN). These were mounted on two float boats (Z1DR, Carmate) and moored. The analyzer used a nondispersive infrared analyzer to monitor both ambient air and air-equilibrated water through a passive equilibrator with a Teflon gas-permeable membrane at 0.3-m depth in the water column. The measurements were conducted in cycles of surface water (109 min), atmosphere (10 min), and drainage (1 min). When the measurement line switched from atmosphere to water, we excluded the first 29 min of measurements from the analysis considering the time lag in gas–liquid equilibrium exchange within the measurement line. The measured xCO₂ values were converted to pCO₂ values using monthly average water temperature and salinity data, assuming an atmospheric pressure of 1 atm.

2.3 Experimental cases

We conducted the (i) control experiment (i.e., CTL case) considering all biogeochemical sources. We then examined additional cases: (ii) CTL-Sed without sedimentary inputs of ammonium and phosphate; (iii, iv) CTL-Riv and CTL-Ocn without inputs of both nutrients and carbonate system components from the Futamata River and ocean, respectively; and (v) the CTL-Bio case, where primary production within the simulation domain was set to zero by setting a parameter of the maximum photosynthetic rate to zero. In these five cases, physical fields such as circulation were calculated to remain identical. For example, in the CTL-Riv case, we set nutrients, DIC, and TA from the Futamata River to 0 μmol/L while maintaining freshwater input to keep the same physical conditions as in the CTL case. In the CTL-Ocn case, boundary conditions for nutrients, DIC, and TA were changed from radiation with nudging to zero-gradient while maintaining the original boundary conditions for temperature and salinity. By maintaining identical physical fields, we focused on the influence of nutrient and carbonate system component inputs from each source on the biogeochemical cycling in Lake Obuchi. Since this study aims to elucidate the seasonal dynamics of biogeochemical cycling in Lake Obuchi throughout the year, we focused our analysis mainly on monthly median values for both simulation results and observational data.

2.4 Empirical orthogonal function analysis

Empirical Orthogonal Function (EOF) analysis was applied to the simulated pCO₂ data to extract the dominant spatiotemporal variability patterns. Prior to analysis, grid cells in the marine water inflow region were excluded to focus on variability within the lake. EOF analysis was conducted using the Eof function from the Python eofs.xarray library, applied to daily-averaged pCO₂ data from the 4-year analysis period. Spatial patterns (eigenvectors) and time series (principal component scores) for each mode were extracted to quantitatively separate the controlling factors of pCO₂ variability.

3.1 Comparison of simulation results with observational data

We compared simulated results for the CTL case with observational data at Station C (Fig. 3). The CTL case well captured observed seasonal temperature variations (Fig. 3a) and generally reproduced the seasonal patterns of salinity gradients (Fig. 3b). However, the observed vertical salinity gradient from summer to autumn exceeded that produced by the model, suggesting potential underestimation of density stratification in the CTL case during summer and autumn.

For Chl-a concentrations, the CTL case successfully reproduced the overall seasonal pattern, although it tended to overestimate concentrations in August and September (Fig. 3c). The CTL case adequately reproduced the general annual nitrate concentrations but failed to capture the characteristic increase in bottom concentrations from August and October (Fig. 3d). For simulated phosphate concentrations, the CTL case produced elevated phosphate concentrations from July to November, reflecting the sedimentary input, while those in April and June were underestimated (Fig. 3e). The observed nitrogen-to-phosphorus (N:P) molar ratios in both layers were observed to fall below 5 from June to September, suggesting nitrogen limitation during this period due to phosphate release from sediments combined with nitrogen depletion. The CTL case also demonstrated this shift to nitrogen-limited conditions between July and October. While the CTL case adequately reproduced the general annual silicate concentrations, it significantly overestimated concentrations in April and underestimated bottom layer concentrations in August (Fig. 3f).

The patterns of DIC and TA were similar, with the CTL case reproducing similar temporal variations generally falling within the range of observed values (Fig. 3g, h). However, the model underestimated the vertical concentration gradients of both DIC and TA in August, and underestimated both surface and bottom layer concentrations in September and October. The CTL case successfully reproduced the summer decrease in dissolved oxygen concentrations (Fig. 3i), but the model showed a general tendency to overestimate concentrations and underestimate vertical gradients compared to the observational data.

Overall, the CTL case successfully captured the general seasonal patterns of the observed variables. However, three primary discrepancies were identified: (1) underestimation of vertical concentration gradients during summer for salinity, DIC, TA, and dissolved oxygen; (2) underestimation of DIC and TA concentrations in September and October; and (3) concurrent underestimation of nitrate and phosphate concentrations and overestimation of silicate concentrations in April.

3.2 Variations in the simulated results

The CTL case successfully reproduced the observed seasonal variations in physical and biogeochemical properties of Lake Obuchi, allowing us to examine the detailed mechanisms controlling biogeochemical cycling using daily-averaged simulation results. The daily average water level at Station C exhibited fluctuations corresponding to spring and neap tides (Fig. 4). Following water level increases in the lake, salinity increased, presumably due to seawater inflow (Fig. 4a). The annual tidal cycle demonstrated minimum water levels in spring and maximum levels in summer. Consequently, spring conditions inhibited seawater inflow to the lake, whereas snowmelt-induced river discharge contributed to the annual minimum salinity values (Fig. 4a). Significant density stratification was observed following major river discharge events, especially during summer (Fig. 4b).

Nitrate concentrations mostly exceeded 1 μmol/L throughout the year, except during the depletion period in summer and autumn. Notable increases in surface nitrate concentrations were observed following discharge events from the Futamata River (Fig. 4c). However, in summer and autumn periods, these discharge-driven nitrate inputs were rapidly assimilated by phytoplankton because of nitrogen limitation, resulting in no observable increase in nitrate concentrations. Phosphate concentrations exhibited particularly elevated values in bottom waters during summer and autumn, which we attribute to sedimentary input (Fig. 4d). The simulated Chl-a concentrations remained at moderate levels from winter through spring, followed by a decrease during the spring–summer transition and a subsequent increase during summer and autumn (Fig. 4e).

To investigate variations in water and exchange of chemical constituents between the lake and ocean, we established a survey line (Fig. 1c) and analyzed simulated inflow and outflow fluxes (Fig. 5). Water flux showed consistently positive values, indicating outflow from the lake (Fig. 5a). The outflow generally reflects the export of freshwater input from the Futamata River, maintaining the long-term water balance of the lake. The outflow exhibited regular oscillations with approximately fortnightly cycles, suggesting a relationship with spring–neap tidal variations. However, this regular pattern was modified during specific periods. From February to early June, when the ocean water level was at its annual minimum, tidal inflow was reduced and these regular oscillations became less distinct. The pattern was also disrupted by river discharge events during summer and autumn, when significant rainfall events (indicated by downward red triangles in Fig. 4) induced marked increases in outflow.

While salinity flux patterns generally followed those of water flux, they showed negative values (indicating salinity inflow to the lake; Fig. 5b) during the periods of low water outflow (Fig. 5a). Although not shown in these daily-averaged data, semidiurnal tidal variations resulted in high-salinity seawater entering the lake when water flowed landward and low-salinity lake water exiting when water flowed seaward. This alternating pattern occasionally resulted in negative daily-averaged salinity flux, indicating net salt transport into the lake.

Nitrate and phosphate fluxes in Lake Obuchi showed contrasting patterns (Fig. 5c, d), reflecting their different biogeochemical cycles. Nitrate flux was always positive, indicating net export from the lake to the ocean (Fig. 5c). This pattern was driven by continuous nitrate input from the Futamata River that exceeded the lake’s consumption capacity. However, from July to November, nitrate flux approached zero. During this period, enhanced biological activity, supported by phosphate release from sediments, led to nitrate utilization within the lake (Fig. 3d, e). In contrast, phosphate flux showed predominantly negative values (Fig. 5d), indicating net import from the ocean to the lake. This pattern occurred because the lake was generally phosphorus-depleted because of high N:P ratios in Futamata River water, while the adjacent subarctic Pacific waters were characterized by relatively high phosphate concentrations. The flux direction reversed to positive values during summer and autumn, when sedimentary phosphate release enriched the lake water and resulted in net export to the ocean.

3.3 Comparison of simulated cases

Results from the CTL-Sed case were nearly identical to those of the CTL case from January to June (Fig. 6). We attribute this similarity to the timing of sediment nutrient supply, which was limited to July through October, thus having no impact during the earlier months. From July to November, simulated Chl-a concentrations were lower compared to those of the CTL case, whereas nitrate concentrations were higher and phosphate concentrations were lower. These patterns in the simulated results suggest that the absence of phosphate supply from sediments resulted in sustained phosphorus limitation from July to November, leading to decreased Chl-a concentrations, and the reduced phytoplankton utilization then led to increased nitrate concentrations. Simulated DIC and TA concentrations remained relatively unchanged from those of the CTL case, despite the significant decrease in Chl-a concentrations from July to November. Simulated dissolved oxygen concentrations were lower than those of the CTL case from July to November, which we attribute to reduced oxygen supply from photosynthetic activity.

Results from the CTL-Riv case exhibited consistently lower Chl-a concentrations throughout the year compared to those of the CTL case (Fig. 6). The simulated nitrate concentrations were nearly depleted in the CTL-Riv case, indicating that the Futamata River functioned as the primary source of nitrate input. The results showed elevated phosphate concentrations, suggesting that intensified nitrogen limitation resulted in excess phosphate concentrations that remained unutilized by phytoplankton. The observed marked decrease in silicate concentrations indicated that the Futamata River also served as the primary source of silicate input. Simulated DIC and TA concentrations in the lake approached marine concentrations because of the cessation of riverine input. Dissolved oxygen concentrations were slightly lower than those of the CTL case throughout the year, suggesting reduced oxygen supply due to decreased photosynthetic activity.

In the CTL-Ocn case, simulated Chl-a concentrations were lower from January to May, while showing minimal deviation from those of the CTL case after June (Fig. 6). Nitrate concentrations remained higher than in the CTL case from January through June. The elimination of phosphate input from the ocean (Fig. 5d) resulted in the excess nitrate that was not utilized by phytoplankton. This case yielded the highest silicate concentrations among all the cases. The increased silicate concentrations can be attributed to the elimination of ocean input, which provides relatively low silicate concentrations. Because oceanic inputs of DIC and TA were not considered in this case, the lake’s DIC and TA concentrations approached those of the Futamata River. Dissolved oxygen concentrations decreased from January to May, which can be attributed to reduced oxygen production due to decreased photosynthetic activity.

In the CTL-Bio case, where primary production was forcibly set to zero, both nitrate and phosphate concentrations increased because of the absence of phytoplankton uptake of nutrients (Fig. 6). Simulated DIC concentrations were slightly higher than those of the CTL case from July to November, while simulated TA showed differed negligible from those of the CTL case. Dissolved oxygen concentrations were lower throughout the year as compared with those of the CTL case, with the most pronounced differences occurring from July to November.

3.4 Variations of pCO₂

We compared the pCO₂ values simulated by the model with those observed at Station C (Fig. 7). The CTL case results fell within the range between observed maximum and minimum values for all months except November and December. Furthermore, for April, June, July, August, and October, the simulated values of the CTL case fell within the range between the 5th and 95th percentiles of the observational data. The September observations showed a bimodal distribution with groups ranging from 100 to 200 μatm and 400 to 600 μatm, with the CTL case results being consistent with the lower group. While the median value of the CTL case was close to the median observed value in July, it was lower in other months. These results indicate that while the CTL case generally reproduced the climatological seasonal variation of pCO₂ in Lake Obuchi, it tended to underestimate the values.

In the CTL-Sed case, Chl-a concentrations decreased from summer to autumn (Fig. 6a) because of the elimination of nutrient supply from sediments, resulting in increased pCO₂ during this period (Fig. 7). In the CTL-Riv case, the absence of riverine nutrient input led to decreased Chl-a concentrations throughout the year, causing an increase in pCO₂ compared to the CTL case. In the CTL-Ocn case, reduced nutrient supply from the ocean during winter and spring resulted in decreased Chl-a concentrations and increased pCO₂. Although the CTL-Ocn case showed similar Chl-a concentrations as in the CTL case in summer and autumn (Fig. 6a), it exhibited a marked decrease in pCO₂. This can be attributed to lower DIC in the CTL-Ocn case. Additionally, the reduced TA in this case decreased the buffering capacity of the carbonate system, amplifying the pCO₂ response to the DIC change. The CTL-Bio case showed increased pCO₂ throughout the year, with values approximately 100–200 μatm higher than the CTL case during all seasons, demonstrating the magnitude of biological influences on pCO₂.

Overall, the CTL case successfully reproduced the seasonal pCO₂ patterns observed in Lake Obuchi, confirming that pCO₂ values remain below atmospheric levels throughout the year. However, the model systematically underestimated pCO₂ values compared to observations, indicating that additional processes not included in the current simulation may contribute to the observed variability.

4.1 Seasonal variation of biogeochemical cycling

The CTL case results successfully reproduced the observed seasonal variations in water temperature, salinity, and biogeochemical variables in Lake Obuchi (Fig. 3). This correspondence suggests that the processes incorporated into the model effectively explain the annual physical and biogeochemical fluctuations in Lake Obuchi. Comparison of simulation cases (Fig. 6) indicated that primary production in Lake Obuchi is supported by various nutrient supply processes.

The Futamata River appears to be the primary source providing the baseline nutrients that sustain primary production in Lake Obuchi. When the nutrient supply from the Futamata River was eliminated in the simulation, Chl-a concentrations in Lake Obuchi decreased by more than half across all seasons, suggesting that nutrient input from the Futamata River serves as a major driver supporting primary production throughout the year. In the Futamata River, the median concentrations were 22.4 μmol/L for nitrate and 0.34 μmol/L for phosphate (N:P ratio, 65.9; Table 1); thus, nitrogen was comparatively abundant in the river water relative to the typical Redfield ratio for phytoplankton. Therefore, the Futamata River appears to be particularly important as a nitrogen source.

Nutrient supply from sediments appears to support primary production in Lake Obuchi particularly during summer and autumn, when primary production reaches its peak. Previous observations indicated that phosphate is released from sediments during summer and autumn, thus increasing concentrations (Ueda et al. 2000). This concentration increase can be interpreted as sediment-derived phosphorus compensating for the insufficient phosphorus supply from the Futamata River, thereby sustaining the high primary production during summer and autumn. In our observational data, despite the high N:P ratio (65.9) in Futamata River water, the N:P ratios in Lake Obuchi were remarkably low from summer to autumn (e.g., 0.35 and 0.19 in surface and bottom waters, respectively, in September). This finding indicates that sufficient phosphorus was supplied from sources other than the river to deplete nitrogen at Station C. Given that phosphate concentrations in the adjacent coastal waters are low during summer and thus insufficient as a source, it is reasonable to conclude that lake sediments serve as the primary source of phosphorus from summer to autumn.

We expect that oceanic nutrient supply support primary production in Lake Obuchi, particularly during winter and spring. In the subarctic North Pacific region facing Lake Obuchi, a mixed layer develops during winter, and nutrient concentrations increase until the spring bloom occurs. Consequently, nutrients are supplied through water exchange between Lake Obuchi and the ocean, especially in winter and spring. In this ocean area, phosphorus is relatively abundant as compared to nitrogen. For instance, the N:P ratio in the data used the lateral boundary condition is 7.1 annually and never exceeds 16 in any month. Therefore, the oceanic phosphorus supply compensates for the lake’s phosphorus deficiency, supporting primary production during winter and spring. As nutrient concentrations in the ocean decrease toward summer, nutrient supply from the ocean weakens from June to July, likely resulting in the decreased primary production.

The CTL case results showed several discrepancies with observational data. The simulation results underestimated the vertical gradients of salinity, DIC, TA, and dissolved oxygen during summer. This is presumably due to the model’s inadequate reproduction of density stratification formed by the inflow of river water and seawater. Increasing the number of vertical layers would improve the reproduction of density stratification.

The CTL case results underestimated nitrate and phosphate concentrations in April, while overestimating silicate concentrations. Because there were insufficient observational data, this study used median annual concentrations for chemical constituents of Futamata River water. The seasonal reproduction of nutrient concentrations in Lake Obuchi could potentially be improved by incorporating seasonal variations in the river’s chemical constituents.

Additionally, the CTL case underestimated DIC and TA concentrations in September and October. During the summer-to-autumn transition, sediments likely experience reducing conditions, leading to elevated DIC and TA concentrations in pore water due to sulfate reduction (Rassmann et al. 2020). The release of these constituents from sediments to the water column may explain the observed high concentrations during this period. Future work should focus on estimating and incorporating sediment-derived DIC and TA fluxes into the model.

4.2 Variation of pCO₂

The monthly variations in pCO₂ simulated in the CTL case generally fell within the 5th to 95th percentile range observed in the measurement data, indicating that the model captured the general trend of seasonal variations in pCO₂ in Lake Obuchi. However, the pCO₂ values in the CTL case tended to be on the lower side compared to the observational data. While the annual mean of the monthly medians of observed pCO₂ was 287 μatm, the annual mean pCO₂ simulated for the CTL case was 165 μatm, calculated from the median values for months for which observational data were available. The September pCO₂ values, for which measurement data were limited, were bimodally distributed with high and low groups, with the median of the observational data falling within the higher group (Fig. 7). Even after excluding the bimodally distributed September data, a substantial difference persisted: the annual means of the monthly median values decreased from 287 to 254 μatm for the observational data and increased from 165 to 167 μatm for the CTL case, maintaining a difference of 87 μatm.

We attribute the discrepancy in pCO₂ to the simulation’s inability to reproduce short-term pCO₂ increases. The observed pCO₂ data includes numerous upper outliers, and the density estimation portion of the violin plot is skewed upward (Fig. 7). This distribution indicates that short-term pCO₂ increases are driving up the median values. Since pCO₂ measurements are taken at 1-min intervals, these short-duration events are captured in the observational data. For instance, strong winds may cause vertical mixing of the lake water, bringing high-pCO₂ bottom water to the surface and resulting in sediment resuspension that mixes high-pCO₂ pore water to the surface, or the decomposition of resuspended organic matter may temporarily elevate surface-water pCO₂ (Abril et al. 2004). Although our model uses sufficiently short time steps to capture such rapid changes, these short-term pCO₂ increases are not reproduced in our simulation because it does not incorporate the effects of sediment resuspension, such as DIC release from pore water and decomposition of resuspended organic matter. This limitation may have resulted in lower monthly pCO₂ median values and their annual averages compared to those of the observational data.

The lack of consideration of particulate organic carbon (POC) input from the Futamata River in our simulations may also have contributed to the underestimation of pCO₂. Frankignoulle et al. (1998) reported that POC decomposition plays a significant role in increasing pCO₂ in estuarine regions. The model employed in this study does not account for horizontal POC transport, and therefore POC input from the Futamata River was not considered. Although the POC concentration in the Futamata River is not especially high at 14 ± 6 μmol/L (Higashi et al., in press), suggesting limited impact, future studies should investigate the lability of POC from the Futamata River and incorporate this information into the model.

During major flood events, increased turbidity inhibits photosynthesis, preventing decreases in pCO₂ (Cloern 1987). The current simulation did not account for such flood-induced increase in turbidity and its subsequent effect on light availability for photosynthesis, likely contributing to pCO₂ underestimation. This exclusion of turbidity effects would be particularly important during summer and autumn, when Lake Obuchi is nitrogen-limited and phytoplankton growth depends heavily on nitrate supply from the Futamata River. During these seasons, the omission from the model of turbidity effects during flood events may result in overestimated primary production and thus underestimated pCO₂. This potential bias is supported by the CTL case results, which show overestimated Chl-a concentrations during summer and autumn (Fig. 3c), highlighting the need to incorporate turbidity effects into the simulation.

The dynamics of pCO₂ are determined by the interaction of various factors, including nutrients, DIC, TA, and biological activity. Therefore, when discussing the ability of the model to simulate pCO₂ values, it is necessary to comprehensively evaluate not only pCO₂ itself but also all the underlying factors that affect it. The variability of observed pCO₂ was large, and since simulation results fall within the range of observational scatter in many months, it is difficult to distinguish the superiority among cases based solely on pCO₂ results. This outcome highlights the importance of examining the underlying components: the CTL-Riv case, while matching the observed pCO₂ in certain months, shows DIC and TA values close to those of seawater (Fig. 6e, f), indicating that it does not reflect realistic carbonate system dynamics in the lake. This finding demonstrates that agreement between the model and simulations in pCO₂ does not necessarily indicate consistency in the underlying factors. In contrast, the CTL case shows good agreement not only in pCO₂ but also in the seasonal variations and concentrations of nutrients and carbonate system components when compared with observational data (Fig. 3), indicating that it better represents the actual lake processes.

The observational data showed that the annual median pCO₂ in Lake Obuchi is 287 μatm, which is lower than atmospheric pCO₂ (approximately 400 μatm). Model simulations estimated that biological activity reduces pCO₂ by 100–200 μatm. Combining these observational and model results suggests that without biological activity, pCO₂ could reach 387–487 μatm, potentially making the lake a seasonal CO₂ source. This significant biological control on pCO₂ is sustained by the relatively long water residence time in Lake Obuchi, which facilitates continuous CO₂ fixation by photosynthetic organisms.

Zostera eelgrasses are known to flourish in summer and autumn in Lake Obuchi (Ueda et al. 2006). Could Zostera eelgrasses photosynthetic CO₂ uptake be an important driver of CO₂ variability in the lake water? The total carbon content of Zostera spp. in Lake Obuchi is estimated at 1400 ± 590 kg C (Ueda et al. 2006), and the average DIC concentration in the lake water is approximately 1500 μmol/L. Given the lake volume of 9.3 ×ばつ 10⁶ m³ (3.7 km² ×ばつ 2.5 m), the total DIC content amounts to 1.7 ×ばつ 10⁵ kg C. Considering that the carbon content of Zostera spp. represents only about 0.8% of the DIC content, coupled with the water residence time of 25 to 30 days (Ueda et al. 2011), it is reasonable to conclude that eelgrasses make only a minimal contribution to CO₂ uptake in Lake Obuchi.

4.3 Spatiotemporal patterns of pCO₂ variability

To analyze the two-dimensional spatial and temporal variability patterns of pCO₂, EOF analysis was applied to the simulated pCO₂ results to clarify and quantitatively separate the dominant spatiotemporal patterns (Fig. 8). The spatial patterns of the first mode (90%) showed uniform with seasonal variation (negative values in winter, positive values in summer), indicating that temperature-driven seasonal patterns dominate pCO₂ variability.

During summer and autumn periods, phosphorus limitation is relieved by sedimentary supply and nitrogen becomes limiting. Since the Futamata River has high N:P ratios, it becomes an effective nitrogen source under these nitrogen-limited conditions. Therefore, discharge variations from the Futamata River can influence spatial variability patterns of pCO₂, which appear as the second and third modes. During other seasons, chronic phosphorus limitation eliminates the connection between Futamata River discharge variations and spatial pCO₂ patterns.

Differences in river discharge magnitude during summer-to-autumn periods generate different pCO₂ response patterns. During normal-scale discharge events, the third mode becomes dominant, with nitrate and Chl-a increasing only on the western side near the Futamata River mouth, leading to localized pCO₂ decline (Fig. 9, Aug. 1). The third mode exhibits an east–west gradient pattern with negative values on the western side and positive values on the eastern side. This occurs because the nitrogen input from normal-scale Futamata River discharge is insufficient to support increased production throughout the entire lake and is consumed locally on the western side.

During large-scale discharge events, the second mode becomes dominant, with Chl-a increasing in a vast area of the lake and widespread pCO₂ decline is observed (Fig. 10, Aug. 11). The second mode shows spatial patterns with negative values in the central lake and positive values in the periphery, representing enhanced CO₂ consumption by primary production throughout the lake.

The temporal changes in the second and third modes suggest dynamic processes in the brackish lake. Following the formation of widespread low pCO₂ conditions during large-scale discharge events (second mode), a phenomenon where the third mode changes in the negative direction was observed (Fig. 8, Aug 11). This potentially represents a restoration process which the lake recovers from the widespread low pCO₂ state back to the typical east–west gradient pattern.

This multi-modal structure quantitatively demonstrates that pCO₂ dynamics in brackish lakes are controlled by the composite action of temperature effects, nutrient supply magnitude, sediment processes, and discharge-driven spatial patterns. However, these spatial response patterns associated with different discharge magnitudes have not been validated through direct field observations and require future verification. The turbidity effects during flood events may also influence the spatial response patterns. The omission of flood-induced turbidity in our model likely resulted in overestimated primary production during discharge events, which could enhance the magnitude of the second and third EOF modes, particularly the second mode representing lake-wide pCO₂ decline.

4.4 Implications for the carbon cycle in estuaries

While estuaries have generally been recognized as CO₂ source regions (Frankignoulle et al. 1998), some estuarine water bodies have been identified as CO₂ sinks depending on their nutrient conditions and the characteristics of inflowing river water (Jiang et al. 2008; Kubo et al. 2022). Previous pCO₂ observations in Lake Obuchi had suggested that it functions as a CO₂ sink (Higashi et al., in press). This study supported the previous observational evidence and further demonstrated, through numerical simulations, that the lake likely maintains its CO₂ sink capacity throughout the annual cycle.

The average CO₂ uptake in Lake Obuchi was maintained by active primary production supported by diverse nutrient supplies, as well as the nonlinear mixing effect of the carbonate system that produces water masses with lower pCO₂ than both river water and seawater end-members (Higashi et al., accepted). While river water is generally known to be phosphorus-limited relative to nitrogen (Meybeck 1982), phosphate release from sediments in a reductive state compensated for this nutrient imbalance. Our results suggest the possibility of a biogeochemical feedback loop: nitrogen-rich river input initiates primary production and, when combined with density stratification, enhances phosphate release from sediments. This additional phosphate supply could further stimulate primary production, potentially creating a self-reinforcing cycle of enhanced biological activity. Further investigation of such biogeochemical feedbacks in estuarine systems would be valuable. Therefore, understanding the CO₂ uptake mechanisms in estuaries requires an integrated approach considering the relationships among physical stratification, nutrient supply from various sources, and biological processes.

Furthermore, the role of oceanic nutrient supply to estuaries through tidal exchange has been underappreciated in previous studies. Different oceanographic regions can significantly influence this process. For example, estuaries adjacent to nutrient-rich subarctic waters receive higher nutrient inputs through tidal exchange compared to those connected to nutrient-poor subtropical waters. These regional differences in oceanic nutrient supply underscore the importance of the characteristics of adjacent marine waters as a controlling factor in estuarine carbon cycling.

This study used a three-dimensional hydrodynamic-ecosystem model to investigate CO₂ dynamics and their controlling factors in Lake Obuchi, a brackish lake in Aomori Prefecture, Japan. Simulation results demonstrated that pCO₂ levels in Lake Obuchi remained below atmospheric levels throughout the year, suggesting that the lake functions as a CO₂ sink. The primary mechanism driving CO₂ uptake in Lake Obuchi was identified as the formation of low-pCO₂ water masses through the mixing of river and seawater. Additionally, primary production, supported by multiple nutrient sources, further contributed to the reduction in pCO₂. The pCO₂ variability is controlled by spatially uniform seasonal patterns driven mainly by temperature, while summer river discharge events generate distinct spatial response patterns depending on river inflow magnitude.

Our findings emphasize that understanding carbon cycling in estuarine systems requires comprehensive analysis of multiple processes, including not only the characteristics and residence time of water masses formed through river–sea mixing, but also biogeochemical feedbacks between riverine nutrient input and sedimentary processes under stratified conditions. The nutrient characteristics of the adjacent coastal waters can also control biological activity. Such an integrated understanding of physical mixing processes, biogeochemical feedbacks, and regional oceanographic settings is essential for evaluating the role of estuarine systems in coastal carbon cycling.

The simulation data are available in https://doi.org/10.5281/zenodo.15093743.

BEC:

Biogeochemical elemental cycling

DIC:

Dissolved inorganic carbon

DOC:

Dissolved organic carbon

pCO₂ :

Partial pressure of CO₂

POC:

Particulate organic carbon

ROMS:

Regional ocean modeling system

TA:

Total alkalinity

N:P ratios:

Nitrogen-to-phosphorus ratios

We thank R. Niwa for preparing input data used in the simulations. We acknowledge Y. Maeda, Y. Kawa, T. Hirayama, J. Furukawa, R. Higashi, C. Ochi, and T. Fukumaki for their support in obtaining the observational data.

This work was partly supported by JSPS KAKENHI Grant Number 21H05056, 23H04824 and 25H01911.

Authors and Affiliations

1. Sustainable System Research Laboratory, Central Research Institute of Electric Power Industry, 1646 Abiko, Abiko, Chiba, 270-1194, Japan

   Kazuhiro Misumi, Takaki Tsubono & Daisuke Tsumune

2. Center for Research in Radiation, Isotopes, and Earth System Sciences, University of Tsukuba, 1-1-1 Tennodai, Tsukuba, Ibaraki, 305-8572, Japan

   Daisuke Tsumune

3. Faculty of Fisheries Sciences, Hokkaido University, Kita 10, Nishi 5, Kita-Ku, Sapporo, Hokkaido, 060-0810, Japan

   Takeshi Yoshimura

4. Department of Radioecology, Institute for Environmental Sciences, Rokkasho, Aomori, 039-3212, Japan

   Shinji Ueda

Contributions

KM planned and conceptualized the overall research, conducted numerical simulations, analyzed observational data, and performed field observations. TT carried out numerical simulations and participated in field observations. DT planned and designed the overall research framework and conducted field observations. TY performed field observations and chemical analyses. SU provided meteorological data and river discharge information. All authors participated in discussions and contributed to the preparation of the final manuscript.

Corresponding author

Correspondence to Kazuhiro Misumi.

Competing interests

The authors declare that they have no competing interest.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions