Acidiﬁcation of the Nordic Seas

. Due to low calcium carbonate saturation states, and deep winter mixing that brings anthropogenic carbon to the deep ocean, the Nordic Seas and their cold-water corals are vulnerable to ocean acidiﬁcation. Here, we present a detailed investigation of changes in pH and aragonite saturation in the Nordic Seas from pre-industrial times to 2100, by using in situ observations, gridded climatological data, and projections for three different future scenarios with the Norwegian Earth System Model 5 (NorESM1-ME). During the period of regular ocean biogeochemistry observations from 1981-2019 the pH decreased with rates of 2-3 10 − 3 yr − 1 in the upper 200 m of the Nordic Seas. In some regions, the pH decrease can be detected down to 2000 m depth. This resulted in a decrease of the aragonite saturation state, which now is close to undersaturation in the depth layer of 1000-2000 m. The model simulations suggest the pH of the Nordic Seas to decrease at an overall faster rate than the global ocean 10 from preindustrial to 2100, bringing the Nordic Seas pH closer to the global average. In the esmRCP8.5 scenario, the whole water column is projected to be undersaturated with respect to aragonite at the end of the 21st century,

There has been extensive research on changes in the carbonate system and pH in the Nordic Seas, facilitated by the many research and monitoring cruises in the area (e.g., Olsen et al., 2006;Ólafsson et al., 2009;Skjelvan et al., 2008;Chierici et al., 2012;Skjelvan et al., 2014;Jones et al., 2020;Skjelvan et al., 2021;Pérez et al., 2021). Between the 1980s and 2010s, the pH has been shown to decrease with rates of -0.0023 to -0.0041 y −1 in surface waters, which is greater than expected from the increase in atmospheric CO 2 alone (Ólafsson et al., 2009;Skjelvan et al., 2014). This is consistent with the many observations 60 that have indicated a weakening of the pCO 2 undersaturation of the Nordic Seas surface waters, i.e., that surface ocean pCO 2 has risen faster than the atmospheric pCO 2 (Olsen et al., 2006;Skjelvan et al., 2008;Ólafsson et al., 2009), over the past decades. The future pH of the Nordic Seas have been assessed with different modelling approaches (Bellerby et al., 2005;Skogen et al., 2014Skogen et al., , 2018. Bellerby et al. (2005) investigated the impact of climate change on the Nordic Seas CO 2 system under a doubling of the atmospheric CO 2 to a value of 735 ppm. It was done by combining observed relationships between carbonate and silicious minerals. Together, this forms the following equilibria: (1) Combined, the concentration of CO 2 , H 2 CO 3 , HCO − 3 , and CO 2− 3 , constitute the concentration of dissolved inorganic carbon (C T ). In seawater, approximately 90% of C T exists in the form of HCO − 3 , 9% as CO 2− 3 and 1% as CO 2 .

90
As seen from Equations 1 -3, the dissolution of CO 2 in seawater results in an increase in H + concentration, which leads to a decrease in pH. On total scale, pH is defined as: where HSO − 4 is sulphate. Apart from C T , pH is influenced by temperature, salinity, and total alkalinity (A T ). A T is mostly determined by HCO − 3 and CO 2− 3 (carbonate alkalinity). Temperature and salinity affect pH by altering the dissociation constants 95 and thus the partitioning of C T between its different constituents. The relation between C T and A T influences pH by affecting the buffer capacity of seawater. The qualitative, direct effects of an increase in each property are shown in Table 1. Note that this Table does not consider indirect effects on pH, for example from the change in air-sea fluxes that will follow from e.g., a temperature driven pCO 2 change (e.g. Jiang et al., 2019;Wu et al., 2019). Equations 1 to 3 show that an increase in anthropogenic CO 2 and C T results in a reduction in CO 2− 3 . This affects the Sarmiento and Gruber, 2006;Orr, 2011). Waters with a higher buffer capacity, i.e., higher concentrations of CO 2− 3 , have the capability of converting a larger fraction of the absorbed CO 2 into bicarbonate. A smaller fraction remains as dissolved CO 2 , implying a smaller increase in the seawater pCO 2 . These waters therefore have the capability of absorbing more CO 2 for any given increase in atmospheric pCO 2 (assuming a uniform increase in pCO 2 between water masses), which also implies a larger decline in CaCO 3 saturation state. pH is, on the contrary, decreasing more in waters with lower buffer capacity as they are less 115 effective in neutralising carbonic acid.

Observational data
As observational data, we used C T , A T , temperature, salinity, phosphate, and silicate data collected between 1981 and 2019 during dedicated research cruises, at two time-series stations, and in the framework of the Norwegian program "Monitoring 120 ocean acidification in Norwegian waters". Sampling locations are shown in Fig. 1.
Time-series data are from Ocean Weather Station M in the Norwegian Sea, and from the Iceland Sea. The data from the 130 Ocean Weather Station M, located at 66 • N and 2 • E, have been described in Skjelvan et al. (2008). At this station, sampling at 12 depth levels between surface and seabed (2100 m) was carried out each month between 2002 and 2009, and 4-6 times each year between 2010 and 2019. For these data, the uncertainty related to the measurements is 0.001 for salinity, 0.7 µmol kg −1 for silicate, 0.06 µmol kg −1 for phosphate, and 2 µmol kg −1 for C T and A T . The time-series station in the Iceland Sea, covering the period of 1985-2019, is situated at 68 • N and 12.67 • W. It is visited approximately 4 times a year and samples are 135 taken at 10-20 depth levels between surface and seabed (1900 m). The uncertainty related to the measurements at this station is 0.005 for salinity, 2% for silicate, 2% for phosphate, and 4 µmol kg −1 for both C T and A T . These data have been described in Ólafsson et al. (2009).
The data from the program "Monitoring ocean acidification in Norwegian waters" were sampled in the full water column along repeated sections in the Nordic Seas in the period 2011-2019 (Chierici et al., 2012(Chierici et al., , 2013(Chierici et al., , 2015(Chierici et al., , 2016 The filled contours illustrate the bathymetry at 250 m intervals. et al., 2018, 2019, 2020). The uncertainties related to the sampled data is 0.005 for salinity, 0.1 for silicate, 0.06 for phosphate, 2 µmol kg −1 for both C T and A T .
Data from the eastern Fram Strait were collected on cruises with RV Helmer Hansen within the CarbonBridge project, and on cruises with RV Lance (Chierici et al., 2019b) organized by the Norwegian Polar Institute.
Analytical methods for C T and A T , in all datasets described above (for GLODAP after the mid 1990s), follow Dickson 145 et al. (2007) and the accuracy and precision is controlled by Certified Reference Materials (CRM), and by participation in international intercomparison studies (e.g. Bockmon and Dickson, 2015).
For estimates of atmospheric CO 2 change, we used the annual mean atmospheric CO 2 mole fraction (xCO 2 ) from the Mauna Loa updated records, downloaded from www.esrl.noaa.gov/gmd/ccgg/trends/. Although the absolute value of atmospheric xCO 2 varies with latitude, the growth rates are the same across the globe.

Gridded climatological data
Climatological distributions of pH and Ω Ar were calculated from C T , A T , temperature, salinity, phosphate and silicate in the mapped GLODAPv2 data product (Lauvset et al., 2016). Pre-industrial pH was determined by subtracting the GLODAPv2 estimate of anthropogenic carbon from the mapped climatology of present (i.e., 2002) C T (Lauvset et al., 2016). We assumed that the changes in the temperature, salinity and A T of the Nordic Seas are of minor importance for changes in pH between 155 pre-industrial times and present-day. The GLODAPv2 estimate of anthropogenic carbon have been calculated using the transit time distribution (TTD) approach. We note that we use the GLODAPv2 estimate of pre-industrial pH only for comparison with the ESM data, specifically in Fig. 4 (5.2 ).

Earth System Model data
For the estimates of past and future ocean acidification and saturation states under various climate scenarios, we primarily used 160 output from the fully coupled Norwegian Earth System Model with interactive atmospheric CO 2 (NorESM1-ME, Bentsen et al., 2013;Tjiputra et al., 2013Tjiputra et al., , 2016. NorESM1-ME includes the dynamical isopycnic vertical coordinate ocean model MICOM (Bleck and Smith, 1990) and the Hamburg Oceanic Carbon Cycle model (HAMOCC5, Maier-Reimer et al., 2005), adapted to the isopycnic ocean model framework. HAMOCC5 simulates lower trophic ecosystem processes up to the zooplankton level, including primary production, remineralization and predation, and full water column inorganic carbon chemistry. For 165 our assessment, we utilised emission-driven historical simulations for the period from 1850 to 2005 and future scenarios simulations for the period from 2006 to 2100, with focus on RCP's 2. 6, 4.5 and 8.5 (RCP2.6, RCP4.5, and RCP8.5;Meinshausen et al., 2011;van Vuuren et al., 2011a). RCP2.6 represents a mitigation scenario, RCP4.5 a stabilization scenario and RCP8.5 a high-emission scenario. For the emission-driven runs used here, the corresponding scenarios are named esmRCP2.6, esm-RCP4.5 and esmRCP8.5. Because the emission-driven scenarios prognostically simulate the atmospheric CO 2 concentration, 170 it normally deviates from the prescribed concentrations in the concentration-driven scenarios (e.g. Friedlingstein et al., 2014). This is most critical for the historical scenario, where the prescribed atmospheric CO 2 follows the observed evolution. Here, deviations in the simulated atmospheric CO 2 might result in pH changes that differ from the actual pH change. The deviation in the simulated atmospheric CO 2 concentration in the emission-driven NorESM1-ME scenarios, from the prescribed one in the concentration-driven scenarios, and its effect on pH, is shown in the supplementary Table S1. Between 1850 and 2005, the 175 model simulates an increase in the atmospheric CO 2 that is 14 ppm too strong, which results in a pH drop that exceeds the expected one by 0.01. This deviation is, however, one order of magnitude smaller than the actual pH change between 1850 and 2005, and has the same order of magnitude as the estimated uncertainty in both observational data (Table 2) and GLODAPv2 pre-industrial pH estimate in the Nordic Seas (Sect. 4.4). The impact of the historical atmospheric CO 2 deviations between emission driven and concentration driven on pH change in our results is therefore negligible. Prior to experiments, NorESM1-180 ME has undergone an extended spin-up procedure (>1000 years). The changes in pH, in all considered depth layers, is minor (more than one order of magnitude less) in the preindustrial control simulation compared to the historical run and the future scenarios, indicating that the impact of model drift on our results is insignificant .
As a means of uncertainty assessment, we use the outputs from an ensemble of emission-driven ESMs that participated in CMIP5 (Taylor et al., 2012). We chose emission-driven rather than concentration-driven scenarios, as they include the 185 feedback of the carbon cycle on the physical climate (Booth et al., 2013) and thus give a more comprehensive estimate of the effect of model-related uncertainties on climate projections, and in particular on atmospheric CO 2 , ocean carbon uptake and ocean acidification. It is well known that the inter-model spread is larger in emission-driven scenarios than in concentrationdriven ones (Booth et al., 2013;Friedlingstein et al., 2014). While NorESM1-ME outputs are available for low to high future emission scenarios, the CMIP5 data-portals only contains emission-driven ESM outputs for the high future emission scenario 190 (esmRCP8.5). Our ESM-ensemble consists of all ESMs that participated in the experiment 'esmrcp85' and RCP8.5, and whose output is publicly available in one of the CMIP5 data portals and contains all variables needed for our analysis. This results in an ensemble of 7 ESMs; 1) CESM1(BGC) (The Community Earth System Model, version 1 -Biogeochemistry, Long et al., 2013), 2) CanESM2 (second-generation Canadian earth system model, Arora et al., 2011)

Cold-water coral positions
To estimate the potential impact of the Nordic Seas acidification on cold-water corals, we used habitat positions in longitude and latitude from EMODnet Seabed Habitats (www.emodnet-seabedhabitats.eu) together with information on depth from ETOPO1 (NOAA National Geophysical Data Center, 2020).

Spatial drivers of pH and saturation states
To identify drivers of observed spatial variability of surface pH and Ω Ar , we calculated pH and Ω Ar by using spatially varying GLODAPv2 climatologies of specific drivers in Table 1, while keeping all other drivers constant (set to the spatial mean value of the Nordic Seas surface waters). Because the relation between C T and A T is a proxy for the buffer capacity, we decided to look at their combined effect on pH, meaning that both changes in C T and A T are included in the calculations. Their combined 210 effect we from now on refer to as C T +A T . First, pH and Ω Ar were calculated with temperature being the only spatially varying climatology (pH(T), Ω Ar (T)). Thereafter, we used spatially varying temperature, C T and A T climatologies to calculate pH(T, C T , A T ) and Ω Ar (T, C T , A T ). Finally, the salinity variability was added to estimate pH(T, C T , A T ,S) and Ω Ar (T, C T , A T , S).
To estimate the constribution of each driver, the pH and Ω Ar fields calculated with the different spatially varying drivers were thereafter correlated with the actual pH and Ω Ar of the Nordic Seas. Measurements of temperature, salinity, C T , and A T (Figs. S1-S4), phosphate, and silicate from the data sets described in Sect.
3.1 were used to calculate pH and Ω Ar , using CO2SYS for MATLAB (Lewis and Wallace, 1998;van Heuven et al., 2011). pH was calculated on total scale at at in situ pressure and temperature. Wherever nutrient data were missing, silicate and phosphate 220 concentrations were set to 5 µmol kg −1 and 1 µmol kg −1 , respectively, which are representative values for the Nordic Seas.
For the CO2SYS calculations, the dissociation constants of Lueker et al. (2000), the bisulfate dissociation constant of Dickson (1990) and the borate-to-salinity ratio of Uppström (1974) were used. This ratio has recently been shown to be suitable for the western Nordic Seas (Ólafsson et al., 2020a). Present-day trends (1981Present-day trends ( -2019 in pH, and Ω Ar were determined for six different regions in the Nordic Seas: the Norwegian 225 Basin (NB), the Lofoten Basin (LB), the Barents Sea Opening (BSO), eastern Fram Strait (FS), the Greenland Sea (GS) and the Iceland Sea (IS) (Fig. 1). These regions were chosen based on the data-availability, being centered around stations and sections where repeated measurements are taken, but also to get a representation of the main surface water masses of the Nordic Seas.
In the surface, the Norwegian Basin, Lofoten Basin, and Barents Sea Opening are influenced by relatively warm and salty northward flowing Atlantic Water, while the Greenland and Iceland Seas are influenced by relatively cold and fresh southward 230 flowing polar waters. As the Fram Strait surface is influenced by Atlantic and polar waters, we constrain the Fram Strait box to the east (hereinafter referred to as eastern Fram Strait) to ensure that it mostly represents Atlantic Waters. The geographical range of each regional box is kept small so that aliasing effects of latitudinal and longitudinal gradients are minimized.
Regional trends were computed from annual means for five different depth intervals (0-200, 200-500, 500-1000, 1000-2000, and 2000-4000 m) using linear regression. The relatively thick upper layer was chosen to keep all depths influenced by the sea-235 sonal cycle in one layer, that is, to minimize the number of layers where the trends may be affected by seasonal undersampling.
As the winter mixed layer reaches approximately 200 m (e.g. Skjelvan et al., 2008), this depth sets the approximate lower limit for the impact of seasonal variations. The significance of the trends (at 95% confidence level), were determined from the pvalue of the t-statistic (as implemented in MATLAB's fitlm function). For the comparison of trends, 95% confidence intervals of the slopes were determined by the use of the Wald method (as implemented in MATLAB's fitlm and coefCI functions).

240
The observed long-term changes in pH were decomposed into contributions from changes in temperature (T), salinity (S), C T and A T (Figs. S1-S4, and Tables S2-S5), following the procedure of Lauvset et al. (2015). First, the effect of each of these processes on the CO 2 fugacity (f CO 2 ) is determined following Takahashi et al. (1993): The long-term mean values for the sensitivities (the f CO 2 partial derivatives) were approximated as in Fröb et al. (2019).

245
Changes in A T and C T are driven by biogeochemical processes, transport, mixing and dilution or concentration by freshwater fluxes, which is in direct proportion to the dilution or concentration of salinity. The freshwater-effect can be separated by introducing salinity-normalized C T (sC T ) and A T (sA T ) (Keeling et al., 2004;Lovenduski et al., 2007): Here we set S 0 to 35 (Friis et al., 2003) and used the intercept of Eq. 6 and 7 in Nondal et al. (2009) as the non-zero 250 freshwater end-member (A 0 and C 0 ). Substituting A T and C T in Eq. 6 by Eqs. 7 yields: Finally, H + in equation 10 was converted to pH by acknowledging that dpH = −([H + ]ln(10)) −1 d[H + ]. Here we consider the sulphate in Eq. 4 to be negligible, and did therefore not include it.
To control whether the observed pH changes are consistent with the changes in atmospheric CO 2 , we additionally determined the pH change that can be expected in seawater where the pCO 2 perfectly tracks the atmospheric pCO 2 (pH perf ) for each 260 region. This was achieved by adding the observed change in atmospheric xCO 2 to the local mean pCO 2 for the first year with observations, and then calculating the pH with CO2SYS with the local temperature, salinity, A T , phosphate and silicate and their respective changes as inputs. We applied no corrections for water vapour and atmospheric pressure as the rates of change for xCO 2 and pCO 2 are proportional. Any deviation between observed pH change and pH perf is a consequence of changes in seawater pCO 2 that are smaller/larger than in the atmosphere, i.e., a change in the air-sea pCO 2 difference. As described in Sect. 2, the total change in pH and saturation states does not only depend on local changes in C T , A T , temperature, salinity, and nutrients, but also on the initial buffer capacity of the seawater. For the calculation of past and future pH changes, we use ESM data, which is usually biased, i.e., there is an offset between modelled fields and reality and this also holds for the buffer capacity. In particular, NorESM1-ME has high A T and low C T relative to observations in deep waters, 270 leading to biased high pH (Fig. S5) and saturation states (not shown). To alleviate this bias in our analysis of past and future pH and Ω Ar , we applied the modelled change of temperature, salinity, C T , A T , phosphate and silicate to the gridded GLODAPv2 climatology. Here, the modelled change between pre-industrial, present-day and future were calculated as differences between 10-year means; i.e., 1850-1859, 1996-2005 and 2090-2099, respectively. We note that we could not center our present-day 10-year mean around the year 2002 to which the GLODAPv2 climatology is normalized as the future scenarios start in 2006.

275
After we obtained past and future states of the properties listed above, we calculated past and future pH, Ω Ar and Ω Ca in CO2SYS. Similar procedures have been used by Orr et al. (2005) and Jiang et al. (2019) to calculate future pH. Additionally, we used these data to calculate the drivers of past-to-present and present-day-to-future pH changes, following the methodology described in the previous section. Here we used a value of zero for the freshwater end members A 0 and C 0 as NorESM1-ME does not include any riverine input of A T and C T .

280
To estimate the impact of acidification on the cold-water corals of the Nordic Seas, we calculated the mean saturation state in our region east of 0 • E, and south of 64 • N, for P.I., present day and for the future under the esmRCP2.6, esmRCP4.5 and esmRCP8.5 scenarios. The exclusion of the western and northern parts was done to constrain the mean to the Atlantic Water where the cold-water corals are located. The saturation horizon was defined as the deepest vertical grid cell where Ω Ar > 1.
In order to facilitate a comparison with other model-based acidification studies, we have chosen to present the past and future

pH or H + change?
In a recent publication, Fassbender et al. (2021) recommend to analyze changes in H + concentrations in addition to changes 290 in pH, when comparing pH trends across water masses with different initial pH. The underlying reason is that a change in pH represents a relative change, and that it is possible to obtain the same pH changes across water masses with different change in H + concentration. We estimated the sensitivity of our results to the choice between pH and H + by plotting the change in H + concentration against the change in pH, for a given change in C T at various initial pH (Fig. 2). The different initial pH were obtained by varying the C T over A T ratio, and the calculations were done with a temperature and salinity of 5 • C linear relationship breaks down, if pH decreases as a result of an increasing C T over A T ratio. The maximum pH change takes place at the buffer minimum, which is close to where C T =A T , approximately at (pK1+pK2)/2 (Frankignoulle, 1994;300 Fassbender et al., 2017;Middelburg et al., 2020), which in our example is at pH 7.6. The linear relationship between the H + and pH change does therefore not hold for pH ranges where relatively low initial pH values are included, as is the case for the examples in Fassbender et al. (2021), as well as for larger C T changes. In these cases it is more appropriate to use H + for diagnosing ocean acidification.

Uncertainty analysis 305
There are several sources of uncertainties (σ) involved in our calculations of pH and Ω: measurement uncertainty (σ mes ), mapping uncertainty (σ map ) for the gridded product, and uncertainties related to dissociation constants (σ Kx ) used in the CO2SYS calculations. To estimate the total uncertainties in our calculations of pH and Ω, we used the error propagation routine in the MATLAB version of CO2SYS (Orr et al., 2018). The uncertainties in the input parameters (A T , C T , temperature, salinity, phosphate and silicate) were set to σ mes for the single measurements, and σ 2 mes + σ 2 map for the mapped product as well as field) from Lauvset et al. (2016), were used, respectively. The correlation between uncertainties in A T , C T were set to 0. This is a reasonable assumption given that C T and A T are measured on different instruments using different analytical methodologies.
In addition, including a positive correlation term would decrease the overall uncertainty and we prefer a potential overestimation. For the dissociation constants the default uncertainties in the errors.m function were used. From here on, the calculated 315 uncertainties will be presented as σ point for discrete data, when σ Kx and σ mes are included, and σ f ield for 3D data, when σ Kx , σ mes and σ map are included.
For the observations described in Section 3.1, the mean, maximum and minimum uncertainties (σ point ) for our calculations of pH, Ω Ar , Ω Ca and pCO 2 are listed in Table 2. Variations in the uncertainties arise from variations in temperature and salinity, which impact the uncertainty of dissociation constants. While systematic uncertainties would tend to cancel out when 320 calculating trends (i.e.,comparing measurements from the same location but from different times), random uncertainties would not (Orr et al., 2018). Therefore, to estimate to what extent these uncertainties could impact our trend estimates, we further investigated whether there is any trend in the uncertainties. This is discussed in Sect. 5.4.
For the GLODAPv2 estimate of pre-industrial C T there is an additional uncertainty coming from the TTD method that was used to calculate the anthropogenic CO 2 . He et al. (2018) published a thorough analysis of the different sources of uncertainty 325 in this method, and concluded that the overall uncertainty is 7.8-13.6%. Combining this with the mapping errors, Lauvset et al.
(2020) estimate that the global ocean anthropogenic carbon inventory calculated from the mapped fields is 167±29 PgC. This results in an uncertainty of 0.02 in the pre-industrial Nordic Seas upper layer pH. In the trends of the uppermost layer (0-200 m), there is also an uncertainty related to seasonal undersampling. Most samples (about 60% in total) from the data sets described in Sect. 3.1 were collected during spring and summer (April-September,

330
Figs. S8 -S13). The uneven sampling frequency of different seasons introduces uncertainty in the annual means of the uppermost ocean layer, and can lead to biases in our trend estimates. Unfortunately, there are not enough data to allow for de-seasonalization in order to remove such potential biases. Therefore, to get an idea of the effect of seasonal undersampling we additionally calculated trends by using annual means containing samples from the productive season only, both for a longer period (April-September) to include both the spring bloom and the summer production, and for a shorter period (June-August),

335
to include only the summer season.
Modelled future projections are uncertain due to incomplete understanding or parameterization of fundamental processes, as well as different and unknown future carbon emission scenarios (Frölicher et al., 2016). Because this study primarily focuses Table 3. Spatial correlation (r) and explained variance (r 2 , in paranthesis) between pH and pH(T), pH(T,CT ,AT ) and pH(T,CT ,AT , S), and that we do not explicitly account for in this study. However, a large part of this variability is eliminated because we use 10 year means for the future and past estimates of pH.

Results and discussion
This Section is organized as follows: we will start to describe the present distribution of pH and Ω Ar and its drivers (Section 5.1). In Section 5.2, we give an overview of pH changes from pre-industrial to 2100. Thereafter we describe regional changes 350 from pre-industrial to present-day (Section 5.3), present-day changes (Section 5.4), changes from present-day to future (Section 5.5) and assess its impacts on cold-water corals (Section 5.6). In Section 5.7 we analyze the drivers of pH change in the different time periods.

Present-day spatial distribution of pH and Ω saturation states
Due to the contrasting properties of Atlantic waters, here defined as waters with salinity > 34.5, (Malmberg and Désert, 1999;355 Nondal et al., 2009) and polar waters (defined as the waters with salinity < 34.5 detached from the Norwegian coast) that meet and mix in the Nordic Seas, there are large spatial gradients in surface (0 m) temperature, salinity and chemical properties ( Fig.   3 and S14). The Atlantic Water, located in the eastern part of the Nordic Seas, is characterized by higher temperature, salinity, and A T , while polar waters are colder and fresher with lower A T . This results in a decrease in temperature, salinity, and A T from southeast to northwest. Within the Atlantic water there is a tendency of increasing C T with decreasing temperature. This 360 is largely as a consequence of the increased CO 2 solubility in colder water, i.e., a cooling of a water mass result in an increase in C T due to an uptake of CO 2 from the atmosphere. In polar waters, C T is lower than in Atlantic waters due to the lower (c,f), temperature and CT + AT (d,g) and temperature, CT + AT and salinity (e,h) in pH and ΩAr, calculated as described in Section 4.1 in Atlantic Water (red) and low salinity waters (blue). Each circle represents a value from a single grid cell.
pCO 2 (Fig. S14)), and also as a result of the large freshwater export from the Arctic Ocean that dilutes not only C T , but also A T and salinity. and the pH calculated with spatially varying temperature only (pH(T)), keeping all other drivers constant, is 0.58. This means that temperature-induced variations (through the thermodynamic effect) are able to explain 34% of the spatial variability in pH ( Fig. 3 and Table 3). Adding C T +A T and salinity contributions explains an additional 55% and 11%, respectively, of the spatial variability in pH. The effect of salinity is largest in the low-salinity regions, i.e., in polar waters and the Norwegian coastal waters.
In contrast to what is suggested by directly correlating pH and C T +A T (Table S9), the results in Table 3 show that 370 C T +A T are important contributors to spatial variations in pH. This indicates that the influence of C T and A T on pH is masked out by temperature variations in Table S9 and Fig. S15, which can be explained by the two cancelling effects that temperature has on pH (Jiang et al., 2019). For example, while the instantaneous, thermodynamic effect of a drop in temperature leads to a pH increase, it also results in a drop in pCO 2 , which subsequently leads to an anomalous CO 2 uptake from the atmosphere.
This increases the C T /A T ratio, which in turn causes a drop in pH that counteracts the initial thermodynamic affect.

375
The saturation state Ω Ar show an opposite pattern to pH, with low saturation states in polar waters, and high saturation states in Atlantic Water. From Fig. 3,d, it becomes clear that the temperature effect on the solubility of Ω Ar (Ω Ar (T)) only explain 11% of the observed Ω Ar range, although it is able to explain 98% of the variability. When adding C T + A T contributions, the observed range in Ω Ar is reproduced, and 100% of the variability is explained. C T +A T strongly influences Ω Ar , because with decreasing C T to A T ratio, the CO 2− 3 concentration decreases as well. The C T to A T ratio itself strongly correlates with 380 temperature as the CO 2 solubility increases with decreasing temperature and vice versa (S9). The strong correlation between Ω Ar and temperature (Table S9)   with a limited spatial and temporal coverage, its representativeness for the entire Nordic Seas is questionable, and we do not expect an exact agreement with the model. For example, the stronger trend obtained from the observational data might be a 400 result of the samples in the beginning of the period being biased to regions with higher pH.
The future evolution of upper layer pH in the Nordic Seas depends strongly on the CO 2 emission scenario (Fig. 4). In the esmRCP2.6 scenario, where the CO 2 emissions are kept within what is needed to limit global warming to 2 • C (van Vuuren et al., 2011b), pH drops by 0.04 from 2020 to 2099, and reaches a value of 8.03±0.01. Note that in this scenario there is a peak and decline, related to the overshoot profile of the atmospheric CO 2 concentration, with a minimum pH value in mid-century.

405
For the esmRCP4.5 scenario, which corresponds roughly to the currently pledged CO 2 emission reductions under the Paris agreement, the surface pH is simulated to drop by about 0.15, reaching an average value of 7.93±0.01 by the end of the 21st century. Under the high-CO 2 esmRCP8.5 scenario, NorESM1-ME simulates the pH to decrease by 0.40 between 2020 and 2099 to an average value of 7.67±0.02. This equals a pH decline of approximately to -5.00 10 −3 yr −1 . The model related uncertainty in the esmRCP8.5 scenario, measured as the inter-model spread of pH in 2099, displays a pH range of 7.59-7.79 410 in the surface layer (Fig. 4, S5). This spread is larger than that observed in the concentration-driven simulations with the same models, 7.69-7.75, as expected from the increased degrees of freedom brought about by the interactive atmospheric CO 2 .
Within the emission-driven model ensemble, the pH-decline from pre-industrial to the end of the 21st century as simulated by NorESM1-ME is among the strongest, which most likely is a result of a simulated stronger increase in atmospheric CO 2 . A full analysis of the reasons behind the inter-model spread is beyond the scope of this paper.

415
The simulated Nordic Seas average upper layer pH is 0.11 higher than the global average in 1850, which is related to the undersaturation of CO 2 in the surface waters of the Nordic Seas (Jiang et al., 2019). Our global average pH is about 0.1 lower than that estimated by, e.g., Jiang et al. (2019) for the surface ocean due to our consideration of a 200 m thick upper layer. The difference between the simulated upper layer pH of the global ocean and the Nordic Seas is decreasing with time. By the end of the 21st century, the Nordic Seas upper layer pH is 0.03, 0.07 and 0.08 higher than the global average for the esmRCP8.5, 420 esmRCP4.5 and esmRCP2.6 scenarios, respectively. This is partially a result of the colder waters of the Nordic Seas, which gives them a lower buffer capacity. Additionally, in esmRCP8.5, there is an increase in the pCO 2 undersaturation of the global ocean that increases the global average pH (Fig. S16). Other factors driving this decreasing pH difference between the global ocean and the Nordic Seas can be differential heating. A quantitative assessment of the drivers is beyond the scope of this paper.

Modelled pH and Ω Ar changes from pre-industrial to present-day
In this Section and the following, we present temporal changes in pH and Ω Ar . Note that results for the modelled changes are referring to the 0 m surface, unlike the 0-200 m depth range that we use for the upper layer in Sect. 5.2 and 5.4.
From pre-industrial to present, the spatial pattern of changes in surface pH and Ω Ar are similar (Fig. 5). The strongest decreases, reaching -0.12 and -0.55, respectively, are found in Atlantic Water along the Norwegian coast both for pH and Ω Ar .

430
The smallest change is found in polar waters (see more in depth discussion in Sect. 5.7.2). The corresponding maps for H + (Fig.   S17 ) show a similar spatial distribution as for pH. Due to the longer ventilation time scales of deeper waters, the pH decrease weakens with depth. As shown in the section across 70 • N (Fig. 6), waters below 2500 m are nearly unaffected. While the entire water column remains saturated with respect to calcite, the saturation horizon (Ω=1) of aragonite shoaled from a mean depth of 2200 m (uncertainty range: 2100-2400 m) during pre-industrial, to a present-day mean depth of 2000 m (uncertainty range: 435 1700-2300 m), across this specific section. Note that these depths were obtained from the contour interpolation when creating Fig. 6, which has a finer vertical resolution than the GLODAPv2 climatology.

Observed present-day changes in pH and Ω Ar
Regional trends in observed seawater pH between 1981 and 2019 for five different depth intervals are presented in Fig. 7 and Table 4. The corresponding trends in H + are shown in Fig. S18 and Table S10. In the upper layer (0-200 m), significant trends . Table 4. pH trends ± standard error (10 −3 yr −1 ) calculated from the data presented in Fig. 7    and -2.40 ± 0.23 10 −3 yr −1 , respectively, agrees well with the trend of -2.3 10 −3 yr −1 that they calculated for both regions.
The non-significant trend we find in the Barents Sea Opening is also in agreement with the results of Skjelvan et al. (2014).
In contrast to their results, we obtained a significant trend in the eastern Fram Strait, which may be a result of the larger time Trends of aragonite saturation states are shown in Fig. 8 and Table 5. As for pH, the rate of change is strongest in the series, able to state that there is a significant decrease in Ω Ar in several regions and at several depth layers.
During the period 1981-2019, we detect trends in the uncertainties of pH and Ω Ar (Figs. S6 and S7 ), reaching -0.04 10 −3 yr −1 and 0.53 10 −3 yr −1 , respectively. These are, however, about two orders of magnitude smaller than the trends in pH and Ω Ar , and they do therefore not significantly impact interpretation of our results.

Modelled pH
and Ω Ar changes from present-day to future In this section we go into regional details of future pH and Ω Ar changes under the esmRCP2.6 and the esmRCP8.5 scenarios.
The results are presented for the surface (0 m), and not for the upper layer 0-200 m as in Sect. 5.2 and 5.4.
In esmRCP2.6, a pH decline of 0.06-0.11 in the surface waters is simulated between present-day (1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005) and future (2090-2099) (Fig. 9c). The largest pH decreases are found in polar waters, leading to a weakening of the present-day zonal pH gradient. Surface Ω Ar is projected to decrease by about 0.2-0.5 under esmRCP2.6, with the largest drops taking place in 475 polar waters. Surface waters remain supersaturated with respect to both calcite and aragonite. Interestingly, the strongest ocean acidification occurs at depths of 1000-2000 m in this scenario (Fig. 10c), which leads to a shoaling of the aragonite saturation horizon to a depth of 1100 m (uncertainty range: 800-1200 m). This is discussed in more detail in Sect. 5.7.2.
Under the esmRCP8.5 scenario, surface pH drops by about 0.4-0.5 between present-day and future (Fig. 11), with the largest decreases in polar waters. Surface Ω Ar drops by around 1.1-1.3. In contrast to esmRCP2.6, the largest decline of Ω Ar take 480 place in the Atlantic Water. The reason behind this is discussed in Sect. 5.7.2. The strong ocean acidification in this scenario leads to a reversal of the pH depth-dependency so that pH increases from surface to depth by the end of the 21st century (Fig. region. For all emission scenarios the spatial distribution of H + and its change (shown in Fig. S19 and Fig. S20) are similar to that of pH.

Implications for cold-water corals
Cold-water corals build their structures out of aragonite, which is the more soluble form of calcium carbonate. These corals can,  (Hennige et al., 2015). Furthermore, dead coral structures, which compose the major part of the reefs, cannot resist corrosive waters and experience increased dissolution rates at Ω Ar <1. Cold-water coral reefs, along with their ecosystems, are consequently likely to collapse if the water they live in becomes undersaturated with respect to aragonite.

495
It has been estimated that about 70% of the deep sea corals globally will be below the aragonite saturation horizon by the end-of-the-century under high-emission-scenarios (Guinotte et al., 2006;Zheng and Cao, 2014). Most of the reef sites that have been identified in the Nordic Seas (321 out of the 324 within the region defined in Fig. 1) are at depths of 0-500 m (Fig. 12, see also Buhl-Mortensen et al. (2015)). The aragonite saturation horizon estimated from the GLODAPv2 climatology for present climate is at 2000 m, with uncertainty range 1750-2500 m. Note that the uncertainty 500 range of the depth of the saturation horizon is not equally distributed around the mean because the uncertainty analysis is done for the saturation state, from which the depth distribution is calculated. From the discrete measurements we also see that the waters in the depth range 1000-2000 m are close to being undersaturated with respect to aragonite (Sect. 5.4). For the time being, the saturation horizon is thus well below the majority of the cold-water corals in the Nordic Seas.
In the esmRCP2.6 scenario, NorESM1-ME projects that the aragonite saturation horizon will shoal to 900 m (uncertainty: 505 800-1100 m), while in the esmRCP4.5 scenario the saturation horizon is projected to shoal to 600 m depth (uncertainty: 400-700 m) by the end of this century. This implies that the deepest observed reefs will be exposed to corrosive waters, and thus experience elevated costs of calcification and dissolution of dead structures. The majority (315 out of 324) of the coral sites in the Nordic Seas are, however, found at shallower depths than the projected saturation horizon with uncertainty, although the margins are small. Also García-Ibáñez et al. (2021) suggested that cold-water corals in the subpolar North Atlantic will be exposed to corrosive waters if the 2-degree goal (which is the aim of RCP2.6) is not met. In the esmRCP8.5 scenario, NorESM1-ME projects the whole water column below 20 m (uncertainty: 10-20 m) to be undersaturated with respect to aragonite at the end of this century, such that all cold-water coral reefs in the Nordic Seas will be exposed to corrosive waters.
For esmRCP8.5 the NorESM1-ME results are consistent with our CMIP5 model ensemble that suggests the future saturation horizon lies in the range of 0 and 100 m. Comparison with the CMIP5 ensemble is not possible for esmRCP2.6 and esmRCP4.5 515 because few of the models have performed emission-driven runs under these scenarios. However, NorESM1-ME simulates one of the stronger pH-declines in all depth layers considered in Fig. S5 (Table S6), and has also been shown to be on the upper end of absorption of anthropogenic carbon in the Arctic Ocean (Terhaar et al., 2020a), suggesting that our estimates of the future saturation horizon lies in the shallower end of possible future states.

Present-day drivers
To understand what has caused the observed pH changes presented in Sect. 5.4, we decompose the trends into their different drivers as described in Sect. 4.2 (Fig. 13). In the upper layer (i.e., 0-200 m) the pH decrease in the period 1981-2019 is in agreement (within 95% confidence) with the pH change expected from the increase in atmospheric CO 2 , except for in the Norwegian Basin and the Iceland Sea where the trends are stronger. This is related to a faster increase in the seawater pCO 2 525 compared with that of the atmosphere (Fig. S21), meaning that the pCO 2 undersaturation of the Norwegian Basin and the Iceland Sea is has decreased. We note that this diminishing undersaturation is sensitive to seasons. In the Norwegian basin there is no significant decrease if using data from only April to September and June to August, respectively. In the Iceland Sea the decreasing undersaturation is absent for April-September, but it becomes stronger than the annual mean if using data only from June-August. The sensitivity to the choice of seasons indicates that the strong positive trend in the air-sea pCO 2 530 difference as seen in our dataset can be a result of seasonal undersampling, and that this should be verified with a larger dataset.
Notwithstanding, diminishing pCO 2 undersaturation has been observed in earlier studies of the North Atlantic (Lefèvre et al., 2004;Olsen et al., 2006;Ólafsson et al., 2009;Metzl et al., 2010;Skjelvan et al., 2014), and could be a result of a change in any of the mechanisms underlying the pCO 2 undersaturation in surface waters of the Nordic Seas (see Sect. 1), including cooling of northward flowing Atlantic waters, primary production and the outflow of pCO 2 undersaturated waters from the 535 Arctic Ocean. One other possible mechanism was suggested in Olsen et al. (2006) and Anderson and Olsen (2002), where they associated the fast increase in seawater pCO 2 with a large advective supply of anthropogenic carbon from the south and corresponding changes in the buffer capacity (see also Terhaar et al. (2020b)).
The main driver of the present-day (1981-2019) pH decrease in the upper layer is increasing C T , which primarily is caused by biogeochemical processes (C T bg), including increasing anthropogenic carbon, along with a small freshwater contribution 540 (C T fw) caused by an increasing salinity (Fig. S2). The increasing salinity also results in an increasing A T (Fig. S4). As seen in Fig. 13, the freshwater components of C T and A T are of equal size but opposite sign, and there is therefore no net effect of freshwater fluxes on the pH change (see Sarmiento and Gruber (2006) for a theoretical explanation). Also the thermodynamic effect of increasing salinity on pH is negligible. This increasing salinity of the Nordic Seas is a result of changes in the inflowing Atlantic Water related to subpolar gyre strength (Holliday et al., 2008;Lauvset et al., 2018). The contribution of 545 the biogeochemical component of A T is generally negligible, except in the Barents Sea Opening where it explains the lack of a significant pH decline (Fig. 7). In our dataset, the effect of changes in temperature on pH in the upper layer is relatively small. In contrast to several studies pointing towards a warming of the Nordic Seas (e.g. Holliday et al., 2008;Blindheim and Østerhus, 2013;Lauvset et al., 2018;Ruiz-Barradas et al., 2018), the Barents Sea Opening, the eastern Fram Strait and the Iceland Sea show no significant change in temperature. This might be an artefact of unequal distribution of sampling over the    seasons. When calculating trends with all available temperature data, not only those accompanying the C T and A T data, we obtain a clear warming signal (not shown).
In deeper layers, there is an overall increase in C T , A T (except in the Iceland Sea), salinity, and temperature. Although the effect of increasing C T bg is reduced away from the surface as a consequence of the gradual isolation of deeper waters from the atmosphere, it remains the main driver of pH change down to 2000 m. The significant trends of C T bg at the 1000-2000 555 m depth level in the Greenland Sea could be a consequence of the deep winter mixing that has been shown to reach down to 1500 m in this region (Brakstad et al., 2019). In the other regions of the Nordic Seas the winter mixed layers have not been documented to reach these depths (Ólafsson, 2003;Skjelvan et al., 2014;Våge et al., 2015, e.g.). However, intermediate water masses of the Greenland Sea has been shown to spread horizontally in the Nordic Seas, which could also explain the significant trends in the Norwegian and Lofoten Basin and in the Icealand Sea (Blindheim, 1990;Blindheim and Rey, 2004;560 Messias et al., 2008;Jeansson et al., 2017). undersampling. One biogeochemical process that could have a potential impact the Barents Sea A T bg trend is the recurrent blooms of calcifying coccolithophorids (Giraudeau et al., 2016), which consumes A T during growth, and releases A T when their shells are decomposed. There are indications of an increase in their presence in the Barents Sea (Giraudeau et al., 2016;Oziel et al., 2020). In which direction this would impact the A T depends on horizontal advection, remineralization and burial, and deserves separate dedicated process studies. The freshwater components of C T and A T are mainly detectable in the upper 570 500 m. As for the surface, the thermodynamic effect of salinity changes on pH are neglibible in the deep water. The warming seen in deep waters, that has a negative contribution on the pH trend, is an additional indication of that the absence of a temperature trend in the upper layer is a result of seasonal undersampling. In deep waters, the warming signal do not only come from local vertical mixing. There is also an indication of decreased deep-water formation in the Greenland Sea, which has caused an increased exchange with warmer Arctic deep waters (e.g. Østerhus and Gammelsrød, 1999;Blindheim and Rey, 575 2004; Karstensen et al., 2005;Somavilla et al., 2013). Below 2000 m, there are barely any detectable changes in the various pH drivers. The water masses at these depths are increasingly dominated by old Arctic deep waters (e.g. Somavilla et al., 2013). With ages exceeding 200 years (Jutterström and Jeansson, 2008;Stöven et al., 2016) they have been isolated from the increasing anthropogenic CO 2 , which explains the weak trends at these depths.

Past and future drivers 580
For past and future changes, the drivers of surface pH change show similar spatial patterns over all time periods, except for temperature (Fig. 14). The main driver is an increase in C T , which is larger in Atlantic Water than in polar waters. This is explained by the dilution of C T in polar waters by the increased freshwater export from the Arctic Ocean (Fig. 15, Shu et al., 2018) that to some degree counteracts the effect of atmospheric CO 2 uptake. A similar freshwater effect has recently been observed also in the Arctic Ocean (Woosley and Millero, 2020). The biogeochemical component of the C T driver (Fig. 15), 585 which is primarily the effect of increasing anthropogenic carbon, is larger in polar waters for the changes from present to future in both the esmRCP2.6 and esmRCP8.5 scenarios, in agreement with is what is expected from their lower buffer capacity (Sect. 2). The effect of A T is most prominent in polar waters, where a reduced A T concentration contribute to a pH decrease that is of the same order of magnitude as that driven by C T (Fig. 14). From the freshwater decomposition in Fig. 15, we see that the A T changes are mainly driven by freshwater fluxes, and that contributions from the biogeochemical component are negligible.

590
A T dilution has also been shown to be important in the future in the Arctic ocean in several CMIP6-models (Terhaar et al., 2021). However, as discussed earlier, the net effect of these freshwater fluxes on pH are minor, as the dilution of A T and C T is similar, but have opposite effects on pH (compare Fig. 15d-f with 15j-k). The increasing freshwater export also results in a dilution of salinity in polar waters that has a positive contribution to the pH trend. The Atlantic Waters show a tendency towards increasing salinity that partly amplifies the decrease in pH. Temperature has an overall negative effect on the pH trend 595 as a result of an overall warming. From pre-industrial to present-day, and present-day to future esmRCP2.6, the temperature increase is almost non-existent in polar waters, indicating that it has been shielded from warming through the presence of sea ice. In some smaller regions there is even a sign of a cooling, which could be a result of an increased presence of polar waters due to the increasing freshwater export.
The combined effect of these drivers explain the zonal gradients in the pH decrease that are described in Sect. 5.3 and 5.5.

600
From past to present-day the largest pH decrease takes place in the Atlantic Water due to a stronger increase of anthropogenic carbon and a stronger warming in these waters. From present-day to future the acidification becomes larger in polar waters compared to Atlantic Water due to the stronger increase of anthropogenic carbon in these waters. The increasing freshwater export from the Arctic that is seen in all time periods is of importance when regarding C T and A T concentrations separately, but their combined effect on pH is negligible. For the changes from past to present-day and present-day to future esmRCP2.6 605 the zonal gradient in Ω Ar trend follows that of pH, showing the importance of the C T driver. It is reinforced by the spatially variations in the warming, i.e., the stronger warming in the Atlantic Water compared polar waters results in a relatively stronger drop in Ω Ar in polar waters. In the esmRCP8.5 future, Ω Ar , in contrast to pH, exhibit a larger drop in the Atlantic Water. This can be explained by the relatively small changes in temperature in this region compared to the rest of the Nordic Seas, which affect Ω Ar in the opposite direction compared to pH.

610
Below the surface layer, C T is also the main driver of past and future pH changes (Fig. 16). The change from pre-industrial to present-day indicates a gradually weaker impact of C T with depth, except for a tongue at about 1000 m depth that connects to the surface in the Iceland sea. This is most likely related to the deep water formation in this region that spreads at depth.
The end-of-the-century C T increase for the esmRCP2.6 scenario is larger in the deep than in the surface layer, resulting in the stronger pH reduction at mid-depths as seen in Fig. 10. This mid-depth layer with a strong acidification is partly a result 615 of the higher atmospheric CO 2 concentrations in the middle of the 21st century in combination with the rapid ventilation of the water column in this area, i.e., when these waters were at surface they were exposed to peak atmospheric CO 2 . However, the large C T increase in deep waters is also partly explained by increased remineralization, as indicated by a ∼1 ml O 2 l −1 increase in the apparent oxygen utilization (AOU) at depths of 1800-2100 m throughout the Nordic Seas in both esmRCP2.6 kg −1 , which results to a pH decrease of ∼0.1 at the alkalinity in question. Impacts of changes in A T , salinity and temperature, are relatively modest at depth.
The residual between the sum of the four drivers and the actual pH change is small (Figs. 14 and 16) and can be attributed to approximations involved in the decomposition, including the approximations of the partial derivatives, the assumption of a linear trend and the use of temporal means (Takahashi et al., 1993;Lenton et al., 2012;Lauvset et al., 2015). Although the 625 absolute numbers related to the drivers should be taken with care, this decomposition still gives a good estimate of the relative importance of temperature, salinity, C T , and A T on pH changes.
In the historical run and all three future projections of NorESM1-ME, the change in surface ocean pCO 2 differs from the change in the atmosphere (Fig. S16). From pre-industrial to present-day, there is an increase in the undersaturation, i.e., the increase in the oceanic pCO 2 lags behind the increase in the atmosphere. This means that the pH decrease is less than 630 that expected from the increase in atmospheric CO 2 . The lag continues into all the future scenarios, but from around 2040 and onward, the oceanic pCO 2 increases faster than that of the atmosphere, resulting in a decreasing undersaturation. In esmRCP2.6 and esmRCP4.5 this causes stronger decreases in pH (from 1996-2005 to 2090-2099) than expected from the rise in atmospheric CO 2 . In esmRCP8.5, however, the difference between the end-of-the century ocean and atmospheric pCO 2 is still larger than the present-day, meaning that the decrease in pH is less than expected. As detailed above there are several 635 mechanisms underlying undersaturation of surface ocean pCO 2 in the Nordic Seas, but further analyses of these, including their potential future changes, is beyond the scope of this paper.

Summary and Conclusions
We have provided a detailed analysis of spatial and temporal variations of past, present-day and future acidification, and its drivers, in the Nordic Seas. We have further assessed the potential impacts of this acidification on aragonite saturation and 640 cold-water coral reefs. This work builds on , who estimated pH trends, and their drivers, for various sub-regions of the Nordic Seas from observational data sampled between 1981 and 2013. Here we have added data from the Iceland Sea and from later years, to obtain the greatest possible temporal and spatial coverage. We have additionally made an analysis of past and future pH changes by the use of the gridded GLODAP climatolgy and ESM-simulations, to put the observed changes into the context of long-term climate change. In contrast to previous studies that have assessed the future pH 645 changes in the Nordic Seas for single scenarios (Bellerby et al., 2005;Skogen et al., 2014Skogen et al., , 2018, we here analyse output from one mitigation scenario, one stabilization scenario and one high-emission scenario. To our knowledge, no previous studies have presented past pH changes in the Nordic Seas.
pH changes and its potential ecosystem impacts From pre-industrial (1850-1860) to present days (1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005), a combination of NorESM1-ME with the GLODAPv2 pre-650 industrial estimate, suggests that the pH of Nordic Seas surface waters has dropped by 0.1. During this period, the aragonite saturation horizon has slightly shallowed, but has remained well below the depths of known cold-water coral habitats. During 1981-2019, when regular sampling of carbon system variables have been made in the region, the pH of the Nordic Seas upper layer has decreased at a rate of -2.79±0.3 10 −3 yr −1 on average, resulting in a pH decline of 0.11. The pH reductions are significant all over the Nordic Seas upper layer (0-200 m), except in the Barents Sea Opening where the lack of significant 655 change is a result of a strong increase in A T . In some regions the acidification is detectable down to 2000 m, which we attribute to the deep water formation and spreading of these water masses at depth. The waters at 1000-2000 m throughout the Nordic Seas are now close to aragonite undersaturation. Our results are in overall agreement with Skjelvan et al. (2014), but the longer timeseries result in statistically significant (p<0.05) trends in even more regions and depth layers. An additional pH drop of 0.1-0.4 in the surface waters is projected until the end of the 21st century, depending on the emission scenario. In the high-emission 660 scenario, esmRCP8.5, all cold-water coral reefs will be exposed to corrosive waters by the end of the 21st century, threatening not only their existence, but also that of their associated ecosystems. This is confirmed by an ensemble of 6 CMIP5 models, who all agree on these consequences. The NorESM1-ME simulations suggest that some cold-water corals will be exposed to undersaturation also under the esmRCP4.5 scenario, and that this only can be avoided by keeping the emissions within the limits prescribed in the esmRCP2.6 scenario. Because NorESM1-ME tends to simulate a relatively strong decline of pH and 665 shallow saturation horizons in comparison to our ESM-ensemble for esmRCP8.5, our estimated aragonite saturation horizons for esmRCP2.6 and esmRCP4.5 should be considered as the shallow, lower bound of possible future states. Our estimates of the future pH and Ω Ar in the Nordic Seas add more possible future states to the ones presented for the A1B and RCP4.5 scenarios by Skogen et al. (2014Skogen et al. ( , 2018.

pH drivers 670
The acidification during the last 39 years is, in all sub-regions, mainly driven by increasing C T in response to the rising anthropogenic carbon concentrations. This is in agreement with the results for the period of 1981-2013 from Skjelvan et al. (2014), who calculated the drivers of pH change for the Norwegian Basin and the Greenland Sea. The effects of increasing C T are slightly opposed by increasing A T . The increasing A T is partly a result of a "salinification" of the Nordic Seas. However, this salinification also results in a decrease in C T , which counteracts the effect of the freshwater-driven increase in A T . The 675 net effect of C T and A T on increasing pH is therefore a result of biogeochemical processes. We find a clear warming signal in deep waters, which has contributed to the decreasing pH. In the upper 200 m, however, there is no clear temperature change.
We find this to be a result of seasonal undersampling, which further complicates a comparison of the changes in sea surface pCO 2 to the atmospheric one. In the Barents Sea Opening, there is an exceptionally strong increase in A T , which we cannot relate to increasing salinity. The reasons behind this strong increase is then either a result of biogeochemical processes, or can 680 also be a result sampling issues. Unfortunately, we cannot pin this down with the dataset we have, and this remains as an open question for future investigations.
For past and future changes, we also find increasing C T to be the main driver of pH change in the Nordic Seas. This is in agreement with Skogen et al. (2014), but we distinguish some regional differences related to different water masses. Increasing temperatures, that amplifiy the effect of increasing C T , have the largest impact in Atlantic Water in changes from pre-industrial 685 to present-day and present-day to the future esmRCP2.6. The absence of a warming signal in polar waters is a result of the shielding effect of sea-ice. In esmRCP8.5, however, the warming is more uniform over the Nordic Seas, which most likely is a result of the significantly reduced sea ice cover. In both past and future scenarios, there is a clear signal of an increasing freshwater export from the Arctic Ocean that dilutes C T , A T , and salinity in polar waters, and there is a tendency to increasing salinity in the Atlantic Water, that also leads to increasing C T and A T . The total effect of this change in freshwater content on 690 pH is negligible as the effect of changing C T and A T oppose each other, and because the thermal effect of salinity is minor in comparison to the other drivers.

695
The data from Ocean Weather station M from 2001-2007 is available in GLODAPv2.2019. Data from the time period 2008-2019 will be available in the next GLODAP version.
The data from the time-series station in the Iceland Sea can be obtained from the NCEI database (Ólafsson, 2012; Ólafsdóttir et al., 2020) The data from the Norwegian ocean acidification monitoring program (

705
Author contributions. AO, FF and FF designed the research. FF, FF, and AO performed the data-analysis with inputs from NG, IS, MC and EJ. FF lead the writing of the manuscript with inputs from all co-authors. JT designed, tested, and performed the NorESM1-ME model simulations.