# Remote sensing the sea surface CO2 of the Baltic Sea using the SOMLO methodology

Gaëlle Parard, Anastase Alexandre Charantonis, A. Rutgerson

► To cite this version:

Gaëlle Parard, Anastase Alexandre Charantonis, A. Rutgerson. Remote sensing the sea surface CO2 of the Baltic Sea using the SOMLO methodology. Biogeosciences, 2015, 12 (11), pp.3369-3384. 10.5194/bg-12-3369-2015. hal-01207654

HAL Id: hal-01207654

https://hal.sorbonne-universite.fr/hal-01207654

Submitted on 1 Oct 2015

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Remote sensing the sea surface CO₂ of the Baltic Sea using the SOMLO methodology

G. Parard<fnref n="1" />, A. A. Charantonis<fnref n="2" />,<fnref n="3" />, and A. Rutgerson<fnref n="1" />

*Correspondence to:* G. Parard (gaelle.parard@geo.uu.se)

Received: 30 June 2014 – Published in Biogeosciences Discuss.: 12 August 2014

Revised: 16 March 2015 – Accepted: 7 May 2015 – Published: 4 June 2015

**Abstract.** Studies of coastal seas in Europe have noted the high variability of the CO₂ system. This high variability, generated by the complex mechanisms driving the CO₂ fluxes, complicates the accurate estimation of these mechanisms. This is particularly pronounced in the Baltic Sea, where the mechanisms driving the fluxes have not been characterized in as much detail as in the open oceans. In addition, the joint availability of in situ measurements of CO₂ and of sea-surface satellite data is limited in the area. In this paper, we used the SOMLO (self-organizing multiple linear output; Sasse et al., 2013) methodology, which combines two existing methods (i.e. self-organizing maps and multiple linear regression) to estimate the ocean surface partial pressure of CO₂ (pCO₂) in the Baltic Sea from the remotely sensed sea surface temperature, chlorophyll, coloured dissolved organic matter, net primary production, and mixed-layer depth. The outputs of this research have a horizontal resolution of 4 km and cover the 1998–2011 period. These outputs give a monthly map of the Baltic Sea at a very fine spatial resolution. The reconstructed pCO₂ values over the validation data set have a correlation of 0.93 with the in situ measurements and a root mean square error of 36 µatm. Removing any of the satellite parameters degraded this reconstructed CO₂ flux, so we chose to supply any missing data using statistical imputation. The pCO₂ maps produced using this method also provide a confidence level of the reconstruction at each grid point. The results obtained are encouraging given the sparsity of available data, and we expect to be able to produce even more accurate reconstructions in coming years, given the predicted acquisition of new data.

## 1 Introduction

The ocean plays an important role in the global carbon budget. It acts as a major carbon sink for anthropogenic carbon dioxide (CO₂) emitted to the atmosphere from fossil fuel burning, cement production, biomass burning, deforestation, and various land use changes. The ocean is currently slowing the rate of climate change, having absorbed approximately 30% of human emissions of CO₂ to the atmosphere since the Industrial Revolution (Stocker et al., 2013). The exchange of CO₂ between coastal environments and the atmosphere is a significant part of the global carbon budget (e.g. Borges et al., 2005; Chen and Borges, 2009; Laruelle et al., 2010). While continental shelves represent only 7% of the oceanic surface area and less than 0.5% of the ocean volume, the estimated overall sink of CO₂ in the continental shelf sea is -0.22 PgC yr⁻<fnref n="1" /> (Laruelle et al., 2010), corresponding to 16% of the open oceanic sink (Takahashi et al., 2009). These estimates are subject to great uncertainty related to sparse data coverage in time and space. Monitoring the oceanic partial pressure of CO₂ (pCO₂) at monthly and seasonal timescales is essential for estimating the regional and global air-sea CO₂ fluxes and reducing this uncertainty. For technical and budgetary reasons, in situ measurements of marine pCO₂ are sparsely distributed in time and space. However, over the last decade, technical improvements and cooperation with the shipping industry have allowed for the installation of several autonomous monitoring systems aboard commercial vessels routinely crossing the ocean basins. These instruments make quasi-continuous measurements, allowing

<footnote n="1">¹Department of Earth Sciences, Uppsala University, Uppsala, Sweden</footnote>

<footnote n="2">²Centre d'études et de recherche en informatique, Conservatoire des Arts et Métiers, Paris, France</footnote>

<footnote n="3">³Laboratoire d'océanographie et du climat: expérimentations et approches numériques, Université Pierre et Marie Curie, Paris, France</footnote>
regional analysis of the highly variable spatial and temporal distributions of pCO₂ (e.g. Lefèvre et al., 2004; Lüger et al., 2004; Corbière et al., 2007; Schneider et al., 2003). The Baltic Sea, a semi-enclosed sea in northern Europe, is relatively well monitored and has been studied for several decades (Meier et al., 2014). Despite the increased number of measurements made in the Baltic Sea, assessing the carbon fluxes in the Baltic Sea remains particularly challenging due to the non-linearity of the emission and absorption system. This non-linearity is complicated by a combination of varying salinity, varying river input of dissolved organic carbon, and large general variability due to the strong seasonal cycle in the region. Using new methodologies could generate additional information from the relatively limited number of existing measurement data. Neural network techniques are empirical statistical tools that somewhat resolve the non-linear and often discontinuous relationships among proxy parameters without any a priori assumptions. In the past decade, several authors have reported the application of a neural network technique to basin-scale pCO₂ sea analysis (Lefèvre et al., 2005; Jamet et al., 2007; Friedrich and Oschlies, 2009; Telszewski et al., 2009; Landschützer et al., 2013; Nakaoka et al., 2013; Schuster et al., 2013), concentrating mainly on the North Atlantic Ocean. Most recently, Telszewski et al. (2009) successfully applied a neural network technique based on a self-organizing map (SOM) to reconstruct the seawater pCO₂ distribution in the North Atlantic (10.5 to 75.5° N, 9.5° E to 75.5° W) for three years (i.e. 2004 to 2006) by examining the non-linear/discontinuous relationship between pCO₂ sea and the ocean parameters of sea surface temperature (SST), mixed-layer depth (MLD), and chlorophyll a concentration (Chl a). In this paper, we applied the SOMLO methodology (self-organizing multiple linear output), which creates an SOM classification of the available explicative oceanic parameters in the Baltic Sea and then calculates multiple linear regression (MLR) parameters for estimating the pCO₂ from the elements belonging to each class separately. The major benefit of this methodology is that it allows the use of a linear model (i.e. the MLR) despite the non-linear relationship between pCO₂ and its explicative parameters, by using the SOM classifier. The SOM classifier can determine the region of the multidimensional data space in which to perform the linear regression. If the classification is fine enough, each region will represent a single type of relationship between the pCO₂ and the explicative data, a type that can be reproduced by an MLR. Due to the temporal and spatial limitations of the in situ pCO₂ data, the satellite data can help estimate the air-sea CO₂ fluxes over the entire Baltic Sea. The satellite data have a higher spatial coverage of the Baltic Sea, allowing estimation of the pCO₂ from in situ data using the SOMLO method. Chierici et al. (2009) demonstrate that, in the North Atlantic Ocean, SST, Chl a, and MLD contributed significantly to the estimation of pCO₂ from a linear relationship. Based on this idea, we applied these parameters to the Baltic Sea, adding two other

parameters, i.e. net primary production (NPP) and coloured dissolved organic matter (CDOM), that provide information about the biological activity occurring in summer. From this, we develop pCO₂ algorithms applicable to the Baltic Sea using in situ pCO₂ values; remotely sensed SST, Chl a, and CDOM; modelled MLD and NPP; and time.

The manuscript is structured in four parts. First, we present a synopsis of the problem studied, including existing studies of pCO₂ reconstruction in other maritime regions. Next we present the available data and briefly describe the methodology used. In the third part of the article we present our results, namely, the topological maps obtained and the reconstructions performed with them. We conclude the article by discussing the results obtained and future possible improvements of the method used.

# 2 Materials and methods

## 2.1 Study area

The Baltic Sea is a semi-enclosed sea with limited exchange with the North Atlantic through the North Sea–Skagerrak system. Previous investigations of the Baltic Proper found large temporal and spatial variability of pCO₂. The amplitude of the annual pCO₂ cycle varies significantly depending on the region, ranging from 400 µatm in the northeastern Baltic Proper to 120 µatm in the transition areas to the North Sea (Schneider and Kaitala, 2006). The Baltic Sea receives significant river run-off from surrounding land (a total of approximately 15 000 m³ s⁻¹ (Bergstrom, 1994) and net precipitation of approximately 1500 m³ s⁻¹ (Omstedt et al., 2004)). This large freshwater addition brings large amounts of nutrients and inorganic and organic carbon to the Baltic Sea basin (Omstedt et al., 2004; Hjalmarsson et al., 2008). The biogeochemical processes in the Baltic Sea marine environment are controlled mainly by the biological production and decomposition of organic matter occurring in the context of the region’s hydrography (Siegel and Gerth, 2012). Physical forcing controls the water transport, stratification, temperature, and salinity in the Baltic Sea; these factors then influence the nutrient and carbon distribution, thereby affecting biogeochemical processes. We divide the Baltic Sea into three basins, i.e. the central part (CP), the Gulf of Bothnia (GB), and Gulf of Finland (GF), as shown in Fig. 1. The Baltic Sea has an average depth of 55 m and a maximum depth of 460 m at the Landsort Deep (Wesslander, 2011).

## 2.2 pCO₂ observations

To compile the pCO₂ maps, we use measured data from three sources.

1. The Östergarnsholm site: this site is located next to the small island of Östergarnsholm in the central Baltic Sea and is further described by Rutgersson et al. (2008) and
<figure><img src="image_3.png" /><figcaption>**Figure 1.** The dashed red lines divide the Baltic Sea into three basins: CP: central part; BB: Gulf of Bothnia; and GF: Gulf of Finland. Monthly data are available from 1998 to 2011 for the Baltic Sea. The colour bar shows the $pCO_2$ values in µatm. The black arrow shows the position of the SAMI sensor.</figcaption></figure>

Norman et al. (2013a). The island is situated 4 km from the east coast of the larger island of Gotland. SST and $pCO_2$ are measured semi-continuously 4 m below the sea surface using a submersible autonomous moored instrument (SAMI) $CO_2$ sensor moored at a buoy 1 km southeast of the tower situated on the island. In addition, SST is also measured using a wave rider buoy (operated by the Finnish Meteorological Institute) at 0.5 m depth situated approximately 4 km southeast of the tower.

2. Cargo ship: this data set derives from continuous measurements of the surface water $pCO_2$ made in the Baltic Sea using a fully automated measurement system deployed on a cargo ship. The Leibniz Institute for Baltic Sea Research, Warnemünde, Germany (IOW – Institut für Ostseeforschung Warnemünde), has made continuous measurements of $pCO_2$ at 5 m depth aboard the cargo vessel *Finnpartner*. This ship crosses between Lübeck and Helsinki at a 2-day interval, alternately crossing the eastern and western Gotland Sea (Schneider and Kaitala, 2006; Schneider et al., 2009). Data from *Finnpartner* were acquired between July 2003 and December 2005.

3. The Swedish Meteorological and Hydrological Institute (SMHI) database Svenskt Havsarkiv (SHARK): pH (measured using the method of Grasshoff et al., 1999) and total alkalinity (TA) (measured using potentiometric titration as described by Grasshoff et al., 1999) are

measured continuously at a monthly or semi-monthly resolution in the Baltic Sea at various stations. All measurements are made at a depth of 5 m depth. The uncertainty of the pH is ±0.03 pH units and of the TA is ±5% (Wesslander et al., 2009). $pCO_2$ is estimated from the pH, TA, salinity, and temperature measurements using the standard CO2SYS program (Lewis and Wallace, 1998) with the equilibrium constant from Weiss (1974) and Merbach et al. (1973) as refitted by Dickson and Millero (1987), as in Wesslander et al. (2009). The $pCO_2$ estimation was compared with the other source of $pCO_2$ data, all the data were compared in time (less than 1 h) and position (around less than 0.2°), the correlation coefficient ($R$) gives 0.98, and the standard deviation (SD) is 9 µatm. The high standard deviation (SD) is explained by the presence of upwelling with high $pCO_2$ present near the coast in the SAMI sensor measurement; when we remove the upwelling event, the SD is 1.5 µatm.

## 2.3 Remote sensing data

The satellite data used in this study are from various sources. We use a monthly temporal resolution and a spatial resolution based on the lowest spatial resolution of our data sets, i.e. based on the lower spatial resolution CDOM data set.

We obtain values for five parameters from various sources:

SST: several data sets are used for SST, and we combine two types of data sets for 2007 and 2011. For 2005–2011, we use data from the Federal Maritime and Hydrographic Agency (BSH), which processed data from the NOAA Advanced Very High Resolution Radiometer (AVHRR), and data from the Group for High-Resolution Sea Surface Temperature (GRHSST) data set for the Baltic Sea, 2007–2011. The spatial resolution is 0.03° at a daily temporal resolution (http://podaac.jpl.nasa.gov/dataset/DMI-L4UHfnd-NSEABALTIC-DMI_OI). For 1998–2004, the data come from a reanalysis of the NOAA/NASA AVHRR data stream conducted by the University of Miami’s Rosenstiel School of Marine and Atmospheric Science (RSMAS) and the NOAA National Oceanographic Data Center (NODC). This data set consists of the monthly average SST (in °C) over the zone, with a spatial resolution of 4 km, extracted from version 5.2 of the AVHRR Pathfinder project (Casey et al., 2010) (http://www.nodc.noaa.gov/SatelliteData/pathfinder4km/). The various SST data sets were compared with measured SMHI temperature data at a monthly resolution from 1998 to 2011, giving a correlation coefficient ($R$) of 0.99 and a mean difference (MD) of 0.05 °C between the two data sets. The difference observed between the measured and satellite data between 1998 and 2004 give a value ($R$ = 0.99 and MD = 0.09 °C) near the difference between 2005 and 2011.
and 2011 (R = 0.99 and MD = 0.14). The two SST data sets used between 2007 and 2011 were also compared with the SAMI sensor data at a daily resolution, giving a good correlation for the BSH data (R = 0.95 and MD = 0.06 °C) and the GRHSST data (R = 0.95 and MD = 0.08 °C).

Chl *a*: this data set consists of monthly averages from the following sensors: Sea-viewing Wide Field-of-view Sensor (SeaWiFS) (September 1998 to December 2002) with 4 km spatial and monthly temporal resolutions and Moderate-Resolution Imaging Spectroradiometer (MODIS-Aqua) (July 2002 to June 2011) with 4 km spatial and monthly temporal resolutions (Casey et al., 2010). A log-normal distribution was assumed for the Chl *a* data. Comparison with SMHI measurement data and in situ data (personal communication from Dr. Tiit Kutser) gives R = 0.67 and MD = 7 mg m⁻³. The chlorophyll levels from the satellite data seem to be overestimated compared with the in situ data. Since we use the same data set over the whole study period, this bias was learned during the classification process and the MLR parameters are calculated considering that any further input will include that bias.

CDOM values come from MODIS-Aqua 4 km monthly average data. The CDOM index quantifies the deviation in the relationship between the CDOM and Chl *a* concentrations, where 1.0 represents the mean relationship for Morel and Gentili (2009) case 1 waters, and values above or below 1.0 indicate an excess or deficit, respectively, in CDOM relative to the mean relationship. The algorithm and its application are fully described by Morel and Gentili (2009). In situ CDOM data (T. Kutser, personal communication, 2014) give a lower correlation coefficient and a low average difference (R = 0.48, MD = 2.3). As for the chlorophyll, the bias is applied to all years so it does not affect the estimation of pCO₂.

NPP values come from two data sources. The first data set comes from the Environmental Marine Information System (EMIS): the EMIS model is depth-integrated but allows for depth-dependent variability in the diffuse attenuation coefficient, which is calculated from a multiple-component semi-analytical inversion algorithm (Lee et al., 2005). The primary production calculation is based on the formulation obtained through dimensional analysis by Platt and Sathyendranath (1993). The photosynthetic parameters are assigned by the combined use of a temperature-dependent relationship for the maximum growth rate Eppley (1972) and a variable formulation to retrieve the C: Chl *a* ratio following the empirical relationship of Cloern et al. (1995). The EMIS data set comprises monthly average values from October 1997 to September 2008. The second data

set, for 2009–2011, uses the Vertically Generalized Production Model (VGPM) of Behrenfeld and Boss (2006) as the standard algorithm. VGPM is a “chlorophyll-based” model that estimates net primary production from chlorophyll using a temperature-dependent description of chlorophyll-specific photosynthetic efficiency. For VGPM, net primary production is a function of chlorophyll, available light, and photosynthetic efficiency. VGPM uses MODIS-Aqua chlorophyll and temperature data, SeaWiFS photosynthetically active radiation (PAR) data, and estimates of the euphotic zone depth from a model developed by Morel and Berthon (1989) and based on chlorophyll concentrations. For the NPP for 2009–2011, the observed maximum value was limited to 10 to be comparable to the data for 1998–2008. Validation of NPP was difficult due to the number of data available in the area. Comparison of the two data sets gives similar values and seasonal cycles. We compared the seasonal cycles between 1998 and 2008 and between 2009 and 2011, obtaining values of the same order of magnitude.

MLD: there are also two sources of the MLD data. Monthly averages from 1998 to 2007 come from a 3-D hydrodynamic model currently used at the Joint Research Centre/Institute for Environment and Sustainability (JRC/IES), i.e. the public-domain General Estuarine Transport Model (GETM; www.getm.eu), which has its roots in developments at the JRC/IES (Burchard and Bolding, 2002). GETM simulates the most important hydrodynamic and thermodynamic processes in coastal and marine waters and includes flexible vertical and horizontal coordinate systems. Different turbulence schemes are incorporated from the General Ocean Turbulence Model (GOTM; www.gotm.net). For 2008–2011, we use data from the carbon-based production model at a monthly resolution (Behrenfeld et al., 2005). The MLD was estimated from SMHI temperature and salinity profile measurements using the density criterion of Boyer Montégut et al. (2004); the comparison between the model estimation and estimated SMHI MLD is good (R = 0.63 and MD = 17 m). Between 1998 and 2007, the correlation coefficient is higher (R = 0.8) than between 2008 and 2011 (R = 0.5). Nevertheless, the MD is lower for the second data source, i.e. 20 m vs. 11 m. This can be explained by the available data coverage: 63% of the in situ data cover the 1998–2004 period. From 2008 to 2011, the maximal value is below 80 m between 1998 and 2007, so for the in situ data estimation, we replace every value above 80 m with “not a number”.

In the Baltic Sea, there are many gaps in the satellite data; this is due to the high proportion of coastal waters, where satellite data are less reliable, and to the frequent large-scale cloud coverage. To increase the total number of our data,
we used a monthly temporal resolution and a method to improve the spatial distribution of the data. For statistical analysis, the irregularly spaced density of the measurements were first uniformly resampled. To this end, Gaussian grinding was used, as described by Greengard and Lee (2004); Dutt and Rokhlin (1995). The data points of the original series are convolved using a Gaussian kernel. As a result, the data points are smeared over their neighbouring equi-spaced points, which are more densely distributed. This method produces more realistic values than does simple interpolation, particularly when there are many data gaps (Schomberg and Timmer, 1995). There is no discontinuity between the different data sets, but NPP, CDOM, and Chl a data are missing for January and December, so no reconstruction can be performed for these months.

## 2.4 Data available

The positions and values of all the in situ pCO₂ data are shown in Fig. 1. We use the spatial resolution of the parameter with the lowest resolution for the final data set chosen (i.e. CDOM). A monthly temporal resolution is used for this study. Rutgersson et al. (2008) demonstrate that the agreement between SAMI sensor data and the ship data from Finnpartner (near the mooring maximum of 23 km) is quite good. This good agreement is confirmed by the comparison between pCO₂ data from the SAMI sensor, the data surrounding the mooring (0.2°), and the other data sets. These analyses give a good correlation factor of 0.98. The in situ data are available mainly for the central basin, but the number of data for the Gulf of Bothnia is very low, coming from two SMHI stations. The in situ pCO₂ data are well distributed over the 12 months (Fig. 2). January is the month for which the number of data is lower (i.e. below 80), but the other months have 110–155 data points each. In our case, each in situ data point is characterized by SST, Chl a, CDOM, NPP, and MLD as well as information on the date the measurements were made. This temporal information was normalized by sine and cosine, as follows:

$$ time(cosine) = \cos \left( \frac{\text{Days} \cdot 2\pi}{365} \right), \qquad (1) $$

$$ time(sine) = \sin \left( \frac{\text{Days} \cdot 2\pi}{365} \right), \qquad (2) $$

where “Days” represents the Julian day.

This definition of time is used to render the values continuous over the course of the year, sidestepping the artificial numerical transition from the last day of one year to the first day of the next, to be able to situate the process in relation to its seasonality.

Although time itself is not affecting the pCO₂, the inclusion of time (time(cosine) and time(sine)) as a parameter is important since, in the database, some situations have similar values for SST, Chl a, CDOM, NPP, and MLD but

<figure><img src="image_4.png" /><figcaption>Figure 2. Histogram showing the number of observations in function of the month between 1998 and 2011.</figcaption></figure>

different pCO₂ values. For example, in May we had a situation (SST = 9.3°C, Chl a = 0.1 mg m⁻³, CDOM = 4.3, NNP = 1.7 mg C m⁻² d⁻¹, and MLD = 10.8 m) whose parameters were very close to those of a situation in November (SST = 9.1°C, Chl a = 0.1 mg m⁻³, CDOM = 4.3, NNP = 1.6 g C m⁻² d⁻¹, and MLD = 10.3 m); however, the pCO₂ values differed greatly, being pCO₂ = 214 and 444 µatm respectively, during the two situations. This dissimilarity in the pCO₂ values informs us that these situations are generated by other drivers than SST, Chl a, CDOM, NPP, and MLD. Since these drivers can be related to seasonal patterns, we included information on the time of the year, as a proxy that allows us to fine-tune our classification. The inclusion of two temporal values instead of only one, even though those are highly correlated, was unavoidable in order to preserve continuity in the values obtained when changing years.

A principal component analysis (PCA) was conducted to highlight the importance of the parameters in the pCO₂ variability (Fig. 3). PCA is a method of analysis which involves finding the linear combination of a set of variables that has maximum variance and removing its effect, repeating this successively (Jolliffe, 2002). The percent of variance explained by each axis of the PCA is shown in Table 1. The results of this PCA indicate that the percentages of variance explained by all axes beyond the first do not differ greatly, with the notable exception of the last two axes, indicating that most parameters can be discriminant in the definition of the SOM states. The need to maintain the totality of these parameters is further demonstrated by the projection of the explicative parameters in the correlation circle, where all values are close to the boundaries of the correlation circle. This informs us that they are all tied to the phenomena explaining 62% of the total variation of our data. While we could reduce the dimension of the problem by using only the projection of these variables on the first few axes, we chose to maintain all the explicative parameters presented when applying SOMLO to estimate the pCO₂ in order to be able to not lose any information when performing the SOM classifications.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Correlation scatter plot with the representation of parameters on axes 1 and 2. The first component plane contains 62% of the total variance of the system under study, and all parameters are close to the exterior of the circle, indicating that they are all important in order to invert the pCO₂.</figcaption></figure>

<table><caption>**Table 1.** Percent of variance explained by each axis of the PCA.</caption><thead><tr><th>Axis number</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th></tr></thead><tbody><tr><th>Percent of variance</th><td>42</td><td>20</td><td>14</td><td>9</td><td>8</td><td>4</td><td>2.5</td><td>0.5</td></tr></tbody></table>

In total, 1445 pCO₂ data are used in this study, after having removed the outliers from our data set. These 170 outliers were beyond 3 standard deviations away from at least one of the explicative parameters. All parameters (i.e. SST, Chl *a*, CDOM, NPP, and MLD) are located around each pCO₂ datum. In winter (i.e. October to March), more data are missing (Table 2, column 1), particularly for Chl *a*, CDOM, and NPP, winter being the period when it is more difficult to measure or estimate these parameters. Between April and September, the number of missing SST, Chl *a*, CDOM, and MLD data is relatively low compared with the total number of data (Table 2, column 2), i.e. fewer than 3% of the total. To increase the number of data available, we completed the data by training the topological map. Further details on this are presented in Sect. 2.5.

## 2.5 Methodology

The relationship between pCO₂ and the environmental parameters is highly non-linear: a slight variation in some of the environmental parameters could correspond to significant variations in pCO₂. We chose to use the SOMLO methodology, which combines two statistical approaches: SOMs (Kohonen, 1990) and linear regression. SOMs are a subfamily of neural network algorithms used to perform multidimensional classification. A defining characteristic of SOMs is that their

<table><caption>**Table 2.** Number of missing values for each parameter for the satellite data for the October–March and April–September periods. The numbers in parentheses indicate the total data points in each period.</caption><thead><tr><th>Parameter</th><th>October–March (685)</th><th>April–September (814)</th></tr></thead><tbody><tr><td>SST</td><td>28</td><td>0</td></tr><tr><td>Chl <i>a</i></td><td>202</td><td>24</td></tr><tr><td>CDOM</td><td>320</td><td>5</td></tr><tr><td>NPP</td><td>468</td><td>571</td></tr><tr><td>MLD</td><td>6</td><td>2</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** The various elements used in training a self-organizing map. On the left we have the data set used to train the SOM, which discretizes it into classes. Each class contains a referent vector containing, for each parameter, the average value of the elements comprising it as well as an index of the class that informs us of its location on the topological map.</figcaption></figure>

classes can represent a Gaussian distribution centred around the typical profile of environmental parameters, if there is high discretization of the training data set (Dreyfus, 2005). We use this hypothesis to classify the environmental parameter data set and then estimate the parameters of a linear regression for each class. In the following section, we present a brief overview of the two statistical algorithms and their application to our data sets.

### 2.5.1 Self-organizing maps

Self-organizing topological maps comprise a clustering method based on neural networks. They cluster a learning data set into a reduced number of subsets, called classes, with common statistical characteristics.

Generating a SOM requires the creation of a training database that contains homogenous vectors. After a training phase, we obtain a SOM. The term “map” corresponds to a 2-D matrix that stores, for each class, its referent vector, $r_i$,
<figure><img src="image_7.png" /><figcaption>**Figure 5.** On top we have a representation of the topological map as a lattice, while on the bottom we have a projection of data in the data space (black circles), as well as the average values of each class (red circles). Adjacent classes on the SOM lattice correspond to adjacent areas in the multidimensional data space. The red lines indicate the connections between neighbouring classes.</figcaption></figure>

which approximates the mean value of the elements belonging to it, and its index, which positions it in the map matrix used to situate it in relation to the other classes (Fig. 4).

SOMs are also called self-organizing topological maps, “topological” indicating that the SOM training algorithm forces a topological ordering of the classes in the matrix, meaning that any two neighbouring classes $C_i$ and $C_j$ on the map matrix have referent vectors $r_i$ and $r_j$ that are close in the Euclidean sense in the data space.

Let us consider a vector **x** that is of the same dimensions and nature as the data used to generate the topological map; we can find the index of the class to which it is classified by choosing index = arg max<sub>i</sub> (||**x** − $r_i$||) (where arg max is the argument maximal), therefore assigning it to the class whose referent is closest to it in the Euclidean sense (Fig. 5). A classified vector **x** will be represented by its class index, $C_{index}$. If we are trying to classify a vector that has some missing values, the comparison is performed between the existing values of **x** and the corresponding values of each $r_i$.

As a version of the expectation-maximization algorithm, the SOM algorithm performs an iterative training. During the early phases of this training, the referent vectors of each class are strongly affected by the changes imparted on their neighbours’ referent vectors in order to capture the shape of the data cloud. Depending on the training parameters of the SOM, in the latter phases of the training, the effect of the neighbouring vectors on the determination of the referent vector can be considered null. In these cases, each referent vector approximates, locally, the mean value of the multidimensional Gaussian random distribution that generated the training data assigned to that class (Dreyfus, 2005).

## 2.5.2 Multiple linear regression

MLR is a modelling method that expresses the value of one response variable, $V$ (in our study $pCO_2$), as a linear function of other explicative variables, i.e. $X = X_1, X_2, ..., X_i$ (in our study SST, Chl $a$, CDOM, NPP, MLD, time$_{sin}$, and time$_{cos}$). An MLR is generally performed either to interpret the relationship between the variable $y$ and each of the other predictive variables $X_i$ or else to predict, from a data set of vectors containing the values of $X$, the corresponding value of $y$. In this paper, we used both aspects of MLR.

However, to perform MLRs in the present case, we had to take into account their limitations and the nature of the problem. Specifically, to perform an MLR we are obliged to assume that the relationship between the predictor variables and the response variable is linear. However, this is not the case in our data sets: $pCO_2$ is not linearly related to the variables presented when considering the entirety of the problem. However, as noted above in Sect. 2.5.1, if we consider the classes created by the SOM, they are very localized regions of the combined explicative and response data space that can be considered to approximate, locally, the mean value of a multidimensional Gaussian random distribution. We therefore assume that, if performed in the reduced neighbourhood of a SOM class, the relationships between $pCO_2$ and the explicative variables are linear.

# 3 Application and results

## 3.1 Statistical imputation

As described in Sect. 2.4, both the satellite and measured data available for the application present missing values. To complete these data sets, we chose to use imputation methods similar to those described by Schafer and Graham (2002) and Malek et al. (2008). The main idea of these methods is to use the classifying abilities of the SOMs to regroup the data in typical situations and replace the missing explicative data values with the corresponding values of the referent vector of the class to which it belongs.

We first selected the database containing SST, Chl $a$, CDOM, NPP, MLD, time$_{sine}$, and time$_{cosine}$. The vectors were sorted according to the number of values missing from each vector and noting the locations of these missing values. We chose all complete data vectors and the first 5% of the sorted vectors containing missing data, and we trained a SOM. We proceeded by replacing the missing values of these first 5% of vectors with the corresponding values of the referent vector of the class to which they each belong. We then included the next 5% of vectors with missing data in a new training data set and created a new SOM. Based on this new SOM, we again filled the new missing values with the corresponding values of the referent vector of the class of the new SOM to which they each belong. In addition, we deleted the
<figure><img src="image_8.png" /><figcaption>**Figure 6.** A schematic of the imputation method used. We initially sort the data depending on the amount of missing values present, then progressively train SOMs on the data set, steadily including more vectors for the training and completing and updating the missing data during each iteration.</figcaption></figure>

values we added to the first 5% of vectors and replaced them with the values of the referent vector of the class of the new SOM to which they each belong. The training parameters of this method, such as the number of classes of the topological map and the number of iterations, were selected by parameter tuning. We then continued iterating this process, updating the previously filled missing values with the values of their corresponding referents belonging to the most recently trained SOM, until all missing values were filled. The updating of the previously filled values allows the method to progressively incorporate information and to rectify the previous errors. A schematic of the imputation method used can be seen in Fig. 6.

After this imputation of the missing data through iterative training, the reconstructed data represent the original data well. Figures 7 and 8 show data for six variables before and after the reconstruction respectively. The main difference is observed for the values of $Chl\ a$, where the peak of over 200 individuals occurs at $0\ mg\ m^{-3}$ because, at the initialization of the imputation process, we decided to replace all null values. The repartitioning of $pCO_2$ (Fig. 8a) is very representative of the data variability with a large range of values. Some very high values occur during local events, such as coastal upwelling. Most of the data range between 180 µatm (value observed in summer) and 550 µatm (observed in winter). The SST (Fig. 8b) is very representative of the variability in the Baltic Sea, with a maximum occurring between July and September in all basins around 18 °C (Siegel and Gerth, 2012). The NPP variability is fairly homogenous, except for the peak at $10\ mg\ C\ m^{-2}\ d^{-1}$. This peak occurs because the first model providing us the NPP values has a set maximum of NPP at $10\ mg\ C\ m^{-2}\ d^{-1}$; therefore, the correction of the NPP satellite data takes this maximum into account.

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Histogram of (a) $pCO_2$ and satellite data, (b) SST, (c) $Chl\ a$, (d) NPP, (e) CDOM, and (f) MLD available for the SOM before reconstruction. Y axis represents the number of data, and the x axis the value of the parameters.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Histogram of (a) $pCO_2$ and satellite data, (b) SST, (c) $Chl\ a$, (d) NPP, (e) CDOM, and (f) MLD available for the SOM after reconstruction. Y axis represents the number of data, and the x axis the value of the parameters.</figcaption></figure>

The variability of chlorophyll results in one subset of the data with a low value and another subset with a value higher than $6\ mg\ m^{-3}$, which can be explained by the fact that the Baltic Sea is a narrow sea, with important coastal regions, and that two blooms take place, in spring and in summer. The chlorophyll value can be very high during these periods, and the reconstruction gives a mean value for this characteristic. A peak at $10\ mg\ m^{-3}$ is observed in the chlorophyll data, not due to the reconstruction but to the maximum value in the satellite data file. The low MLD occurs in summer, and in the model the minimum is 10 m deep, which appears to be around the minimum value observed in Fig. 8f. Absorption by CDOM decreases with increased distance from the riverine sources, reaching a relatively stable absorption background in the open sea. Most of our CDOM data capture open-sea conditions, so the values are quite low.
## 3.2 pCO₂ estimation

### 3.2.1 Topological map

We classified the explicative variables (i.e. SST, Chl *a*, CDOM, NPP, MLD, time<sub>sin</sub>, and time<sub>cos</sub>) into classes that share similar characteristics. In order to optimize the SOM map size for the method and to calculate the method's performances, we randomly sampled 90% of our completed data set (1300 vectors) to be used for the training phase, keeping 10% (144 vectors) to be used for computing the performances of the method. We iterated this process, selecting many different random samplings for each map size, and selected the map size with the best average reconstruction when applying the SOMLO methodology.

At the end of our optimization, we selected a SOM consisting of 77 classes. The number of observations captured by each class ranges from 0 to 38 (Fig. 9). The order of magnitude of the number of observations is constant throughout the SOM, and we can regard the classes as having spread in multidimensional space in order to accurately represent the data space of the explanatory parameters. The presence of classes that did not capture any elements can be justified as preventive: they preserve the topological aspect of the SOM by preventing classes that are not similar enough from becoming neighbours.

To estimate the average concentration of pCO₂ in each class, the measurements of pCO₂ associated with vectors consisting of SST, CDOM, NPP, MLD, and CHL *a* components were presented to the already-trained SOM as input data (Fig. 10). The average value computed for the vectors belonging to each class corresponds to the average value of pCO₂ for that class.

In the final map, the distribution of pCO₂ is strongly dependent on the SST distribution, with low values of pCO₂ correlating with high values of SST (Fig. 10). This is in agreement with the seasonal pCO₂ cycle, which is characterized by a large amplitude, ranging from a high value in winter (≈ 500 µatm) to a low value in summer (≈ 150 µatm), described, for example, by Wesslander (2011). According to Schneider and Kaitala (2006), the high winter value of pCO₂ is a consequence of mixing with a deeper water layer enriched in CO₂, which is in agreement with the distribution of the MLD (Fig. 10h), with the higher value in winter and autumn correlating with the high value of pCO₂. High values can also be explained by the mineralization, which exceeds production in winter (Wesslander, 2011). Biological production starts in spring when sunlight and nutrients are sufficient. The chlorophyll begins to increase in March–April due to the spring phytoplankton bloom, which reduces the pCO₂ level during this period. The more intensive decrease occurs in April and May, which is consistent with the higher value of NPP (Fig. 10). Studies in the central Baltic Sea identify two summer minima, the first in April–May and the second in July–August, resulting from a second production period.

<figure><img src="image_11.png" /><figcaption>**Figure 9.** (a) Distribution of the neural map of (a) the average difference in pCO₂ for each neuron represented by the colour bar. (b) The standard deviation (SD) in pCO₂ for each neuron represented by the colour bar. For (a) and (b), the numbers inside the neurons correspond to the number of data that construct the characteristics of the neuron.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Distribution of each parameter in the neural map: (a) pCO₂ in µatm, (b) SST in °C, (c) and (e) time cosine and sine respectively, (d) Chl *a* in mg m⁻³, (g) NPP in mg C m⁻², and (h) MLD in m.</figcaption></figure>

Higher variability is observed during this period, with a standard deviation between 39 and 50 µatm for different regions (Wesslander, 2011; Schneider and Kaitala, 2006).

### 3.2.2 Linear regression in the neurons

To perform an MLR, we must assume that the relationship between the predictor variables and the response variable is linear. We could take this to be a valid hypothesis only when performing the MLR in the reduced neighbourhood of a SOM class, where the relationship between pCO₂ and the explicative variables can be assumed to be linear.

For each class *j* a separate training data set was created containing all the vectors assigned to that class and to all its adjacent classes. Based on that data set, we computed the linear regression coefficient parameters for every explicative parameter and for a constant value.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Distribution of the neural map of the coefficient from linear regression for each parameter: (a) SST in °C, (b) and (d) time cosine and sine respectively, (c) Chl *a* in mg m⁻³, (f) NPP in mg C m⁻², and (g) MLD in m.</figcaption></figure>

The calculated linear regression coefficient parameter values for each class are shown in Fig. 11. Note that all parameters are important in specific regions of the SOM, having both positive and negative correlations in different classes.

More importantly, the fact that each parameter has a significantly varying range of values over the different classes demonstrates that each parameter is important in reconstructing the *pCO₂* in the Baltic Sea, even though a parameter may be highly significant in some classes and relatively stable in other regions of the topological map.

The addition of vectors belonging to adjacent classes did not generally perturb the estimation of the coefficient parameters because, as seen in Fig. 10, the values of all parameters are generally organized coherently on the map. The assumption that they are close in the data space is not as robust as it would have been had we solely considered the vectors belonging to each class, but, given the limited number of data available for modelling this highly non-linear and complex system, we would not have sufficient elements to correctly estimate the linear regression coefficients. Given the projected increase in available data in coming years, further applications of this approach will limit themselves to the elements belonging to each class.

### 3.2.3 Validation of reconstruction

To validate our results, we calculated the difference and SD between the value of *pCO₂* reconstructed in each neuron and the observation defining that neuron (Fig. 9). On average, the SD is approximately 38 µatm and the difference observed is 25–30 µatm. By cross-validating a data set (divided by 10 sequences), we obtained a mean SD of 48 µatm with a variation of 32–57 µatm and *R* equal to 0.9 with a variation of 0.86–0.96. Nevertheless some points indicating higher values can be identified (shown in red in Fig. 9). These values are explained by the positions of these points, which are at the edges of the cloud and therefore more likely to include

<table><caption>**Table 3.** Maximum and mean value observed in the difference between the data used for the trainee and the value in the neurons.</caption><thead><tr><td>Parameter</td><td>Maximum</td><td>Average</td></tr></thead><tbody><tr><td>pCO₂ (µatm)</td><td>56,4</td><td>29.15</td></tr><tr><td>SST (°C)</td><td>1.9</td><td>0.98</td></tr><tr><td>time(cos)</td><td>0.33</td><td>0.07</td></tr><tr><td>Chl *a* (mg m⁻³)</td><td>0.14</td><td>0.06</td></tr><tr><td>time(sin)</td><td>0.4</td><td>0.06</td></tr><tr><td>CDOM</td><td>5.15</td><td>0.06</td></tr><tr><td>NPP (mg m⁻³)</td><td>1.19</td><td>0.25</td></tr><tr><td>MLD (m)</td><td>9.7</td><td>3.2</td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** *pCO₂* reconstructed in function of the measured. The red line represents the linear relationship. The data represent 10% of the total of data used for the validation data set.</figcaption></figure>

outliers that disturb the estimation of the MLR coefficients. For the reconstruction of *pCO₂*, with this identifiable point, it is quite easy to organize a flag system. The flag can give information about the quality of the reconstructions based on the root mean square errors (RMSEs) of the neurons used for the reconstruction. The difference obtained for *pCO₂* in each neuron ranges from 0 to 56 µatm (Table 3), but 58% of the values observed are under 30 µatm. The difference can be quite high for a parameter such as SST, with a maximum value of 1.9 °C, but most of the values are lower than 1 °C and CDOM ranges from 0 to 5.15 (Table 3). The other parameters have quite low variability, such as MLD, which ranges from 0 to 9.7 m. The average is 2 to 3 times lower than the maximum value observed, which gives a low value for all the satellite parameters.

The *pCO₂* validation data set gives a quite good correlation (*R* = 0.93) with the results of the reconstruction method (Fig. 12), the root mean square (rms) being 36.7 µatm: 12% of the validation data have a value higher than 20 µatm, and 45% have a value between 20 and 30 µatm (Fig. 13). The characteristics of time, SST, MLD, CDOM, Chl *a*, and NPP do not explain the difference observed in the reconstruction.
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Difference between *pCO₂* reconstructed and measured in function of time. The red crosses and black dotted lines represent differences greater than 20 µatm. The data represent 10% of the total of data used for the validation data set.</figcaption></figure>

A reconstruction has been done using the satellite data from 1998 to 2011. The seasonal cycle of *pCO₂* is well reproduced and in agreement with the results of other studies. The maximum is observed in winter, with a *pCO₂* of 437 µatm on average, while the level is 274 µatm in summer. These values are comparable to the averages estimated in the central Baltic Sea of 500 µatm in summer and 150 µatm in winter (Wesslander, 2011). The *pCO₂* decreases in April due to the biological activity and increases slowly in September (Fig. 14).

We also evaluate these results by comparing them with modelling results. The model output used in the present study is from a process-oriented biogeochemical ocean model in which the Baltic Sea is divided into 13 natural sub-basins (e.g. Omstedt et al., 2009; Norman et al., 2013b). The properties of each sub-basin are horizontally averaged and vertically resolved, and the various sub-basins are horizontally coupled to each other using strait flow models. The model is forced by meteorological gridded data with a 3 h temporal resolution and by river run-off and net precipitation data with a monthly resolution (Omstedt et al., 2005). To compare the output of the model with our results, we couple the 13 basins in optic to reduce the number to three basins, corresponding to the three basins defined in Fig. 1. The modelled and estimated *pCO₂* are compared for the entire Baltic Sea and for the three basins from 1998 to 2009 (Fig. 14). The seasonal cycle for the entire Baltic Sea is well reproduced with a quite good correlation (*R* = 0.7) between the modelled and estimated *pCO₂* values (Table 4), whose standard deviations differ by 74 µatm. The modelled and estimated *pCO₂* values for the gulfs of Finland and Bothnia are not as correlated (*R* = 0.6), while the order of magnitude of the variability of *pCO₂* in the Gulf of Finland as calculated with SOMLO is closer to the model estimate (122 µatm for the modelled and 142 µatm for the estimated *pCO₂*). This lower correlation could be due to the lower number of data in this region available for these basins. The central basin is well reproduced, but the amplitude of the seasonal *pCO₂* cycle is lower

<table><caption>**Table 4.** Coefficient of correlation (*R*) between the modelled *pCO₂* and the *pCO₂* data estimated and the SD for the model and the data.</caption><thead><tr><th>Basin</th><th><i>R</i></th><th>SD model (µatm)</th><th>SD data (µatm)</th></tr></thead><tbody><tr><td>Baltic Sea</td><td>0.7</td><td>48</td><td>96</td></tr><tr><td>Central basin</td><td>0.7</td><td>69</td><td>92</td></tr><tr><td>Bothnian Basin</td><td>0.6</td><td>44</td><td>101</td></tr><tr><td>Gulf of Finland</td><td>0.6</td><td>122</td><td>144</td></tr></tbody></table>

in the simulation. In the southwest part of the Baltic Sea, SOMLO underestimates the *pCO₂* concentration by 60 µatm compared with the model. In the eastern and western parts of the basin, SOMLO produces good estimates of *pCO₂* compared with the model, with an average difference of 20 µatm. In Omstedt et al. (2009), the simulated *pCO₂* agrees quite well with the calculated values based on observations in the eastern Gotland Basin.

A simple flag was constructed to monitor the reconstruction quality and give an idea of the confidence in the estimated *pCO₂*. The difference between the estimated and neutral values is computed. The flag equals 1 for classes in which the average difference is less than 20 µatm, equals 2 for an average difference of 20–30 µatm, and equals 3 for higher average differences. In the example shown here, the flag values are high (i.e. 3), so the confidence in the reconstruction is low, but some points have flag values of 1 or 2 (Fig. 15d, e, and f) so the reconstruction is more reliable. On the geographic map (Fig. 15d, e, and f), the values of 4 correspond to the presence of ice, which is estimated using the satellite data of the National Snow and Ice Data Center based on NOAA level 3 data Njoku (2007). The flag gives confidence in our reconstruction; for example, in March 2010 (Fig. 15a), the southern portion of the map (i.e. the Bornholm and Arkona basins) shows lower *pCO₂* values than does the northern portion and than in February (not shown here). In March 2010, this region corresponds to a flag value of 2, which was described medium confidence. In July 2010, the flag value is quite good and the variability of *pCO₂* seems to be in line with the monthly variability (Fig. 15b and e). In September 2010, the value of *pCO₂* has a good order of magnitude when the flag is 2 but seems slightly too high when there is poor confidence (i.e. a flag value of 3) (Fig. 15c and f).

In conclusion, the reconstruction of *pCO₂* needs to be improved to increase the confidence in the reconstruction data, particularly in the gulfs.

### 3.2.4 Sensitivity analysis

In order to estimate the sensitivity of the reconstruction of the *pCO₂* to noisy data, white noise was added on all incomplete data before the reconstruction. We performed three tests by adding white noise to each parameter, which was related to the SD of that parameter. The tested configurations were 1 × sigma, 2 × sigma, and 0.5 × sigma.
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Comparison between modelled *pCO₂* (dotted lines) and *pCO₂* estimated using the SOM linear method (solid lines) for (a) the Baltic Sea (BS, blue), (b) the Gulf of Finland (GF, grey), (c) the Gulf of Bothnia (BG, red), and (d) the central basin (C, green).</figcaption></figure>

<table><caption>**Table 5.** Coefficient of correlation (*R*) between the measured *pCO₂* and the *pCO₂* data estimated with SOMLO and the rms in function of the sigma apply for the noise 0.5, 1 and 2 times the sigma.</caption><thead><tr><td>X SIGMA</td><td>R</td><td>rms (μatm)</td><td>Numbers of tests</td></tr></thead><tbody><tr><td>0.5</td><td>0.85</td><td>61</td><td>701</td></tr><tr><td>1</td><td>0.80</td><td>77</td><td>733</td></tr><tr><td>2</td><td>0.59</td><td>134</td><td>709</td></tr></tbody></table>

<figure><img src="image_17.png" /><figcaption>**Figure 15.** (a, b, and c) reconstruction of the *pCO₂* map. (d, e, and f) the flag for each map. (a, d) March 2010, (b, e) July 2010, (c, f) September 2010. The flag value 1 corresponds to high confidence, 2 to medium confidence, and 3 to low confidence.</figcaption></figure>

The results as seen in Table 5 imply that the method is sensitive to noisy data, since when increasing the noise we progressively get worse reconstructions. It is important to note that the values obtained when using 1 × sigma give an average rms of 77 and an average *R* coefficient of 0.80, which indicate that the method does not degrade too much.

Another area where errors might appear is the completion of missing data. In the data set there were vectors with missing parameters. We kept only those missing up to three parameters. The difference observed between the *pCO₂* measured and reconstructed in function of month and of number of missing parameters can be seen in Fig. 16. We chose to present the data this way to also give a better understanding of the seasonal variation of the reconstruction. As seen in Table 6 the imputation method does introduce errors. When reconstructing the *pCO₂* with data that have no missing data, we obtain a correlation of 0.96 and an rms of 25.7 μatm; when reconstructing the *pCO₂* with data missing three parameters, the results are less reliable, with a correlation of 0.81 and an rms of 51.4 μatm. The rms we obtain in the case of the reconstruction using data that contain one missing value is higher than the one obtained when using data containing two missing values, while retaining a higher cor-
<figure><img src="image_18.png" /><figcaption>**Figure 16.** Absolute difference between *pCO₂* reconstructed and *pCO₂* measured in function of the month. The colour bar corresponds to the number of empty parameter before the reconstruction (blue: 0 missing data; light blue: 1 missing datum; yellow: 2 missing data; red: 3 missing data).</figcaption></figure>

<table><caption>**Table 6.** Coefficient of correlation (*R*) between the measured *pCO₂* and the *pCO₂* data estimated and the SD for the model and the data in function of the number of missing parameters.</caption><thead><tr><td>Number of parameter missing</td><td>R</td><td>rms (μatm)</td></tr></thead><tbody><tr><td>0</td><td>0.96</td><td>25.7</td></tr><tr><td>1</td><td>0.93</td><td>39.5</td></tr><tr><td>2</td><td>0.86</td><td>31.9</td></tr><tr><td>3</td><td>0.81</td><td>51.4</td></tr></tbody></table>

relation. This could be due to the higher number of vectors containing only one missing value.

## 4 Discussion and conclusions

In this paper, we used the SOMLO methodology to reconstruct the *pCO₂* from satellite data for the Baltic Sea. SOMLO was used to accommodate the non-linearity of the mechanics driving the *pCO₂*. It uses artificial neural networks to classify data into situations and then performs a reconstruction by using an MLR in each class. The process involves classifying the explicative parameters (i.e. SST, CDOM, Chl *a*, time, NPP, and MLD) and then using the linear regression coefficients corresponding to that class in order to reconstruct the *pCO₂*. The satellite data used were also completed using an iterative process of SOM training.

We also performed a statistical analysis of the reconstructions obtained, which allowed us to add a flag to each class, informing us of the quality of the reconstruction obtained. This could both influence further numerical modelling of other phenomena depending on the *pCO₂* and allow for an informed interpretation of the reconstructions obtained.

The current results obtained using this method, based on 1445 vectors, gave a high correlation coefficient of 0.93% and an rms of 36 μatm. These results are promising given the

conditions under which we obtained them since, in addition to having a limited number of in situ *pCO₂* measurements, the co-localized satellite data were frequently incomplete.

In comparison, existing studies were performed over the North Atlantic and North Pacific, based on a minimum of 10 000 data points (which take into account all the data from SOCAT) to a maximum of 800 000 data points (e.g. Friedrich and Oschlies, 2009; Hales et al., 2012; Landschützer et al., 2013). Friedrich and Oschlies (2009) obtained an RMSE of 19 μatm. A similar study over the totality of the Atlantic Ocean obtained an RMSE of 17 μatm for independent time series (Landschützer et al., 2013). Hales et al. (2012) obtained an RMSE of 20 μatm with a correlation coefficient of 0.81. The RMSE we obtained here was higher than that obtained in a previous study of the Atlantic Ocean but, taking into account the much smaller number of data available, and a stronger spatial *pCO₂* variability might in the Baltic Sea compared to the open ocean. the results presented are promising.

The organization of the values of various MLR coefficients over each class indicated that all the satellite data parameters are important for reconstructing the *pCO₂* in the Baltic Sea, even if only in certain cases. Improved satellite data availability could therefore also improve the performance of our reconstruction.

This study could be further developed so as to reconstruct the spatial fields of *pCO₂*. Specifically, one could imagine a Bayesian approach that would select which class to use for the MLR by also taking into account the potential classes attributed to the neighbouring grid points of a geographic study area. This, however, remains dependent on the acquisition of additional in situ measurements to allow for the robust estimation of such Bayesian probabilities.

Many programs exist for the acquisition of new data. Data from the Östergarnsholm site are still being acquired; 2012 did not yield much data from this site, and the data from 2013 and 2014 still need to be validated. In time, the SMHI station could also supply additional data. The cargo ship transect data are not yet available for 2012–2014, but these measurements will continue, and some data will soon be available. Data are also being gathered from ferries sailing the Gothenburg–Kemi–Oulu–Lübeck–Gothenburg route. This Gothenburg transect is weekly (see http://www.hzg.de/imperia/md/content/ferryboxusergroup/presentations/fb-ws2011_karlson.pdf). The first tests of these data were conducted in 2010 and 2011, so some data should soon be available. In addition, new measurements of *pCO₂* began in 2012 at the Utö Atmospheric and Marine Research Station (see http://en.ilmatieteenlaitos.fi/GHG-measurement-sites#Uto).

Given the amount of new data soon to be available, we remain optimistic that comprehension and statistical modelling of *pCO₂* in the Baltic Sea will continue to improve in coming years.
*Acknowledgements.* We wish to thank the two reviewers for their useful comments. The authors would like to thank Tiit Kutser and Melissa Chierici for their valuable help with data collection. We would like to thank Sylvie Thiria for her advice on the SOM classifications and Marion Leduc-Leballeur for her help with the method of filling the gaps in satellite data. We wish to thank Maria Norman for providing us with the model output. We would also like to thank Bernd Schneider and the BONUS/Baltic-C programs for cargo ship observations. We would like to thank the Swedish National Space Board (number 120/11:3) for financing the project.

Edited by: C. Garbe

## References

Behrenfeld, M. J. and Boss, E.: Beam attenuation and chlorophyll concentration as alternative optical indices of phytoplankton biomass, J. Mar. Res., 64, 431–451, 2006.

Behrenfeld, M. J., Boss, E., Siegel, D. A., and Shea, D. M.: Carbon-based ocean productivity and phytoplankton physiology from space, Global Biogeochem. Cy., 19, 1–25, 2005.

Bergstrom, S.: River runoff to the Baltic Sea: 1950–1990, Ambio, 23, 280–287, 1994.

Borges, A. V., Delille, B., and Frankignoulle, M.: Budgetting sinks and sources of CO₂ in the coastal ocean: Diversity of ecosystems counts, Geophys. Res. Lett., 32, L14601, doi:10.1029/2005GL023053, 2005.

Boyer Montégut, C., Madec, G., Fischer, A. S., Lazar, A., and Ludi-cone, D.: Mixed layer depth over the global ocean: An examination of profile data and a profile-based climatology, J. Geophys. Res., 109, C12003, doi:10.1029/2004JC002378, 2004.

Burchard, H. and Bolding, K.: GETM: A General Estuarine Transport Model, Scientific Documentation, European Commission, Joint Research Centre, Institute for Environment and Sustainability, 2002.

Casey, K. S., Brandon, T. B., Cornillon, P., and Evans, R.: The past, present, and future of the AVHRR Pathfinder SST program, in: Oceanography from Space, Springer, the Netherlands, 273–287, 2010.

Chen, C.-T. A. and Borges, A. V.: Reconciling opposing views on carbon cycling in the coastal ocean: Continental shelves as sinks and near-shore ecosystems as sources of atmospheric CO₂, Deep-Sea Res. Pt. II, 56, 578–590, doi:10.1016/j.dsr2.2009.01.001, 2009.

Chierici, M., Olsen, A., Johannessen, T., Trinañes, J., and Wan-ninkhof, R.: Algorithms to estimate the carbon dioxide uptake in the northern North Atlantic using shipboard observations, satellite and ocean analysis data, Deep-Sea Res. Pt. II, 56, 630–639, 2009.

Cloern, J. E., Grenz, C., and Vidergar-Lucas, L.: An empirical model of the phytoplankton chlorophyll: carbon ratio-the conversion factor between productivity and growth rate, Limnol. Oceanogr., 40, 1313–1321, 1995.

Corbière, A., Metzl, N., Reverdin, G., Brunet, C., and Takahashi, T.: Interannual and decadal variability of the oceanic carbon sink in the North Atlantic subpolar gyre, Tellus B, 59, 168–178, 2007.

Dickson, A. and Millero, F.: A comparison of the equilibrium constants for the dissociation of carbonic acid in seawater media, Deep-Sea Res. A, 34, 1733–1743, 1987.

Dreyfus, G.: Neural networks: methodology and applications, Springer, Berlin, New York, 2005.

Dutt, A. and Rokhlin, V.: Fast Fourier transforms for nonequispaced data, II, Appl. Comput. Harmon. A, 2, 85–100, 1995.

Eppley, R. W.: Temperature and phytoplankton growth in the sea, Fish. Bull., 70, 1063–1085, 1972.

Friedrich, T. and Oschlies, A.: Neural network-based estimates of North Atlantic surface pCO₂ from satellite data: A methodological study, J. Geophys. Res.-Oceans, 114, 1–12, C03020, doi:10.1029/2007JC004646, 2009.

Grasshoff, K., Kremling, K., and Ehrhardt, M.: Methods of seawater analysis, third Edn., VCH Publishers, Weinheim, New York, Chichester, Brisbane, Singapore and Toronto, 1999.

Greengard, L. and Lee, J.-Y.: Accelerating the nonuniform fast Fourier transform, SIAM rev., 46, 443–454, 2004.

Hales, B., Strutton, P. G., Saraceno, M., Letelier, R., Takahashi, T., Feely, R., Sabine, C., and Chavez, F.: Satellite-based prediction of pCO₂ in coastal waters of the eastern North Pacific, Prog. Oceanogr., 103, 1–15, doi:10.1016/j.pocean.2012.03.001, 2012.

Hjalmarsson, S., Wesslander, K., Anderson, L. G., Omstedt, A., Perttilä, M., and Mintrop, L.: Distribution, long-term development and mass balance calculation of total alkalinity in the Baltic Sea, Cont. Shelf Res., 28, 593–601, 2008.

Jamet, C., Moulin, C., and Lefèvre, N.: Estimation of the oceanic pCO₂ in the North Atlantic from VOS lines in-situ measurements: parameters needed to generate seasonally mean maps, Ann. Geophys., 25, 2247–2257, doi:10.5194/angeo-25-2247-2007, 2007.

Jolliffe, I. T.: Principal component analysis, Springer, New York, 2002.

Kohonen, T.: The self-organizing map, P. IEEE, 78, 1464–1480, 1990.

Landschützer, P., Gruber, N., Bakker, D. C. E., Schuster, U., Nakaoka, S., Payne, M. R., Sasse, T. P., and Zeng, J.: A neural network-based estimate of the seasonal to inter-annual variability of the Atlantic Ocean carbon sink, Biogeosciences, 10, 7793–7815, doi:10.5194/bg-10-7793-2013, 2013.

Laruelle, G. G., Dürr, H. H., Slomp, C. P., and Borges, A. V.: Evaluation of sinks and sources of CO₂ in the global coastal ocean using a spatially-explicit typology of estuaries and continental shelves, Geophys. Res. Lett., 37, L15607, doi:10.1029/2010GL043691, 2010.

Lee, Z.-P., Darecki, M., Carder, K. L., Davis, C. O., Stramski, D., and Rhea, W. J.: Diffuse attenuation coefficient of downwelling irradiance: An evaluation of remote sensing methods, J. Geophys. Res.-Oceans, 110, C02017, doi:10.1029/2004JC002573, 2005.

Lefèvre, N., Watson, A. J., Olsen, A., Ríos, A. F., Pérez, F. F., and Johannessen, T.: A decrease in the sink for atmospheric CO₂ in the North Atlantic, Geophys. Res. Lett., 31, L07306, doi:10.1029/2003GL018957, 2004.

Lefèvre, N., Watson, A. J., and Watson, A. R.: A comparison of multiple regression and neural network techniques for mapping in situ pCO₂ data, Tellus B, 57, 375–384, doi:10.1111/j.1600-0889.2005.00164.x, 2005.
Lewis, E. and Wallace, D. W. R.: Program developed for CO₂ system calculations, Tech. rep., ORNL/CDIAC, US Department of Energy, Oak Ridge, Tennessee, 1998.

Lüger, H., Wallace, D. W., Körtzinger, A., and Nojiri, Y.: The *pCO₂* variability in the midlatitude North Atlantic Ocean during a full annual cycle, Global Biogeochem. Cy., 18, GB3023, doi:10.1029/2003GB002200, 2004.

Malek, M. A., Harun, S., Shamsuddin, S. M., and Mohamad, I.: Imputation of time series data via Kohonen self organizing maps in the presence of missing data, Engineering and Technology, 41, 501–506, 2008.

Merbach, C., Culberson, C. H., and Hawley, J. E.: Measurements of the apparent dissociation constants of carbonic acid in seawater of atmospheric pressure, Limnol. Oceanogr., 18, 897–907, 1973.

Morel, A. and Berthon, J.-F.: Surface pigments, algal biomass profiles, and potential production of the euphotic layer: Relationships reinvestigated in view of remote-sensing applications, Limnol. Oceanogr, 34, 1545–1562, 1989.

Morel, A. and Gentili, B.: A simple band ratio technique to quantify the colored dissolved and detrital organic material from ocean color remotely sensed data, Remote Sens. Environ., 113, 998–1011, 2009.

Nakaoka, S., Telszewski, M., Nojiri, Y., Yasunaka, S., Miyazaki, C., Mukai, H., and Usui, N.: Estimating temporal and spatial variation of ocean surface *pCO₂* in the North Pacific using a self-organizing map neural network technique, Biogeosciences, 10, 6093–6106, doi:10.5194/bg-10-6093-2013, 2013.

Njoku, E.: AMSR-E/Aqua L2B Surface Soil Moisture, Ancillary Parms, & QC EASE-Grids v002, Información digital. Actualizado diariamente, Boulder, Colorado EEUU: National Snow and Ice Data Center, 2007.

Norman, M., Raj Parampil, S., Rutgersson, A., and Sahlée, E.: Influence of coastal upwelling on the air-sea gas exchange of CO₂ in a Baltic Sea Basin, Tellus B, 65, 21831, doi:10.3402/tellusb.v65i0.21831, 2013a.

Norman, M., Rutgersson, A., and Sahlée, E.: Impact of improved air-sea gas transfer velocity on fluxes and water chemistry in a Baltic Sea model, J. Mar. Syst., 111, 175–188, doi:10.1016/j.jmarsys.2012.10.013, 2013b.

Omstedt, A., Elken, J., Lehmann, A., and Piechura, J.: Knowledge of the Baltic Sea physics gained during the BALTEX and related programmes, Prog. Oceanogr., 63, 1–28, 2004.

Omstedt, A., Chen, Y., and Wesslander, K.: A comparison between the ERA40 and the SMHI gridded meteorological databases as applied to Baltic Sea modelling, Nord. Hydrol., 36, 369–380, 2005.

Omstedt, A., Gustafsson, E., and Wesslander, K.: Modelling the uptake and release of carbon dioxide in the Baltic Sea surface water, Cont. Shelf Res., 29, 870–885, 2009.

Platt, T. and Sathyendranath, S.: Estimators of primary production for interpretation of remotely sensed data on ocean color, J. Geophys. Res.-Oceans, 98, 14561–14576, 1993.

Rutgersson, A., Norman, M., Schneider, B., Pettersson, H., and Sahlée, E.: The annual cycle of carbon dioxide and parameters influencing the air–sea carbon exchange in the Baltic Proper, J. Mar. Syst., 74, 381–394, 2008.

Sasse, T. P., McNeil, B. I., and Abramowitz, G.: A novel method for diagnosing seasonal to inter-annual surface ocean carbon dynamics from bottle data using neural networks, Biogeosciences, 10, 4319–4340, doi:10.5194/bg-10-4319-2013, 2013.

Schafer, J. L. and Graham, J. W.: Missing data: our view of the state of the art, Psychol. Methods, 7, 147–77, 2002.

Schneider, B. and Kaitala, S.: Identification and quantification of plankton bloom events in the Baltic Sea by continuous *pCO₂* and chlorophyll a measurements on a cargo ship, J. Mar. Syst., 59, 238–248, 2006.

Schneider, B., Nausch, G., Nagel, K., and Wasmund, N.: The surface water CO₂ budget for the Baltic Proper: a new way to determine nitrogen fixation, J. Mar. Syst., 42, 53–64, 2003.

Schneider, B., Kaitala, S., Raateoja, M., and Sadkowiak, B.: A nitrogen fixation estimate for the Baltic Sea based on continuous *pCO₂* measurements on a cargo ship and total nitrogen data, Cont. Shelf Res., 29, 1535–1540, 2009.

Schomberg, H. and Timmer, J.: The gridding method for image re-construction by Fourier transformation., IEEE T. Med. Imaging, 14, 596–607, doi:10.1109/42.414625, 1995.

Schuster, U., McKinley, G. A., Bates, N., Chevallier, F., Doney, S. C., Fay, A. R., González-Dávila, M., Gruber, N., Jones, S., Krijnen, J., Landschützer, P., Lefèvre, N., Manizza, M., Mathis, J., Metzl, N., Olsen, A., Rios, A. F., Rödenbeck, C., Santana-Casiano, J. M., Takahashi, T., Wanninkhof, R., and Watson, A. J.: An assessment of the Atlantic and Arctic sea-air CO₂ fluxes, 1990–2009, Biogeosciences, 10, 607–627, doi:10.5194/bg-10-607-2013, 2013.

Siegel, H. and Gerth, M.: Baltic Sea environment fact sheet Sea Surface Temperature in the Baltic Sea in 2011, HELCOM Baltic Sea Environment Fact Sheets, http://www.helcom.fi/baltic-sea-trends/data-maps/ (last access: 1 June 2015), 2012.

Stocker, T. F., Qin, D., Plattner, G.-K., K., Tignor, M., Allen, S. K., Boschung, J., Nauels, A., Xia, Y., Bex, V., and Midgley, P. M.: Climate Change 2013, The Physical Science Basis, Working Group I Contribution to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change-Abstract for decision-makers, Tech. rep., Groupe d'experts intergouvernemental sur l'evolution du climat/Intergovernmental Panel on Climate Change-IPCC, C/O World Meteorological Organization, 7 bis Avenue de la Paix, CP 2300 CH-1211 Geneva 2 (Switzerland), 2013.

Takahashi, T., Sutherland, S. C., Wanninkhof, R., Sweeney, C., Feely, R. A., Chipman, D. W., Hales, B., Friederich, G., Chavez, F., Sabine, C., Watson, A., Bakker, D. C., Schuster, U., Metzl, N., Yoshikawa-Inoue, H., Ishii, M., Midorikawa, T., Nojiri, Y., Körtzinger, A., Steinhoff, T., Hoppema, M., Olafsson, J., Arnarson, T. S., Tilbrook, B., Johannessen, T., Olsen, A. O., Tilbrook, B., Bellerby, R., Wong, C. S., Delille, B., Bates, N. R., and De Baar, H. J. W.: Climatological mean and decadal change in surface ocean *pCO₂*, and net sea–air CO₂ flux over the global oceans, Deep-Sea Res. Pt. II, 56, 554–577, doi:10.1016/j.dsr2.2008.12.009, 2009.

Telszewski, M., Chazottes, A., Schuster, U., Watson, A. J., Moulin, C., Bakker, D. C. E., González-Dávila, M., Johannessen, T., Körtzinger, A., Lüger, H., Olsen, A., Omar, A., Padin, X. A., Ríos, A. F., Steinhoff, T., Santana-Casiano, M., Wallace, D. W. R., and Wanninkhof, R.: Estimating the monthly *pCO₂* distribution in the North Atlantic using a self-organizing neural network, Biogeosciences, 6, 1405–1421, doi:10.5194/bg-6-1405-2009, 2009.
Weiss, R. F.: CO₂ in water and seawater: The solubility of a non ideal gas, Mar. Chem., 2, 203–215, 1974.

Wesslander, K.: The carbon dioxide system in the Baltic Sea surface waters, PhD thesis, University of Gothenburg, Göteborg, 2011.

Wesslander, K., Omstedt, A., and Schneider, B.: Inter-annual variation of the air-sea CO₂ balance in the southern Baltic Sea and the Kattegat, in: EGU General Assembly Conference Abstracts, vol. 11, p. 8629, 2009.