Drivers of the spatiotemporal distribution of dissolved nitrous oxide and air–sea exchange in a coastal Mediterranean area

Among the well-known greenhouse gases (GHGs), nitrous oxide (N₂O) is the third most impactful, possessing a global warming potential approximately 300 times greater than that of carbon dioxide (CO₂) over a century. The distribution of N₂O in aquatic environments exhibits notable spatial and temporal variations, and emissions remain inadequately constrained and underrepresented in global N₂O emission inventories, particularly for coastal zones. This study focuses on N₂O levels and air–sea fluxes in the coastal waters of the Balearic Islands Archipelago in the Western Mediterranean basin. Data were gathered between 2018 and 2023 at three coastal monitoring stations: two on the densely populated island of Mallorca and the third in the well-preserved National Park of the Cabrera Archipelago. Seawater N₂O concentrations varied from 6.5 to 9.9 nmol L⁻¹, with no significant differences being detected across the sites. When these sink–source strengths are integrated on an annual basis, the Balearic Sea is close to equilibrium with atmospheric N₂O, resulting in a neutral atmosphere–ocean exchange (0.1 ± 0.2 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). A consistent seasonal pattern was noted during the study period. Machine learning analysis indicated that seawater temperature was the primary factor influencing N₂O concentrations, with lesser contributions from chlorophyll levels and salinity.

1 Introduction

Nitrous oxide (N₂O) is a potent greenhouse gas (GHG) with 300 times higher warming potential per mole than carbon dioxide (CO₂) on a 100-year timescale (Solomon et al., 2007). Atmospheric N₂O levels have risen by more than 18 % since preindustrial times and increased by 332 ppb between 2011 and 2019 (Masson-Delmotte et al., 2021). The ocean N₂O budget is highly sensitive to climate change and influences the climate system significantly (Ravishankara et al., 2009). Variations in temperature, ocean circulation, and biological activity can alter the production and release of N₂O from the oceans. Consequently, the concentration of N₂O in the atmosphere influences global warming and climate patterns, creating a feedback loop between oceanic processes and the climate system. Marine N₂O sources account for one-third of the natural emissions to the atmosphere, yielding a net source of 3.5 (2.5–4.7) Tg N yr⁻¹, excluding coastal contributions (Tian et al., 2024). N₂O is primarily produced through nitrification and denitrification pathways (Freing e  al., 2012. Nitrification is an aerobic process in which ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$) is oxidized to nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$), producing N₂O as a byproduct, favoured by low-oxygen conditions. In contrast, denitrification occurs under suboxic to anoxic conditions and involves the stepwise reduction of nitrate to dinitrogen gas (N₂), also releasing N₂O as an intermediate. These microbial processes are particularly sensitive to environmental gradients in oxygen, nutrients, and the availability of organic matter. The nitrogen cycle is one of the most complex regulating factors of primary production, highly dependent on dissolved oxygen concentration and the prevailing redox conditions (Codispoti, 2010). In coastal environments, significant variability exists in the nitrogen cycle, where land-derived nutrient inputs, coastal upwelling events, and complex biogeochemical processes play crucial roles (Doney, 2010), contributing significantly to the spatiotemporal variability of the N₂O concentrations (Nevison et al., 1995). In estuarine and coastal waters, the effects of climate change may be more pronounced, such as in ocean acidification (OA, Carstensen and Duarte, 2019), which can enhance the generation of N₂O (Wan et al., 2023; Zhou et al., 2023). Despite the necessity of a better understanding of atmospheric and oceanic inventories of non-CO₂ GHGs to provide realistic and accurate models for future scenarios, there are limited open-ocean and coastal monitoring time series networks compared to for CO₂ (Bakker et al., 2014; de la Paz et al., 2015; Farías et al., 2007; Ma et al., 2019; Wilson et al., 2017).

The Mediterranean Sea is a semi-enclosed basin surrounded by highly sensitive coastal zones that are particularly vulnerable to human activities. Factors such as high population density, widespread urbanization, and intensive agriculture have escalated the risks of pollution and habitat degradation in the region (UNEP/MAP, 2012). Due to its distinct biogeochemical and hydrodynamic features, this basin has been identified as a “hotspot” for climate change research (Giorgi, 2006). The impacts of global warming and extreme weather events are expected to be more severe in the Mediterranean compared to in other oceanic regions (Giorgi, 2006; Giorgi and Lionello, 2008; Masson-Delmotte et al., 2021). Despite representing just 0.82 % of the global ocean surface, the Mediterranean hosts 4 %–18 % of the world's marine biodiversity, including numerous endemic species (Bianchi and Morri, 2000; Mouillot et al., 2011). Rising temperatures and OA threaten the biodiversity of the region (Micheli et al., 2013). Additionally, anthropogenic pressures along the Mediterranean coast have intensified due to rapid population growth and economic activities. In the Western Mediterranean, high tourism and coastal development levels have left only a small fraction of the coastline in a natural state, with even fewer areas under protection (EEA, 1999).

The Balearic Islands Archipelago, located in the Western Mediterranean, comprises the islands of Mallorca, Menorca, Ibiza, and Formentera, with a combined coastline of 1723 km. Renowned as a major European tourist destination, tourism accounts for approximately 45 % of the total gross domestic product of the archipelago. Visitor numbers have surged dramatically over the past century, reaching nearly 18 million in 2023 (Institut d'Estadística de les Illes Balears, Spain) compared to a resident population of around 1.2 million. Coastal ecosystems in the Balearic Islands are vital for the local economy. Meadows of the endemic seagrass Posidonia oceanica extend across depths of up to 45 m in the Balearic Sea (del Valle Villalonga et al., 2023), providing critical ecosystem services such as carbon sequestration (Duarte et al., 2005), oxygen production (Hendriks et al., 2022), biodiversity support, coastal erosion prevention, sediment stabilization, and water transparency (Barbier et al., 2011). However, these ecosystems face increasing pressure from recreational activities and other anthropogenic impacts.

Given the increasing threats to these ecosystems, understanding the relationship between anthropogenic pressures and GHG emissions, particularly N₂O, has become urgent. The absence of long-term N₂O datasets for the Mediterranean Sea, along with the uncertainties surrounding current emission estimates, underscores the importance of assessing coastal areas with varying human impacts. These evaluations are essential for improving coastal N₂O emission estimates and refining global ocean N₂O budgets. In this study, we evaluate the spatial and temporal N₂O concentrations in surface waters and the air–sea exchange in the coastal area of the Balearic Islands Archipelago and estimate the potential drivers of the observed variability in N₂O concentrations. We focused on three different sites in the coastal zone: a highly impacted site and a medium-impacted site, both located near the island of Mallorca, and a pristine site in the Cabrera National Park Archipelago.

Figure 1 (a) Map of the stations location in the Western Mediterranean basin and (b) detailed locations of the Bay of Palma (PB; blue dot), Cabrera National Park (CA, red dot), and Cape Ses Salines (CS; green dot) study sites and Palma Airport wind station (black star). Dashed lines represent bathymetry levels. Maps were developed with the Python software version 3.12.3.

2 Methods

2.1 Study area

We collected physicochemical and biogeochemical parameters from three stations in the Balearic Sea within the Western Mediterranean basin (Fig. 1a); these are part of the Balearic Ocean Acidification Time Series (BOATS; Flecha et al., 2022). The data presented here were collected from 2018 to 2023; however, the data availability varies among the different stations.

Two sampling sites are fixed monitoring stations equipped with autonomous sensors. The first station, established in December 2018, is located in the Bay of Palma (PB: 39.492848° N, 2.700405° E; at a depth of approximately 30 m; Fig. 1b) and is part of the fixed monitoring network of the Balearic Islands Coastal Observing and Forecasting System (Tintoré et al., 2019; 2013-SOCIB; https://www.socib.es/, last access: December 2024). Temperature (°C) and salinity (PSU) data were obtained from the SOCIB buoy sensors; see Tintoré and Casas (2022) for details about the sensors. Additionally, a MiniDot sensor (PME Inc.©) recorded dissolved oxygen (DO) hourly, with a manufacturer accuracy of ± 5 %. The sensor location and water sampling depth were set at 1 m.

The second fixed monitoring station is located in the Bay of Santa Maria (Fig. 1b) within the Cabrera Archipelago National Park (CA: 39.151395° N, 2.950823° E; ∼ 8 m depth) and was installed in November 2019. It is situated in an area under governmental protection and is regarded to be a pristine site with no apparent human influence. Temperature (°C), salinity (PSU), and DO data were collected hourly using SBE37-SMP-ODO (Sea-Bird Scientific Electronics^®) and a MiniDot. Both sensors were secured to a mooring line at approximately 4 m depth, and water samples were taken at the same depth. The manufacturer accuracy of measurements was ± 0.002 °C, ± 0.003 mS cm⁻¹, and ± 5 % for temperature, conductivity, and oxygen sensors, respectively.

At PB and CA stations, samples of dissolved N₂O; DO; dissolved organic carbon (DOC); Chlorophyll a (Chl a); and inorganic nutrients, including nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$), nitrite (${\mathrm{NO}}_{\mathrm{2}}^{-}$), phosphate (${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$), silicate (Si(OH)₄), ammonia (NH₄), and total nitrogen (TN) were collected monthly.

The third sampling point is located in the coastal area near the Cape Ses Salines lighthouse (CS: 39.2649° N, 3.0535° E; Fig. 1b). At this site, with a total bottom depth of 2 m, data were collected biweekly from surface water directly off the coast at approximately 0.5 m depth. Sampling began in August 2018. Temperature (°C), salinity (PSU), N₂O, DOC, Chl a, and inorganic nutrients were collected from the same volume of surface water. DO data were gathered using a MiniDot sensor starting in August 2018. Validation of DO sensor data was performed with DO water samples collected from PB and CA stations, as detailed in Agueda-Aramburu et al. (2024).

The three sites differ in terms of their physical setting and benthic characteristics. PB is located in an open bay with moderate anthropogenic influence and limited seagrass coverage within the sampling area. CA, on the other hand, is located in a semi-enclosed bay within a protected marine park and is surrounded by extensive Posidonia oceanica meadows, which are known to influence nutrient and gas dynamics. CS is a shallow nearshore site with intermittent seagrass patches and high exposure to coastal land–sea interactions. These contrasts are relevant to the interpretation of local biogeochemical processes, including N₂O production and air–sea fluxes.

2.2 Data collection and analysis

2.2.1 Biogeochemical variables

To determine N₂O levels, samples were collected in duplicate using 120 mL serum vials sealed with grey-butyl rubber stoppers and aluminium crimps. After being sealed, the samples were preserved with HgCl2 and stored upside-down until analysis. N₂O concentrations were analysed at the AQUANITROMET laboratory (https://www.iim.csic.es/en/about-iim/organization/aquanitromet-analysis-greenhouse-effect-gases-natural-waters, last access: December 2024) of the Instituto de Investigaciones Marinas (IIM-CSIC, Vigo, Spain), employing a static headspace equilibration technique combined with gas chromatography (GC) equipped with electron capture detection, following the methodology detailed by De la Paz et al. (2015). To create the headspace, 20 mL of nitrogen gas from a Tedlar bag at atmospheric pressure was introduced into the vials while simultaneously extracting the same volume of water sample using a double-needle setup. The vials were shaken and left to equilibrate for at least 12 h in a temperature-controlled environment. For injection into the GC, a brine solution was added through one needle to displace the headspace gas into the GC via the second needle.

The GC, an Agilent 7890 GC, was calibrated using three standard gas mixtures: a NOAA-certified primary standard resembling atmospheric air composition and two additional N₂O-in-N₂ mixtures supplied by Air Liquide (De la Paz et al., 2015). While participating in the first large-scale international inter-laboratory comparison experiment for seawater N₂O measurements (Wilson et al., 2018), an additional certified standard from the Scientific Committee for Oceanographic Research (SCOR) was utilized. The precision of the analysis was determined to be 0.5 %, calculated based on the average coefficient of variation across 400 replicate measurements.

To determine the DO concentration, samples were collected in five borosilicate bottles and then analysed using the Winkler method modified by Benson and Krause (1984), employing potentiometric titration with a Metrohm 808 Titrando. The precision of the DO analysis was estimated to be less than ± 2 µmol kg⁻¹.

Chl a samples were collected in glass bottles, and three replicates were filtered in the laboratory using a Whatman GF/F glass fibre filter. Chl a extraction was performed with 90 % acetone for 24 h under dark conditions, and the samples were measured with a fluorometer (Turner Designs Instrument, Model 7200-00). The fluorometer was calibrated with a pure Chl a standard derived from Anacystis nidulans algae (Sigma Chemical Company) (Knap et al., 1996).

DOC samples were collected in two pre-combusted borosilicate bottles and then filtered through pre-combusted Whatman GF/F glass fibre filters and stored in two pre-combusted borosilicate vials containing 25 µL of orthophosphoric acid (H₃PO₄). The samples were analysed with a Shimadzu TOC-L analyser following the method described by Álvarez-Salgado and Miller (1998), which is based on catalytic oxidation at a high temperature of 680 °C.

Inorganic nutrient samples were collected in four 10 mL Falcom vials at the same time as the DOC sample processing. Concentrations of ${\mathrm{NO}}_{\mathrm{3}}^{-}$, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$, Si(OH)₄, NH₄, and total nitrogen were obtained from the analysis laboratory at the Mediterranean Center for Marine and Environmental Research (CMIMA, Barcelona, Spain) using the Autoanalyser AA3 HR (Seal Analytical, United Kingdom) through continuous flow analysis. The precision, estimated from the coefficient of variation based on replicate analyses of the same water samples (n = 10), ranged from 0.13 % to 0.5 %.

2.2.2 Meteorological and atmospheric data

Wind speed at a height of 10 m was provided by the Agencia Estatal de Meteorología (AEMET) from the Sant Joan Airport station in Palma, Spain (Fig. 1b). Monthly averaged data of atmospheric N₂O molar fraction were obtained from the monitoring station at Lampedusa (LMP), Italy, as part of the NOAA (National Oceanic and Atmospheric Administration) monitoring network (http://www.esrl.noaa.gov/gmd/dv/site/, last access: December 2024; Lan et al. 2024).

2.3 Flux estimation and other calculations

To calculate the water–atmosphere N₂O fluxes ($\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), the following equation was used:

$$\begin{array}{}\text{(1)}& F=k\cdot (C_{\mathrm{w}}-C^{\ast }),\end{array}$$

where k (cm h⁻¹) is the gas transfer velocity, C_w is the concentration of N₂O dissolved in water samples (mol⁻¹), and C^∗ is the gas saturation concentration calculated as the product of the atmospheric fraction of N₂O and the solubility coefficient proposed by Weiss et al. (1980). To compute the fluxes of N₂O, we used the monthly mean atmospheric molar N₂O obtained from NOAA.

The most suitable parameterization for gas transfer (k) in coastal areas with seagrass ecosystems, which are characterized by limited wind fetch, representing the study area in the Balearic Sea, was utilized, the equation for which was described by Dobashi and Ho (2023) as follows:

$$\begin{array}{}\text{(2)}& k=\mathrm{0.143}{U}_{\mathrm{10}}^{\mathrm{2}}.\end{array}$$

Emissions to the atmosphere are indicated by positive values. Additionally, due to the bottom depth differences between stations and to compare with other studies, the gas transfer parameterizations described by Wanninkhof (2014) and Cole and Caraco (1998) were also used to determine water–atmosphere N₂O fluxes, and differences are presented in the obtained annual N₂O fluxes.

Annual N₂O fluxes were calculated by integrating daily flux estimates over time for each station and year using the trapezoidal rule. For each group of samples, the area under the flux–time curve was computed based on the sampling dates, providing an estimate of the total annual air–sea flux. This method accounts for the uneven temporal distribution of data and offers a reliable approach to derive annual fluxes from discrete observations.

Furthermore, saturation levels for N₂O (Sat % N₂O) are expressed as a percentage, determined by the ratio of the observed N₂O concentration to the calculated equilibrium concentration. ΔN₂O represents the difference between the observed N₂O concentrations. The apparent oxygen utilization (AOU) is calculated using the DO values along with the measured solubility of oxygen in seawater (Benson and Krause, 1989).

2.4 Data analysis

To test for differences in N₂O concentrations and fluxes across regions and over temporal scales, simple multifactorial general linear model analyses were conducted using MATLAB version 9.10.0 (R20211, MathWorks Inc.).

Prior to model fitting, we performed a preliminary screening of all available environmental variables to identify those with the highest potential predictive value for N₂O dynamics. We tested the influence of environmental drivers on N₂O levels using the supervised machine learning (ML) method known as gradient-boosting machine (GBM), which is based on decision tree models using a cross-validated boosting (CVB). GBM utilizes a specific implementation of an ensemble method that combines sequenced base weak models to create a stronger one by applying gradient descent to minimize the model loss function. This process included correlation analysis, variance inflation factors, and variable importance metrics derived from preliminary GBM models. In each iteration, a new model is adapted to the residuals of the previously combined model. The combination of all of the base models is achieved by summing their predictions, and, to prevent overestimation, a weighted learning factor is applied. The CVB technique allows for more precise model evaluation by repeatedly using different data sheets to train and test the model. Both GBM and CVB were implemented in Python version 3.12.3 using the XGBoost library (Chen and Guestrin, 2016). Although ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ were included in the initial evaluation, they were excluded from the final model due to either lower relative importance or missing data across sites, which would have compromised model consistency.

The three most influential variables identified (temperature, Chl a, and salinity) were further analysed to explore their individual and combined effects using ordinary least squares regression (OLS) and generalized additive models (GAMs). The correlation between the response of N₂O and the resulting dominant environmental variables (relative importance > 10 %) obtained from the CVB was assessed through OLS using the Statsmodels library (Seabold and Perktold, 2010) in Python version 3.12.3 and by employing a GAM with the pyGAM library (Servén and Brummitt, 2020) in Python version 3.12.3 to address the nonlinear dependency. This multi-step approach enabled robust variable selection and the identification of mechanistic relationships driving N₂O spatiotemporal patterns.

To determine the simple linear correlation between environmental variables and FN₂O, Pearson correlation coefficients and p values were calculated for each station using the SciPy library (Virtanen et al., 2020) in Python version 3.12.3. Furthermore, a multiple linear regression analysis was performed to assess the impact of various predictor variables on FN₂O by applying an OLS.

Table 1 Seasonally averaged values of temperature (°C), salinity (PSU), Chlorophyll a (Chl a; µg L⁻¹), dissolved oxygen (DO; µmol kg⁻¹), apparent oxygen utilization (AOU; µmol kg⁻¹), nitrate (µM), nitrite (µM), silicate (µM), dissolved organic carbon (DOC; µM), N₂O (nanomol L⁻¹), and percentage saturation of N₂O measured in samples collected from PB (Palma Bay), CA (Cabrera National Park), and CS (Cape Ses Salines). Variability is represented as ± standard deviation.

^∗ DO data are only available from August 2022 at the CS station. N₂O saturation percentages differ significantly from 100 %, as evaluated by a t test (p > 0.005).

Figure 2 Time series of (a) daily averaged temperature (°C) for PB (blue dots) and CA (red dots) stations, as well as instantaneous values for the CS site (green dots). (b) Daily averaged salinity (PSU) for the PB (blue dots) and CA (red triangles) stations, along with instantaneous values for the CS site (green squares). (c) Daily averaged wind speed (m s⁻¹) at 10 m height from the Palma station (Spain). Data for PB, CA, and CS were collected at depths of 1, 4, and 0.5 m, respectively. Figures were created using Python software version 3.12.3.

Figure 3 Time series of (a) nitrous oxide (N₂O) concentrations in nmol L⁻¹, (b) N₂O saturation percentage (Sat % N₂O, %), (c) temperature (°C), (d) Chlorophyll a (Chl a) in µg L⁻¹, (e) dissolved oxygen (DO) in µmol kg⁻¹, (f) dissolved organic carbon (DOC) in µM, (g) nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$) in µM, and (h) nitrite (${\mathrm{NO}}_{\mathrm{2}}^{-}$) in µM as measured in the samples collected from PB (blue dots), CA (red triangles), and CS (green squares). Data for PB, CA, and CS were collected at depths of 1, 4, and 0.5 m, respectively. Figures were developed with the Python software version 3.12.3.

3 Results and discussion

3.1 Environmental variables description

The observed seawater temperature patterns during the study period reflect the climatic characteristics of the Mediterranean region, with peak temperatures being recorded in summer (June–August), notably reaching their maximum in August 2022 at 30.3, 29.9, and 29 °C at the PB, CA, and CS sites, respectively (Figs. 2a and 3c and Table 1). Conversely, minimum temperatures recorded during winter (December–February) were 13.7 °C in 2023 for PB, 13.04 °C in 2023 for CA, and 13 °C in 2019 for CS. Over the entire study period, average (± standard deviation) temperature values stood at 20.9 ± 5.1, 20.7 ± 4.9, and 20.4 ± 5.0 °C for PB, CA, and CS, respectively, with no statistical differences between sites (p > 0.05, Table 1), but variability was denoted between months and years (p < 0.005, Table 1).

Surface salinity levels showed significant differences between sites (p < 0.005 with CS values), months (p < 0.005), and years (p < 0.05), maintaining average values of 37.5 ± 0.2 practical salinity units (PSU) for PB, 37.4 ± 0.2 PSU for CA, and 37.7 ± 0.3 PSU for CS. The highest salinity levels were recorded during the summer of 2018 at PB (38.2 PSU), during the winter of 2023 at CA (38.0 PSU), and during both the summer of 2018 and the winter of 2023 at CS (38.2 PSU) (Fig. 2b, Table 1). The differences in the factory precision of the sensors that measure conductivity – and, thus, salinity – could explain the variations observed between the CS, PB, and CA stations in terms of salinity.

Wind speed measurements taken at a height of 10 m displayed a noticeable seasonal trend (Fig. 2c). However, no statistical differences were found between years and months (p > 0.05). The highest average values occurred during spring (March–May), averaging around 3.7 ± 0.4 m s⁻¹, and decreased during winter to 2.9 ± 0.5 m s⁻¹ over the study period. Maximum peak values were recorded in December and January, reaching levels of up to 19.9 m s⁻¹.

DO levels showed distinct seasonal patterns across all stations, with the highest values being recorded in spring (up to 300 µmol kg⁻¹) and the lowest being recorded during late summer, particularly at CA, reaching around 195 µmol kg⁻¹. Average DO concentrations were similar at PB and CA at 233.5 ± 22.0 and 231.8 ± 25.3 µmol kg⁻¹, respectively, while limited data were available for CS. AOU values were primarily negative, indicating net photosynthetic production. Nutrient concentrations displayed spatial and temporal variability. ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations ranged from 0.1 to 6.3 µM, with higher values being observed at CS (average 1.5 ± 1.5 µM) compared to at PB and CA. ${\mathrm{NO}}_{\mathrm{2}}^{-}$ levels were generally low (0.1–0.3 µM) but exhibited significant interannual differences. DOC showed seasonal variation, with values ranging from 52 to 116 µM, with the highest values in summer and lower concentrations in winter (Table 1).

3.2 Nitrous oxide and related biogeochemical variables

Over the monitored period from 2018 to 2023, dissolved N₂O concentrations exhibited a seasonal pattern (p < 0.005, Fig. 3a), ranging from 6.5 to 9.9 nmol L⁻¹, with no significant differences between stations and years (p > 0.05). The highest concentrations were recorded during winter, averaging 9.0 ± 0.2 nmol L⁻¹ (with a peak in February), followed closely by autumn (September–November) at 8.2 ± 0.6 nmol L⁻¹ and spring at 7.9 ± 0.7 nmol L⁻¹ (Table 1). In contrast, the lowest values were observed in summer, averaging 6.9 ± 0.2 nmol L⁻¹, with the minimum recorded in August 2021 at PB. The opposite trends were noted in Sat % N₂O percentages, with maximum values in summer and autumn and minimum values in winter, ranging from 93 % to 116.7 % (Fig. 3b).

Chl a followed a distinct seasonal pattern (Fig. 3d, Table 1), with significant differences being observed between sites, years, and months (p < 0.005). In January, we recorded the highest productivity, with a mean Chl a value of 0.4 ± 0.2 µg L⁻¹, while May showed the lowest productivity at 0.2 ± 0.05 µg L⁻¹. PB exhibited the highest levels at 0.3 ± 0.2 µg L⁻¹, whereas CA (0.2 ± 0.2 µg L⁻¹) and CS (0.2 ± 0.1 µg L⁻¹) displayed comparable concentrations. In 2019, all stations showed the strongest Chl a signal (Fig. 3d).

DO presents significant monthly differences observed across all stations (p < 0.005, Fig. 3e and Table 1). During the spring, DO levels peaked at up to 300 µmol kg⁻¹, while at the end of the summer, minimum values around 160 µmol kg⁻¹ were recorded, particularly at the CA site. The average levels from 2018 to 2023 were consistent between PB and CA at 233.5 ± 22.0 and 231.8 ± 25.3 µmol kg⁻¹, respectively. In contrast, only limited sensor data were available for CS (Fig. 3e, green dots). However, it is well-known that relevant oxygen production in the area is closely related to the presence of seagrass meadows, dominated by the endemic species Posidonia oceanica, as previously noted by Agueda-Aramburu et al. (2024) for the precise CA station location of this study.

Figure 4 The relative importance of the leading environmental and biogeochemical variables driving the variability of nitrous oxide concentration variability analysed through cross-valued boosting (CVB). These are, in decreasing order, temperature, Chlorophyll a (Chl a), salinity, dissolved organic carbon (DOC), nitrate (NO₃), and nitrite (NO₂). The figure was developed using Python software version 3.12.3.

AOU values were predominantly negative (Table 1), with average values of −14.1 ± 12.0, −11.9 ± 15.3, and −15.1 ± 4.9 µmol kg⁻¹ for the PB, CA, and CS sites, respectively. The negative AOU levels indicate that the DO produced by the excess of photosynthesis compared to respiration, which aligns well with the organic matter content observed from the DOC data (Fig. 3f). In the coastal Balearic Sea, DOC showed similarities between sites and years (p > 0.05), with significant differences among months (p < 0.05). During the summer season, the highest average levels of 83.9 ± 0.9 µM were recorded, while winter exhibited the lowest at 73.8 ± 1.5 µM (Fig. 3f, Table 1). Variability in DOC data was notable in spring and autumn, with average values of 76.0 ± 5.9 and 76.6 ± 3.1 µM, respectively.

${\mathrm{NO}}_{\mathrm{3}}^{-}$ levels ranged from 0.1 to 6.3 µM during the study period, showing no significant differences between months but significant differences between years (p < 0.005) and sites (p < 0.005). CS exhibited the highest ${\mathrm{NO}}_{\mathrm{3}}^{-}$ levels, with an average of 1.5 ± 1.5 µM, followed by PB with 0.3 ± 0.4 µM, while CA had the lowest at 0.1 ± 0.2 µM (Fig. 3g, Table 1). ${\mathrm{NO}}_{\mathrm{2}}^{-}$ concentrations were significantly different between years (p < 0.005; Fig. 3h, Table 1) and between stations (p < 0.05 for CA; Table 1), with total average values of 0.03 ± 0.02, 0.05 ± 0.05, and 0.07 ± 0.06 µM for CA, PB, and CS, respectively.

The analysis of the effects of environmental and biogeochemical variables on N₂O concentrations through the CVB revealed only seven main parameters associated with the variability of N₂O concentrations (Fig. 4). These parameters, in order of importance, were temperature (65 %), Chl a (17 %), salinity (7 %), DOC (5 %), ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (4 %), and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ (2 %).

Temperature values, the first variable modelling N₂O concentrations, exhibited a strong linear decreasing correlation with N₂O levels, with p < 0.005 and an R² value of 0.94 for all evaluated stations (not shown). The solubility of N₂O declines in warmer waters, which explains the observed maximum and minimum seasonal patterns of N₂O, primarily governed by the thermodynamic effects of temperature on N₂O concentrations. This strong correlation reflects the thermodynamic influence of temperature on N₂O solubility and its broader role in regulating stratification and microbial activity. The Sat % N₂O used here as a proxy for the biological consumption and production of N₂O shows a positive linear correlation with temperature (p < 0.005 and an R² value of 0.59; Fig. 5). This positive correlation is anticipated since the microbial activities of nitrification and denitrification are expected to increase with rising temperatures, leading to greater N₂O production (Wu et al., 2018). However, our results indicate that the seasonality control of N₂O by thermodynamic solubility surpassed the enhancement of microbial activity. The dominant role of temperature observed here aligns with previous large-scale coastal assessments (e.g. Yang et al., 2020), which identified temperature as a key variable influencing both gas solubility and microbial turnover, particularly in surface waters affected by seasonal thermal stratification.

Figure 5 Linear relation saturation percentage of N₂O (Sat % N₂O) versus temperature (°C) from the data obtained at PB (blue dots), CA (red triangles), and CS (green squares). Linear equations are represented in blue for PB, red for CA, and green for CS stations. Figures were developed with the Python software version 3.12.3.

The DO availability and the predominance of negative AOU levels, with no correlation to Sat % N₂O (p > 0.05, R² = 0.002), suggest that in situ N₂O production makes a minor contribution to the observed temporal variability of N₂O concentration. Overall, the low significance of predictors in the CVB analysis, such as Chl a, DO, DOC, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, and ${\mathrm{NO}}_{\mathrm{2}}^{-}$, is consistent with limited N₂O production from nitrification, which is associated with the seasonal remineralization cycle of organic matter, particularly enhanced under high-productivity conditions. This interpretation aligns with previous studies in coastal and shelf environments that also reported low N₂O accumulation under oxygenated and nutrient-limited conditions (e.g. Babbin et al., 2015; Lin et al., 2016). These conditions can suppress canonical N₂O production through nitrification and denitrification, emphasizing the sensitivity of microbial N₂O cycling to ambient redox and substrate gradients.

DOC and Chl a also contribute to the availability of organic substrates and oxygen fluctuations, shaping microbial N₂O dynamics, especially under stratified or high-productivity conditions. In fact, DOC enrichment has been linked to elevated heterotrophic activity and subsequent shifts in nitrogen cycling pathways in productive coastal areas (Herbert, 1999), potentially influencing the balance between N₂O production and consumption. Microbial production by nitrification and denitrification are generally considered to be the dominant pathways of N₂O biological production in the shelf sea and open ocean (Burgos et al., 2017; Chen et al., 2021; de la Paz et al., 2024; Sierra et al., 2020). However, in the oxygenated coastal Balearic Sea waters, near-equilibrium values for Sat % N₂O and low ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ concentrations (Figs. 3b, g, and h and 5, Table 1) suggest a negligible contribution from pelagic nitrification. Nevertheless, ${\mathrm{NO}}_{\mathrm{3}}^{-}$ variability may still indicate diffuse anthropogenic inputs, potentially enhancing local nitrification under favourable conditions. While benthic N₂O fluxes were not measured during this study due to logistical constraints, their contribution cannot be dismissed. A follow-up study is being conducted to evaluate potential sediment-derived N₂O inputs at the same study sites.

In addition, the CVB analysis identifies Chl a as the second driver explaining the observed variability in N₂O concentrations. The seasonal cycle of Chl a is closely coupled to that of N₂O, with maximum values in winter and minimum values in summer. A hypothetic mechanism that could be contributing to the N₂O concentration variability at our study site is the poorly described production of N₂O by epipelagic photosynthetic organisms in the light through NO reduction. Given the low ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ concentrations, the lack of correlation between AOU and N₂O, and the seasonal coupling between Chl a and N₂O peaks, we propose that an additional, non-canonical pathway may contribute to N₂O production in our system. Recent studies have identified a light-dependent mechanism in which marine photosynthetic organisms, including green algae and cyanobacteria, reduce intracellular nitric oxide (NO) to N₂O, potentially as a detoxification strategy under fluctuating nitrogen conditions or during photosynthetic stress. This NO reduction pathway involves the activity of flavodiiron proteins (FLVs) or other reductases and may become relevant in oxygenated, nutrient-depleted surface waters with active phytoplankton populations (Burlacot et al., 2020). While our dataset does not allow for direct confirmation of this process, the observed seasonal patterns and weak traditional N₂O signals warrant further consideration of this alternative mechanism, especially in coastal systems dominated by epipelagic productivity.

Finally, significant unexpected negative linear correlations (p < 0.05) were found between N₂O concentrations and silicate only at the PB and CA sites (Table 1), with correlation coefficients of −0.41 and −0.34, respectively. This relationship and the differences in salinity values between PB and CA compared to CS may indicate the presence or absence of groundwater input pulses in these areas (Basterretxea et al., 2010; Sospedra et al., 2018), which can affect N₂O concentrations in specific coastal regions (Calvo-Martin et al., 2024). Since this relationship was only observed at two of the three stations considered, this parameter was excluded from the CVB analysis. Together, these factors illustrate a multifactorial control of N₂O production and accumulation in the coastal Balearic Sea, where solubility, productivity, and nutrient inputs act in combination.

An OLS regression was conducted to evaluate the relationship among the three variables with the highest relative importance for N₂O concentrations (importance > 5 %), accounting for 89 % of the variability (temperature, Chl a, and salinity; Fig. 4) identified in the CVB analysis. The resulting equation is as follows: N₂O = 7.42 − 0.18 ⋅ temperature + 0.28 ⋅ Chl a + 0.11 ⋅ salinity (R² = 0.94), with a p value of < 0.005 for temperature and p values of 0.07 for Chl a and 0.152 for salinity.

To resolve the nonlinear dependency of Chl a and salinity on N₂O concentrations, a GAM was applied, explaining 96 % of the response variable, with a p value of < 0.005 for temperature, < 0.05 for Chl a, and 0.1390 for salinity (see Fig. S2b–c in the Supplement). The non-significance of salinity may be due to the fact that GAM smooth terms are designed to capture broader, nonlinear effects. In this case, salinity may not have a robust and smooth effect, or its impact may be masked by other variables, such as temperature. Salinity also influences stratification dynamics and can indirectly affect oxygen availability and microbial processes linked to N₂O. The GAM residuals ranged from −0.6 to 0.6 nmol L⁻¹ (see the figures in the Supplement).

When compared to other temperate and subtropical coastal regions, our results support the emerging view that N₂O dynamics in seagrass-dominated systems may diverge substantially from classic estuarine and shelf models, with solubility and light-driven processes playing an increasingly relevant role (Rosentreter et al., 2023; Al-Haj et al., 2021). Placing the obtained data into a global ocean context, specifically among studied coastal areas, the Balearic Sea presents N₂O concentrations similar to those found in 2008 in the Bohai Sea, ranging from 7.14 to 8.32 nmol L⁻¹ (Gu et al., 2022), but these are low in comparison to the values of upwelling and estuarine coastal areas. The maximum N₂O concentrations found in this study represent approximately baseline levels for the other coastal and open-ocean areas evaluated previously (Bange, 2006; Bange et al., 1996; de la Paz et al., 2024; Yang et al., 2009; Zhang et al., 2008).

Figure 6 Bar chart of monthly averaged values of air–sea nitrous oxide (N₂O) transfer (FN₂O) in $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ for all study years (left axis) and averaged temperature (°C) values (right axis; grey line) for the PB (a), CA (b), and CS (c) sites. Error bars represent ± standard deviation. Analysis of individual FN₂O flux measurements revealed that none differed significantly from zero (p > 0.5). Figures were developed using Python software version 3.12.3.

3.3 Air–sea N₂O exchange

FN₂O values ranged from −0.3 to 0.6 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ during the study period, indicating that the coastal Balearic Sea functions as a minor sink or is nearly in equilibrium in winter and spring and acts as a light source in summer and autumn (Fig. 6; p < 0.005), closely following the seasonal temperature variations (Fig. 6, grey line). No differences between years were observed. However, significant differences were noted between sites, with PB exhibiting the highest FN₂O levels and CA showing the lowest (p < 0.05).

At the PB site, positive linear correlations were found between FN₂O and temperature (R² = 0.67; p < 0.005) and between FN₂O and ΔN₂O (R² = 0.96; p < 0.005). In the OLS analysis, wind speed also appears to exert a positive influence on FN₂O values, detectable only in PB data (OLS, R² = 0.94, p < 0.05). This effect is likely to be attributable to PB site geographical characteristics as it represents an open bay, which is in contrast to the enclosed CA station in the Santa Maria Bay sampling site and the shallow coastal waters of CS.

FN₂O at the most pristine station (CA) showed positive correlations with temperature (R² = 0.79; p < 0.005), ΔN₂O (R² = 0.97; p < 0.005), and Sat % N₂O (R² = 0.97; p < 0.005).

The CS shallow site showed positive linear correlations with temperature (R² = 0.73; p < 0.005), ΔN₂O (R² = 0.98; p < 0.005), and Sat % N₂O (R² = 0.97; p < 0.005). Salinity was also a predictive variable included in the OLS analysis (OLS, R² = 0.96, p < 0.05) and was negatively related. This characteristic may be associated with the differences in CS salinity values compared to the other stations.

Table 2 Annual N₂O flux ($\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) per station since 2019 for the PB (Buoy of Palma) and CS (Cape Ses Salines) stations and from 2020 for the CA (Cabrera National Park) site, obtained using the Dobashi and Ho (2023), Wanninkhof (2014), and Cole and Caraco (1998) gas transfer parameterizations. Variability is represented as ± propagated error.

Annual average fluxes for the three sites were calculated; when monthly data were unavailable, linear interpolation was applied (Table 2). The different parameterizations for the gas transfer velocity indicate higher values using the approach of Cole and Caraco (1998), ranging from 8.7 to 204.2 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, followed by results obtained with the Wanninkhof (2014) equation, which range from 9.9 to 190.1 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, and the lowest values defined by Dobashi and Ho (2023), ranging from 5.4 to 107.4 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (Table 2). These differences are clearly influenced by the environment since the Dobashi and Ho (2023) equation was derived from studies conducted in shallow areas at approximately 3.5 m depth, while Wanninkhof (2014) and Cole and Caraco (1998) focused on deeper ocean areas and a lake, respectively. However, since the benthic zone of this study is heavily covered with P. oceanica seagrass meadows and because the fetch area is limited, we opted to retain the equation of Dobashi and Ho (2023) for the primary results.

Similar low air–sea N₂O fluxes have also been reported in shallow vegetated environments, such as seagrass-dominated coastal lagoons in Australia (Rosentreter et al., 2023), where negative or near-zero fluxes reflect the combined effects of high primary production, oxygenation, and limited nitrogen availability. In contrast, higher fluxes, often exceeding 1 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, are typical of eutrophic estuaries or upwelling regions (Bange, 2006; Zhang et al., 2008), highlighting the relatively pristine nature of the Balearic coastal zone.

The results showed that all the stations are significantly neutral with regard to N₂O emissions to the atmosphere. However, notably higher values were detected at the PB station, as seen in the daily FN₂O values (Fig. 6). PB exhibited significant interannual variability during the study period (Table 2). The CS and CA sites followed the same trend as the PB station but had considerably lower values, particularly for CA. Observed annual N₂O fluxes in the Balearic Sea are much lower than the expected ranges for this ocean region, as described by Resplandy et al. (2024) based on global reconstruction products from 1985 to 2018 (Yang et al., 2020). These differences may relate to the substantial variability in local N₂O fluxes in coastal areas with vegetation. Coastal regions with seagrass meadows are expected to exhibit the lowest N₂O flux ranges (Rosentreter et al., 2023). These areas are considered to be negligible sources (Al-Haj et al., 2021) or sinks in near-pristine seagrasses, leading to significantly biased coastal global estimations (Chen et al., 2022). These low fluxes can be partly attributed to the biogeochemical effects of seagrass meadows. Furthermore, elevated rates of primary production reduce the availability of inorganic nitrogen in the water column, thereby limiting N₂O generation. Additionally, the dense root systems and organic exudates of seagrass meadows support microbial communities that further facilitate N₂O consumption under oxic conditions (Rosentreter et al., 2023; Al-Haj et al., 2021). Collectively, these mechanisms contribute to the role of vegetated coastal systems as potential N₂O sinks.

These observations place the Balearic Sea at the lower end of the global N₂O flux spectrum for coastal systems. For example, in the East China Sea and Pearl River Estuary coastal areas, fluxes ranged from 5 to 17 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ (Zhang et al., 2008; Lin et al., 2016), often driven by hypoxia and strong anthropogenic nitrogen inputs – conditions absent in our study area.

Considering the fact that European seagrasses represent approximately 6 % of global seagrass cover and that most of this is found in the Mediterranean Sea (Jordà et al., 2012), it is crucial to understand more about the drivers of N₂O formation, consumption seasonality, and transport pathways. Additionally, increasing the number of observations is vital since N₂O data in the Mediterranean Sea are limited. Furthermore, we must consider the ongoing coastal eutrophication linked to anthropogenic inputs from estuaries, sewage discharge from densely populated coastal areas, and industrial effluents that may significantly alter seagrass habitats, along with current emissions of N₂O from coastal zones (Bakker et al., 2014).

The low N₂O fluxes observed in our study are in line with recent assessments that emphasize the underestimated climate mitigation potential of seagrass meadows and vegetated coastal habitats (Al-Haj et al., 2021; Chen et al., 2022). Therefore, understanding the ecological conditions that suppress N₂O emissions in these systems is crucial for refining regional and global GHG budgets.

As coastal N₂O emissions are assumed to offset 30 %–58 % of the net CO₂ coastal uptake radiative effect (Resplandy et al., 2024), estimating the water N₂O emission-based global warming potential (GWP) can provide a clear picture by comparing it with existing data on CO₂ and methane fluxes in the area (Flecha et al., 2023). By following the Intergovernmental Panel on Climate Change (IPCC) Assessment Report 6 (Arias et al., 2021) updated 100-year GWP, we obtained the ${\mathrm{GWP}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}$ (i.e. the 100-year time-integrated radiative forcing from the instant release of 1 kg of N₂O is 273 times larger than the forcing of 1 kg of CO₂).

In the coastal Balearic Sea, considering an area 1 km offshore and 1428 km in coastal longitude, the average ${\mathrm{GWP}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}$ obtained in this study is 8.1 × 10⁷ ± 5.7 × 10⁷ Pg CO₂-eq yr⁻¹ for all sites from 2020 to 2022, which was the period with the most available data. This area estimate was used as a first-order approximation of the active N₂O exchange zone in the Balearic region, following a similar approach to other regional-scale studies (e.g. Resplandy et al., 2024; Rosentreter et al., 2023). It encompasses both vegetated and non-vegetated shallow coastal habitats and is not intended to represent specific benthic types. Therefore, it reflects a conservative estimate of the total coastal footprint, acknowledging that further habitat-specific refinement would improve spatial resolution. The associated uncertainty in the GWP estimate incorporates the observed variability in fluxes across sites and seasons.

4 Conclusions

N₂O concentrations in seawater during the study period ranged from 6.5 to 9.9 nmol L⁻¹ without significant differences between the three sampling sites. Several drivers dominated the variability of N₂O concentration, with temperature being the most essential factor, followed by the less critical Chlorophyll a and salinity. Even with the possible biological implications in N₂O formation, atmospheric forcing may control the surface concentrations in this area. Averaged estimated N₂O fluxes oscillated between −0.3 and 0.6 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, with the most influenced stations showing the highest N₂O fluxes. Some sites appear to be weak sources of N₂O, although many fluxes are close to equilibrium and are not significantly different from zero, following a robust seasonal pattern throughout the sampling period. Reducing anthropogenic pressures in coastal regions is crucial to preserving marine ecosystems as increasing human impacts, such as pollution, overfishing, and habitat destruction, can lead to irreversible damage, loss of biodiversity, and the degradation of vital ecosystem services. If these pressures continue to escalate, future projections indicate more frequent and severe environmental crises, including increased coastal erosion, declining fish stocks, and heightened vulnerability to climate change, threatening marine life and human coastal communities. There is a strong need to enhance observations of the GHG N₂O in coastal areas to better understand dominant drivers and to make more accurate predictions of future consequences under additional anthropogenic impacts. Moreover, the values presented for the coastal Balearic Sea will aid in improving global N₂O emission budgets in coastal vegetated areas overall.