Optical model for the Baltic Sea with an explicit CDOM state variable: a case study with Model ERGOM (version 1.2)

Colored dissolved organic matter (CDOM) in marine environments impacts primary production due to its absorption effect on the photosynthetically active radiation. In coastal seas, CDOM originates from terrestrial sources predominantly and causes spatial and temporal changing patterns of light absorption which should be considered in marine biogeochemical models. We propose a model approach in which Earth Observation (EO) products are used to define boundary conditions of CDOM concentrations in an ecosystem model of the Baltic Sea. CDOM concentrations in riverine water derived from EO products serve as forcing for the ecosystem model. For this reason, we introduced an explicit CDOM state variable in the model. We show that the light absorption by CDOM in the model can be improved considerably in comparison to approaches where CDOM is estimated from salinity. The model performance increases especially with respect to spatial CDOM patterns due to the consideration of single river properties. A prerequisite is high-quality CDOM data with sufficiently high spatial resolution which can be provided by the new generation of ESA satellite sensor systems (Sentinel 2 MSI and Sentinel 3 OLCI). Such data are essential, especially when local differences in riverine CDOM concentrations exist.


Introduction
Colored dissolved organic matter (CDOM) is a major light absorption constituent in the marine environment and especially in coastal seas. The spectral absorption characteristic of CDOM follows an exponential function with highest absorption towards shorter wavelengths (e.g., Nelson and Siegel, 2002). By modifying the underwater light climate, CDOM has an impact on primary productivity; e.g., in clear water sufficient light intensity to enable phytoplankton growth is available down to greater depths than in turbid waters (Dutkiewicz et al., 2015). Water temperature is affected by CDOM absorption as well. In turbid water, the shortwave light absorption is located in the upper water column increasing the temperature, while in clear water a thicker layer is warmed but to a lesser degree (Jolliff and Smith, 2014). This process impacts the sea surface temperature (SST) in particular. Model studies show an SST increase up to 2 K in coastal regions when colored organic materials are considered (Gnanadesikan et al., 2019).
Jerlov (1976) developed a classification for different water masses based on specific optical properties. This classification is widely used in numerical ocean models (e.g., Griffies, 2004). For global models, this parametrization works reasonably well. However, coastal ecosystems with substantial terrestrial runoff require more detailed parametrization of light penetration. In particular, variable light attenuation in river Table 1. Constants of the optical model.

Constant
Value Unit k w 0.027 m −1 k c 0.029 m 2 (mg Chl) −1 k det 0.0039 m 2 (mg N) −1 k don 0.0009 m 2 (mg N) −1 k cdom 0.221 1 dr 0 8.75 × 10 −5 d −1 r 0.5 1 plumes and their environs affect the hydrodynamic and ecological response (Cahill et al., 2008). Marine CDOM comprises humid substances (yellow substances) of terrestrial origin, and autochthonously produced CDOM. Its degradation is governed by photochemical bleaching and bacterial activity. In freshwater-dominated systems, like the Baltic Sea, terrestrial CDOM dominates (e.g., Stedmon et al., 2010, and references herein). In such systems, salinity and light absorption due to CDOM show a robust relationship (Neumann et al., 2015). The relationship is hyperbolic indicating the degradation processes.
For coupled physical-biogeochemical models of freshwater-influenced coastal seas, a spatially resolved CDOM concentration is important for realistic light climate estimates. Based on the nearly conservative character of CDOM, statistical models have been developed which estimate absorption from salinity and, for example, chlorophyll (Kowalczuk et al., 2006;Neumann et al., 2015). These models usually deliver reasonable results. However, two distinct disadvantages are prominent: (i) uncertainties in model salinity propagate into the biogeochemical model and (ii) variability in CDOM riverine load (Skoog et al., 2011;Asmala et al., 2013) cannot be resolved. These disadvantages can be eliminated by introducing an independent CDOM state variable into the biogeochemical model (Dutkiewicz et al., 2015). A necessary prerequisite is boundary data for riverine CDOM loads of sufficiently good quality, which is not commonly available for most rivers (Pefanis et al., 2020).
In this study, we present the implementation of a CDOM state variable in the biogeochemical model ERGOM (Ecological ReGional Ocean Model, Leibniz Institute for Baltic Sea Research, 2015) and the generation of CDOM boundary data with the aid of satellite imagery, and we discuss the effect of the proposed model extension on the Baltic Sea ecosystem. . The green line is the coastline of the 3 nm model. We refer to the labeled river later in the text. The map was created using the software package GrADS 2.1.1.b0 (http://cola.gmu.edu/grads/, last access: 6 August 2021), using published bathymetry data (Seifert et al., 2008).

Model development and description
We start presenting the optical model including an explicit CDOM state variable and its implementations as part of a biogeochemical model, and we briefly introduce the circulation and biogeochemical model used for this study.

Fundamentals of the optical model
Starting point of the development is the optical model as in the study by Neumann et al. (2015). The photosynthetical active radiation (PAR) follows an exponential decay with depth z:  The map was created using the software package GrADS 2.1.1.b0 (http://cola.gmu.edu/grads/, last access: 6 August 2021), using published bathymetry data (Seifert et al., 2008).
K PAR is the underwater bulk light attenuation and is described by five components: where k c , k det , and k don are material-specific constants and Chl, DET, and DON are concentrations of chlorophyll, detritus, and dissolved organic nitrogen, respectively. These concentrations are state variables of the ecosystem model ERGOM or, in the case of chlorophyll, can be estimated from model phytoplankton (Sect. 2.1.2). k w is the attenuation coefficient for pure water. In Neumann et al. (2015), K CDOM (S) is a statistical relationship between in situ salinity and CDOM absorption for the Baltic Sea derived from observations (Eq. A1). For the new approach, we use the additional state variable CDOM: The PAR attenuation now reads as follows: Terrestrial CDOM behaves nearly conservatively in the ocean. An indication is the linear salinity-CDOM relationship in the northern Baltic (Harvey et al., 2015;Neumann et al., 2020). This is due to the fact of high freshwater supply with high CDOM concentrations. However, this relation does not apply to the central Baltic. In this region with longer residence time, the effects of CDOM degradation processes become more pronounced and observable (Skoog et al., 2011).
Two processes control CDOM degradation, photobleaching and biological degradation. Moran et al. (2000) studied the degradation of terrestrial CDOM in the coastal ocean and found that photobleaching accounts for 80 % of the degradation. Furthermore, they showed that CDOM decay closely follows simple first-order kinetics. In accordance with these findings, we implement CDOM as with the degradation rate dr = dr 0 · PAR(z).
PAR(z) is the ambient PAR at depth z, dr 0 is a constant, and PAR(z) can be estimated as

5052
T. Neumann et al.: CDOM in marine coastal models I 0 is the solar radiation at sea surface depending on sun zenith angle, which is a function of latitude and time. Factor r is the fraction of I 0 available as PAR (spectral range 400 to 700 nm). The integral considers the depth dependence of model concentrations of, for example, chlorophyll in Eq. (4).

Implementation of the optical model
A comprehensive overview of CDOM in the ocean is given in Nelson and Siegel (2002). CDOM concentration is usually given by a proxy, the light absorption for a specific wavelength, e.g., a CDOM (440) for 440 nm. The spectral distribution can be parameterized by an exponential decline with wavelength.
s is the exponential slope parameter and varies between 0.015-0.025 nm −1 . For a given slope s and an absorption a CDOM (λ 0 ), any a CDOM (λ) can be estimated for wavelengths longer than 320 nm. We use a reference wavelength of 440 nm. The biogeochemical framework for implementing the CDOM state variable is the model ERGOM. In this model, state variables are given as concentrations of an element, e.g., mol carbon per m 3 , because these models primarily describe cycles of elements (carbon, phosphorus, nitrogen etc.). In order to model CDOM, a relationship between CDOM absorption and the concentration is required. Following the Lambert-Beer law, a linear relation exists. Neumann et al. (2020) derived a relationship based on optical measurements and measurements with a calibrated CDOM sensor, which we used to convert CDOM absorption into concentration and vice versa. We have to note that the accuracy of the conversion does not impact the performance of the optical model because both the satellite-derived CDOM in freshwater and CDOM in the optical model is given as a CDOM (440). CDOM loads are estimated from concentration times runoff which is available from Gustafsson et al. (2012).
Most model constants used in the model are provided by Neumann et al. (2015). We use 50 % of total solar radiation for PAR (e.g., Stigebrandt and Wulff, 1987) since the invisible, longwave part is absorbed at the water surface (factor r, Eq. 6). dr 0 (Eq. 6) has been estimated with a series of calibration simulations. The aim of the calibration was to find an optimal match of observed and simulated absorption values a CDOM (440). The constants used are listed in Table 1. Model state variables in Eq. (4) have to be converted into appropriate units before entering the optical model. We use a volumebased concentration. CDOM should be given as absorption because of uncertainties in the absorption-concentration relationship.
The technical implementation is done by an automatic code generation. Fundamentals are a set of text files describing the biogeochemistry independently of computer language and the host system. Code templates describe physi-cal and numerical aspects and are specific for a certain host, e.g., a circulation model. All necessary ingredients (the code generation tool, text files, and templates for several systems) can be downloaded from http://www.ergom.net (last access: 22 September 2020). The same technique is used for example in Radtke et al. (2019).

Circulation and biogeochemical model
For model testing, we have used a similar model system to that in Neumann et al. (2015). The circulation model is MOM5.1 (Griffies, 2004) adapted for the Baltic Sea. The horizontal resolution is 3 nautical miles. Vertically, the model is resolved into 152 layers with a layer thickness of 0.5 m at the surface and gradually increasing with depth up to 2 m. The circulation model is coupled with a sea ice model (Winton, 2000) accounting for ice formation and drift.
Coupled with the circulation model is the biogeochemical model ERGOM. It describes a marine nitrogen and phosphorus cycle. Primary production, forced by PAR, is provided by three functional phytoplankton groups (large cells, small cells, and cyanobacteria). Chlorophyll concentration can be estimated from the phytoplankton groups, which is used in the optical model. Dead particles accumulate in the detritus state variable which is another compartment in the optical model. A bulk zooplankton grazes on phytoplankton and constitutes the uppermost trophic level in the model. The metabolism of phytoplankton and zooplankton produces DON which has only little impact on light absorption. Phytoplankton and detritus can sink down in the water column and accumulate in a sediment layer. In the water column and the sediment, detritus is mineralized into dissolved inorganic nitrogen and phosphorus. Mineralization is controlled by temperature and oxygen. Oxygen is produced by primary production and consumed due to all other processes, e.g., metabolism and mineralization. Coupled to the nitrogen and phosphorus cycle is a carbon cycle as described in Kuznetsov and Neumann (2013). A schematic of the model structure is provided in Appendix B. The estimated shortwave absorption (Sect. 2.1.2) feeds back into the physical part of the model and hence impacts the temperature distribution.
The new CDOM variable in the current model development state is not involved in the biogeochemical processes. This is justified by the fact that CDOM is relatively refractory and has a long residence time, and autochthonous CDOM produced by for example phytoplankton is a small fraction. In later developments, it will be included in the carbon cycle. For this purpose, it is essential to realize CDOM as a carbon-based concentration. If this will be not the case, the CDOM state variable could be implemented as an absorption and thus conversions between absorption and concentration could be prevented.
The model has been forced by meteorological data from the coastDat-2 dataset (Geyer and Rockel, 2013). We run the model from 1948-2019. A first run was used to spin up the new CDOM tracer. In a second run, CDOM was initialized with data from the first run. In addition to the 3-nautical-mile resolution, we use a 1-nautical-mile resolution for the period 2017-2019. The model has been successfully used in several applications (e.g., Neumann, 2010;Neumann et al., 2015).

CDOM boundary data from Earth observation products
Aim of the development is to improve the simulated light climate by a more realistic representation of CDOM concentration compared to available statistical models (Sect. 1). However, in situ data of CDOM loads are not available in sufficient spatial and temporal resolution; i.e. in the Baltic Sea, CDOM is not a parameter of the HELCOM (http://www. helcom.fi, last access: 22 September 2020) monitoring program. New instruments and technologies in Earth observation (EO), now in operation, are ideal tools to overcome these limitations.

Characteristics of satellite data
In order to estimate the load of CDOM coming from a river, it is necessary to derive it from observations within the river or as close to the discharge point as possible. The rivers in the Baltic Sea are usually small and thus high-resolution (HR) instruments such as Sentinel-2 Multi Spectral Instrument (S2-MSI) are required. The MSI has a spatial resolution of 10-60 m depending on the central wavelength of the band. Water quality products are usually generated in 60 m resolution in order to reduce noise. This is sufficient for estimating the CDOM absorption of most rivers. Two Sentinel-2 satellites are currently in orbit: S2A was launched on 23 June 2015 and S2B on 7 March 2017. Together they provide a global revisit time of 5 d next to the Equator. Due to the high latitude of the Baltic Sea, the revisit time amounts to 2-3 d in this area. Despite the frequent cloud cover, sufficient observations can be gathered to monitor river CDOM throughout the open water season (typically from March-April to October in the northern Baltic Sea).

Earth Observation processor for CDOM absorption estimation
Earth Observation (EO) processors are a set of algorithms designed to convert the radiance signal acquired by the satellites into values of geophysical parameters. The estimation is based on the scattering and absorption features of the material suspended or dissolved in water. In addition to CDOM, these materials include phytoplankton cells, represented by Chlorophyll a (Chl-a), and suspended particulate matter. One commonly used water quality processor is Case 2 Regional Coast Color (C2RCC) (Brockmann et al., 2016). It utilizes one artificial neural network (ANN) first to remove the effects of the atmosphere from the signal (atmospheric cor-rection) and then another to estimate inherent optical properties (IOPs) of water from the marine reflectance.
For estimation of CDOM absorption (a CDOM ), we utilized the C2RCC (version 1) output called a dg (combined absorption by detritus and yellow substances), which was calibrated to a CDOM (440) values (absorption coefficient by CDOM at 440 nm) using this equation: The equation is based on in situ sampling made with a flowthrough device (ac-9) (Lindfors et al., 2005;Koponen et al., 2007) during two coastal estuary measurement campaigns. These data are not yet published but a similar local calibration method has provided good results with other water quality parameters in earlier studies such as Attila et al. (2013). The data processing and extraction are done in a Calvalus massive parallel processing system (http://www. brockmann-consult.de/calvalus, last access: 22 September 2020). Data extraction areas are manually defined in the vicinity of the mouths of 69 rivers that represent ERGOM input locations (Fig. 1).
The areas are designed so that islands, mixed pixels, and shallow areas are excluded. All valid pixels (not masked as land or cloud by the pixel classification processor Idepix) within each area and image are collected and analyzed, and the 75th percentile value is chosen to represent the river a CDOM . The cases in which the number of valid pixels is less than 50 % of all available pixels from an area are removed from the analysis. Presumably, these represent cases with partial cloud cover and they are discarded to keep only estimates with highest quality and low uncertainty. The arithmetic means of the 75th percentile pixel values of all valid days within each calendar month during years 2017-2019 are then computed for each extraction area.
We are aiming to provide the ecosystem model ERGOM with an annual cycle of CDOM loads based on monthly data. Since optical EO methods cannot provide a CDOM estimates in darkness and throughout times with ice coverage, the values for the winter months have been interpolated. As a result, the dataset contains a a CDOM value for each month for each of the areas under investigation (69 extraction areas in total). Figure 2 shows four examples of the annual CDOM absorption cycle. The behavior of the data follows the annual cycle well: spring values are high due to the terrestrial matter brought into the coastal water by melting snow, summer values are low due to lower influx, and the fall values are higher due to increasing rainfall.

Results and discussion
In this section, we show the improved model CDOM representation and its impact on the simulation results. Owing to the changed shortwave distribution in the water column, an effect on the biogeochemistry in particular is ex-  Absorption estimates based on simulated CDOM are shown in black, based on a conversion from simulated salinity in green, and red diamonds are observations. The map was created using the software package GrADS 2.1.1.b0 (http://cola.gmu.edu/grads/, last access: 6 August 2021), using published bathymetry data (Seifert et al., 2008). pected. All CDOM data presented are converted into absorption at 440 nm (a CDOM (440)). Especially for observations at different wavelengths, we use Eq. (8) with a slope s of 0.018 nm −1 (Kratzer and Moore, 2018). For comparison, we show a CDOM (440) values estimated from the CDOM state variable and from model salinity. The models differ only in the estimation of PAR which becomes evident when we show the impact on biogeochemistry.

CDOM absorption
In Fig. 3a, we show the simulated CDOM absorption at the sea surface. The snapshot clearly illustrates the spatial patterns. Strong absorption is visible in the northern Baltic and the river mouths. The difference to the salinity-based estimate (Eq. 2) is depicted in panel (b). The strongest differences appear in the Gulf of Bothnia and the Gulf of Finland while in the central Baltic differences are small. Strong Figure 7. 2018 annual mean profiles at station 4 (Fig. 6d). Blue curve shows data from the model with an explicit CDOM state variable, and the red curve is from the model with CDOM absorption estimates based on its relation to salinity. For oxygen (d), we show the whole water column.
differences are also pronounced in river estuaries. Owing to the low salinity, the salt-CDOM relationship overestimates CDOM content in these areas. Furthermore, the new, EObased method considers individual CDOM concentrations of different rivers. Rivers of the northern catchment area carry higher CDOM loads compared to rivers of the southeastern catchment area due to a high fraction of peat land.
The range of a CDOM values is demonstrated in Fig. 4. Both datasets were compared against in situ data collected from monitoring stations in coastal waters of Finland and from the northern coast of Sweden. As shown in Fig. 5, the improvement becomes obvious. With the salinity method, the correlation is low (R 2 = 0.15) and there are some clear overestimates (difference between a CDOM (440) from salinity and the one-to-one line more than 2 m −1 ) while most data points are underestimated. With the EO CDOM method, the correlation improves significantly (R 2 = 0.61). There are no large overestimates and the data points move closer to the oneto-one line. Large in situ values (a CDOM (440) > 3 m −1 ) are still underestimated with the new EO method. This underestimation is most likely caused by the following two major inaccuracies in the present version of the model input: -The coast has many small rivers. Not all of them are yet included in this version.
-In river estuaries with low bottom depth or complex morphological structure, the shapes and formulations of the extraction areas do not sufficiently capture the incoming CDOM loading from the river. In order to avoid EO observations contaminated by bottom reflectance it was necessary to use pixels that are sufficiently far from the shoreline. Therefore, pixels of the extraction area may not represent river water as it has been already mixed with sea water. In some cases, this leads to lower concentrations, especially during the low runoff season. Figure 6 demonstrates the different CDOM absorption estimates. Shown are time series of surface CDOM absorption at six stations (Fig. 6d). The green curve is the salinity-based estimate and the black curve the estimate from model CDOM concentration. Red diamonds are in situ observations. At stations 1-4, absorption values from simulated CDOM are much closer to observations compared with salt-based estimates. The absorption difference at stations 5 and 6 is less pronounced, which is also evident from Fig. 3. The seasonal variability is stronger for absorption derived from the model CDOM variable (compared to salt-based absorption) and reflects the observed variability. The reason is the annual cycle of riverine CDOM concentration in addition to the runoff . The shaded area is the range between 10th and 90th percentile of the black line. Simulated data are from the 3-nautical-mile model version. The map was created using the software package GrADS 2.1.1.b0 (http://cola.gmu.edu/grads/, last access: 9 August 2021), using published bathymetry data (Seifert et al., 2008). cycle. In the central Baltic Sea, both methods overestimate CDOM absorption.

Impact on biogeochemistry
As an example of the changed light absorption impact on the biogeochemistry, we show annual mean profiles of selected variables in Fig. 7. We have chosen station 4 from Fig. 3 since at this station the CDOM absorption is considerably increased due to the new, EO-based, optical model and it is located in the center of the Bothnian Bay. Owing to the increased CDOM absorption, PAR is reduced in the EO model approach (Fig. 7a) as expected. Consequently, primary production (PP) is reduced (Fig. 7b). However, in the uppermost layer, PP is increased. As a result of reduced PP, phytoplankton concentration shows lower values (Fig. 7c). An integrated response is the increased bottom oxygen concentration (Fig. 7d). Less net PP results in less accumulation of organic matter in the deep water of the basin and subsequently reduced oxygen consumption. The impact on water temperature is small (order of 0.01 K, not shown). The effect is a temperature increase in the surface layer and a lower temperature below.
We demonstrate changes in the biogeochemistry with a climatology of surface nutrient concentrations at three stations in Fig. 8. For this analysis, we use data from the 3nautical-mile model version because of the longer simulation period. Shown are data from the EO model (black) and the previous (salinity-based) model (green) together with observations (red). In the Gulf of Finland at station KAS-11, the spring-bloom-related nutrient depletion is delayed by 2 weeks (Fig. 8a and b). Sufficient PAR intensity, initiating a bloom, is available later in the season. The winter nutrient concentrations are elevated compared to the salinity model version. At the Bothnian Bay station ( Fig. 8c and d), the spring bloom delay is less pronounced. In this area, the longer sea ice coverage dominates the PAR in spring. At station BY15 in the Baltic Proper ( Fig. 8e and f), the difference between both model versions is small due to similar CDOM absorption (Fig. 6).

Summary and conclusions
In this study, we propose an approach for considering light absorption due to terrestrial CDOM in a marine ecosystem model for the Baltic Sea. An explicit consideration is necessary if large amounts of terrestrial CDOM enter the marine system and strong coastal-sea gradients develop. In such cases, a uniform light absorption due to CDOM cannot account for the in situ light climate in a sufficient way. An often-applied approach uses CDOM-salinity relationships for CDOM absorption estimates (Kowalczuk et al., 2006(Kowalczuk et al., , 2010Neumann et al., 2015) but with distinct disadvantages (see Sect. 1).
Our approach uses an explicit CDOM state variable as part of the biogeochemical model. In order to improve the simulated absorption compared to the salt approximation, a high quality dataset of riverine CDOM loads is necessary. This has been accomplished by using Earth observation data from Sentinel-2 MSI. The high spatial resolution (10-60 m) allows the river mouths to be observed directly. A difficulty in the regions of higher latitude like the Baltic Sea area is the insolation, the occurrence of sea ice, and the frequent cloud cover in winter. Continuous observations are not possible during this time. We have used a linear interpolation to bridge the winter data gap. This could be validated by ground truth measurements in winter possibly guiding for a non-linear interpolation.
The results (Sect. 3) show that the proposed approach clearly improves the ability of the model to estimate CDOM and thus light absorption, especially in the northern parts of the Baltic Sea where the impacts of terrestrial CDOM are large. This underlines the performance of the combined approach to increase the predictive capability of ecosystem models. The method can be further improved by adding more rivers to the model and improving the quality of CDOM data from Sentinel 2 MSI.
For the model CDOM, we have applied a light-sensitive degradation. Although this is the dominating degradation process for terrestrial CDOM (Moran et al., 2000), bacterial breakdown contributes to the degradation as well. Technically, such a process can easily be implemented. However, to our knowledge comprehensive process studies in the Baltic Sea are not done yet. Therefore, we have decided that bacterial breakdown is subject to later developments.
We consider only terrestrial CDOM in our model. In regions with high runoff, like the Baltic Sea, terrestrial CDOM is the dominating fraction (Harvey et al., 2015;Stedmon and Markager, 2003). However, a further step toward a more sophisticated model could be the inclusion of autochthonous CDOM as for example in Dutkiewicz et al. (2015). Photobleaching accounts for slow decomposition of CDOM. Although the CDOM decomposition is slow compared to decomposition of in situ detritus for example, in water bodies with longer residence time it becomes important. For example, in the salt-CDOM absorption relationship, decomposition is considered by a hyperbolic function (Neumann et al., 2015, Eq. 10): We demonstrate the effect of the implemented photobleaching on the model CDOM concentration by comparing a simulation without photobleaching with the control run. The 3-nautical-mile setup was used for this study. The photobleaching was switched off in 1980 and then the simulation was continued until 2019. Figure A1 shows the differences developing due to the lack of a degradation process. Absorption values are far away from observations and even after 40 simulation years, a new steady state is not achieved (Fig. A1a). Spatial patterns after 39 simulation years show that the largest differences occur in the central basins. Differences are smaller in freshwater-dominated regions (Fig. A1b) like river estuaries. Small differences are also in the vicinity of the open boundary toward the North Sea. Figure A1. The effect of photobleaching on CDOM absorption. The black line in (a) is from the model without photobleaching at station BY15 (Fig. 8). In (b) the difference in 2018 of the surface CDOM absorption is shown (without photobleaching minus control run). Simulated CDOM concentrations have been converted into absorption. The map was created using the software package GrADS 2.1.1.b0 (http://cola. gmu.edu/grads/, last access: 6 August 2021), using published bathymetry data (Seifert et al., 2008).

5060
T. Neumann et al.: CDOM in marine coastal models