Modelling CO2 budget of mussel farms across the Mediterranean Sea

The role of bivalve aquaculture as a carbon sink is highly debated, without a general consensus on the components to include in the budget. This study proposes to estimate the terms of the budget using a scope-for-growth-based model. The model was applied at 12 Mediterranean sites, with environmental forcings provided by operational oceanography data spanning over 12 years. Mussels were found to be slight sinks, with a limited variability across sites, if all components of the budget, i.e. accumulation in soft tissue, emissions associated with calcification and respiration, are included. The differences found among stations concerning the calcification and soft tissue contributions to the budget were found to be related to site-specific productivity and water chemistry parameters. This led to the identification of a set of meta-models, which could be used for relating the budget to local conditions, at a screening level, rendering them useful for optimal site selection.


INTRODUCTION
'Restorative aquaculture' has been defined as commercial or subsistence aquaculture that provides direct ecological benefits to the environment, with the potential to generate net-positive environmental outcomes (The Nature Conservancy 2021).This role is mostly provided by non-fed, low trophic, aquaculture (Barrett et al. 2022).Within this, shellfish are increasingly being looked at as sustainable food sources that can provide additional ecosystem services, such as biodiversity enhancement, nitrogen removal and carbon sequestration (van der Schatte Olivier et al. 2020; Theuerkauf et al. 2021).Low GHG emissions are associated with bivalve cultivation (e.g.GHG emissions by oyster culture are below the 0.5% of those of land based meats, Ray et al. 2019), and they are also often considered as a carbon sink (Lee et al. 2020) granting them the definition of 'climate positive' foods (Jones et al. 2022).Theuerkauf et al. (2019) proposed a global analysis aiming at identifying the regions where marine aquaculture (shellfish and seaweed) would provide more benefits for nature and people, showing a particularly high opportunity value in temperate regions.
Molluscs represent 16% of Mediterranean mariculture by volume: of these, mussels represent 11% of the share, being the fourth most cultivated species overall, however they represent only 3% of the total production by value (Carvalho and Guillen 2021).Mussels are considered to have a high potential to become a major source of protein for human consumption, as well as playing a role for the sustainability of seas and oceans, if best management practices are adopted (Suplicy 2020).Additional finances to the sector could come via a payment for ecosystem services, for example with the generation of tradable credits, tax-payer funding, subsidies or increasing value to consumers (van den Burg et al. 2022).
Proper quantification of these ecosystem services and indicators is important for the implementation of financial mechanisms.Specifically, with reference to the carbon sequestration, there is still a debate on whether mussel farms are actually a carbon sink or a source (e.g.Munari et al. 2013;Aubin et al. 2018;Filgueira et al. 2019;Tamburini et al. 2022).While it is now widely accepted that the release of CO 2 associated with the shell calcification process and respiration should be accounted for, views range from considering CO 2 release associated with calcification as the only process to include in the calcification budget (A ´lvarez-Salgado et al. 2022), to potentially ignoring it as being part of the 'short' cycle together with the CO 2 from respiration, given the likelihood of it being up-taken by phytoplankton, and therefore account only the CO 2 stored in the shell in the form of carbonate (Frankignoulle et al. 1994;Tamburini et al. 2022).Another discordance arises as to whether to include soft tissue, as mussels are cultivated for food purposes, or only the shell as a by-product (Filgueira et al. 2015), depending also on its potential fate (Aubin et al. 2018).The methodologies used for quantifying the budget may also lead to different estimates: the carbon budget is often quantified using indirect methods, e.g.estimating 'shell contribution to total harvest' (Alonso et al. 2021) but to make generalised assumptions it should be realised that growth (thus shell making and calcification) has a strong seasonal component and so do the physicochemical parameters that drive the emission factor linked to calcification (e.g.Bertolini et al. 2021), which are particularly sensitive to temperature (Frankignoulle et al. 1994;Morris and Humphreys 2019).Respiration has also been quantified mostly from indirect methods (Hily et al. 2013;Munari et al. 2013), which may result in inaccurate estimations.
In this study, we attempted a quantification of the mussel carbon budgets across the Mediterranean Sea basin, focusing on areas where aquaculture of this species is practised, to obtain the estimates from areas characterised by different temperature, food availability and water carbonate chemistry (pH, alkalinity).Removal and emissions from the two compartments of meat production and shell (by-product) were kept separated, in accordance with Filgueira et al. (2015), to focus on site-specific differences and avoid arbitrary choices concerning the processes to be included in the carbon budget, considering the following processes: (1) CO 2 accumulation in the soft tissue; (2) CO 2 emissions associated with calcification and respiration.Ingested and egested carbon was also reported for completeness but investigating the fate of egested material and rejection through pseudo-faeces was beyond the scope of this study.

Model
The carbon budget associated with mussels was simulated using a dynamic bioenergetic model, based on the Scope for Growth (Brigolin et al. 2009), which explicitly includes shell calcification (Bertolini et al. 2021).The model requires as forcings time series of water temperature, chlorophyll-a, particulate organic carbon POC, salinity, pH and alkalinity: the last three variables are involved in the computation of the CO 2 emissions related to shell formation (Humphreys et al. 2018;Morris and Humphreys 2019).
Most parameters of the bioenergetic model were taken from the literature (Brigolin et al. 2009).In this paper, we focused on the estimation of the parameter 'K sh ', which represents the fraction of the Scope for Growth invested in shell formation.This parameter was estimated by fitting the model to a time series of shell weight and soft tissue dry weights collected in a mussel farm located in the southern part of the Venice Lagoon in 2020-2021 (coordinates 45.23,12.27)at five time points ('26-05-2020', '29-07-2020', '29-10-2020', '28-05-2021').The function chosen for the minimisation was the standardised residual sum of squares error (Eq.1), in which: Wb i d Wbi and represent, respectively, the observed and simulated weight of soft tissue at time ''i'', and r Wb the data set standard deviation, while Wsi and d Wsi and r Ws represent, respectively, the same quantities concerning shell weight This function was minimised using the fminbnd package Pracma (v 2.2.9) in R (version 4.2.1).Minimisation was repeated 100 times, in order to estimate the possible range of k SH , bound between 0.1 and 0.6 (10-60% of the energy).The bound was set to leave at least 10% of energy to somatic tissue, given that reproduction investment was set as 30%.A Montecarlo methodology was used, with each simulation starting with state variables at t 0 extracted from a normal distribution defined by the means and standard deviations of the sampled population at t 0 using rnorm (W b = 0.06 ± 0.025; W r = 0.01 ± 0.005; W s = 0.4 ± 0.14).As already done by Brigolin et al. (2009), each individual was also assigned a specific maximum clearance rate extracted from a normal distribution assuming the mean value and standard deviation of the parameter, between 107 ± 50 l day 1 .The comparison between model predictions and empirical data was made at the 5 temporal steps coinciding with sampling (times 0, 63,153,267,365).

Applying the model to the whole Mediterranean basin
In order to quantify carbon fluxes for the main farming area in the Mediterranean basin, 12 representative musselfarming sites (Fig. 1; Table S1) were selected.These sites are characterised by remarkable differences in the level of development of the mussel-farming industry, ranging from locations at which the activity has a long-term tradition and a remarkable extension (e.g.Cattolica, Italy) to those at which mussel farms started only recently (e.g.Aı ¨n Taya, Algeria).
The input data were estimated from products downloaded from the Copernicus Marine Service portal,1 which provides both earth observations and operational oceanography products.These included the following: Twelve yearly simulations per site were performed, with a standard cycle starting on September 1st (Brigolin et al. 2009(Brigolin et al. , 2018)).The initial conditions were set in accordance with standard seeding practices in Italy: the initial shell length was extracted from a normal distribution with mean 20 mm and standard deviation 2 mm and converted into dry weight using the allometric relationship from Brigolin et al. (2009).The shell to soft tissue ratio was based upon the sampling described above, and it was assumed that the reproductive tissue weight at the beginning of the simulation was nearly negligible given the young age (W b = 0.02 ± 0.002, W r = 0.001 ± 0.0001, W s = 0.15 ± 0.015).The maximum Clearance Rate (CR max ) was again imposed to vary among individuals as described previously, while k sh was randomly extracted based on the mean and standard deviation of the parameter as estimated above, in both cases assuming a normal distribution of the parameters.

Carbon fluxes
The ingested carbon was quantified based on the absorbed/ digested particulate organic carbon (AE*CR(t)*POC(t)), while the remaining part, not absorbed, was considered as egested ((1-AE)*CR(t)*POC(t)), where AE is the absorption efficiency and CR the clearance rate.The egested carbon was not transformed into CO 2 equivalent (CO 2eq ), as its fate is highly dependent upon local environmental All other carbon fluxes included in the budget were expressed as CO 2eq : emissions from calcification and respiration are expressed as CO 2, whereas the carbon content of soft tissue and shell [g] are converted into CO 2eq by multiplying them by 3.664 (ratio between the molar weights 44/12).Fluxes associated with shell formation were calculated according to Morris and Humphreys (2019), who pointed out that calcification entails the removal of DIC and a decrease in total alkalinity, which results in a net flux of CO 2 to the atmosphere In this study, this term of the budget was estimated stepwise: (1) the calcification rate [g CaCO 3 /day] was calculated at each time step by multiplying the growth rate by CaCO 3 percentage in shell, 0.95.(2) The rate of CO 2 released to the atmosphere was then obtained by multiplying the calcification rate by the function U(t), estimated by means of the R package seacarb, and the CO 2 /CaCO 3 molar ratio.In accordance with A ´lvarez-Salgado et al. ( 2022), the calcification therefore represents always a CO 2 emission, as the function U (t) already includes an account of the carbon incorporated in the shell.The energy derived from the respiration was allocated to the maintenance of somatic tissues, reproductive tissues and shell, similar to Filgueira et al. (2019).Fluxes associated with the soft tissue were calculated based on the daily weight changes.The gross value considered that carbon makes up to 40.8% of the tissue (Brigolin et al. 2009).To obtain the net value, the % of respiration, according to the % of energy invested towards soft tissue growth, was then subtracted.

Meta-model for estimating C budget
The model outputs were post-processed, to identify a set of statistical meta-models which could be used for a preliminary estimation of C budget of farmed mussels.Yearly median values were computed for the three terms of the C budget and for the 5 environmental input variables.The terms of the C budget were first aggregated as follows: Equations 2 and 3 consider a shell-soft tissues partitioning approach, in accordance with Filgueira et al (2015), but using k SH as partitioning coefficient.An analysis of variance was used to address whether significant differences existed between stations or years.6 statistic meta-models (generalised additive models and linear models) were identified, by considering as independent variables the environmental variables and as dependent variables the terms of the C budget (Eqs.2-4).Meta-models construction was carried out for each of the three dependent variables (CO 2 shell net, CO 2 soft tissues net, and CO 2 total).Generalised additive models (function gam in package mcgv in R, Wood 2006) were constructed, considering the three subsets of the 5 environmental variables producing an effect.Linearity was therefore tested, and linear models were identified.An analysis of variance was performed on residuals with respect to stations, in order to exclude the presence of a location effect.Models were compared based on the % of deviance explained and the Akaike Information Criterion (AIC).

Model parameterisation-Venice
From the estimation of K sh parameter, it was found that 27 ± 0.9% of the energy budget is invested in shell making (Fig. 2).The model could predict well the growth data collected in the Venice lagoon, aside from a little overestimation of soft tissue at day 150 (Fig. 2).

Model application-Mediterranean basin
Figure 3 shows the median annual values of the environmental forcings at the 12 representative sites: the minimum values of Chla, 0.05-0.07lg L -1 , were reached in station 2 (Mersin Bay, Turkey) and station 11 (Gulf of Tunis, Tunisia) while the maximum values, 2.47 lg L -1 , were reached in station 7 (Cattolica, Italy).This pattern was followed by POC (lg L -), reaching values of 20-22 lg L -1 in stations 2 and 11 up to 546 lg L -1 in station 7. The stations with the two lowest trophic status also had the two highest annual median temperature values (station 2, up to 23.5 °C followed by station 11, up to 20.7 °C), while station 7 and 4 (Montpellier, France) had the lowest values (15.2 °C).Annual median pH had a relatively small range from 8.05 in station 2 up to 8.2 in station 5 (Gulf of Orfani, Greece).Total alkalinity ranged from 2.5 mM in station 3 (Cala Iris, Morocco) to 2.8 mM in station 7.
In terms of growth (Fig. 4), the site with the lowest total final dry weight was again station 2 (5.3 ± 0.5 g) whereas the highest values of 9.6 ± 0.5 and 9 ± 0.6 g were reached in station 3 and 12 (Ain Chrob, Algeria).Despite the high primary production in station 7, mussels in this site were only second in size to mussels in station 2 (5.6 g in the worst growing year but averaging 6.8 ± 0.6 g through the whole period).Station 2 had also the lowest values of ingested and egested carbon (488 and 209 g year -1 ind -1 , respectively) whereas the highest values (of an order of magnitude greater) are observable in station 7 (up to 6745 and 2890 g year -1 ind -1 , Online Resource 1 S2).

CO 2 fluxes
The total budgets at each site, see (Fig. 5), did not show any significant effect of station or year (Table S2).On the other hand, fluxes concerning shell and soft tissue were found to be different across sites (F 1,141 : 8.6, p \ 0.01; F 1,141 : 11.4, p \ 0.001).Station 2 (Mersin Bay, Turkey) and 7 (Cattolica, Italy) had the lowest emissions from calcification, but also the lowest amount of CO 2eq associated with the soft tissue, station 3 (Cala Iris, Morocco) showed the opposite pattern with the highest emissions from calcification and respiration but also the largest amount of CO 2eq associated with the soft tissue.
When linking these outputs to the environmental forcings (Fig. 6), it can be observed that the model explaining the largest fraction of variability is the one considering the shell CO 2 as a dependent variable (Table 1), with total alkalinity having the highest significance.In this case, the linear model presents a slightly higher deviance explained and lower AIC with respect to the GAM.When looking at soft tissue CO 2 , temperature and Chla had, respectively, the highest and the lowest significance.Models considering CO 2 total as the dependent variable provided the lowest explanation.In these two latter cases, linear models provided lower deviance explained and higher AIC values with respect to GAMs, with linearity for Chla and POC not verified.

DISCUSSION AND CONCLUSION
This study took a modelling approach to quantify carbon fluxes associated with mussels aquaculture at the scale of the Mediterranean basin, in areas characterised by different temperatures, food availability and water carbonate chemistry.The findings highlighted that the mussels could be considered a slight sink in all regions and the mean net value during 1 year (approximately corresponding to a culture cycle) was estimated to be around -0.55 g CO 2eq ind -1 , varying between different years and sites, resulting in occasions as a source.
The results obtained in this study align with A ´lvarez-Salgado et al. ( 2022), with CO 2 emissions linked with the calcification process ranging between 0.79 and 2.21 g CO 2eq ind -1 vs 0.94-3.2g CO 2eq ind -1 , although the accumulated amount in soft tissue in this case reached comparable yet slightly higher values in this study (1.85-5.15g CO 2eq ind -1 ) compared to theirs 1.93 g CO 2eq ind -1 , likely due to different methodologies.The values of this study also generally align with values per individual calculated by Filgueira et al. (2019) for Norwegian Fig. 4 Growth curves of total dry weight (somatic tissue, gonads and shell) at each of the 12 stations.Full line represents mean and dotted lines represent standard deviations mussels, with regard to net soft tissue of -7.7 to -3.2 g CO 2eq ind -1 which was slightly higher than our results averaging at the different sites -3.07 to -1.45-g CO 2eq ind -1 (min-max).Within this study, fluxes were calculated using a dynamic model, which considered seasonal variability (e.g.periods of slow growth or even tissue loss).Following Filgueira et al. (2019), a partitioning rule was used to calculate the amount of respiration linked to shell compared to soft tissue, but the amount of shell calcified per day was considered starting from energetic investment parameterised for the species, and while this value may be slightly site dependent (e.g.salinity, predator presence (Bertolini et al. 2021)), our results were in line with Duarte et al. (2010) which estimates the investment at 20-28%.The site specificity of the K sh parameter may lead to slightly increased or decreased emissions associated with the calcification process, both with the actual amount of CaCO 3 deposition and the allocated fraction of respiration.Shell growth is not often accessible in the literature as empirical data: condition index of mussels can be used as an indicator of shell vs soft tissue; however, it should be considered that changes in its value can be dependent on loss of soft tissue and not only on shell growth.Condition index was demonstrated to be dependent upon local water temperatures and productivity (higher when colder and more productive, Martinez et al. (2018)); however, smallscale habitat characteristics may also play a role, and therefore, further empirical quantifications would be a useful addition to this research.
As found in this study, the extent of accumulation depends upon site-specific productivity, with lower values of both accumulation and emissions in areas where mussels are smaller (e.g.stations 2 and 7).This may seem an obvious effect of the setup, given that accumulation is directly linked to the amount of soft tissue, but larger individuals also have a higher respiration rate and would emit more CO 2 due to more shell being produced.It has been already shown by A ´lvarez-Salgado et al. ( 2022) how the metabolic carbon footprint is influenced by management strategies, particularly the seeding time and harvesting size.
The shell growth process is also influenced by the physicochemical parameters in the area, i.e. temperature, pH, TA, and salinity.Primary productivity appears as an important driver.Temperature is also important, promoting the highest net sink at intermediate values, likely due to the effects of water temperature on growth.When considering only shells, the emissions associated with calcification are influenced by water chemistry, in particular total alkalinity, but should also be lower in areas with less growth and higher in areas with more growth.Emissions were, in fact, highly associated with both total alkalinity and temperature, with lower emissions at higher temperatures, as could be expected from both having lower values of U (t) and lower growth, and lower emissions associated with higher TA values, also linked with lower values of U (t). Lowtemperature extreme was also associated with lower emission values, likely due to slower growth and less respiration.As found in this study, lowest outputs of biomass and thus carbon accumulation are obtained in both station 2 (Turkey) and 7 (Cattolica, Italy), sites that, however, are characterised by different environmental forcings, with the first being characterised by low food availability in terms of Chla and POC but high annual median temperatures (22 °C), whereas the latter has the opposite situation with high food availability but low median temperatures.Both food availability and temperature have high predicting powers in terms of variance explained, with optimal values at intermediate temperatures and low to intermediate POC.In this case, site 7 was also an 'outlier' in terms of both POC and Chla, suggesting that in this site, temperature may have been the greater discriminating factor.Indicators of food availability, POC and Chla, were extracted from satellite data.While Chla from satellites has been widely used (Gernez et al. 2017;Palmer et al. 2020), the POC algorithm was validated in oceanic context (Stramski et al. 2008;Le et al. 2017;Moutier et al. 2019).While the values reported should therefore be used with caution, they are useful to make the comparisons within the study, especially given the same approach was taken for all sites.As an example, median annual egested carbon ranged between 600 and 3800 g per individual (Supplementary S2), The two sites with the highest POC (stations 4 and 7) were the sites with the greatest absorption and egestion of carbon.This consideration is particularly relevant for station 7, given the limited growth resulting in low levels of CO 2 accumulation by the mussels themselves.Different authors agree on the potential relevance of biodeposit contribution on the overall budget: Jansen and van den Bogaart (2020) reported fluxes due to biodeposition larger Fig. 6 CO 2eq budgets calculated following Eqs.2-4 (empty dots: CO 2 total; full black dots: CO 2 soft tissues net; full grey dots: CO 2 shell net), plotted against environmental forcings (from top to bottom: Chlorophyll-a, water temperature, POC, pH, total alkalinity) than the respiration ones under gradients of environmental conditions and Filgueira et al. (2019) shown a relevance for these fluxes in low seston areas.The fate of the egested and rejected carbon will be dependent on the local conditions in terms of currents and physical factors driving the sinking rate, substrate type which will determine the burial rates, oxygen, bacterial community and other factors that contribute to remineralisation.As shown in Brigolin et al. (2018), the combination of temporal variability in deposition fluxes, with sediment heterogenous conditions, can drive to differences in early-diagenesis pathways, which are in turn expected to affect benthic-pelagic fluxes of mineralised products.This focus on ecosystem-scale processes affecting the budget was beyond the scope of this study, and something that should be investigated at the relevant spatial scale.However, some research suggests that a large fraction of the deposited carbon is remineralised (Labrie et al. 2022); therefore, caution should be taken in considering this a form of long-term storage.
The inclusions of different aspects in the carbon budget will depend upon the needs, end goals and management practices within farms (e.g.selling of whole vs processed product, fate of shells, Alonso et al. 2021); therefore, the approach proposed in this study can be used as a useful tool for budget calculations also as support for financial actions (e.g.carbon credits).The total value based upon total harvest for each country should be quantified, and the losses should be included in the budget (Jansen and van den Bogaart, 2020), since they remain in the system with remineralisation of the organic carbon accumulated leading to emissions on the assessed time scale of the cultivation year and shell dissolution likely leading to some sequestration on longer time scales.In this context, mass mortalities events, for example as a result of heatwaves, may sway the carbon budget suggested in this study and should be given careful consideration.Lastly, the approach developed in this study made extensive use of freely accessible data from operational oceanography, a relevant (i) daily reprocessed data on sea surface temperature (L4) at 5 km resolution from reprocessing of different sources of L3 climate data records from the 'Mediterranean Sea-High Resolution L4 Sea Surface Temperature Reprocessed'; (ii) monthly data on chlorophyll-a concentration (L4) at 1 km resolution derived from multiple sensors (SeaWiF, MODIS, MERIS, VIIRS-SNPP & JPSS1, OLCI-S3A) from the 'Mediterranean Sea Ocean Colour Plankton MY L4 daily gap free observations and climatology and monthly observations'; and (iii) monthly average data on pH and total alkalinity at surface obtained from 'Mediterranean Sea Biogeochemistry Reanalysis' (MedBFM3 model system).POC concentration was not available from Copernicus Marine Service and was therefore based on L3 ModisAqua data (Nasa Oceandata).The spatial time series of daily SST values and monthly ones of the other input variables cover 13 years, starting from 2007 up to 2019.Monthly data were interpolated using splines to obtain daily model forcings.Salinity was kept as the average Mediterranean salinity of 38 ppm, although it can reach values closer to 39 in the Easternmost part and to 36 in the Westernmost part (Christopoulos et al. 2018).

Fig. 1
Fig. 1 Location of the 12 representative sites considered within this work.Coordinates of each site are provided in Table S1 (Online Resource 1 S1)

Fig. 2
Fig. 2 Model simulations (grey lines) for dry soft (top) and shell (bottom) tissue dry weights (in grams), compared to means (black dots) and standard deviations (error bars) sampled from the Venice lagoon

Fig. 3
Fig. 3 Boxplots of median annual values from 2007 to 2019 of environmental model forcings (Chlorophyll-a, POC, temperature, total alkalinity, pH) at each of the 12 stations

Fig. 5
Fig.5Boxplots of annual CO 2 fluxes at the 12 study sites (positive values refer to emissions, negative to removal).A total budget (calculated as described in Eq. 4), B amount stored in the soft tissue, C amount emitted through calcification, D amount emitted through respiration.All fluxes are reported as g CO 2eq ind -1 year -1 .Box edges indicate the 1st and 3rd quartiles, the central line marks the median, while whiskers indicate extreme (max-min) data points not considered outliers (plotted individually using the ''°'' symbol)

Table 1
Meta-models relating CO 2eq budgets to environmental forcing feature to consider when designing tools for management support.