The total electricity generation worldwide estimated from the dataset with geolocation information was 19,163 TWh in 2017, which is 78.6% of the 24,394 TWh reported in the database provided by the EIA (https://www.eia.gov/international/overview/world; lastly accessed on 2 January 2024). This discrepancy may arise from the scarcity of geolocation data for wind and solar power plants (Byers et al., 2021). On the other hand, geolocation data for conventional thermal power plants are more widely available (~100%). The electric generation within 20 km of the coastline (6,838 TWh) accounts for 35.7% of the total electric generation. Regional electric generation differs notably between latitudinal bands ( Figure 1 ). In particular, along the coastline, 28.1% and 37.9% of worldwide electricity generation are observed in the 20°N–30°N and 30°N–40°N latitude bands, respectively.

The water use intensity and thermal efficiency vary seasonally depending on the difference in temperature between the internal heat source and the external environment (De La Guardia et al., 2022). In South Korea, the outlet temperature exhibited seasonal variations, ranging from 0 to 15°C higher than the inlet temperature, with an average difference of 6.6°C between 2018 and 2020 ( Supplementary Figure 1 ). To the best of our knowledge, there is no publicly available data for the temperature information, including the temperature of the source water and the difference between the inlet and outlet of the cooling water, that can be applied to global estimates. Therefore, the water use intensity was calculated based on the CWD and electricity generation using monthly data from 22 fossil fuel and 5 nuclear power plants in South Korea. An annual electricity generation of 375 TWh and a CWD of 58.5 km³ has been reported from 2018 to 2020, from which the corresponding water use intensity was estimated to be 0.156 km³ TWh⁻¹, which is consistent with the range reported in previous regional studies (0.136–0.207 km³ TWh⁻¹) for once-through cooling systems in fossil fuel and nuclear power plants (Peer & Sanders, 2016; Pan et al., 2018; Lin et al., 2021; De La Guardia et al., 2022) ( Figure 2 ). The average water use intensity is 0.137 km³ TWh⁻¹ for fossil fuel power plants and 0.184 km³ TWh⁻¹ for nuclear power plants, which likely reflects the differences in the thermal efficiency between the two fuel sources ( Figure 2 ). Although individual nuclear and fossil fuel plants can vary widely in their efficiency, the average thermal efficiency of fossil fuel plants (40−60%) is greater than that of nuclear plants (34−36%) (Wibisono & Shwageraus, 2016). As shown in Figure 2 , coal plants in South Korea also have higher thermal efficiency and require less cooling water than nuclear power plants.

The global CWD was then calculated using annual electric generation (6,838 TWh) for coastal areas and the water use intensity estimated from the South Korea data (0.156 km³ TWh⁻¹). Overall, 1,066 km³ of cooling water was estimated to be discharged yearly, which is comparable to the global wastewater discharge (~1,000 km³; Lu et al., 2018) and the discharge from Changjiang River (900 km³; Yang et al., 2015). However, our estimate is three times higher than the reported seawater withdrawal rate (300 km³) for coastal areas based on satellite imagery (Lohrmann et al., 2019). One reason for this difference could be that we assume that only once-through cooling systems are used in coastal areas, which can result in the overestimation of the CWD due to the higher water use intensity of once-through cooling systems compared to other cooling technologies (e.g., recirculating cooling towers and dry cooling systems). Another possible explanation is that estuaries are often found close to large rivers and oceans, which can be an important source of fresh (or brackish) water for cooling; however, our estimates did not take into account freshwater use. Despite this, it should be noted that once-though cooling systems are the most common cooling technology for the power plants built along the coastline, and differences in CWD estimates are apparent in country-level comparisons (De La Guardia et al., 2022).

To verify our methodology for estimating the CWD, our CWD estimates for China, Japan, the United States, and South Korea were compared with those from previous research. Our estimates for the CWD are 222.3, 138.7, 109.0, and 79.8 km³ for China, Japan, the United States, and South Korea, respectively, which are higher than the 14.4, 29.4, 1.8, and 28.6 km³ reported by Lohrmann et al. (2019). However, previously reported seawater withdrawal volumes of 108–201 km³ for China (Zhang et al., 2016; Lin et al., 2021) are comparable to our estimates. The seawater withdrawal data reported by Zhang et al. (2016) for China in 2011 was approximately scaled for 2017 based on the annual electric generation of China between 2011 and 2017. In addition, for the United States, the reported seawater withdrawal of 31.8–42.8 km³ (Harris & Diehl, 2019) is 29–39% of our estimates. The variation in the estimates for different countries can be attributed to the relative proportion of freshwater- and seawater-based cooling systems in coastal areas. Nevertheless, the order of magnitude difference in the estimated CWD at the country and regional level suggests that additional assessments are needed using a more extensive database for CWD based on estimates from various methods.

When a specific power generation capacity is given, the temperature difference (ΔT) between the inlet and outlet is inversely proportional to the amount of CWD. Therefore, this study utilized the ΔT of ~6.6°C, corresponding to the water use intensity (0.156 km³ TWh⁻¹). This ΔT was estimated to lead to a variation of 70–170 μatm in pCO₂, which is comparable to the seasonal variation in pCO₂ observed for the temperate open ocean (<90 μatm; Takahashi et al., 2014). To examine the potential impact of the CWD on the absorption of atmospheric CO₂, we calculated the rise in the pCO₂ of the discharged water and the change in the DIC (ΔDIC) resulting from perturbations in the air–sea CO₂ flux. It should be noted that the calculation for assessing the ΔDIC was conducted using pCO₂–TA pairs, and a coefficient (∂ln pCO₂/∂SST) of 0.0423°C⁻¹ (Takahashi et al., 1993; Takahashi et al., 2002) was not directly applied. This is because the influence of temperature on pCO₂ exhibited values different from 0.0423°C⁻¹, especially under low-temperature or salinity conditions (Jiang et al., 2008; Calleja et al., 2013; Gallego et al., 2018).

The upper limit of the ΔDIC was calculated assuming that the elevated pCO₂ returns to the original condition through the air–sea CO₂ exchange without considering further biophysical processes. It was also assumed that the relaxation back to initial conditions would occur independent of oceanic and meteorological variables. The meridional variation in the ΔDIC indicates that warm and salty waters in the subtropics (−57 ± 3 μmol kg⁻¹) are more sensitive to pCO₂ perturbations compared to cold polar waters with low salinity (−47 ± 6 μmol kg⁻¹; Figure 3A ). For initially oversaturated seawater with respect to the atmosphere, a negative value of the ΔDIC represents the expected value from the greater release of CO₂ due to warming. Thus, a negative value for the ΔDIC for initial and post-discharge undersaturated seawater represents the probable decline in CO₂ uptake potential. Based on the expected ΔDIC and the relative contribution of CWD in each latitudinal band, we evaluated the impact of CWD on the average surface DIC (ΔDIC = −55 ± 4 μmol kg⁻¹) in global coastal waters. Overall, the estimated warming caused by CWD (1,066 km³) could potentially lead to a reduction in the oceanic CO₂ uptake of 0.69 Tg C yr⁻¹ (~2.5 Tg CO₂ yr⁻¹).

The effect of seawater pCO₂ perturbation depends on how rapidly the perturbed seawater is restored to its initial state in terms of pCO₂ and the SST. Thermal discharge to a net CO₂ source area would quickly lead to a return to the initial pCO₂ due to the enhanced sea-to-air CO₂ flux (i.e., outgassing). The thermal stratification and surface residence time associated with the seasonal SST cycle can also regulate the response of discharged water. When air–sea heat exchange (with heat loss reducing oceanic pCO₂ in accordance with ∂ln pCO₂/∂SST = 0.027–0.0423°C⁻¹) is faster than air–sea CO₂ exchange, the pCO₂ of cooling water would recover quickly. In this case, the DIC deficit of the thermal plume would be smaller than that expected solely from CO₂ outgassing. Therefore, our estimate of ~0.69 Tg C yr⁻¹ is likely to represent an upper theoretical limit based on the current data.

In addition to the upper limit estimation, we calculated ΔDIC while considering varying gas exchange rates to investigate the realistic response of discharged water associated with instantaneous warming at different wind speeds and air–sea CO₂ gradients. In the actual environment, discharged cooling water is mixed with surrounding (source) waters. Assuming a 10-fold dilution (leading to the ΔT of 0.66°C), the ΔDIC is proportionally reduced to a factor of one-tenth. Since it is challenging to observe the dilution signal with the gridded product used in this study, the assumption of one-tenth dilution was made with the expectation that ΔDIC aligns with a level comparable with accuracy and precision ( ± 2–3 μmol kg⁻¹) of the DIC measurement ( Supplementary Figure 2 ). It should be noted that taking into account the dilution factor, the overall change in DIC of CWD remains relatively consistent. Initially, we calculated the ΔDIC of the source water associated with the air–sea CO₂ gas exchange using its original condition, which corresponds to an initial pCO₂ ranging from 257 to 557 μatm, and the average value of TA of 2248 μmol kg⁻¹, S = 33.7, and T = 17.1°C of the coastal ocean datasets ( Supplementary Figure 2 ). Subsequently, we calculated the ΔDIC for the cooling water, assuming a 10-fold dilution, resulting in a temperature increase of 0.66°C. The difference in ΔDIC between source water and cooling water represents the impact of temperature elevation on air–sea CO₂ exchange and the DIC inventory.

These calculations of the realistic response of carbonate parameters associated with instantaneous warming were explicitly performed for the 6-month warming period (spring to summer) when SST increased. It was assumed that the sea-to-air heat flux was minimized during this period. Considering the average wind speed of 4 m s⁻¹ between 30–40°N of our study region and an assumed mixed layer depth of 15 m, the final DIC concentration of discharged cooling water was ~27 μmol kg⁻¹ lower than that of source water ( Figure 4B ). This decrease in DIC is primarily attributed to the enhanced CO₂ outgassing or weakened CO₂ uptake depending on the air–sea CO₂ gradient ( Supplementary Figure 2C ). The magnitude of the ΔDIC is ~49% of our upper limit of ~55 ± 4 μmol kg⁻¹ expected from the complete relaxation of pCO₂ perturbation. Since only a 6-month period was considered, the ΔDIC was applied to half of the annual discharged water volume. Additionally, taking into account the average duration of the discharged water during the warming period, the CWD can potentially impact the surrounding water for approximately 3 months. Overall, it was estimated that the upper limit (~0.69 Tg C yr⁻¹) of the impact of CWD has reduced to ~13% of its original magnitude (~0.09 Tg C yr⁻¹).

Our estimates are sensitive to variations in wind speed, mixed layer depth, and air–sea CO₂ gradient. It is difficult to obtain precise environmental information for individual coastal areas. Furthermore, due to the non-linear response of the ΔDIC to variations in individual factors, our exploration has been limited to simple sensitivity information within a narrow range. Deviations in wind speed by 1 m s⁻¹ from the assumed 4 m s⁻¹ resulted in approximately 30% change in the magnitude of ΔDIC, indicating a positive correlation ( Figure 4B ). Utilizing mixed layer depths of 10 and 20 m led to a 20% increase and decrease in the magnitude of ΔDIC, respectively, compared to employing a mixed layer depth of 15 m. The uncertainty associated with using gridded products to assess the biogeochemical properties of the complex coastal areas was challenging to quantify accurately. Instead, we conducted a sensitivity analysis using different air–sea CO₂ gradients (seawater pCO₂ – atmospheric pCO₂; −150 to 150 μatm). The assumed ΔpCO₂ range of ±150 μatm corresponds to 3σ ( ± 125 μatm) of the gridded pCO₂ database and the seasonal amplitude of pCO₂ measured at coastal areas (Torres et al., 2021). These variations of ±150 μatm contributed to an 11% difference in our ΔDIC estimate ( Supplementary Figure 2C ). The sensitivity of ΔDIC to variations in the pCO₂ is also significantly influenced by the Revelle factor, which will be further discussed in the following section. Due to the significant variations in ΔDIC associated with each variable, accurate quantification requires precise measurement of CO₂ and environmental parameters.

The effect of the CWD-associated carbon flux is minor compared to that of other carbon fluxes, which are several orders of magnitude higher, in the global ocean. However, for regions with significant coastal thermal pollution, this impact should be taken into account during the assessment of the regional CO₂ budget, which includes the oceanic CO₂ flux with the atmosphere and terrestrial environment (via river runoff) on a national scale. For example, China, Japan, and South Korea contribute to 43.2% of global thermoelectric generation and subsequent CWD in coastal regions. In the coastal areas of China, a CWD of ~222.3 km³ yr⁻¹ is produced by coal and nuclear power plants, which is ~25% of the discharge from Changjiang River (~900 km³ yr⁻¹), while a CWD of 79.8 km³ yr⁻¹ estimated for South Korea is greater than the discharge from its five main rivers (~43 km³ yr⁻¹; Lee E. et al., 2021). Of the carbon fluxes in the coastal areas of South Korea, our upper estimate for the potential decrease in the CO₂ influx (~58 Gg C yr⁻¹) is comparable to the amount of organic carbon sequestrated by a tidal flat (~71 Gg C yr⁻¹; Lee et al., 2021) and ~10% of the DIC release from the five main rivers (~450 Gg C yr⁻¹; Lee E. et al., 2021). Our estimate (~58 Gg C yr⁻¹) could decrease to ~13% when considering other physical factors such as a mixed layer depth and warming period. Globally, the estimated ΔDIC due to CWD (0.09–0.69 Tg C yr⁻¹) is relatively minor (<11%) compared to the carbon sequestration rate of the top sediment (6.8 Tg C yr⁻¹; Chen & Lee, 2022) in global tidal flats (127,921 km²).

Currently, available data for physical factors (e.g., hydrodynamics, mixed layer depth, and air–sea heat flux) are insufficient for the accurate estimation of the gas exchange timescale for the restoration of perturbed oceanic CO₂. Given the difficulty in modeling the relatively slow process of gas exchange (i.e., over weeks to months) for the CO₂, we focus on variation and heterogeneity in coastal systems with a variety of salinity and carbonate parameters. This analysis can serve as a roadmap for future attempts to measure the impact of the CWD on more specific regions (e.g., estuaries).

Buffer factors determine the relative sensitivity of carbon parameters (e.g., pCO₂, DIC, pH, and calcium carbonate saturation state) to perturbations from normal conditions. Because we assessed the potential impact of the CWD-induced change in pCO₂, the associated ΔDIC is proportional to the inverse of the conventional buffer factor (i.e., the sensitivity of pCO₂ to changes in the DIC) and varied regionally. In coastal oceans, the TA is generally higher than the DIC, with a DIC/TA ratio ranging from 0.82 to 1.02. The lowest DIC/TA ratios are observed in the warm subtropical regions with the most robust buffering capacity (Revelle factor of ~8). The lowest buffering capacity (Revelle factor of ~17) is observed in subpolar regions, where the DIC/TA ratio approaches unity ( Figure 3 ).

Near-shore areas and estuaries are generally known as net sources of CO₂ emitted to the atmosphere due to the microbial oxidation of organic carbon and the input of poorly buffered freshwater (Chen & Borges, 2009; Cai et al., 2021). Due to the lower buffering capacities (i.e., a higher Revelle factor) in estuarine and low-salinity coastal waters, a fractional change in pCO₂ results in a weaker fractional change in the DIC compared with seawater. As the salinity approaches zero, the decrease in the DIC concentration with a lower buffer factor reduces the ΔDIC associated with CWD ( Figure 3 ). However, when the riverine TA is high (> 1,800 μmol kg⁻¹) with a strong buffering capacity even within a lower salinity range (e.g., the Mississippi River and Changjiang River; Cai et al., 2021), ΔDIC may exceed than expected value from the ΔDIC/salinity trend derived from the endmember (TA = ~548 μmol kg⁻¹ at zero salinity) obtained from coastal ocean dataset (OceanSODA-ETHZ; Gregor & Gruber, 2021). Additionally, the ΔDIC/salinity trend can be significantly disrupted when freshwater mixes with seawater due to changes in biological production. Enhanced biological production, for example, can lead to an increase in buffer capacity. Consequently, it is essential to assess regional variations in the sensitivity to thermal stress to comprehensively evaluate the physical effect of thermal discharge on CO₂ dynamics.