Tracking freshwater browning and coastal water darkening from boreal forests to the Arctic Ocean

The forest cover of Northern Europe has been steadily expanding during the last 120 years. More terrestrial vegetation and carbon fixation leads to more export to surface waters. This may cause freshwater browning, as more degraded plant‐litter ends up as chromophoric (colored) dissolved organic matter. Although most freshwater ultimately drains to coastal waters, the link between freshwater browning and coastal water darkening is poorly understood. Here, we explore this relationship through a combination of centennial records of forest cover and coastal water clarity, contemporary optical measurements in lakes and coastal waters, as well as an ocean drift model. We suggest a link between forest cover in Northern Europe and coastal water clarity in the Baltic, Kattegat, and Skagerrak Sea and show how brown‐colored freshwater from Northern European catchments can dictate coastal water clarity across thousands of kilometers, from the Baltic lakes to the Barents Sea.

growing season (Barichivich et al. 2013), more vegetation (Myneni et al. 1997) and forest growth (Thomas et al. 2010;Høgberg et al. 2017). For instance, the growing stock of Fennoscandian forests (Norway, Sweden, and Finland) has more than doubled during the last century (Claesson et al. 2015;Breidenbach et al. 2020;Anon 2022), partly also due to large forest planting and management programs started in the 1950s (Östlund et al. 1997).
Higher terrestrial plant biomass and primary production adds dissolved organic carbon (DOC) to water, either directly by leachate from litterfall or roots with associated mycorrhiza, or indirectly via bogs and soils. This is the dominant pool of carbon in most freshwaters in the boreal zone (Cole et al. 2007), where it is visible as a yellow or brown color, often referred to as chromophoric (colored) dissolved organic matter (CDOM). The greening of catchments therefore has a profound browning effect on recipient freshwaters (Meyer-Jacob et al. 2015;Finstad et al. 2016;Kritzberg 2017;Skerlep et al. 2020), which in turn drain to marine coastal waters. Changed hydrology (de Wit et al. 2016), peatland ditching (Nieminen et al. 2021), and reduced acid deposition (Evans et al. 2006;Monteith et al. 2007) will also contribute to this browning.
In a study of 77 Swedish and Norwegian freshwater lakes, CDOM was found to be the main driver of light attenuation, even in highly eutrophic lakes (Thrane et al. 2014). In the Baltic Sea, which supplies around 70% of the freshwater to the southern part (Skagerrak) of the Norwegian Coastal Water (NCW; Aure et al. 1998), CDOM from freshwater sources is an important driver of light attenuation (Kowalczuk et al. 2006). Here, changes in chlorophyll a (Chl a) concentration typically explains maximum 10-17% of the variation in water clarity (Fleming-Lehtinen and Laamanen 2012;Harvey et al. 2019), while CDOM can account for up to 70-90% of that variation (Kowalczuk et al. 2006;Harvey et al. 2019). The importance of Chl a concentration for light attenuation is highest during periods of low freshwater discharge and high biological activity (Kowalczuk et al. 2006;Fleming-Lehtinen and Laamanen 2012).
For the North Sea, Skagerrak Sea, and Baltic Sea, several studies have found that light attenuation has been increasing throughout the 20 th century (Sandén and Håkansson 1996;Aarup 2002;Dupont and Aksnes 2013;Capuzzo et al. 2015;Opdal et al. 2019;Kahru et al. 2022), inferred from a gradual decreasing Secchi disk depth-a phenomenon often referred to as coastal water darkening. Interestingly, the term darkening originally referred to the freshening of coastal waters, and associated increase in light attenuation (Aksnes et al. 2009), but is more commonly used as any nonbiotic driver of reduced light penetration in coastal waters, including both dissolved and suspended matter (Frigstad et al. 2013;McGovern et al. 2019;Blain et al. 2021;Konik et al. 2021).
There is little evidence that increasing Chl a concentration has been an important driver of coastal water darkening in these areas (Fleming-Lehtinen and Laamanen 2012;Opdal et al. 2019;Kahru et al. 2022), but see Sandén and Håkansson (1996) for a general discussion. Rather, increasing CDOM concentration from freshwater sources has been suggested as a likely cause (Fleming-Lehtinen and Laamanen 2012;Dupont and Aksnes 2013;Opdal et al. 2019;Kahru et al. 2022), and possibly also wind driven resuspension of sediments-at least for the most recent decades and in the shallower parts of the North Sea (Capuzzo et al. 2015;Wilson and Heath 2019).
In this study, we investigate the terrestrial and freshwater drivers of 120 years of coastal water darkening in the Skagerrak Sea and explore its extended geographical impact on light attenuation along 3000 km of Norwegian coastline.

Materials and methods
See Table 1 for a list of abbreviations, units and explanations, and Opdal (2022) for all raw data used in this study.

Estimation of nonphytoplankton light attenuation in the Skagerrak Sea
Nonphytoplankton light attenuation (K NON ) from substances such as suspended and dissolved organic matter was estimated based on a combination of in situ measurements of Secchi disk depth (S, m, N = 256,536), salinity (SAL, psu, N = 714,400), and Chl a (CHL, mg m À3 , N = 943,548, includes 6944 estimates from cell counts and phytoplankton color index; see table S1) taken in the Baltic, Kattegat, and Skagerrak seas between 1903 and 2021 (Supporting Information Table S1; Fig. S2). As the focus of this study is largely on the freshwater endmember of coastal waters, we represent salinity (SAL) as the freshwater fraction, FWf = 1 À SAL/35.2, where 35.2 is the salinity of the assumed ocean endmember, North Atlantic Water (Aksnes 2015).
Due to variable data coverage, we applied two separate approaches (A1970 and A1903) for estimating K NON . Both approaches rely on the mechanistically derived relationship between Secchi disk depth (S) and the attenuation coefficient of downwelling irradiance, K S = 1.48/S (Lee et al. 2018). For a general discussion, see Turner et al. (2022).
The first approach (A1970) is based on a single generalized additive model that incorporates 48,143 simultaneous (same day and location) measurements of Secchi disk depth (S), salinity (SAL) and Chl a (CHL), such that where s 1 and s 2 are smoothing terms and ε is the error term. Low spatiotemporal coverage before 1970 limited this approach to the years 1970-2021. However, it allowed us to disentangle the effects of salinity (SAL) and Chl a (CHL) on light attenuation (K S ) as well as predicting the annual (1970-2021) change in K NON for a fixed value of SAL and CHL. See the Supporting Information for model details.
The second approach (A1903) covered the period 1903-2021 and follows the methodology described in Opdal et al. (2019). For photosynthetically active radiation (400-700 nm) we assumed that light attenuation from sources other than phytoplankton and water itself (K NON ) can be derived from a quasi-inherent composite light attenuation, such that where K PHY is the light attenuation from phytoplankton (K PHY = 0.121 CHL 0.428 , m À1 ; Morel 1988) and K W is the contribution of water itself (K W = 0.0089 m À1 , Morel and Maritorena 2001). Based on 84,719 pairs of Secchi disk and salinity measurements and 943,548 Chl a measurements and estimates (Supporting Information Fig. S2G-I), we constructed two separate generalized additive models (Supporting Information Eqs. S1, S2) that allowed us to obtain annual estimates of Secchi disk depth (S, Supporting Information Eq. S1) and Chl a (CHL, Supporting Information Eq. S2), and their respective light attenuation coefficients K S and K PHY , from 1903 to 2021 in the southern part of the NCW in the Skagerrak Sea (LAT = 58 N, LON = 8.5 E, star in Fig. 1) in January (MONTH = 1), when the effect of Chl a is at its lowest.
Connecting coastal water darkening to freshwater browning Coastal water darkening in the southern part of the NCW, as expressed by increased K NON , was linked to freshwater browning through the nonchlorophyll light attenuation of the freshwater endmember of the coastal water (K FW NON ).
In the first approach (A1970, YEAR = 1970-2021), a general estimate of K FW NON was obtained from Eq. 1 by setting Chl a and salinity both to zero, such that K FW NON (YEAR) = K S (YEAR, SAL = 0, CHL = 0). In the second approach (A1903, YEAR = 1903-2021), K FW NON was derived for the southern NCW in January from Eq. 2, where K S is obtained from Secchi disk depth (S) at zero salinity (Supporting Information Eq. S1), such that K FW NON (YEAR) = K S (YEAR, LAT = 58, LON = 8.5, MONTH = 1, SAL = 0) -K PHY (YEAR, LAT = 58, LON = 8.5, MONTH = 1) À K W . An increase in K FW NON over time indicates that the freshwater entering the coastal water has undergone browning. The latter estimate of K FW NON was then compared to the change in forest cover (see Supporting Information) in the relevant drainage area (Fig. 1).

Reduced K FW NON from the Skagerrak to the Barents Sea
The freshwater endmember of the southern NCW originates from rivers draining into the Kattegat and Baltic Sea ($ 15,000 m 3 s À1 ), North Sea ($ 4500 m 3 s À1 ), and Skagerrak ($ 2500 m 3 s À1 ), making up a total of ca 22,000 m 3 freshwater per second (Aure et al. 1998). From Skagerrak, the NCW is transported northwards by the Norwegian Coastal Current (NCC; Saetre 2007).
A reduction in K FW NON in the NCW downstream from Skagerrak, is expected according to two processes. The first is CDOM degradation caused by bacterial consumption and Opdal et al. Tracking freshwater browning photo oxidation/degradation (Fransner et al. 2016;Massicotte et al. 2017). The second is freshwater drainage to the NCW from additional sources along the western Norwegian coast. Because western Norwegian drainage areas are less vegetated and more alpine than those further east, the light attenuation from these freshwaters will be lower (cf. Thrane et al. 2014). Hence, we expect that the NCW freshwater endmember will gradually be more influenced by the clearer freshwater of western Norwegian origin than that of Baltic origin, as distance from the Skagerrak increases (referred to as replacement effect below). Although the effect of CDOM degradation rate on K FW NON (deg, m À1 km À1 ) downstream the NCC is unknown, we estimate the freshwater replacement effect (fwr, m À1 km À1 ) from information on the discharge (m 3 s À1 ) and light attenuation (m À1 ) of the added freshwater along the coast (see below). In a previous study (Aksnes 2015; see Supporting Information), the total reduction in K FW NON with distance was estimated, tot = 0.058 m À1 km À1 . Thus, the effect from CDOM degradation, deg, can be inferred from the difference, deg = tot À fwr. Knowing deg and the NCC transportation speed (see below) of a CDOM molecule (km d À1 ) provides an estimate of the CDOM degradation rate (m À1 d À1 ).
Freshwater replacement effect on K FW NON The mean winter (December-February, 1990-2020) freshwater discharge (FW, m 3 s À1 ) into the NCW was retrieved for six drainage regions (A 1 to A 6 ) between Skagerrak and the Barents Sea (see Supporting Information). The corresponding nonchlorophyll light attenuation of this freshwater was drawn from measurements taken in 2011 in 13 lakes draining to the Western Norwegian coastline (Supporting Information Fig. S6; Thrane et al. 2014). Starting in Skagerrak (A 0 ) we then obtained updated estimates of K FW NON for each additional drainage area (A i ) using the fraction (f i ) between the freshwater of Skagerrak (Baltic) origin (FW(A 0 ) = 22,000 m 3 s À1 ) and the added freshwater from the respective drainage areas along the western Norwegian coast (FW(A i ), m 3 s À1 ) where Transport and retention of a CDOM molecule The drift trajectory and speed of 100,000 neutrally buoyant virtual CDOM molecules from the Skagerrak to the Barents Sea, was estimated according to advection and turbulent diffusivity, using a Lagrangian particle tracking model in the Python-based framework OpenDrift (Dagestad et al. 2018). See Supporting Information for details.

Results and discussion
Overall, we found that variation in Secchi disk depth (S) was adequately described by a statistical model including Chl a, salinity, and year (A1970, Eq. 1; Fig. 2; Supporting Information Fig. S3). The relationship between the light attenuation coefficient and Chl a concentration (Fig. 2C), is close to the general function K PHY = 0.121ÁCHL 0.428 derived by Morel (1988). This adds confidence to approach A1903 (Eq. 2), where K PHY is derived directly from CHL (Supporting Information Eq. S2) using this function. Also, the relationship between light attenuation and the freshwater fraction (Fig. 2D), resembles the functional form described by Højerslev et al. (1996) for CDOM light absorption a(380), suggesting the model captures some fundamental mechanisms of light attenuation in the region. An apparent increase in light attenuation (0.8-1.1 m À1 ) of the nonphytoplankton freshwater endmember (Fig. 2E) from 1970 to 2020, could indicate increasing CDOM concentrations in freshwater sources.
During the period 1900 to 2020 the net forest cover in the region increased by around 200,000 km 2 , or by ca 25% (Fig. 3A; Supporting Information Fig. S1). Most of this increase ($ 70%) stem from Fennoscandia (Supporting Information Fig. S1). Over the same period (1903-2021, A1903), we found an increase in the light attenuation, K S , of the southern NCW (Supporting Information Fig. S4A), unrelated to Chl a concentration (Supporting Information Fig. S4A), even at low freshwater fractions. Model A1903 thus suggests a centennial increase in the nonphytoplankton light attenuation of the freshwater endmember in the southern NCW, K FW NON , from 0.6 to 1.4 m À1 (Fig. 3B). This increase was significantly larger than the seasonal variations in light attenuation (Supporting Information Fig. S4D,E). The model-estimate of K FW NON for the year 2008 (1.38 AE 0.06 95% CI) was close to that estimated (for 440 nm) by Aksnes (2015) in the same year (1.47 AE 0.24 95% CI), supporting our model estimate. The increase in K FW NON (Fig. 3B) and its apparent association with forest cover (Fig. 3A) corresponds well to a series of studies suggesting that terrestrial primary production and vegetation coverage (greening) drive freshwater browning through increased concentrations of CDOM (Kritzberg 2017;Creed et al. 2018;Skerlep et al. 2020).
While ca 90% of the centennial forest increase in Europe stem from the conversion from grassland and cropland (Fuchs et al. 2015), some afforestation is also facilitated through peatland ditching which cause organic carbon leaching and freshwater browning (Nieminen et al. 2021;Williamson et al. 2021;Härkönen et al. 2023). Ditching also increases freshwater Fe concentrations (Estlander et al. 2021), further enhancing light absorption (Xiao et al. 2013;Weyhenmeyer et al. 2014). However, dissolved Fe is likely to have little impact on light attenuation in coastal waters, as around 95% The increase in the forest cover (1000 km 2 ) of catchment areas draining to the Baltic Sea, North Sea, and Skagerrak in the period 1900-2020. Colors correspond to the geographical areas in the map insert (Fig. 1). (B) The estimate (and 95% CI) of nonphytoplankton light attenuation of the freshwater fraction (K FW NON ) of the Skagerrak Sea (star in map insert) for the period 1903-2021 using the approach A1903 (Eq. 2). A simple linear regression between forest cover in the Baltic Sea and North Sea drainage area (FC) and K FW NON in the southern NCW, K FW $ α FC + β, gives an R 2 = 0.90 (α = 0.0034 m À1 1,000 km À2 AE 3 Â 10 À4 95% CI, N = 13).  Morel (1988) and used in A1903 (Eq. 2). A simple linear regression between our modeled K S (CHL) and the K PHY from Morel (1988), K PHY $ K S (CHL) α + β, gives an R 2 = 0.98 (α = 1.03 AE 0.03 95% CI and β = 0.01 AE 0.006 95% CI). (D) The influence of freshwater (FWf) on K S at zero CHL in 2021; (E) the change in the nonphytoplankton light attenuation of the freshwater endmember, K FW NON (CHL = 0, FWf = 1). See also Supporting Information Fig. S4 for a detailed comparison between the model predictions and observations. is removed through flocculation during estuarine mixing (Sholkovitz 1976;Sholkovitz 1978). In comparison, flocculation only removes 3-16% of the DOC (Sholkovitz 1976;Asmala et al. 2014). Additional drivers of freshwater browning are changed hydrology (de Wit et al. 2016), and the reduction in atmospheric sulfur deposition (Evans et al. 2006;Monteith et al. 2007).
While K FW NON is clearly influenced by CDOM, wind driven resuspension of bottom sediments has also been suggested to decrease water clarity (Capuzzo et al. 2015;Wilson and Heath 2019). By selecting our location in the southern NCW, more than 20 km offshore and at a depth of 520 m, the effect of riverine borne or bottom resuspended particulate matter should be minimized (Opdal et al. 2019;Thewes et al. 2021). We also note that our models account for only 39% (Eq. 1; Supporting Information Fig. S3), 47% (Supporting Information Eq. S1; Fig. S4) and 36% (Supporting Information Eq. S2; Fig. S4) of the deviance in the observations, leaving ample room for additional sources of variation (see Supporting Information).
Having established a centennial change in the freshwater endmember in the southern NCW, mainly of Baltic origin (Aure et al. 1998), we then investigated how far this Baltic freshwater persist in the NCW. Results from two previous studies (Thrane et al. 2014;Aksnes 2015) suggest that the K FW NON of the NCW (Fig. 4C) was higher than the average of 13 Western Norwegian lakes ( Fig. 4C; Supporting Information   Fig. S6). Even though measurements of light attenuation in the lakes are not directly transferable to those at the river outlet, the steep terrain, small watersheds, and short transit times in Western Norwegian catchments (Thrane et al. 2014), will likely minimize the discrepancy compared to other landscapes.
Along the Norwegian coast we found that a total decrease in K FW NON , (0.058 m À1 1000 km À1 , Aksnes 2015) could be explained by the freshwater replacement effect (0.0575 m À1 1000 km À1 ; Eq. 3). This suggests little or no degradation of CDOM between the southern NCW and Lofoten, and that the reduction in K FW NON over distance is due to increased influence of relatively clear freshwater from western Norwegian water sheds (Thrane et al. 2014). Most of the degradation of CDOM from microbial or photochemical oxidation happens in lakes and rivers, and within the first few months after entering the coastal waters (Vähätalo and Wetzel 2004;Timko et al. 2015;Allesson et al. 2021). This leaves more recalcitrant CDOM further away from the freshwater sources (Tranvik et al. 2009;Fransner et al. 2016;Massicotte et al. 2017), which in our case are mainly located around the Baltic Sea. A recent estimate suggests that 80% of the terrestrial DOC in the Baltic Sea is degraded before it reaches the Skagerrak (Fransner et al. 2016). In addition, we can assume an in situ production of new CDOM by marine phytoplankton accounting for up to 25% of the total CDOM (Yamashita and Tanoue (2015), the total reduction (tot) of K FW NON from Skagerrak to the Barents Sea was 0.058 m À1 1,000 km À1 at 440 nm. This estimate was based on measurements from 40 locations downstream the NCC (A, squares). The freshwater replacement model (see Eq. 3) calculates the expected reduction (fwr), in K FW NON from the Skagerrak Sea (A 0 , orange color) and northwards (arrows in A) based on the addition of new (clearer) freshwater from Western Norwegian drainage areas (A 1 -A 6 , A,B). Light attenuation of this added freshwater (shaded area in C) is the average of 13 Western Norwegian lakes (circles in A) measured by Thrane et al. (2014), and was below that of the NCW (squares in C). The freshwater replacement model (Eq. 3) predicts a reduction in K FW NON of 0.0575 m À1 1,000 km À1 (at 440 nm, shown as circles in D) which is very close to the total reduction estimated by Aksnes (2015). This suggests minor effect of CDOM degradation (deg) in the NCW. Also shown in the top x-axis in D is the mean transportation time of a CDOM molecule from the Skagerrak to the Barents Sea, based on the drift simulation model shown in Supporting Information Fig. S7. phytoplankton concentrations in the NCW appeared unchanged (Supporting Information Fig. S4A), this should not affect the overall patterns described.
Moreover, we found that the transportation speed of a virtual CDOM molecule from Skagerrak to the Barents Sea is measured on the order of months (Supporting Information Fig. S7), rather than years and decades as in the Baltic Sea (Fransner et al. 2016). The median horizontal speed of a CDOM molecule was 0.11 m s À1 , which is close to the reported average speed of the NCC (Saetre 2007). After 150 and 350 d, 50% of the molecules had reached the Lofoten and Barents Sea, respectively. Hence, the absence of a CDOM degradation effect could thus either be due to a recalcitrant fraction of CDOM (Fransner et al. 2016), or simply due to the short transportation time, or both.
In summary, we propose that a centennial increase in vegetation cover in Northern Europe has led to browner freshwaters due to increased CDOM concentrations. This browning and subsequent decrease in light attenuation is detectable downstream in the freshwater endmember of the NCW across several thousand kilometers. This connects the Barents Sea to the Baltic Sea, and possibly to the surrounding lakes and forests, pointing to an ecosystem connectivity from land to sea over great distances.