Sea-Air Exchange of Methane in Shallow Inshore Areas of the Baltic Sea

We report sea-air fluxes of methane in physically and biologically distinct inshore habitats of the Baltic Sea with the goal to establish empirical relationships that allow upscaling of local site-specific flux measurements. Flux measurements were conducted using floating chambers with and without bubble shields, and by using a boundary layer gas transfer model before, during, and after an annually occurring algal bloom from June to October 2019. Water and air temperature, salinity, wind, sediment organic content, and organic content of floating algal biomass were found to successfully discriminate the different habitats in terms of methane flux, both over periods of days and over a season. Multivariate statistical analysis was used to establish the relative environmental forcing of methane emissions over one growth season for each flux method. Floating algal biomass carbon and sediment organic content were identified as the most important controlling factors for methane emissions based on flux chamber measurements over a period of days to weeks, whereas water and air temperature and wind velocity were the most important factors based on the gas transfer model on these time scales. Over the season, water and air temperature were the most important controlling factors with both methods. We present a first attempt how our observations can be extrapolated to determine the coastal methane emission along the coastline.


INTRODUCTION
The trace gas methane contributes about 23% to the radiative balance of Earth's atmosphere and its increasing emissions contribute significantly to ongoing global warming (Mende et al., 2019). Atmospheric methane burdens have increased over pre-industrial levels 2.5 times making present and future methane levels highly relevant in discussions about climate forcing and effective mitigation measures (Nisbet et al., 2019). Aquatic ecosystems are responsible for about 53% of the total global methane emissions from anthropogenic and natural sources . Methane has anthropogenic sources from fertilizer application, animal food stocks, industrial and agricultural emissions, and fossil fuel combustion (Seitzinger et al., 2006;Saunois et al., 2016), but also natural sources such as thawing permafrost, wetlands, natural gas seeps, and emissions from lakes and the coastal ocean that react sensitively to climate change and human activities (Nisbet et al., 2019). A fundamental concern in budgeting these different sources are the discrepancies that arise between top-down or bottom-up approaches (Nisbet and Weiss, 2010;Crill and Thornton, 2017), and in the large range of global estimates for bottom-up approaches that vary between 594 and 880 Tg y −1 . Although emissions from lakes and terrestrial wetlands comprise the bulk of the global aquatic methane emissions, shallow-water coastal environments dominate the emissions from aquatic marine environments by area (Weber et al., 2019;Bižić et al., 2020;. Coastal ecosystems, and in particular the inshore habitats, are of interest for global methane emissions, because they are most directly affected by anthropogenic disturbances of the land-ocean interface (Smith and Hollibaugh, 1993;Gattuso et al., 1998;Battin et al., 2008;Torres-Pulliza et al., 2020). Inshore coastal habitats have also been identified as potentially important blue carbon repositories to mitigate CO 2 accumulation in the atmosphere (Macreadie et al., 2019), but the efficiency of this carbon sink is partially offset by methane emissions (Rosentreter et al., 2018).
Part of the problem with inshore bottom-up estimates of methane emissions lies in the challenge to scale emissions in habitats that vary strongly over both spatial and temporal scales. Inshore coastal habitats worldwide are very diverse and range from highly productive coastal ecosystems such as coastal wetlands to bare rocky cliffs largely devoid of significant sediment organic carbon accumulation or dense vegetation coverage. In inshore habitats, environmental factors that lead to large spatial and temporal variations in methane production and emission include vegetation density, organic richness of the bottom substrate, air and water temperature, wind velocity and direction relative to the coastline, and wave activity (Jeffrey et al., 2019). An additional factor is the release of methane as bubbles from the sediment, which can escape to the atmosphere before significant oxidation (Keller and Stallard, 1994;Joyce and Jewell, 2003;Jeffrey et al., 2019;Weber et al., 2019). Bubble emissions comprising up to 95% of methane emissions have been suggested (Casper et al., 2000;Walter et al., 2006;Bastviken et al., 2011;Delwiche et al., 2015), but these emissions are highly localized. Lastly, inshore methane emission also takes place via plant vascular tissue transport, but the impact from plants varies by season (Jeffrey et al., 2019) and life cycle phase (Van Der Nat and Middelburg, 1998;Van der Nat and Middelburg, 2000). For example, hornbeam (Typha albida) and reed (Phragmites australis) can reduce methane emissions by oxidation during their growth phase, but when the plants have matured the oxidation is significantly reduced.
A combination of biological, chemical, and physical properties of the inshore coastal environment can be used to evaluate correlations that weigh in the different environmental forcing factors for methane emission. Quantification of methane emissions from coastal habitats in which such steering variables are analyzed may also allow testing which parameters exert a dominating influence on methane fluxes in different coastal habitats. This parametric approach could ultimately be used for empirical scaling approaches in which the geographic and temporal variability of the steering parameters are used as guiding variables to estimate methane emissions over larger coastal areas. At present, however, we do not know the methane emissions of many coastal ecotypes and the variation of these emissions throughout the year -data that are key for a more accurate bottom-up upscaling of the coastal methane emission.
The Baltic Sea coastline has a large diversity of inshore ecotypes such as rocky shorelines, shallow bays, small fjärds, wetlands, small deltas, and coastal lagoons with bottom types ranging from organic-rich soft muds to unvegetated bare rocks (Bartley et al., 2001;Niemelä et al., 2015). Baltic Sea coastal wetlands have some of the highest methane emissions in the Baltic Sea (Heyer and Berger, 2000), but are most directly affected by land-side anthropogenic perturbances. Inshore water temperatures range from up to 27 • C in the summer to below freezing in the winter, along with air temperature variations of over 40 • C during the year. Light and temperature variability make productivity of the coastal habitat highly seasonal and locally variable. In addition, algal blooms increase the amount of organic matter transferred to sediment, which cause hypoxia and enhanced methanogenesis at shallow sediment depth (Bange et al., 2010;Maltby et al., 2018). The current work had as its purpose to evaluate the forcing of trace CH 4 gas emission rates by different environmental factors and to contribute directly measured near-shore trace gas exchange rates. This is a new approach for establishing scaled-up CH 4 fluxes in the littoral coastal zone that explicitly takes near-shore vegetation variations, bubble emissions, and seasonal variability of physical environmental forcing factors in these inshore habitats of the Baltic Sea into account. As the main objective the present study therefore lays the ground for scaling scenarios to arrive at empirically constrained trace gas emissions from inshore habitats that can be incorporated into comparisons of bottom-up and top-down emission budgets.

Study Location and General Approach
The island of Askö is located in the Trosa archipelago on the east coast of Sweden in the central Baltic Sea (Figure 1). Askö extends in NW-SE direction and is about 10 km long and 1 km wide and covers an area of 272 ha. The whole island is a nature reserve, covered by forest and open land and managed by a resident landholder, with cows and sheep who graze the land to help maintain the open landscape. The SW facing side of the island is dominated rocky cliffs and shallow embayments and is relatively open to the Baltic Sea, whereas the NE side of the island faces the archipelago sea with calmer wave conditions. The dominant wind directions in the area are S to SW and secondarily NW to NE (Station Landsort A 1 ).
We categorized five sampling areas of the island shoreline based on salinity, wind velocity, water temperature, air temperature, water depth, organic content of floating algal biomass, shore vegetation, and sediment organic carbon content into four different habitats (Table 1). Discriminators were indicator species for eutrophication (see descriptions for habitats A and B below), algal mat occurrence, occurrence of bubbles  during sampling, indicators of human activities such as cut reeds, and the occurrence of hydrogen sulfide. Habitat A was characterized by reeds and soft organic-rich clays with a water depth up to 0.5 m in small, relatively wind-protected bights. Deciduous tree shorelines characterize this habitat. The shoreline has dense reed areas (P. australis), hornbeam (T. albida), seagrasses (Carex pendula and Sagittaria sagittifolia), and is covered by thick layers of algae mats (Archaeplastida or Cyanophyceae) during July and August, when water lilies (Nymphaea alba) also occupy the surface. Visible bubble emissions appear between June and August in combination with a strong smell of hydrogen sulfide. Human activities such as reed cutting and agriculture affect this habitat.
Habitat B also has soft, organic-rich sediment, a generally vegetated shoreline with deciduous trees and water depths between 0.5 and 1.0 m. There are dense areas of sessile aquatic vegetation with different species of brown algae (Phaeophyceae) and reed (P. australis), hornbeam (T. albida) and dense summer coverage of bladder wrack (Fucus vesiculosus). Benthic algae cover about half of the area. Visible bubble emissions appear in July and August, but hydrogen sulfide is generally absent. Human activities and agriculture generally do not affect this habitat.
In habitat C the sea bottom consists of sandy sediment and unvegetated rocky hard ground. It is located in more open areas with higher currents and wave activity and water depth-up to 2.0 m. There is generally no sessile aquatic vegetation, algal blooms, or any significantly smell of hydrogen sulfide.
Habitat D is characterized by its unvegetated hard ground and strong water currents. The surrounding area of this habitat has rocky cliffs with sparse plants, sessile aquatic vegetation, and summer algal accumulations, and water depths between 0.5 and 1.5 m. There is also no smell of hydrogen sulfide in this habitat.

Field Methods
Methane fluxes were determined using round, anchored, ca. 8 L large floating chambers with and without bubble shield between June and October 2019 following the design described in Schilder et al. (2016). Chamber fluxes were determined with and without a bubble shield at the same station in order to quantify the contribution of ebullition to the total flux. The bubble shield consisted of a 0.12 m 2 large transparent PVC plate that was mounted underneath the chamber such that it was positioned 15 cm below the floating chamber. Time intervals for the flux measurements varied between 15 and 78 h, with 89% of the measurements were conducted with less than 32 h intervals. Only sampling intervals with linear increase in CH 4 concentration with r 2 above 0.78 were considered as significant gradients, and 85% had r 2 above 0.90. At stations, where intermittent bubbling rapidly increased the methane concentration over an interval of less than 15 h, only that time interval was used in the final calculation and the other time intervals were ignored. At several stations, the flux measurement was reset to conduct multiple flux measurements over a period of up to 6 days per monthly measuring campaign. Air samples were collected from the chamber headspace every 24 h over a period of 3-6 days. Example time series evolution of concentrations in Supplementary Figures 3, 4. Sampling was done according to the method of Keller and Stallard (1994). To collect samples, the sampling syringe was connected via a three-way stopcock to the chamber and was flushed at least 10 times before filling the syringe with 60 ml of air from the chamber. The first 20 ml of the air sample was emptied to the atmosphere to flush the transfer needle. Subsequently, the remaining air in the syringe was transferred through a butyl rubber stopper into a 20 ml glass vial filled, without headspace, with saturated salt solution.
Dissolved methane was collected with the headspace equilibration method described in Bastviken et al. (2003). A 60 ml syringe was filled, flushed and re-filled with seawater three times. A total of 40 ml of the water was equilibrated in the syringe with 20 ml air for 1 min and exactly 20 ml air was then transferred to a glass vial filled, without headspace, with saturated salt solution by replacing the salt solution with the sample air. The glass vials were kept upside-down until all vials were analyzed.
Three to five locations were selected in each habitat totaling 23 sampling stations. On each sampling occasion and before collecting each sample the following data were recorded: air and water concentration of CH 4 , date, time, salinity in water, air temperature, O 2 concentration in water, water temperature and GPS location. Wind data and air temperature were obtained from the database of the Swedish Meteorological and Hydrological Institute (SMHI) (SMHI, 2010). O 2 , water temperature, and salinity data were obtained with a handheld WTW Salinometer.

Laboratory CH 4 Analysis
Water and air samples were analyzed by gas chromatography on a Shimadzu GC-8A with FID detector and a 1 ml loop injection system. The instrument was calibrated using the average of multiple gas standards that were analyzed during each analysis session with concentrations of 11, 100, and 10,000 ppmv. Replicate standard analyses were performed with a coefficient of variation between 1 and 3%. For the analysis of air and airequilibrated water samples, 2 ml of a sample were injected on the loop of the Gas Chromatograph (GC) compensating for the withdrawn volume in the air-equilibrated samples with an equivalent volume of salt solution of the same salinity as the sample (6 ). Gas concentrations in (ppm) were converted to (nmol/L) according to Johnson et al. (1990) where C • W is the concentration of the dissolved gas in the sample (mol/L), n g are the moles in headspace after equilibration, n w are the moles in the aqueous phase after equilibration, and V w is the volume of the water sample. To calculate ng, the ideal gas was used where P is the air pressure in Pa, p is the partial pressure of the gas in Pa, V g is the volume of the headspace (L), T a is air temperature (K) and R is the universal gas constant 8314.46 (L Pa mol −1 K −1 ). n w was calculated using Henry's law constant and the volume of water V w according to Johnson et al. (1990)

Flux Calculations
The flux calculations reported in this study are explained and discussed in detail in Supplementary Material and are briefly summarized here: The flux (F) of methane in the floating chambers, henceforth described as the "chamber method" was calculated according to: where F is the flux (mg/m 2 /day), M is the molar mass for the component CH 4 (g/mol), C t n+1 , C t n are concentrations (nmol/L) of component in air sample at sampling times t n+1 and t n , respectively, t n+1−n is the time difference between sampling times t n+1 and t n in hours. V ch is the volume of the flux chamber (L) and A ch is the water surface area of the flux chamber (m 2 ). The volume of the floating flux chamber was 7.5 L with a surface area of 0.08 m 2 and followed the general design described by Schilder et al. (2016). Fluxes determined from dissolved methane concentrations and wind velocity, henceforth termed "gas transfer model" were calculated according to the general boundary layer flux equation: where F is the flux (mg/m 2 /day), M is the molar mass for the component CH 4 (g/mol), k is the gas transfer velocity (cm/h), C w is the concentration of the dissolved gas in surface water (nmol/L), C a is the concentration of the gas in the air at the surface (nmol/L), α is the Ostwald solubility coefficient (dimensionless). From the Ostwald solubility coefficient (α) for methane from the Bunsen solubility coefficient (β) and water temperature (T w ) (Battino, 1984) where β is the dimensionless Bunsen solubility coefficient of the dry mole fraction where T w is the water temperature (K), S is the salinity ( ). A 1 , A 2 , A 3 , B 1 , B 2 , B 3 are constants specific for CH 4 (Wanninkhof, 2014)). A number of different models are used to determine the gas transfer velocity k and for this context the equation used by Wanninkhof (2014) was used: where k is the gas transfer velocity (cm/hour), U 10 is the wind velocity at 10 m height (m/s), and Sch is the Schmidt number for the gas component (dimensionless). We have also used the k models of Ho et al. (2006); Nightingale et al. (2000b), Wanninkhof et al. (2009), andWanninkhof (1992) and calculated averages from the different k models (Figure 2). The average k-values were close to those of the Wanninkhof (2014) relationship for the range of wind velocities encountered in our study and we have therefore decided to use this relationship. The k model for CH 4 used by Rosentreter et al. (2017) and Jeffrey et al. (2018) for estuaries (mangroves) was considered, and in these studies, tidal currents were a significant factor. However, the central Baltic Sea region in the study does not have any significant permanent or tidal currents, and estuarine models may not provide better results for the setting of this study. 2 Wind velocity measurements were taken from the coastal measuring buoy at Askö (Askö Kustmätsystem). Wind velocity was adjusted by using measurements at 1.5 m height to get the corresponding wind velocity at 10 m height according to Amorocho and DeVries (1980).
U z is the wind velocity at z meters height (m/s), U 10 is the wind velocity at 10 m height (m/s), c 10 is the surface drag coefficient for shallow water depth c 10 = 1.3 × 10 −3 , κ is the von Kármán constant κ = 0.41. The Schmidt number for Baltic Sea brackish water (salinity ∼ 6 ) with measured salinity (S balticsea ) was calculated by interpolation of the Schmidt number for fresh water (salinity 2 https://www.smhi.se/en/theme/tides-1.11272 FIGURE 2 | Relationship between gas transfer velocity and wind velocity at 10 m height for five different gas transfer models and the corresponding mean. The relationship by Wanninkhof (2014) was used in this study.
where t w is the water temperature ( • C), A, B, C, D, and E are constants which depends on the gas component (Wanninkhof, 2014). The constants A,B,C, D, and E in Equation 9 to calculate the Schmidt number were used for fresh water and sea water (Wanninkhof, 2014). A detailed description with examples for the flux calculations can be found in Supplementary Material.

Organic Matter
The concentration of organic carbon in sediment and in algal mats was determined by the loss on ignition method (LOI) according to Hoogsteen et al. (2015). About 3-5 g of surface sediment (top 2 cm) or floating algal biomass covering approximately 20-30 cm 2 was oven-dried at 105 • C and subsequently combusted at 450 • C. Water and organic matter content were determined by difference and expressed in weight %. The analytical uncertainty of the method is about 1-2%.

Hierarchical Clustering Analysis
The parameters CH 4 flux, water and air temperature, salinity, wind velocity, sediment organic content, and algal biomass carbon were used in a hierarchical cluster analysis (Wilks, 2011) to determine the relative influence of the measured environmental forcing factors on CH 4 fluxes for either the gas transfer model or the chamber method. The different forcing parameters have different units and value ranges. To ensure that parameters had equal weight in the cluster analysis (unweighted), the value range for each parameter was scaled to a value between 0 and 1, where 0 represents the minimum value of that parameter and 1 the maximum value. For flux and sediment organic content parameters, a log transformation was applied prior to this normalization to reduce skewness. All parameters were applied to each flux result and each parameter was given equal weight in the analysis. The parameter water depth was categorized into three depth ranges. The first range was up to 0.5 m, the next range was between 0.5 and 1.0 m and the last range was between 1.0 and 2.0 m, which were then represented by the numbers 2, 1, and 0, respectively. If a parameter had identical values for all observations, the chosen value was 0.5 to represent the middle range. This was used in some exceptional cases in an analyzed sampling period with few values for a particular parameter. The hierarchical cluster analysis used Euclidian distance to calculate distance between observations and a bottom-up method to determine which clusters were closest to each other at each stage, using unweighted pair group method with arithmetic mean (UPGMA). Two types of clustering categories were used -category one used grouping of stations based on flux and environmental forcing factors and category two used grouping of flux and environmental forcing factors across all stations. The hierarchical cluster analysis was performed to the point where only a single cluster remained for each analyzed data set. Dendrograms were used for graphical representation of the cluster analysis results, with hierarchical clustering distance representing the calculated distance between two clusters, using UPGMA at the junction point for a pairwise grouping of clusters. Correlation between different parameters are shown by early clustering (low distance) if there was a strong correlation between parameters, while weak correlation make the clusters late (high distance) in different groups. Dendrograms only show proximity (hierarchical clustering distance) between the clusters that are closest to each other at each clustering stage -not between arbitrary data points. Nodes in a dendrograms can be rotated at its junction points without changing the topology and meaning.

Pairwise Correlation Coefficients (Kendall)
Given the limitations of dendrograms in terms of indicating correlations between all of the different parameters, we extended the analysis to include pairwise correlation calculation (Felipe-Lucia et al., 2020) of environmental forcing factors and flux. Pairwise correlations between CH 4 flux gas transfer model and environmental forcing factors were calculated for all stations for all sampling periods, using Kendall rank correlation coefficient. Possible value range is between −1 and 1. Using Cohen's standard (Sawilowsky, 2009;Chen et al., 2010) for evaluation of the relationships, a weak association is between 0.10 and 0.29, medium association between 0.30 and 0.49 and strong association for 0.50-1.00. The Kendall rank correlation coefficient was chosen as it does not impose any requirements on distribution of analyzed parameter values, other than they can be ranked and provided more robust values than a Spearman rank correlation coefficient or the Pearson correlation coefficient (Croux and Dehon, 2010), since these required that values adhered to normal distribution, a linear relationship between variables, and equal distribution along a regression line. This was not true for all parameters and hence Pearson was not used.

Multiple Factor Analysis
Our data consists primarily of quantitative data, but also qualitative data such as water depth categories. We applied multiple factor analysis (MFA) in order to combine all these data series in the analysis and to allow for grouping of different variables to compress and simplify the data set and explore the structure of the observations. MFA (de Tayrac et al., 2009;Abdi et al., 2013) is a generalization of principal component analysis (Abdi and Williams, 2010;Kherif and Latypova, 2019), which allows grouping of data sets in terms of weight and also can combine quantitative and qualitative data in the analysis. We used this method to get a better understanding of the variances of fluxes and environmental forcing factors and their correlations. Multiple factor analysis calculates new variables referred to as principal components, which are calculated as linear combinations of the original variables. The contribution to the total variance for the principal components was investigated for all data for chamber method, and gas transfer model, as well as different groupings per sampling period and habitat. The variance for all different combinations was explained to 68-96% with the first four principal components, which we considered to explain the variance good enough. Thus our analysis focused on these four components. The contribution of the original variables to the different principal components was then used to study correlation between these variables. In our analysis, we focused on the two to three principal component dimensions. As a graphical depiction of the analysis, the variables are represented by the length and direction of an arrow in a plot of the principal components. Arrows close to each other in the same quadrant are correlated in both dimensions. If variables have similar values on one axis they are correlated for that principal component. Variables in opposing quadrants are negatively correlated in one or two dimensions.

Variability in Dissolved Methane Concentrations
The concentration of dissolved methane differed widely between the studied seasons (spring-summer and fall) and the four investigated habitats (Figures 3A-E). Considering all sampling periods, methane concentrations ranged between 26 and 6596 nmol/L in habitat A, between 42 and 1167 nmol/L in habitat B, 47 and 1201 nmol/L in habitat C, and 69 and 874 nmol/L in habitat D. All methane samples were oversaturated between 908% and 83,865% relative to the dissolved equilibrium concentrations of 2.7-3.6 nM for salinities between 5.4 and 6.7 and temperatures between 8.5 and 28.6 • C at the respective ambient air pressures for the different sampling periods. A Wilcoxon-Mann-Whitney test of all habitats showed that the concentration ranges were significantly different between the habitats except C and D (Supplementary Table 1). Habitat A had the largest proportion of high CH 4 concentrations, likely due to the location of the sampling stations in protected bays, which had significantly more inshore and floating vegetation compared to the other habitats. Average concentrations were calculated by excluding major outliers, thus defining a valid concentration range as where IQR is the interquartile range, defined as IQR = (75th percentile − 25th percentile). The averaged dissolved CH 4 concentrations during all sampling periods in all habitats showed a clear seasonal trend and increased from 210 nmol/L (n = 31) in June through the summer to highest concentrations of 329 nmol/L (n = 89) in July, from where they decreased in the fall to the lowest seasonal concentrations of 75 nmol/L (n = 22) in October (Supplementary Table 2). Median values were in almost all cases lower than average values, indicating that the concentration distribution was skewed toward the high end.

CH 4 Fluxes
Methane fluxes calculated with the gas transfer model ranged from 0.6 to 8.3 mg m −2 d −1 in habitat A, 1.5-5.6 mg m −2 d −1 in habitat B, 0.4-6.6 mg m −2 d −1 in habitat C, and 0.6-6.6 mg m −2 d −1 in habitat D ( Figure 4A). Monthly averaged methane fluxes were calculated by excluding outliers in analogy to the procedure used for dissolved methane. There was a good first order relationship between methane fluxes calculated with the gas transfer model and gas transfer velocity k for all habitats during all whole sampling periods, but the data show significant variability of the fluxes for a given gas transfer velocity (Figure 5). This variability is best explained by high variability in methane concentrations and their influence on the calculated flux. It cannot be excluded that some of the variability is due to the presence of microbubbles in some of the samples. High fluxes in June and August were observed during periods with high wind speeds when rising gas bubbles were observed at the water surface, indicating destabilization of trapped shallow gas in sediments of the near-shore area.
Fluxes calculated with the chamber method in habitat A were in the range 0.6-162.6 mg m −2 d −1 , 0.3-22.9 mg m −2 d −1 in habitat B, 0.4-10.4 mg m −2 d −1 in habitat C, and 0.7-2.0 mg m −2 d −1 in habitat D ( Figure 4B). Figures 6A,B show CH 4 concentrations and fluxes for those sites, at which chambers were deployed with and without bubble shields covering the whole sampling period. From June to August, unshielded fluxes in habitats A, B, and C were higher by factor 5-10 compared to shielded fluxes indicating that bubbles accounted for the majority of the total flux to the atmosphere at all stations. In October, there were no significant difference between fluxes from chambers with and without bubble shields indicating that diffusive fluxes dominated the transport to the atmosphere (Additional variance and mean flux data in Supplementary Table 3). The Wilcoxon-Mann-Whitney test (Fay and Proschan, 2010)

Habitat Categorization of Methane Fluxes and Their Environmental Forcing
Water temperatures less than 15 • C generally showed fluxes lower than 5 mg m −2 d −1 . Wind velocities varied between 1 and 8 m/s during all sampling periods, but there were large variations between sampling periods with low winds in July and October and during periods with higher winds in June, August, and September. There was an expected significant positive relationship between methane flux and wind velocity due to the dependency of k on wind velocity, and a slightly negative relationship between flux and salinity for fluxes calculated with the gas transfer model. There was no correlation between flux and wind velocity for the floating chamber method, likely due to contribution from ebullition for CH 4 , but the correlation between CH 4 flux and water and air temperature was weak, but positive (Kendall's τ < 0.29). Algal biomass and sediment organic contents also showed a significant positive correlation with fluxes, when each habitat was considered individually, but the overall correlation, when all habitats were combined, was weak. The biomass of algae including floating algal mats varied seasonally with highest accumulations in August and a virtual absence of floating algal mats in October. Accumulation in the inshore water also depended on the prevailing wind direction relative to the shoreline orientation, where south to southwesterly winds pushed algal mats into the shore region. Both algal biomass carbon and sediment organic content have a similar correlation with flux for the chamber method and gas transfer model.

Flux Variation in Different Habitats
Flux data from our different sampling stations and environmental parameters were used to test the statistical significance of our station classification scheme into four different habitats. This classification was then used to identify relationships between flux and environmental forcing factors in different habitats and their habitat-specific and seasonal variation. Objective hierarchical statistical clustering was used to determine whether the chosen environmental forcing factors and gas transfer mode-derived methane fluxes discriminated between the four habitat types ( Figure 7A). The analysis showed that the statistics-based grouping reproduced the a priori choice of different habitats for the different stations relatively well. For the fluxes determined with the gas transfer model, habitats A and B group closely together because they have similarities in several environmental parameters. It should be noted that stations in area 5 only had two observations for two sampling periods, compared to stations from the other areas that were sampled 4-5 times, which may have affected the results of the cluster analysis. Habitat C and D are grouped together in the cluster analysis because they have more similarity with each other than the other two habitats. Habitat C is the most distant because it has deeper water than the other habitats. Hierarchical clustering was also performed with the fluxes determined with the floating chamber method (Figure 7B). These results differed somewhat from the gas transfer modelbased results because the calculated fluxes represent the integral of a 12-72-h time interval, and the values for water and air temperature, wind velocity, and salinity were taken at the end of the flux calculation interval, whereas the values used for the gas transfer model were taken at the beginning of each flux experiment. In addition, only fluxes from chambers without bubble shield were used in this analysis. The chamber fluxes include both the diffusive and ebullitive flux, whereas the gas transfer model, at least nominally, only accounts for the diffusive flux. The hierarchical cluster distance of stations according to habitats based on the floating chamber method was larger than for the gas transfer model, likely reflecting the temporal asynchroneity between the flux calculation and some of the environmental parameters, foremost wind velocity and air temperature, suggesting that the analysis for the gas transfer model was statistically more reliable than the floating chamberbased cluster analysis.
The cluster analysis for the four different habitats indicated early clustering between fluxes, organic content of floating algal biomass, and sediment organic content on timescale of a few days and between fluxes, wind velocity, and water temperature on timescale covering the whole sampling period of this study. A relatively late clustering was found between flux, salinity, and water depth, likely because all our sampling areas had about the same water depth and only small variations in salinity (Figures 8A,B). Both dendrograms in Figures 8A,B form clusters with flux, algal biomass, sediment organic content, salinity, and wind velocity, although with slightly different ordering of individual elements. The main difference was a tighter clustering between flux and sediment organic content instead of algal mats for the chamber method, which may indicate the added contribution of sediment gas ebullition. Water and air temperature are expected to have a more significant impact across longer time periods (seasonal time ranges). When the hierarchical cluster analysis was broken down for each month in the different habitats, the same pattern resulted for habitats A, C, and D with exception of the months of July and October, but the clustering results were different for Habitat B, likely because of the small number of observations for this habitat.
In order to further evaluate the strength of the correlation between individual parameters and methane fluxes, we used Kendall rank correlation coefficient for each pair of flux and environmental forcing factors ( Table 2). For the gas transfer model, the strongest association was found between flux and wind velocity. For the chamber method, sediment organic content and algal biomass carbon provided the strongest associations with methane flux. Overall, with the exception of salinity, most correlations were positive. The associations between parameters other than flux were similar for both gas transfer model and chamber method, with strong positive associations between air and water temperature and medium negative associations for both wind velocity and temperatures (air and water) with salinity. Most other associations were weaker. When the analysis was broken down into individual habitats focusing on correlations with methane flux, wind velocity still had the strongest correlation in the gas transfer model, in particular for habitats B, C, and D, and a medium correlation for habitat A. For the chamber method, methane flux showed medium correlations in individual habitats with air and water temperature and a weak correlation between methane flux and air and water temperature when all habitats were pooled together.
Variations and correlation between flux and environmental forcing factors in the different habitats were further analyzed using MFA by reducing the multiple environmental variables into smaller sets of dimensions. Figures 9A,B shows the quantitative variable correlation in all habitats and sampling periods for the two primary principal components (dimensions) of the MFA for fluxes derived with the gas transfer model and chamber method. Altogether, the MFA confirmed the results of the pairwise correlation and the hierarchical cluster analysis indicating wind velocity as the major driving factor for the CH 4 flux in the gas transfer model and sediment organic content and algal biomass carbon as the main drivers for the CH 4 flux in the chamber method. There was no association between flux and wind velocity in the chamber method and only a weak negative correlation between flux and salinity. The secondary factor for both methods were air and water temperature. The MFA analysis showed variable association patterns between flux and environmental forcing factors when conducted for individual habitats, which is probably due to the different number of observations for some of the habitats; habitats B and D have less than one third of the number of observations compared to habitats A and B. The contribution from the two primary principal components varied between 41 and 71%. The inclusion of a third dimension increased the total contribution from the three primary principal components to between 56 and 93%. For the gas transfer model in habitats B, C, and D the methane flux showed the strongest correlation with wind velocity in dimension 3, while in habitat A the strongest correlation existed between CH 4 flux and temperature. In the chamber method a strong correlation was found between CH 4 flux and sediment organic content in habitat A, while habitat B and D had stronger correlations between flux and wind velocity. In habitat C, dimension 3 did not enhance the correlation indicated in dimensions 1 and 2. The flux correlation deviated during July sampling period for gas transfer model and in June sampling period for chamber method. The supplementary material contains additional quantitative variable correlation for the two primary principal components (dimensions) of the MFA per habitat (Supplementary Figure 1) and per sampling period (Supplementary Figure 2).

Method Comparison and Implications for Environmental Forcings
An important observation of this study has been that the statistical associations of the environmental forcing factors varied with the flux method. Other studies that compared a gas transfer model-derived with chamber-based flux measurements including day-and nighttime measurements found higher methane emissions with the chamber method than the gas transfer model (Erkkilä et al., 2018). Our measurements with the chamber method generally occurred over a 12-24-h period several times over a week thereby integrating diurnal flux variations over this time period. This is not the case for the gas transfer model, which represents daytime methane measurements and air temperature and wind velocity measurements that were retrieved from measurements from nearby weather stations and do not consider local and short-term temporal variations.
The underlying assumption of the gas transfer model is that transport through the sea surface boundary layer takes place via diffusion of methane in one phase. It cannot be said with certainty  that the high concentrations of the water samples represent truly dissolved methane, a precondition for the application of the gas transfer model. It is possible that some of the water samples contained microbubbles that were not detected visually and interpreted as dissolved methane due to high concentration of methane in water. This raises the question whether the gas transfer model represented the best suited equation for all sampling periods. In addition, as expected, the gas transfer model fluxes were most closely related to variations in wind velocity, however, the strong correlation does not reflect the overall importance of wind as regulating factor considering the dominance of bubble transport in the shallow inshore areas. The chamber-based fluxes make no assumptions about the transport mechanisms diffusion, ebullition and plant-mediated transport, but here shielding effects of the chamber regarding wind can have affected the forcing by wind.
In addition, the seasonally varying presence of algal biofilms at the water surface may have affected the gas transfer velocity and thus the flux of methane when using a gas transfer model. In July and August, thick algal mats covered the water surface of most of the sampling stations of habitats A and B. The algal mat thickness was between 0.5 and 2 cm. Goldman et al. (1988) noted that natural surfactant films in algal mats can affect the oxygen exchange coefficient and reduced gas transfer velocity to 63-85% of a control value for gas transfer for distilled water. This is due to changes in the oxygen exchange coefficient that is sensitive to small changes in carbohydrate concentration, since they affect the viscosity and thereby the Schmidt number (Frew et al., 1990). Reductions of k by up to 50% were reported in near coastal estuarine areas (Goldman et al., 1988) contrary to open ocean water (Nightingale et al., 2000a). Cyanobacterial biofilms can have high concentrations of organic calcium and magnesium complexes (Arp et al., 2001;Dupraz et al., 2009). These may also cause a decrease in k by up to a factor 2 (Frew et al., 1990). Due to the variability in the presence and thickness of the biofilms at different stations and during the different sampling periods, we would only have been able to make qualitative adjustments to the gas exchange coefficient, but we note that the gas transfer flux calculations in our study at certain times and in areas with thick biofilms may have overestimated the actual gas transfer model flux. A surface algal biofilm can also increase the retention of rising bubbles in the viscous biofilm and may therefore reduce the bubble flux released to the overlying chamber. This would imply that methane fluxes determined with the chamber method in areas of thick algal mats may actually be lowered due to the biofilm presence.
In summary, the chamber-based flux method appeared to be superior in capturing variations in flux as a function of the different environmental forcing factors. The chamber method also provided a more stable environmental forcing pattern across the different habitats during the whole season. We therefore conclude that this approach may be more useful than the gas transfer model for spatial scaling of flux measurements provided that a good categorization of different habitats and their prevalence can be performed. Practically, however, the chamber method presents greater technical challenges in comparison to the gas transfer model, with which measurements generally are easier to obtain. Thus, in practice a combination of methods will be needed and also compared to eddy correlationbased flux methods.

Magnitude of Methane Fluxes
The monthly averaged fluxes of our study were in the same order of magnitude as fjord-type environments and coastal lagoons, whereas fluxes in estuarine waters tend to be up to an order of magnitude lower (Borges and Abril, 2012). Our highest fluxes also compare well with those reported by Heyer and Berger (2000) for coastal wetlands of the southwestern Baltic Sea. Shallow inshore areas generally have well-mixed waters compared to the more stratified offshore areas in deeper water. Compared to fluxes in the offshore Baltic Sea, the range of the spatial and temporal variability of fluxes in our study was almost two orders of magnitude higher. The lowest fluxes in our study were 0.1 mg m −2 d −1 , which is similar to the lowest fluxes reported in Gülzow et al. (2013), but the highest values in our study were 162.6 mg m −2 d −1 , compared to the highest fluxes of 1.5 mg m −2 d −1 reported by Gülzow et al. (2013). A similar distinct onshore-offshore difference was observed by Borges et al. (2016), who also observed the highest inshore fluxes in the summer season and attributed these to a positive feedback from marine seepage in near-shore shallow areas with gassy sediment. However, the combined effects of seasonal variability and onshore-offshore patterns were different from our study, since Borges et al. (2016) found the lowest fluxes in the offshore areas in late spring and summer, and the highest fluxes in late winter, whereas our study observed the highest fluxes in the summer and the lowest in autumn. Borges et al. (2016) observed similar trends in relation to salinity, where higher fluxes correlated with lower salinity, although the range of salinity was an order of magnitude larger than in our study and also related to greater range in water depth. These comparisons indicate that the environmental effects of the forcing parameters are different for inshore and offshore areas, which make simple extrapolations of our study to greater water depths difficult.

Environmental Forcing
Three main paths have been identified for atmospheric methane emissions from inshore shallow areas: ebullition, diffusion, and plant mediated transport (Jeffrey et al., 2019). In shallow, littoral waters, most methane is produced in the underlying sediment by methanogenic archaea and methylotrophs, from where it diffuses to the atmosphere by diffusion or is emitted as gas bubbles (Crill and Martens, 1986;Chanton et al., 1989;Madigan et al., 2003). The formation of methane bubbles in sediment or floating algal mats reflects the prevailing physicochemical conditions that determine its solubility at a given hydrostatic pressure and temperature, the rate of methane production and oxidation (Martens et al., 1998;Wever et al., 2006;Mogollón et al., 2011). Methane bubble fluxes between 1.6 and 277 mg m −2 d −1 have been reported in various coastal settings (Chanton et al., 1989;Chanton and Dacey, 1991;Leifer and Patro, 2002;Borges et al., 2011), which agrees well with measurements in our study.
Transport of methane through plants ( Van der Nat and Middelburg, 2000) and the shallow water depth favor direct transfer of gas bubbles to the atmosphere minimizing oxidation (Van der Nat and Middelburg, 2000). In certain types of habitats, plant-mediated methane transport may account for the majority of methane emissions of these three paths annually (Jeffrey et al., 2019).
In addition, sudden air pressure changes and periods with high winds also enhance the release of methane (Lohrberg et al., 2020). It has been suggested that at least 45% of the initial methane gas content in sediment may escape to the atmosphere during a storm event (McGinnis et al., 2006;Lohrberg et al., 2020). Our study was conducted the sampling over a period of 5 months, with samples taken during a period of a few days each calendar month. Although the weather was relatively stable during this period, there were periods with high winds and heavy rain, which occurred both between and on sampling days during June and August, when it was not practically possible to obtain sampling data. Weather data for the region during our sampling season indicate high-wind events during eight occasions with wind velocities higher than 10 m/s. The effects of these high-wind events in the middle of July and end of August/beginning of September can be noted in Figure 4 and show higher fluxes than those measured during the other sampling months. It is also important to note that these events improve the correlations between methane fluxes and wind velocities in the statistical analyses.
Our data show that ebullition was related to habitat type and season. Wind-and wave-protected embayments with organicrich fine-grained sediment covered by dense mats of floating algae (habitats A and B) had the highest ebullition rates, but this mode of transport only dominated during the late spring and summer, and ceased in the fall. Both temperature and salinity affect the solubility of methane, but the observed salinity in our study (5.4-6.7 ) had a negligible effect on the solubility itself. Temperatures measured during the sampling days in each period were generally considered representative of each sampling period representing a whole month with the exception of July 2019, which had an unusually large temperature range from 9 • C in the beginning to 33 • C at the end of the month. The full consequences of this large temperature increase are insufficiently captured with the discrete week-long measurements of our study. Taking these considerations into account, the most significant environmental forcing factor on seasonal timescales were temperature and organic content. This is in line with other studies in the Baltic Sea, which found that the primary ecological driver for methane emission was the seasonal variation of the sediment organic matter content, while temperature controlled daytime and nighttime as well as longterm seasonal long-term variability (Heyer and Berger, 2000;Gülzow et al., 2013).
Habitats A and B with their dense reed beds and floating algal cover in only 0.5 m water depth had the highest fluxes in our study emphasizing the role of plant-mediated methane transport as additional contribution. This was also concluded by Van der Nat and Middelburg (2000), who found that transport through reed and bulrush represented over 85% of the total net methane emission from such habitats. In situ production in floating algal mats contributes additional methane, but has generally been considered secondary compared to the sediment source. However, the high fluxes in habitat C, where boulders covered the seafloor and methanogenic sediment was absent, indicate that in certain habitats floating algal mats are the sole source of the methane flux to the atmosphere during the summer.
Hierarchical clustering analysis helped to validate the categorization of our sampling stations into the four different habitats, but the clustering results showed different group constellations of parameters depending on whether methane fluxes were calculated with the gas transfer model or using floating chambers. There were also different forcing constellations depending on habitat. These results suggest that the calculated hierarchy of forcing parameters for methane fluxes was dependent on the field method and habitat. In the case of the chamber method, this is due to the fact that the methane flux was more strongly forced by methane bubble transport than by wind velocity. Since the bubblemediated flux far exceeded the diffusive flux for all months except October, we conclude that the environmental forcing constellations associated with the chamber measurements are more realistic for habitats A and B, whereas for the vegetationpoor and open areas habitats C and D, where wind forcing plays the dominant role, the forcing constellations of either gas transfer model or chamber-based measurements yield realistic constellations.
The basic forcing constellations derived from the hierarchical clustering were confirmed by the pairwise correlation analysis using Kendall correlation coefficients for each of the flux measurement methods. As expected, the Kendall correlation coefficients provided very good correlation coefficients between wind velocity and methane flux for gas transfer model-derived data, but very weak correlations of the same parameters for the chamber-based measurements. However, since wind velocity only had a minor influence on the actual fluxes in habitats A and B, where vegetation density and organic richness played a much stronger role, these high correlation coefficients can be misleading for understanding the environmental multiparameter forcing. These problems were overcome with the MFA, which provided a measure of the correlation quality and the relative weight of different forcing parameters in the different habitats. The primary control of vegetation and sediment organic content followed by air/water temperature as secondary and, finally, wind velocity as tertiary control provided a realistic assessment of the hierarchy of environmental forcing factors in habitat A. The weighting shifted in the more open habitat B, with the primary control parameters to be the same as in habitat A, but with wind velocity and air/water temperature switching positions as secondary and tertiary controls. Finally, the statistical selection of wind velocity as the primary control followed by water depth as secondary, and air/water temperature as tertiary control, for the barer and deeper habitats C and D is a realistic reflection of the environmental drivers in these more open and rocky bottom habitats.

Spatial Extrapolation
Multiple factor analysis provided the first steps toward a scaling approach to predict the methane fluxes in areas represented by our four habitats, for which no direct measurements were taken. At present, the flux database per habitat is too limited to develop algorithms based on algal biomass density, sediment organic content, temperature, and wind velocity for these shallow-water habitats. As an alternative, we chose to use the monthly averages  for our habitats as guideline values and partitioned the coastline of the island of Askö into the habitat groups by calculating areas extending 10 m from the shore using aerial photographs. This made it possible to extrapolate the methane emissions of our investigated areas with known habitats to the coast line of the entire island ( Table 3). The total length of the Askö coastline was estimated to be 34 km, and the total calculated coastal area within 10 m of the shore is 0.34 km 2 . Habitat A comprises 0.08 km 2 , habitat B 0.05 km 2 , habitat C 0.13 km 2 , and habitat D 0.09 km 2 , representing 24, 14, 37, and 25% of the total area (Figure 10), and the presence of habitat categories that were not included in our study such as a sandy beach, which covers 500 m of the northern shoreline of the island. Here, habitat D was selected as the closest analog. The estimated emission of methane from the inshore habitat of the island of Askö varied from 8.8 kg/day in August to 0.17 kg/day in October for the chamber-based flux. Notably, habitat A-and B-type areas contributed 89% to the total flux in August, but no more than 50% in October, when ebullition had ceased and the emissions were more evenly distributed across the habitats. Extrapolation of the emission over a full year assuming similar emissions for October, the winter, and early spring months yielded a very high estimate of 594 kg per year for the whole coast of the island. Based on the seasonal trend in emissions, we estimate that about 90% of the emission flux occurred during July and August and dominated by habitats A (72%) and B (7%). These daily emissions only provide firstorder estimates of emissions and probably overestimate the actual emission. Uncertainties of our approach arise from the resolution of published aerial photographs and small-scale variability of fluxes within a defined habitat area, the local inshore wind field on the seaward and landward side of the island, and the variable shoreline shape and geographic orientation. The data require comparisons from ongoing, supplemental, continuous eddy correlation, and continuous water sampling measurement campaigns.

CONCLUSION
The current study demonstrates that a habitat-based approach to upscale local, time-limited sea-air exchange measurements of methane fluxes can only be successful, if a weighting of the multiple primary drivers for methane emissions in shallow water coastal habitats is carried out, because no single environmental driver dominates fluxes. In addition, the forcing mechanisms between environmental drivers and methane flux are not linear, because the onset of ebullition is characterized by the passing of temperature and primary productivity thresholds above which the emission type switches from diffusion to ebullition thereby increasing the emission rates more than tenfold. It is likely that such 'switches' in the dominant transport process do not occur synchronously along shorelines. In rocky bottom habitats that lack sufficient organic richness or surface algal accumulations, transport mode changes likely do not occur at all and the exchanges are likely controlled by wind forcing and temperature only. The study also demonstrates that apparent good correlations between environmental forcing variables such as wind velocity and temperature only provide good correlations for boundary layer models, but not when ebullition is explicitly accounted for in the measurement. There is a need for complementary approaches to determine the exchange fluxes, because the spatial and temporal extrapolations have many uncertainties. One possible approach may be to use flux measurements carried with eddy correlation towers in the different habitats. This would overcome the problem of discontinuous measurements and could account for the habitatspecific variability in fluxes within an eddy footprint.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) are: Bolin Centre Database, https://doi.org/ 10.17043/lundevall-zara-2021-asko-methane-1.

AUTHOR CONTRIBUTIONS
ML-Z: main researcher, planning and execution of fieldwork, collection of samples, laboratory analysis of samples, data analysis, and visualization of data. EL-Z: software development for statistical analysis and visualization of data. VB: supervisor and scientific review and editing. All authors contributed to the article and approved the submitted version.