Estimating the seasonal impact of optically significant water constituents on surface heating rates in the western Baltic Sea

Heating rates induced by optically significant water constituents (OSCs), e.g. phytoplankton and coloured dissolved organic matter (CDOM), contribute to the seasonal modulation of thermal energy fluxes across the ocean–atmosphere interface in coastal and regional shelf seas. This is investigated in the western Baltic Sea, a marginal sea characterised by considerable inputs of freshwater carrying nutrients and CDOM and by complex bio-optical and hydrodynamic processes. Using a coupled bio-optical ocean model (ROMS–Bio-Optic), the inherent optical properties of different OSCs are modelled under varying environmental conditions, and the underwater light field is spectrally resolved in a dynamic ocean. We estimate the relative contribution of these OSCs to the divergence of the heat flux and heating rates and find that, while phytoplankton and CDOM both contribute to surface heating in summer, phytoplankton dominates the OSC contribution to heating in spring, and CDOM dominates the OSC contribution to heating in autumn. The study shows that seasonal and spatial changes in OSCs in the western Baltic Sea have a small but noticeable impact on radiative heating in surface waters and consequences for the exchange of energy fluxes across the air–sea interface and the distribution of heat within the water column. In the Pomeranian Bight, where riverine influx of CDOM is strongest, water-constituent-induced heating rates in surface waters in 2018 are estimated to be between 0.8 and 0.9 K m⁻¹ d⁻¹ in spring and summer, predominantly as a result of increased absorption by phytoplankton and CDOM. Further offshore, OSC-induced heating rates during the same periods are estimated to be between 0.4 and 0.8 K m⁻¹ d⁻¹. Warmer surface waters are balanced by cooler subsurface waters. Surface heat fluxes (latent, sensible and longwave) respond to warmer sea surface temperatures, with a small increase in heat loss to the atmosphere of 5 W m⁻² during the period April to September. We find relatively good agreement between our modelled water constituent absorption and in situ and satellite observations. More rigorous co-located heating-rate calculations using an atmosphere–ocean radiative transfer model provide evidence of the suitability of the ROMS–Bio-Optic model for estimating heating rates.

1 Introduction

Radiant energy fluxes impact biological production in the ocean and are modulated in turn as a result of biological production. This has fundamental consequences for upper-ocean physics, surface nutrient supply, net primary and export production, and the exchange of soluble gases across the air–sea interface into the marine atmospheric boundary layer. The contribution of optically significant water constituents (OSCs) to heating rates in the upper ocean is connected to net primary and export production through the direct effect of temperature on metabolic rates of marine plankton and increased stratification and reduced vertical exchange of nutrients. This plays an important role in controlling the flow of carbon and energy through pelagic systems (Wohlers et al., 2009; Taucher and Oschlies, 2011), particularly the partitioning between particulate and dissolved organic carbon, the transfer of primary produced organic matter to higher trophic levels, the efficiency of the biological carbon pump and the exchange of CO₂ across the air–sea interface. Shelf seas and coastal waters are often characterised by a highly variable presence of inorganic suspended particulate matter and coloured dissolved organic matter (CDOM). CDOM is the fraction of dissolved organic matter (DOM) that absorbs light in natural waters in parts of the ultraviolet and visible spectral ranges (ca. 200–550 nm). It is present throughout the world oceans, in both open and deep waters and in coastal and shelf seas. It contributes significantly to the attenuation of light in natural waters and thereby impacts ocean heat content, particularly in coastal and shelf seas (Soppa et al., 2019; Gnanadesikan et al., 2019; Kim et al., 2015, 2016, 2018; Hill, 2008). In the Baltic Sea, CDOM is prevalent and displays strong seasonal and spatial variability (Kowalczuk, 1999; Kowalczuk et al., 2006). Sources of CDOM and changes to its composition through non-conservative processes are tightly coupled to the underwater light field. These will vary with environmental conditions and phytoplankton community structure. Moreover, heterogeneity in phytoplankton pigments and other water constituents will have implications for sub-mesoscale vertical mixing and advective fluxes and thus for water temperature and density and the supply of nutrients to the surface. Understanding how the variable presence of water constituents impacts energy fluxes in the upper ocean and across the air–sea interface and the accumulative effect on the upper-ocean heat budget in shelf seas and coastal waters is of particular importance for our capacity to adequately model regional ocean climate.

1.1 Ocean radiant heating and biological production

For studies of heat transfer modulated by biological production in the upper ocean, it is important to accurately prescribe the shortwave solar radiation in the upper water column. Downward solar radiation penetrating into the upper ocean can be partitioned into the following three spectral domains: visible (UV–VIS; ∼ 0.30–0.75 µm), near-infrared (NIR; ∼ 0.75–1.3 µm) and shortwave infrared (SWIR; ∼ 1.3–3.5 µm). SWIR radiant energy plays an important role in the surface thermal structure of the water column; however, its attenuation can be considered to be invariable to changes of water constituents (Morel and Antoine, 1994) as it is almost completely dominated by water absorption and is fully attenuated very close to the sea surface. NIR radiant energy penetrates a bit deeper into the ocean but is still almost entirely absorbed within the topmost 1 m layer due to the still-strong absorption of pure sea water at these wavelengths. In contrast to that, the (spectral) attenuation of UV–VIS radiant energy within the water body is strongly dependent on the presence of water constituents and may therefore vary considerably horizontally and vertically. More specifically, the variability of UV–VIS radiant energy in the water column is determined by the absorption and scattering of optically significant water constituents, e.g. phytoplankton, detritus, CDOM and inorganic suspended sediment (Sathyendranath et al., 1989). The properties of the individual constituents determine how they absorb and scatter light in different parts of the visible spectrum; CDOM preferentially absorbs light in the blue end of the spectrum, while phytoplankton absorb light in the blue or green and red part of the spectrum – exactly how will depend on the pigment composition of the functional group (Fig. 1).

Figure 1 Spectral absorption coefficients for (a) water, relict and labile CDOM (Bissett et al., 1999b; Kowalczuk et al., 2005b) and for (b) phytoplankton pigments (Bidigare et al., 1990) used in the Bio-Optic model.

A number of feedback mechanisms determine the biogeochemical dynamics in the upper-ocean layer. Absorbed solar radiation is mostly transformed into heat and thus directly controls heating rates and subsequently impacts the vertical stratification of the euphotic layer. A portion of the light absorbed by autotrophic protists is used for photosynthesis and consequently contributes to biomass production. The vertical distribution of absorbing material may be altered significantly due to biogenic (and in coastal areas, non-biogenic) processes (e.g. by the development of a subsurface algae bloom or increased turbidity arising from sediment transport by river plumes), which in turn leads to a significant change of the depth range at which heating occurs (e.g. increased heating within the algae or turbid layer) and the availability of light (e.g. strongly reduced light availability below the algae or turbid layer).

Biogeochemical dynamics are especially complex in shelf and coastal waters where organic and inorganic particulate matter, as well as CDOM, may be present in individually highly varying concentration ranges, e.g. caused by riverine inputs or sediment resuspension from the seafloor. For example, accounting for the highly variable light attenuation in turbid river plumes is critical if nearshore physics are to be resolved correctly (Cahill et al., 2008; Kim et al., 2020). Changes in surface temperature and buoyancy-driven circulation have important consequences for the development, transport and fate of phytoplankton biomass. The resulting carbon fluxes across the air–sea interface, exported to the benthos or advected off the shelf system, are key to understanding the carbon budgets of shelf systems and the open ocean.

1.2 Biogeochemical ocean models

A number of studies in productive open-ocean waters elegantly demonstrate how upper-ocean chlorophyll concentrations regulate radiant energy transmission and heating rates in the mixed layer (Simpson and Dickey, 1981; Lewis et al., 1990; Morel and Antoine, 1994; Ohlmann et al., 1996, 1998, 2000; Ohlmann and Siegel, 2000; Dickey and Falkowski, 2002; Murtugudde et al., 2002; Oschlies, 2004; Manizza et al., 2005, 2008). Enhanced near-surface stratification can have a positive feedback on phytoplankton growth by restricting phytoplankton within shallower mixed layers with more available light, which in turn increases near-surface local heating (Dickey and Falkowski, 2002). A 10 W m⁻³ change in the solar radiation absorbed within a 10 m layer can represent a temperature change within that layer of more than 0.6 ^∘C month⁻¹ (Simpson and Dickey, 1981). However, as light limitation is replaced by nutrient limitation, increased stratification will inhibit the exchange of deeper nutrient-rich water with the surface and limit phytoplankton growth. Ohlmann et al. (2000) demonstrated that an increase in chlorophyll concentration from 0.03 to 3 mg m⁻³ in the upper 10 m of the water column can decrease the solar flux in the waters below by as much as 35 W m⁻².

A few studies have tried to explore the full biophysical feedbacks using coupled physical–biological ocean models (Oschlies, 2004; Manizza et al., 2005, 2008) and fully coupled atmosphere–bio-physical ocean models (Jolliff and Smith, 2014; Wetzel et al., 2006). Notably, results from Oschlies (2004) include a net cooling of the North Atlantic by biota of about 1 W m⁻², with enhanced upper-ocean stratification in summer and deeper winter mixed-layer depths (> 100 m) in parts of the subpolar gyre. Coastal upwelling and associated nutrient supply are reduced, especially in coastal upwelling regions of West Africa. Overall, there is a negative feedback of biotically induced radiative heating on chlorophyll a concentrations, except in parts of the subpolar North Atlantic where intensification of the spring bloom results in increased annual mean chlorophyll a concentrations. Wetzel et al. (2006) further highlighted the importance of marine biology to the radiative budget of the upper ocean and found that positive feedbacks with the climate system cause a global shift of the seasonal cycle, with the onset of spring occurring about 2 weeks earlier. Increased wind stress and changes in the shortwave radiation led to significant warming in the middle latitudes in summer and to seasonal modifications of the overall warming in the equatorial Pacific. Jolliff and Smith (2014) demonstrated a regional example of biological modulation of upper-ocean physics in Monterey Bay, California, and showed how the spatiotemporal pattern of a phytoplankton bloom can persist due to enhanced thermal stratification promoting vertical stability and more efficient use of macronutrients. Furthermore, biothermal warming of surface waters modifies the local surface pressure gradient and modulates wind stress patterns.

More recent studies which investigate the role of OSCs and surface heating highlight the role of CDOM in Arctic amplification (e.g. Soppa et al., 2019; Pefanis et al., 2020) and the impact of CDOM on the annual cycle of sea surface temperature in coastal and northern subpolar regions (Gnanadesikan et al., 2019; Kim et al., 2015, 2016, 2018). Soppa et al. (2019) found that a CDOM absorption at 443 nm of 1.77 m⁻¹ contributed to an increased radiative heating of 0.6 ^∘C d⁻¹ in the upper 2 m in the Laptev Sea shelf waters, implying increased sea ice melt rates and changes in the surface heat fluxes to the atmosphere. Pefanis et al. (2020) confirm that increases in CDOM in the Arctic amplify surface warming by increasing surface temperatures in summer and decreasing sea ice concentrations. They also show that summertime surface warming associated with increases in CDOM induces more heat loss to the atmosphere, primarily through latent and sensible heat fluxes. Gnanadesikan et al. (2019) demonstrate that the presence of CDOM leads to an increase in the amplitude of the seasonal cycle of sea surface temperature (SST) over coastal and northern subpolar regions, with potential implications for extreme ocean temperatures. Importantly, they find the size and sign of the change in amplitude are controlled by the interplay between enhanced surface shortwave heating, shading and cooling of the subsurface and the extent to which these are connected by vertical mixing. They show that the interplay between heat term balances varies regionally. In the central Baltic Sea (58^∘ N, 19.5^∘ E), changes in the seasonal cycle of the heat budget are explained by a 1D balance between the penetration of shortwave radiation and vertical mixing (see Fig. 3a in Gnanadesikan et al., 2019), with advective and diffusive terms being relatively small. In other regions around the world, the heat term balance is represented by a more complicated interplay between the penetration of shortwave radiation, vertical and horizontal mixing, and advection (see Fig. 3b, c, d in Gnanadesikan et al., 2019). Löptien and Meier (2011) show that increased water turbidity affects the summer sea surface temperature trends in the Baltic Sea significantly. While Skákala et al. (2022) demonstrate a significant impact of biogeochemistry on physics in the northwest European shelf, with the light attenuation by chlorophyll being responsible for a 1 ^∘C warming in the upper 20 m of the ocean, with comparable cooling taking place between 20 and 200 m. They also show that accounting for this water-constituent-induced heating improves the timing of the simulated phytoplankton bloom in the region.

Despite these findings, coupled ecosystem–circulation models rarely share the same parameterisation or source of radiative forcing to drive the hydrodynamics and fuel photosynthesis even though their requirements for information on light and heat overlap. This is in part due to the fact that, historically, circulation and ecosystem models have evolved independently, and it is only in the last 10 to 15 years that coupling between the two has made significant advances. It is typical that the ecosystem model is plugged into a circulation model and that communication between the two is in one direction only: state variables (such as temperature) computed in the circulation model are communicated to the biological model at each time step; however, any change to the radiative fluxes as a consequence of biological activity is not necessarily accounted for or communicated back to the circulation model so that potentially available information related to heat transfer in the upper ocean and across the ocean–atmosphere interface is not used.

Many parameterisations of the subsurface vertical distribution of shortwave solar radiation in ocean models have evolved over the last years (e.g. Paulson and Simpson, 1977; Zaneveld and Spinrad, 1980; Simpson and Dickey, 1981; Morel, 1988; Morel and Antoine, 1994; Ohlmann and Siegel, 2000; Manizza et al., 2008). For photosynthesis purposes, one of the more simple parameterisations of light attenuation is based on the surface photosynthetically available radiation (PAR) computed as a fraction of the net surface solar flux (typically 43 %) and then attenuated through the water column as a function of chlorophyll concentration (e.g. Fasham et al., 1990; Fennel et al., 2006). Zielinski et al. (2002) compared the effect of some different light parameterisations in biogeochemical models on primary production and phytoplankton evolution in the subtropical North Atlantic and showed that there can be significant changes in the vertical distribution of simulated phytoplankton depending on how the underwater light field is treated.

Chlorophyll-based approaches to underwater light attenuation are reasonably accurate for the open ocean, where phytoplankton dominates the inherent optical properties of the water constituents (Morel and Prieur, 1977); however, they are inadequate in shelf and coastal oceans as they neglect important contributions from CDOM, detritus and suspended sediments. Neumann et al. (2015) showed that, in the Baltic Sea, including more water constituents in the estimation of light attenuation in their model yields a more realistic representation of the light climate and improved estimates of primary productivity, Secchi disc depth and oxygen concentrations. They estimated light attenuation by explicitly accounting for modelled phytoplankton biomass, detritus and dissolved organic matter due to metabolism and degradation processes and by parameterising CDOM as a function of salinity. More recently, Neumann et al. (2021) showed that explicitly considering light absorption due to terrestrial CDOM in their ecosystem model of the Baltic Sea significantly improved CDOM estimates, particularly in the northern parts of the Baltic Sea, where the impacts of terrestrial CDOM are large. Including directional and spectral light in coupled biogeochemical–circulation–radiative models has been shown to be important for ocean biology, especially for studies of community structure and succession (Gregg and Rousseaux, 2016). It is also important for regional studies which examine the role of other optical constituents such as CDOM and detritus in carbon cycling (Bissett et al., 1999a, b).

1.3 Estimating the impact of optically significant water constituents on surface heating in the western Baltic Sea

In this work, we use a spectrally resolved underwater light field to explore the relationship between OSCs – in particular, CDOM, phytoplankton and detritus – and heating rates in the western Baltic Sea. High concentrations of CDOM optically distinguish the Baltic Sea from other coastal seas (Simis et al., 2017), making it an interesting study site for this application. CDOM also exhibits strong seasonal and spatial variability in the region, which is dependent on sources of CDOM and physics, e.g. periods of intensive mixing and high riverine discharge versus periods of thermal stratification, reduced riverine discharge, enhanced biological production and production of CDOM (Kowalczuk, 1999; Kowalczuk et al., 2005a). This interplay between physics and OSCs is examined using a coupled bio-optical ocean model which incorporates the optical properties of key water constituents and explicitly resolves sources of both terrestrial and autochthonous CDOM as a state variable in a 4D ocean state. We model the inherent optical properties of different water constituents under varying environmental conditions and spectrally resolve the underwater light field in a dynamic ocean. From this, we estimate the contribution of key water constituents to surface heating rates and feedbacks with the marine atmospheric boundary layer heat fluxes. Modelled inherent and apparent optical properties are evaluated with in situ and satellite observations, and estimates of surface heating rates are compared with those derived from an ocean–atmosphere radiative transfer model which accounts for both the directionality and spectral dependence of the underwater light field.

2 Methods

2.1 Study site

Kowalczuk et al. (2006) have shown that there are three pools of CDOM in the waters of the southern Baltic Sea: a riverine pool, an aged marine pool and a pool primarily produced in offshore waters. They explored the seasonal dependence between the light absorption coefficient of CDOM at 375 nm, aCDOM(375), and salinity and chlorophyll a concentrations in the southern Baltic Sea and found a seasonal dependence between physical processes and the source of CDOM. In March, April and November, months of intensive mixing and high riverine discharge, most of the variability in aCDOM(375) values could be explained by dilution of terrestrially derived CDOM alone. In February, May and September, months of thermal stratification, reduced riverine discharge and enhanced biological activity, autochthonous production of CDOM was found to be a significant source of CDOM in the southern Baltic Sea. Changes in the values of spectral slope coefficients are regarded as an indicator of compositional changes in CDOM. These changes can be a result of either conservative mixing processes, i.e. mixing, or non-conservative processes, e.g. production, degradation or flocculation (Kowalczuk et al., 2006).

Figure 2 Western Baltic Sea model domain bathymetry (m) with location of model output analysis stations Darß Sill (DS), Arkona Sea (AS), Oder Bank (OB) and Bornholm Basin (BB) (blue dots) and with in situ CDOM and NAP (non-algal particle) absorption measurements from the Institute of Oceanology of the Polish Academy of Sciences, IOPAN (red dots).

Our study site in the western Baltic Sea (Fig. 2) includes the Bornholm Basin, where we expect the seasonal cycle to be explained by a 1D balance between the penetration of shortwave radiation and vertical mixing (Gnanadesikan et al., 2019), and the Darß Sill, Arkona Sea and Oder Bank, where advection and diffusion will also contribute to the seasonal heat balance, making for an interesting contrast between local regimes. At the Bornholm Basin, we expect to find marine CDOM; at the Darß Sill and Arkona Sea, we expect to find a mixture of riverine and marine CDOM, depending on the season, while at the Oder Bank, we expect the CDOM pool to be dominated by riverine sources from multiple inlets and rivers connecting the Oder River outlet through Szczecin Lagoon with the Greifswalder Bodden and the coastal Baltic Sea (Kowalczuk et al., 1999).

2.2 Model system

The coupled modelling system has two components: the Regional Ocean Modelling System, ROMS, which drives the physics and the advection and diffusion of tracers, and Ecosim/Bio-Optic, which drives the ecosystem and underwater light field. These components interact as shown in Fig. 3 and are described in more detail below.

Figure 3 Model system components and how they interact.

Light penetrating a water body can be described as consisting of three streams (Aas, 1987; Ackleson et al., 1994; Gregg, 2002; Dutkiewicz et al., 2015). These are the downward direct irradiance, E_dir; the downward diffuse irradiance, E_diff; and the upward diffuse irradiance, E_u. E_dir+E_diff is commonly referred to as downward irradiance, E_d. For studies of heat transfer and photosynthesis, we need to know the scalar irradiance, E₀, which describes the light field integrated over a sphere and is thus independent of direction. All of these irradiance quantities (E_dir, E_diff, E_u and E₀) are a function of wavelength and depth.

Following Morel (1988), the rate of radiant energy converted into heat can be estimated as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=-{\displaystyle \frac{\mathrm{d}\left(E_{\mathrm{d}}-E_{\mathrm{u}}\right)}{\mathrm{d}z}}{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }{C}_p}},\end{array}$$

where the term on the right-hand side is the heat flux; E_d and E_u are the downward and upward irradiances, respectively; ρ is the in situ density; and C_p is the specific heat capacity of water. In a horizontally homogeneous water body, the divergence of the radiative flux can be approximated as follows (Morel, 1988):

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\mathrm{d}\left(E_{\mathrm{d}}-E_{\mathrm{u}}\right)}{\mathrm{d}z}}\cong -a{E}_{\mathrm{0}}\approx K_{\mathrm{d}}{E}_{\mathrm{d}},\end{array}$$

where a is the local absorption coefficient, E₀ is the scalar irradiance at the depth in question, and K_d is the downward diffuse attenuation coefficient for downwelling irradiance. These quantities are all dependent on depth, concentrations of OSCs (e.g. phytoplankton pigments, CDOM, detritus) and wavelength. Thus,

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=-{\displaystyle \frac{{\int }_{\mathrm{400}}^{\mathrm{700}}\left[E_{\mathrm{d}}\left(\mathit{\lambda },z\right)K_{\mathrm{d}}\left(\mathit{\lambda },z\right)\right]\mathrm{d}\mathit{\lambda }}{\mathit{\rho }{C}_p}}.\end{array}$$

K_d varies with both absorption, a, and scattering, b, as well as with the angular distribution of the incoming light field. It can be calculated from E_d as follows (Gordon et al., 1980):

$$\begin{array}{}\text{(4)}& K_{\mathrm{d}}={\displaystyle \frac{-\mathrm{d}\ln E_{\mathrm{d}}\left(\mathit{\lambda },z\right)}{\mathrm{d}z}}={\displaystyle \frac{-\mathrm{1}}{E_{\mathrm{d}}\left(\mathit{\lambda },z\right)}}{\displaystyle \frac{\mathrm{d}{E}_{\mathrm{d}}\left(\mathit{\lambda },z\right)}{\mathrm{d}z}}.\end{array}$$

Biogeochemical–optical relationships vary significantly over different regions and/or seasons; therefore, regional and temporal relationships have been adopted to cope with such variations when information concerning the directionality of the underwater light field is limited. For example, in open-ocean waters, where attenuation of underwater light is primarily a function of chlorophyll concentration, Sathyendranath and Platt (1988) parameterise K_d as follows:

$$\begin{array}{}\text{(5)}& K_{\mathrm{d}}={\displaystyle \frac{a+b}{{\mathit{\mu }}_{\mathrm{0}}}},\end{array}$$

where a is the absorption and b is the total scattering (forward and backscatter) of OSCs, while μ₀ is the average cosine, which tells you how much the light field differs from isotropic conditions.

In more complex coastal waters, Lee et al. (2005) have derived an empirical algorithm to parameterise K_d as follows:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{K}_{\mathrm{d}}\left(\mathit{\lambda },z\right)& =\left(\mathrm{1}+\mathrm{0.005}\mathit{\theta }\right)a\left(\mathit{\lambda },z\right)\\ & +\mathrm{4.18}\left(\mathrm{1}-\mathrm{0.52}{e}^{-\mathrm{10.8}a\left(\mathit{\lambda },z\right)}\right)b_{\mathrm{b}}\left(\mathit{\lambda },z\right),\end{array}\end{array}$$

where θ is the solar zenith angle in degrees, and b_b is the backscatter coefficient.

If the absorption and scattering properties of different water constituents are known, K_d can be estimated using Eq. (5) or (6), and E_d can then be calculated using Eq. (7).

$$\begin{array}{}\text{(7)}& E_{\mathrm{d}}=E_{\mathrm{d}}\left(\mathrm{0}\right)e^{-K_{\mathrm{d}}z}\end{array}$$

Thus, the heat balance relationship described in Eq. (3), can be used to estimate heating rates.

2.2.1 Regional Ocean Modelling System, ROMS, and Ecosim/Bio-Optic

The ocean model component, ROMS, is widely used for shelf circulation (e.g. Haidvogel et al., 2008; Wilkin et al., 2011) and coupled physical–biological applications (e.g. Cahill et al., 2008, 2016; Fennel et al., 2006, 2008; Fennel and Wilkin, 2009). The ROMS computational kernel (Shchepetkin and McWilliams, 2005) produces an accurate evolution of tracer fields, which is a particularly attractive feature for biogeochemical modelling because it facilitates the correct interaction among tracers and the accounting of total nutrient and carbon budgets. ROMS is coupled to Ecosim, the carbon-based, ecological–optical modelling system (Bissett et al., 1999a, b) which was developed for simulations of carbon cycling and biological productivity. Ecosim simulates up to four phytoplankton functional groups, with a characteristic pigment suite which varies with the group's carbon-to-chlorophyll-a ratio (C : Chl a). The properties of each functional group evolve over time as a function of light and nutrient conditions (i.e. NO₃, NH₄, PO₄, SiO and FeO). Marine and riverine sources of dissolved organic carbon (DOC and CDOC) are accounted for and explicitly resolved into labile (e.g. available for biological and photo-degradation) and relict (e.g. available for photo-degradation) forms. Dissolved inorganic carbon is also accounted for. Riverine sources of carbon and nutrients are introduced via point sources. The underwater light field is spectrally resolved between 400 and 700 nm, which allows for differential growth of different phytoplankton groups that have unique pigment complements. The interaction between Ecosim's components describes autotrophic growth of and competition between phytoplankton groups, differential carbon and nitrogen cycling, nitrogen fixation, and grazing. Coupled ROMS–Ecosim applications include a deployment in the New York–New Jersey sea bight, which demonstrates the feedback of turbid buoyant plumes originating from the Hudson River on near-shore biogeochemistry and physics (Cahill et al., 2008).

Ecosim contains a daylight module which is central to this work. Light energy just beneath the sea surface is calculated using a derivative of the RADTRAN code described in Gregg and Carder (1990) as a function of the model's meteorological forcing (i.e. wind speed, relative humidity, air temperature and pressure) and cloud cover, atmospheric gases (i.e. water vapour, ozone and oxygen), marine aerosols and the surface roughness and reflectance at the ocean–atmosphere interface. A constant fraction of 0.3 cloud cover is assumed for clouds, while 1.5 cm precipitable water is assumed for water vapour. The underlying algorithms used to compute ozone, water vapour and oxygen absorption coefficients are described in detail in Gregg and Carder (1990). Marine aerosols are computed according to the simplified version of the Navy marine aerosol model, also described in detail in Gregg and Carder (1990). The surface solar downwelling spectral irradiance, $E_{\mathrm{d}}(\mathit{\lambda },\mathrm{0}-)$ (which is the sum of the direct and diffuse irradiance), and the average cosine zenith angle, ${\mathit{\mu }}_{\mathrm{0}}(\mathit{\lambda },\mathrm{0}-)$, are provided at 5 nm wavelength intervals between 400 and 700 nm and are used as inputs to Ecosim's daylight module.

The spectrally resolved downward light stream, E_d(λ,z), is calculated according to Eq. (10) and is attenuated by absorption, a, and scattering (forward, b, and backward, b_b), of the OSCs. Phytoplankton and detritus both absorb and scatter light. Phytoplankton absorption is calculated for the four functional groups as a function of biomass, weight-specific pigment absorption coefficients (Fig. 1b, Bidigare et al., 1990) and packaging effect (Bissett et al., 1999b; Kirk, 2011). Detrital absorption is calculated as an exponential function of wavelength (Gallegos et al., 2011). Phytoplankton and detrital scattering and backscattering are accounted for as total particulate scattering and backscattering according to Morel (1991, 1988), respectively (see Eqs. 16 and 17 in Bissett et al., 1999b). CDOM only absorbs light and is calculated as a function of CDOM concentration and the weight-specific absorption coefficients adapted from Kowalczuk et al. (2005b) (Fig. 1a). The average cosine is modified with depth as a function of absorption and backscattering. This is simplified as a linear function of the optical depth between two levels (see Eq. 22 in Bissett et al., 1999b). The total scalar irradiance, E₀(λ,z), which is the light available to phytoplankton, is calculated following Eq. (5) after Morel (1988).

Bio-Optic is a new option within Ecosim's daylight module which adds some diagnostics and functionality. These are

• the explicit output of inherent optical property diagnostics (absorption, scatter and backscatter) of each of the OSCs (i.e. phytoplankton, detritus and CDOM) and apparent optical property diagnostics (downward attenuation; downward and scalar irradiance fields; surface solar downwelling spectral irradiance, $E_{\mathrm{d}}(\mathit{\lambda },\mathrm{0}-)$; and the average cosine zenith angle, ${\mathit{\mu }}_{\mathrm{0}}(\mathit{\lambda },\mathrm{0}-)$)

• an option to calculate a downwelling irradiance attenuation coefficient, K_d, which accounts for some of the optical complexity found in coastal waters, according to Lee et al. (2005)

• an option to couple the bio-optically calculated downward-irradiance term back into the hydrodynamic solution.

Bio-Optic is activated as an option within Ecosim during compilation.

The explicit calculation of in-water spectrally resolved absorption, scattering and backscattering coefficients, average cosine, and downwelling irradiance attenuation coefficient, K_d, in addition to the scalar, E₀, and downward, E_d, irradiance fields, has important implications. The spectrally resolved underwater light field drives the evolution of OSCs in the ecosystem model, while the OSCs in turn determine the evolution of the light field in each layer by absorption and scattering of the light. This means that the OSCs' contribution to the divergence of the heat flux (Morel, 1988) can be accounted for within the full hydrodynamic solution. Furthermore, water-constituent-induced heating rates can be assessed, and their impact on the ocean sea surface temperature can be communicated to the bulk flux formulation of the atmosphere in the modelling system. While this still represents a very simplified treatment of radiative transfer within the water column, it does permit a direct evaluation of the optical terms and heating rates with those derived from a full solution of the radiative transfer equations and provides a means to improve the parameterisation of water-constituent-based heat flux algorithms in ocean models. For this purpose, we use the vector radiative transfer model MOMO (Matrix Operator Model, described below) to evaluate the more approximate solution provided by ROMS–Bio-Optic.

2.2.2 Vector radiative transfer model, MOMO

A more rigorous treatment of the vertical structure of the light field is provided by atmosphere–ocean radiative transfer models, such as MOMO (Fell and Fischer, 2001), which simulate the light field in the stratified atmosphere–ocean system for the VIS and NIR spectral ranges. MOMO uses the matrix operator method to calculate zenithally and azimuthally resolved light fields for different types and concentrations of optically active components in the ocean and atmosphere; thus, the full directionality of the light field is accounted for. The main advantage of the matrix operator method is its efficiency in simulating light propagation in optically dense media. It is therefore particularly suited to being used in the development of remote sensing algorithms for the retrieval of water constituents. It has most recently been described in Hollstein and Fischer (2012) and is based on previous work by Fischer and Grassl (1984) and Fell and Fischer (2001). It has been successfully applied to the analysis of hyperspectral ocean colour data to derive surface fluorescence signals (Guanter et al., 2010), the analysis of ocean colour data from MERIS measurements (Zhang et al., 2003) and a new retrieval of sun-induced chlorophyll fluorescence in water from ocean colour measurements (Kritten et al., 2020). For our purposes, the most pertinent elements of MOMO include the calculation of the spectrally resolved downward surface irradiance for the VIS and NIR ranges, the direct and diffuse downwelling, and the diffuse-upwelling components of the underwater light field.

2.3 Experimental setup

The ROMS Ecosim/Bio-Optic modelling system was configured for the western Baltic Sea (Fig. 2) with a horizontal resolution of ∼ 1.8km (285×169 grid points) and 30 terrain-following σ levels in the vertical. A bulk flux atmosphere was forced with DWD-ICON output (Zängl et al., 2015), and river forcing, including runoff and biogeochemistry, was derived from HELCOM (Helsinki Commission) PLC (pollution load compilation) data (Thomas Neumann, personal communication, 2020). Open boundaries to the north and east were forced with output from GETM (General Estuarine Transport Model) physics using a combination of Chapman and Flather conditions for u and v velocities and transports and radiation + nudging for temperature and salinity. This 3D setup is based on an existing GETM physics setup which has been previously evaluated and published (Gräwe et al., 2015a, b). It captures the annual cycle of temperature and salinity in the western Baltic Sea and episodic inflows of saline, oxygen-rich North Sea water which control the salinity content and stratification in the Baltic Sea and are important for ventilating the deeper basins of the Baltic Sea (Omstedt et al., 2004; Meier, 2007).

Ecosim was configured with four phytoplankton functional groups representative of small and large diatoms, large dinoflagellates, and cyanobacteria. Two experiments covering the period 1 January 2018 to 31 December 2018 were carried out as follows:

1. three-dimensional western Baltic Sea – feedback of constituent-induced heating to hydrodynamic solution (herein referred to as biofeed)

2. three-dimensional western Baltic Sea – no feedback of constituent-induced heating to hydrodynamic solution (herein referred to as nobiofeed).

MOMO simulations were performed at relatively high angular resolution (27 angles in the atmosphere between 0 and 88^∘ plus 9 additional angles in the ocean to cover the angular domain of total internal reflection) to allow for an accurate calculation of the in-water light field. Up to 120 terms were used for the Fourier expansion of the azimuth dependence of the light field. The oceanic vertical structure in MOMO has been chosen to be identical to the ROMS–Bio-Optic vertical structure; i.e. the light field has been calculated at the 30 ROMS–Bio-Optic layer boundaries located between 0 and ca. 90 m. Absorption and scattering coefficients for phytoplankton, CDOM and detritus are taken directly from ROMS–Bio-Optic output. Spectral resolution was done in steps of 5 nm between 400 and 700 nm. Two Fournier–Forand phase functions (Fournier and Forand, 1994; Freda and Piskozub, 2007) with differing backscattering-to-scattering ratios have been applied to phytoplankton ($b_{\mathrm{b}}/b=\mathrm{0.001}$) and detrital material ($b_{\mathrm{b}}/b=\mathrm{0.1}$), in line with phase functions measured by Siegel et al. (2005) for various Baltic Sea coastal waters. Seasonal heating rates were derived from MOMO simulations at the Bornholm Basin location and compared to the corresponding fluxes from ROMS–Bio-Optic in order to assess the suitability of the simplified treatment of radiative transfer in the latter and the implications of not resolving the full directionality of the light field therein. MOMO results are presented for the 38^∘ solar incident zenith angle, representative of late spring to mid-summer in the western Baltic Sea (Fig. 11).

2.4 Model evaluation strategy and supporting data

Evaluation of our model output was carried out primarily at the Oder Bank, Darß Sill, Arkona Sea and Bornholm Basin sites within our model domain. These have been previously discussed in Sect. 2.1 and are shown as blue dots in Fig. 2.

The following three aspects of our model results were examined:

1. The seasonal cycle of modelled temperature versus observations at four locations was considered. Darß Sill and Arkona Sea mooring data, shown in Fig. 4, middle panel, were obtained from the BSH (Bundesamt für Seeschifffahrt und Hydrographie) MARNET mooring database. SST data, shown in the right panel of Fig. 4, were obtained from the NOAA OI SST v2 high-resolution dataset (Huang et al., 2021).

2. Model surface chlorophyll a, phytoplankton and non-algal particulate absorption at 443 nm and the diffuse attenuation coefficient at 490 nm are compared with the Sentinel 3 Ocean and Land Colour Instrument, OLCI, level-3 300 m data products (https://doi.org/10.48670/moi-00294) on 2 consecutive clear days in May 2018, when a bloom event occurred. Modelled monthly mean CDOM absorption is compared with MERIS-derived and in-situ-measurement-derived seasonal climatologies (see the Supplement for details). Seasonal phytoplankton and non-algal particle absorption (CDOM + detritus) at 440 and 442 nm are compared with seasonal estimates from Meler et al. (2016a).

3. Heating-rate estimates at Bornholm Basin derived from ROMS–Ecosim/BioOptic diagnostic calculations are compared with heating-rate estimates derived from comparable full radiative transfer calculations using MOMO.

3 Results

In Sect. 3.1, we show the results from the biofeed experiment, which includes the feedback from OSC-induced heating to the hydrodynamic solution. In Sect. 3.2, we show the difference between the biofeed experiment and the nobiofeed experiment, where no feedback from OSC-induced heating is included in the hydrodynamic solution.

3.1 Seasonal cycle of temperature at Oder Bank, Darß Sill, Arkona Sea and Bornholm Basin in western Baltic Sea

The modelled versus observed annual cycle of temperature at the different locations is shown in Fig. 4. High-resolution temporal and vertically resolved observations for 2018 were only available at the Darß Sill and Arkona Sea sites (middle plots, Fig. 4). Oder Bank and Darß Sill are shallow, well-mixed locations, where seasonal warming and cooling of the whole water column takes place between May and October. At the deeper Arkona Sea and Bornholm Basin locations, the onset of seasonal stratification occurs in early May and starts to break down in September. Intense summertime warming in late July and early August (SST ∼ 25 ^∘C) leads to a deepening of the thermocline from ca. 20 m to the seafloor in the Arkona Sea and to ca. 38 m in the Bornholm Basin. In the Arkona Sea, the model captures observed summertime baroclinic inflows between 15 and 30 m depth. These inflows are intrusions of deep, saltier, cool water which are pushed over the Drogen and Darß sills into the deeper Arkona Sea. Due to the estuarine nature of the Baltic Sea circulation, these inflows are not unusual in the western Baltic Sea (Fennel and Sturm, 1992). Overall, there is very good agreement between the modelled biofeed results and observed temperature fields at all locations, especially in terms of the sea surface temperature (see Table 1 for r², RMSE and bias statistics). This is especially important as 2018 was a year where two significant marine heat waves (defined as periods where the surface temperature exceeds the 90th percentile of the 30-year local mean for longer than 5 d) took place in May–June (38 d) and July–August (17 d). This result confirms the importance of accounting for the contribution of OSCs to the transfer of light energy.

Figure 4 Modelled (left) versus observed (middle; note that the white triangles in the Arkona Sea observation plot indicate periods where observations are missing from the time series) annual cycle of temperature and sea surface temperature (right) in 2018 at the Oder Bank, Darß Sill, Arkona Sea and Bornholm Basin sites (Legend abbreviations: ROMS – model output; CLIM – 30-year climatological mean calculated from OI SST dataset; 90th – 90th percentile of the 30-year climatological mean (CLIM); OI SST – 2018 daily optimum interpolation sea surface temperature (Huang et al., 2021).)

Table 1 Model versus observed sea surface temperature (^∘C) statistics.

3.2 Inherent and apparent optical properties of OSCs at Oder Bank, Darß Sill, Arkona Sea and Bornholm Basin in western Baltic Sea

OLCI level-3 products of chlorophyll a, phytoplankton and non-algal particle absorption at 443 nm and the diffuse attenuation coefficient at 490 nm at 300 m resolution were used to evaluate our modelled equivalents. We chose 2 days in May 2018 where full satellite data coverage was available and which coincided with peak OSC-induced heating rates found in our model results. Figure 5 shows modelled chlorophyll a, phytoplankton and non-algal particle absorption at 443 nm and the diffuse attenuation coefficient at 490 nm and related RMSE values. The white cross marks on the plots represent the position of the different analysis locations where matchups between the OLCI data and our model output have been extracted.

Figure 5 Modelled mean (29 and 30 May 2018) chlorophyll a (a), phytoplankton absorption at 443 nm (c), non-algal particle absorption at 443 nm (e) and diffuse attenuation coefficient at 490 nm (K_d490) (g) and related RMSE (b, d , f, h).

Table 2 OLCI versus model matchup mean values (29 and 30 May 2018) for Chl a, phytoplankton (aPhy) and non-algal particle (aNAP) absorption at 443 nm and the diffuse attenuation coefficient at 490 nm (K_d490).

The matchups (Table 2) highlight how we can only reasonably compare OLCI and model output at the Darß Sill and Bornholm Basin locations as the bloom event evident in the OLCI data in the Arkona Sea and Oder Bank (Fig. 5) is not fully captured in the model. At these locations, Chl a and NAP absorption are underestimated by the model by as much as 6 mg m⁻³ and 0.27 m⁻¹, respectively. Phytoplankton absorption is slightly overestimated in the model at all locations, but the values are in better agreement with the OLCI data (within 0.1 m⁻¹ difference range), as are the modelled non-algal particle absorption values at the Darß Sill and Bornholm Basin sites (within 0.03 m⁻¹ difference range). Modelled K_d490 also compares reasonably well with the OLCI data at all locations (within 0.2 m⁻¹ difference range). We do not expect the model to capture the dynamic bloom event observed by OLCI without further tuning or data assimilation. As it stands, there is good agreement of the model and OLCI data with the background values at Darß Sill and, especially, at Bornholm Basin, which give us confidence in the model performance and supports the selection of Bornholm Basin for further evaluation of the heating rates and air–sea fluxes.

We also compared modelled monthly mean CDOM absorption with MERIS-derived and in-situ-derived climatologies, as well as seasonal phytoplankton and non-algal particle absorption with seasonal estimates from Meler et al. (2016a). Modelled monthly mean surface CDOM absorption is underestimated compared to the MERIS-derived climatological CDOM absorption (Fig. 6b) (with r² ranging from 0.35 to 0.66 and RMSE ranging from 0.19 to 0.1 at Oder Bank and Bornholm Basin, respectively) but is in better agreement with the seasonal observed estimates of Meler et al. (2016a) (Fig. 6c) (r²=0.7 and 0.64 and RMSE = 0.05 and 0.1 for non-algal particle absorption and phytoplankton absorption, respectively) (Fig. 6d).

Figure 6 (a) MERIS and in situ monthly climatology of surface CDOM absorption (mean value calculated over western Baltic Sea region shown in Fig. 2); (b) mean monthly surface CDOM absorption at model stations and matching MERIS locations; seasonal mean surface non-algal particle absorption (CDOM + detritus) (c) and phytoplankton absorption (d) at model stations compared with similar water type values found in Meler et al. (2016a).

Modelled spectrally resolved surface phytoplankton, CDOM and detritus absorption at the Oder Bank, Darß Sill, Arkona Sea and the Bornholm Basin sites (Fig. 7) show typical absorption characteristics for the individual constituents. CDOM and detritus have high absorption values at the blue end of the spectrum, while phytoplankton shows two maxima, one between 440 and 490 nm and a smaller one around 670 nm. There is a clear seasonal pattern for each of the constituents, with spring and summer being peak seasons for phytoplankton blooms and summer and autumn favouring increased CDOM and detrital absorption. Considerable variability in absorption characteristics is evident between the locations. The highest absorption for all the constituents is seen at the coastal Oder Bank location, which is strongly influenced by riverine inputs from the Oder River. There is a decreasing gradient, especially in CDOM and detrital absorption, moving from the coastal zone to the offshore regions. The summer phytoplankton bloom in the Arkona Sea has a higher peak than the Darß Sill. CDOM-, detritus- and phytoplankton-specific absorption curves intersect around 442 nm, making this an interesting wavelength to explore further with respect to the impact of these constituents on the vertical distribution of absorption and the downward attenuation and irradiance fields. The vertical profiles of phytoplankton, CDOM and detrital absorption at 442 nm (Fig. 8) show the vertical extent of water constituent absorption to be the full water column at Oder Bank and Darß Sill and between 15 and 20 m depth in the Arkona Sea and Bornholm Basin. In spring and especially in summer, phytoplankton dominate sub-surface absorption at all locations, followed by CDOM and then detrital absorption.

The spectrally resolved surface downward attenuation (K_d) and downward irradiance (E_d) at each of the locations shown in Fig. 9 reflect the seasonal impact of the water constituent absorption and solar irradiance. Irradiance at the surface peaks in summer and is at its lowest in winter, as expected. The slight modification of downwelling irradiance intensity in the Baltic Sea depends on atmospheric conditions. Results of direct measurements and local parameterisations of radiative transfer models summarised by Dera and Woźniak (2010) (and initially reported by Rozwadowska and Isemer, 1998, and Isemer and Rozwadowska, 1999) indicate that observed monthly averaged solar irradiance intensities at the sea surface in the Baltic Sea are always lower than model estimates based on the clear-sky assumption. Atmospheric conditions have a regional and seasonal impact on observed solar irradiance entities – e.g. in the southern Baltic Proper and western Baltic Sea, the long-term monthly average for E_d at the surface in May is lower only by 4.8 and 1.8 W m⁻², respectively, than the E_d intensity observed in June in both regions. This is caused by the much lower cloud cover over the Baltic Sea observed in May than in June. Our monthly mean modelled surface irradiances converge with those reported in Dera and Woźniak (2010) (see Fig. S3 in the Supplement). We applied a constant fraction of 0.3 cloud cover, while in Dera and Woźniak (2010), the clear-sky assumption was applied. This would explain why our irradiances are lower than those of Dera and Woźniak (2010), especially in May, June and July.

Figure 7 Surface spectral phytoplankton, CDOM and detrital absorption at the Oder Bank, Darß Sill, Arkona Sea and Bornholm Basin sites in 2018 from the ROMS–Bio-Optic 3D western Baltic Sea model experiment.

Figure 8 Vertical structure of phytoplankton, CDOM and detrital absorption at 442 nm at the Darß Sill, Arkona Sea, Oder Bank and Bornholm Basin sites in 2018 from ROMS–Bio-Optic 3D western Baltic Sea model experiment.

Figure 9 Surface spectral downward diffuse light attenuation and downward irradiance at the Oder Bank, Darß Sill, Arkona Sea and Bornholm Basin sites in 2018 from ROMS–Bio-Optic 3D western Baltic Sea model experiment.

Variability in the surface layer attenuation is greatest between 400 and 550 nm, especially during the stratified spring, summer and autumn seasons, reflecting the seasonal dynamics of phytoplankton, CDOM and detritus. Vertical profiles of K_d and E_d at 442 nm (Fig. 10) show light penetrating deeper in winter, indicating relatively well-mixed (clear) waters, contrasted by seasonally stratified waters in spring, summer and autumn. Variability between the locations is also much higher during these seasons, revealing the different influence of constituents at these locations, for example, the impact of the spring and summer phytoplankton blooms in the Oder Bank and Arkona Sea on attenuation. High attenuation values at the red end of the spectrum are mostly related to the absorption of pure water itself.

Figure 10 Vertical structure of downward diffuse light attenuation and downward irradiance at 442 nm at the Oder Bank, Darß Sill, Arkona Sea and Bornholm Basin sites in 2018 from ROMS–Bio-Optic 3D western Baltic Sea model experiment.

It should be noted that seasonal and spatial variabilities in the concentration of optically significant water constituents not only impact the penetration of solar energy into the water column but also influence the spectral properties of the underwater light field. Elevated absorption by CDOM and phytoplankton pigments in the spring and summer at the Oder Bank, Darß Sill and Arkona Sea sites causes a red shift in the solar irradiance maximum transmission waveband to 570 nm from the 500 nm estimated for the Bornholm Basin (Fig. 9). This is consistent with observations reported by Kowalczuk et al. (2005a), who reported a shift in the solar irradiance maximum transmission waveband from 550 nm in the Baltic Proper to 575 nm in Pomeranian Bay and the Gulf of Gdansk. An even bigger shift in the solar irradiance maximum transmission waveband was observed between Atlantic Ocean coastal water off the west coast of Ireland (maximum solar irradiance transmission at 490 nm) and the Baltic Sea in the Gulf of Gdansk (maximum solar irradiance transmission at 570 nm). This shift was attributed to elevated CDOM absorption, which was ca. 2 times higher in the Baltic Sea compared to in the coastal Atlantic Ocean, while the chlorophyll a concentration was at a similar level in both regions (Darecki et al., 2003).

3.3 Heating rates and surface heat fluxes

The vertical and temporal evolutions of water-constituent-induced heating rates at each of the locations are shown in Fig. 11. Maximum heating rates occur in late spring and mid-summer and are between 0.8 and 0.9 K m⁻¹ d⁻¹ at Oder Bank and between 0.4 and 0.8 K m⁻¹ d⁻¹ at the other locations. Vertical profiles of two heating-rate maxima in May and July indicate that approximately 70 % of the water-constituent-induced heating is contained within the top 5 m, and decreases exponentially to zero by 10 to 15 m depth. We compared the Bio-Optic heating-rate estimates at Bornholm Basin with a comparable full radiative transfer calculation by MOMO for the two heating-rate maxima events in May and June (Fig. 11, bottom right). Bornholm Basin is chosen as the evaluation site for the heating-rate calculations because the seasonal cycle of the heat balance there can be approximated as a one-dimensional balance between the penetration of solar radiation and vertical mixing (Gnanadesikan et al., 2019), and advective and diffusive terms will be relatively small. The main difference between the two calculations, Bio-Optic and MOMO, is that the MOMO takes into account the full directionality of the light field while Bio-Optic does not. There are differences in the seasonal heating-rate results between the two approaches, but they are not so large. At the surface, the Bio-Optic estimates are 0.3 K m⁻¹ d⁻¹ smaller in spring and 0.25 K m⁻¹ d⁻¹ smaller in summer compared to the MOMO estimates. In the MOMO calculations, most of the water-constituent-induced heating (ca. 80 %) is contained within the top 2 m, and this decreases exponentially more rapidly than Bio-Optic to zero by 5 m depth. We find that, by accounting for the full directionality of the light field, as shown by the case investigated by MOMO, the impact water constituents have on the heating rates is contained within the top 2 to 3 m, consistent with the findings of Soppa et al. (2019). However, MOMO may be overestimating the actual magnitude of water-constituent-induced surface heating rates as none of the other physics (i.e. advection, diffusion) and environmental forcing represented in the Bio-Optic experiments are taken into account in MOMO. It could also be that the algorithm used to calculate K_d in Bio-Optic (Lee et al., 2005) is not optimal for the conditions in the Baltic Sea (we elaborate upon this point further in the discussion).

Figure 11 Surface heating rates (left panel) and vertical profiles of two heating-rate maxima in May and July 2018 (right panel) for the Oder Bank, Darß Sill, Arkona Sea and Bornholm Basin sites.

Figure 12 shows the temperature and chlorophyll a anomalies (biofeed minus nobiofeed experiments) for selected days during the productive period at Bornholm Basin. Accounting for the feedback of OSC-induced heating in the hydrodynamic solution has the effect of increasing the surface layer (ca. top 10 m) water temperature by between 0.1 and 0.2 ^∘C in spring and late summer and by as much as 0.5 ^∘C in mid-summer. Below the thermocline, the water temperature is cooler by 0.1 to 0.2 ^∘C. Differences in the thermal structure when the feedback is accounted for impact the development, transport and fate of phytoplankton biomass. This consequence is seen in differences in the chlorophyll a structure at different times during the productive period. The increase in light in spring supports phytoplankton growth and increases the surface temperature (due to both water and phytoplankton absorption) in the surface layer. Thus, the availability of light below the algae layer is strongly reduced, and phytoplankton are restricted within the shallow mixed layer with more availability of light, which will in turn increase surface heating. The net effect is more biomass production in the surface layer at the beginning of the spring bloom in biofeed compared to in nobiofeed. As nutrients become depleted in the surface layer and the supply of nutrients from deeper waters is inhibited by the stronger thermocline in mid-summer, the net effect is less biomass production in the surface layer in mid-summer in biofeed compared to in nobiofeed. As the water column becomes less stable in late August and as nutrients are mixed back into the surface, biomass production is larger again in biofeed compared to n nobiofeed.

Figure 12 Hovmöller plots of temperature and chlorophyll a anomalies (biofeed minus nobiofeed experiments) in 2018 at Bornholm Basin.

The impact this has on surface heat fluxes during the productive period at Bornholm Basin is shown in Fig. 13. The increase in OSC-induced surface temperature captured in spring and summer lead to an increase in heat loss to the atmosphere, with the average difference for the period April to September being on the order of 5.2 W m⁻². This is primarily a result of latent (2.6 W m⁻²) and sensible (1.7 W m⁻²) heat fluxes. Putting this into context with modelled estimates by Omstedt and Nohr (2004) of between 5 and 18 W m⁻² for the net annual heat losses in the Baltic Sea indicates that it may be important to consider OSC-induced heating rates in regional heat balance budgets.

Figure 13 Surface heat fluxes for both biofeed and nobiofeed experiments during the entire productive period of April to September (left panel) and a zoom-in on the period where the difference in surface heat fluxes between experiments is greatest (area shown in rectangular box in top left panel) at Bornholm Basin.

4 Discussion

Modelled seasonal and spatial changes in OSCs in the western Baltic Sea have a small but noticeable impact on radiative heating in surface waters, especially in spring and summer as a consequence of increased absorption of light by phytoplankton and CDOM. Our modelled estimates for 2018 show phytoplankton dominating absorption in spring and summer as a result of a succession of phytoplankton blooms and CDOM dominating absorption in summer and autumn. Simis et al. (2017) found that phytoplankton pigment visibly influences K_d(675) in spring and summer, while absorption by CDOM at 412 nm can account for 38 %–70 % of the total OSC absorption in the area influenced by the Oder River in autumn. First-order variability in CDOM absorption in the Baltic Sea is driven by terrestrial sources. Second-order variability is driven by autochthonous DOM production during phytoplankton blooms and photodegradation. The spatial and temporal variabilities in our modelled OSC absorption at the different locations, especially in spring, summer and autumn, are in good agreement with seasonal observations for different water types in the southern Baltic Sea reported by Meler et al. (2016a) (Fig. 6c, d). This is also bolstered by good agreement between the model and OLCI data match-ups with the background values at Darß Sill and Bornholm Basin, which gives us confidence in the model performance. This is encouraging for future modelling studies of this nature, as more consistent, long-term time series of the optical properties of the Baltic Sea are realised, e.g. using automated measurement systems such as Bio-Argo floats equipped with a simple spectral radiometer. Such a strategy has been applied with significant success in the Mediterranean Sea (Terzić et al., 2019, 2021a, b). We also find it encouraging that the (simplified) Bio-Optic and (full) MOMO radiative transfer heating-rate estimates were somewhat comparable and informative. The directionality of the light field appears to be important for understanding the depth of the influence of water-constituent-induced heating rates, while accounting for the spatial and temporal variability in the physics of the environment is important in determining the magnitude of the heating rates. However, we think further work is needed to optimise the Bio-Optic diffuse attenuation coefficient (K_d) algorithm for the Baltic Sea.

K_d, which describes the transfer of light energy through the water column, also reflects the seasonal variability of water types, i.e. winter (well-mixed) versus spring, summer and autumn (seasonally stratified), and the influence of constituents in different water types during stratified seasons (i.e. spatial variability). Our results show a gradient in K_d and in heating rates which decreases as you move offshore. In late spring, at the Oder Bank, the water constituent contribution to surface heating can be as much as 0.9 K m⁻¹ d⁻¹, while at the Darß Sill, Arkona Sea and Bornholm Basin sites, the water constituent contribution to surface heating in spring and summer is less, between 0.4 and 0.8 K m⁻¹ d⁻¹. Reports on the spectral properties and the temporal and spatial variability of the diffuse attenuation coefficient in the Baltic Sea based on field observations are limited and date back to the early 2000s (Kratzer et al., 2003; Lund-Hansen, 2004; Darecki and Stramski, 2004; Kowalczuk et al., 2005a; Lee et al., 2005). Darecki and Stramski (2004) have assessed that locally optimised satellite remote sensing algorithms for estimating K_d(490) based on MODIS data yield the least uncertainty compared to other variables, e.g. chlorophyll a. However, information on the full K_d spectrum is needed to assess the individual impact of the most significant optical seawater constituents on surface heating rates. Until recently, the only solution was empirical or semi-analytical modelling based on either remote sensing data (Lee et al., 2005; Löptien and Meier, 2011; Alikas et al., 2015) or in situ measurements of apparent or inherent optical measurements (Gonçalves-Araujo and Markager, 2020). The most accurate estimation of K_d could be achieved by using the semi-analytical model; however, uncertainty in those estimates heavily depends on the local parameterisation of the specific inherent optical properties which, in the Baltic Sea regions, have contrasting and highly variable seasonal cycles (Simis et al., 2017). Kratzer and Moore (2018) concluded that the correct choice of the volume-scattering phase function in the Baltic Sea determines the accuracy of the prediction of inherent and apparent optical properties in the Baltic Sea region. CDOM and suspended particles are the most significant optical constituents controlling water transparency. CDOM absorption is regulated mostly by riverine discharge, especially in coastal waters; however, under certain conditions, CDOM absorption in the Baltic Sea is statistically correlated with phytoplankton biomass (Kowalczuk et al., 2006; Meler et al., 2016b). Particulate absorption and scattering is significantly correlated with phytoplankton biomass, which has a well-defined seasonal and spatial pattern in the Baltic Sea (Meler et al., 2016a, 2017). By including a spectrally resolved underwater light field in our model and by diagnosing inherent and apparent optical properties, we are able to resolve the full K_d spectrum and better understand the role different OSCs play in determining the temporal and spatial variability in K_d and the impact on heating rates. Further optimisation of the Bio-Optic K_d algorithm for the Baltic Sea is currently in progress.

Climate change scenarios for central Europe predict significant change in the precipitation regime, which will be manifested in a shift in the seasonal distribution of precipitation: increased rainfall and a decline in snowfall in winter and persistent droughts in summer with episodic intensive thunderstorms (IPCC, 2019). Changes in the precipitation regime coupled with an increase in mean temperatures will significantly impact the outflow of freshwater from the Baltic Sea catchment into the marine basin itself (Meier et al., 2022). We could anticipate that the flux of terrestrial CDOM would be affected most because currently observed climatic changes in the southern part of the Baltic Sea catchment have caused mild winters with a reduced number of frost days and an almost-total reduction in snow fall. As a result, CDOM that was previously immobilised in the frosted ground, streams and rivers is now being transported to the sea in late winter and spring. In the summer, a deepening minima of flows in rivers reduces CDOM input to the Baltic Sea. Recent results by Zabłocka (2017) indicate that the monthly averaged Vistula River flow maximum during the period 1993 to 1998 occurred in April, while from 2008 to 2010, this maximum shifted to March. As the Baltic Sea is warming at a rate of up to 4 times the global mean warming rate (Belkin, 2009), we can expect this trend in earlier river flow maxima to continue and a higher contribution of CDOM to the absorption budget in winter and spring as the chlorophyll a concentration (phototrophic protist biomass proxy) maximum still occurs in April (Stoń-Egiert and Ostrowska, 2022).

Changes in the hydrological regime and a reduction in mineral nutrient input (Łysiak-Pastuszak et al., 2004) have noticeably impacted both phototrophic protist biomass and functional structure. Stoń-Egiert and Ostrowska (2022) have reported a statistically significant decreasing trend of 2.11 % yr⁻¹ of the total chlorophyll a concentrations over the last 2 decades (1999 to 2018), with decreasing pigment markers for such protist groups as diatoms, dinoflagellates, cryptophytes and green algae and an increase of cyanobacteria. As a consequence, primary production in the southern Baltic Sea also declined in the period from 1993 to 2018 compared to its maximum in the late 1980s (Zdun et al., 2021). Kahru et al. (2016) have also reported on changes in the seasonality in the Baltic Sea environment: the cumulative sum of 30 000 W m⁻² d⁻¹ of surface incoming shortwave irradiance (SIS) was reached 23 d earlier in 2014 compared to 3 decades earlier; the period of the year when the sea surface temperature is at least 17 ^∘C has almost doubled (from 29 d in 1982 to 56 d in 2014); the period when K_d(490) is over 0.4 m⁻¹ increased from about 60 d in 1998 to 240 d in 2013 (quadrupled); the period when satellite-estimated chlorophyll is at least 3 mg m⁻³ has doubled from 110 d in 1998 to 220 d in 2013; and the timing of both the phytoplankton spring and summer blooms has advanced, with the annual chlorophyll maximum that, in the 1980s, corresponded to the spring diatom bloom in May has now shifted to the summer cyanobacteria bloom in July. It is interesting to note that we found two OSC-induced heating-rate maxima in May and July in our model results, which coincide with two observed marine heatwave events. At the Darß Sill and Arkona Sea sites, these heating-rate maxima were larger in May by 0.18 and 0.35 K m⁻¹ d⁻¹, respectively, compared to in July, while at Oder Bank, the heating rate maxima were larger in July by 0.1 K m⁻¹ d⁻¹.

5 Conclusions

Heating rates due to the absorption of short wave radiation (UV–VIS) in the western Baltic Sea are controlled by the combined effects of the seasonal solar cycle and the concentration and distribution of OSCs. The intensity of radiative energy reaching the sea surface is locally modified by radiative transfer through the atmosphere, which is mostly controlled by cloudiness, of which the long-term climatology minimum is observed in May (Dera and Woźniak, 2010). Further modulation of heating rates in the western Baltic Sea in UV and VIS spectral domains is dependent on water transparency, which is a complex function of the magnitude and seasonal cycles of inherent optical properties and the directionality of the light field. Our study found that, in 2018, the combined effect of CDOM and particulate absorption on surface heating rates in the western Baltic Sea could reach up to 0.4 to 0.8 K d⁻¹ during the productive period of April to September and is relevant from the surface down to 2–5 m depth. Moreover, this modelled OSC-induced surface warming results in a mean loss of heat (ca. 5 W m⁻²) from the sea to the atmosphere, primarily in the form of latent and sensible heat fluxes, which may be significant for regional heat balance budgets. Two-way coupling with the atmosphere is not included in our experiment, but we expect that this would modulate (decrease) the magnitude of the net loss of heat to the atmosphere.

Anticipated and recently observed changes in phytoplankton functional types and their seasonal pattern and CDOM terrestrial input patterns due to global warming will further modulate the spatial and temporal pattern of heating rates in the Baltic Sea. Observed changes in the quantity and quality of CDOM, the composition and concentration of phytoplankton functional types and photosynthetic pigments, and thus changes to the optical properties of the Baltic Sea need to be communicated to coupled hydrodynamic–biogeochemical models such that the consequences of radiative feedbacks can be better understood and better predictions of the future Baltic Sea environment can be made. Further improvements to coupled hydrodynamic and ecological models are heavily dependent on the correct parameterisation of the downwelling irradiance diffuse attenuation coefficient, K_d, which requires a proper understanding of the seasonal and spatial variability of the optical properties in different water types. This work highlights the importance of K_d as a bio-optical driver: K_d provides a pathway to estimate heating rates and connects biological activity with energy fluxes.

Appendix A: Western Baltic Sea model setup

Table A1 Model configurations. ERGOM: Ecological Regional Ocean Model; DT: baroclinic time step; NDTFAST: number of barotropic time steps between each baroclinic time step.