Lesley A Salt
email: lesley.salt@sb-roscoff.fr.

Helmuth Thomas
email: helmuth.thomas@dal.ca

Yann Bozec
email: bozec@sb-roscoff.fr

Alberto V Borges
email: alberto.borges@ulg.ac.be

Hein J W De Baar
email: hein.de.baar@nioz.nl

The internal consistency of the North Sea carbonate system

Keywords: carbon dioxide, North Sea, internal consistency, total alkalinity 1

In 2002 (February) and 2005 (August) the full suite of carbonate system parameters (Total Alkalinity (A T ), Dissolved Inorganic Carbon (DIC), pH, and partial pressure of CO 2 (pCO 2 ) were measured on two re-occupations of the entire North Sea basin, with three parameters (A T , DIC, pCO 2 ) measured on four additional re-occupations, covering all four seasons, allowing an assessment of the internal consistency of the carbonate system. For most of the year there is a similar level of internal consistency, with A T being calculated to within ±6 mol kg -1 using DIC and pH, DIC to ±6 mol kg -1 using A T and pH, pH to ±0.008 using A T and pCO 2 , and pCO 2 to ±8 atm using

Introduction

Increasing atmospheric CO 2 has been partly mitigated by oceanic uptake, which currently accounts for 25-30% of the total CO 2 emissions [START_REF] Sabine | The Oceanic Sink for Anthropogenic CO 2[END_REF][START_REF] Canadell | Contributions to accelerating atmospheric CO 2 growth from economic activity, carbon intensity, and efficiency of natural sinks[END_REF][START_REF] Khatiwala | Reconstruction of the history of anthropogenic CO2 concentrations in the ocean[END_REF][START_REF] Quéré | Global carbon budget 2014[END_REF]. The uptake of CO 2 in surface waters causes a shift in the chemistry of seawater, notably increases in [H + ] and a decrease in [CO 3 2-] leading to the phenomenon known as ocean acidification (OA) (Caldeira andWickett, 2003, 2005;[START_REF] Raven | Ocean acidification due to increasing atmospheric carbon dioxide[END_REF]. Recent findings of Wootton et al. (2008), [START_REF] Provost | Seasonal and long-term changes in pH in the Dutch coastal zone[END_REF] and [START_REF] Borges | Carbonate chemistry in he coastal zone responds more strongly to eutrophication than to ocean acidification[END_REF] indicate that ocean acidification is proceeding in coastal areas more rapidly than has been predicted. The associated changes related to OA have been shown to have a variety of effects on marine phytoplankton [START_REF] Fabry | Impacts of ocean acidification on marine fauna and ecosystem processes[END_REF]. The coastal zone is responsible for a disproportionately large amount of primary production per surface area as compared to the open ocean (Wollast, 1998) and is responsible for a similarly large proportion of carbon export [START_REF] Dunne | A synthesis of global particle export from the surface ocean and cycling through the ocean interior and on the seafloor[END_REF]. As such, understanding how these waters respond to acidification is vital for accurately predicting the future consequences of continued CO 2 increase.

The North Sea has been identified as an effective continental shelf-pump of CO 2 [START_REF] Thomas | Enhanced Open Ocean Storage of CO 2 from Shelf Sea Pumping[END_REF] and has thus undergone intensive study with respect to the CO 2

A C C E P T E D

M A N U S C R I P T

ACCEPTED MANUSCRIPT

system. High temporally resolved datasets are necessary to understand the large, observed seasonal variability of the carbonate system in the North Sea [START_REF] Frankignoulle | European continental shelf as a significant sink for atmospheric carbon dioxide[END_REF][START_REF] Bozec | Assessment of the processes controlling seasonal variations of dissolved inorganic carbon in the North Sea[END_REF][START_REF] Prowe | Mechanisms controlling the air-sea CO 2 flux in the North Sea[END_REF]. The use of Voluntary

Observing Ships (VOS) to collect underway pCO 2 data have been used extensively in the North Sea [START_REF] Omar | Spatiotemporal variations of fCO 2 in the North Sea[END_REF], complemented by measurements of one other of the measurable carbonate variables (total alkalinity (A T ), dissolved inorganic carbon (DIC), and pH). From at least two of the measured variables (a combination pair of A T , DIC, pH or pCO 2 ), the carbonate system can then be fully solved and quantified using thermodynamic relationships (Park, 1969), which can be integrated in computer programs, such as CO2SYS [START_REF] Lewis | Program developed for CO 2 system calculations[END_REF]. These calculated parameters, however, are prone to errors, which have been well documented [START_REF] Roy | Determination of the ionization constants of carbonic acid in seawater[END_REF][START_REF] Lee | The recommended dissociation constants for carbonic acid in seawater[END_REF][START_REF] Koeve | Potential impact of DOC accumulation on fCO 2 and carbonate ion computations in ocean acidification experiments[END_REF][START_REF] Hoppe | Implications of observed inconsistencies in carbonate chemistry measurements for ocean acidification studies[END_REF]. These errors are due to the accuracy of the two measured variables, and to the reliability of the first and second dissociation constants of carbonic acid (K 1 , K 2 ). The assessment of the reliability of dissociations constants can be achieved with internal consistency studies.

The first internal consistency studies took place in the open ocean and focused on the reliability of the various sets of carbonic acid dissociation constants, which are used to calculate the distribution of inorganic carbon species. Acid dissociation constants describe the interaction between acids (HA), their dissociated proton (H+) and the conjugate base (A-), shown in equation1.

K = {H+}{A-}/{HA}, ( equation 1)
where {} denote activities. The aforementioned determinations of K 1 and K 2 of carbonic acid were performed analytically on artificial seawater (ASW) [START_REF] Roy | Determination of the ionization constants of carbonic acid in seawater[END_REF][START_REF] Goyet | New determination of carbonic acid dissociation constants in seawater as a function of temperature and salinity[END_REF][START_REF] Hansson | A new set of acidity constants for carbonic acid and boric acid in sea water[END_REF] and natural seawater (NSW) [START_REF] Mehrbach | Measurement of the apparent dissociation constants of carbonic acid in seawater at atmospheric pressure[END_REF] and the data were later refitted using improved computergenerated models [START_REF] Dickson | A comparison of the equilibrium constants for the dissociation of carbonic acid in seawater media[END_REF]. These constants are technically apparent dissociation constants (denoted by K'), as they are based on concentrations (indicated by []), not activities, so that

K' = {H + }[A -]/[HA], ( equation 2)
This approach also means that the K' values must be determined across the range of

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

salinities, temperatures and pressures at which they are to be applied and that the measurements take place in the same medium as the solution to which the K' value is to be applied. Assessments and comparisons of each set of constants have taken place in several locations, including the northwest European shelf [START_REF] Ribas-Ribas | Intercomparison of carbonate chemistry measurements on a cruise in northwestern European shelf seas[END_REF], equatorial Pacific [START_REF] Roy | Determination of the ionization constants of carbonic acid in seawater[END_REF] and the North Atlantic [START_REF] Lee | The carbon dioxide system in the Atlantic Ocean[END_REF], and spanning across several ocean basins [START_REF] Wanninkhof | The optimal carbonate dissociation constants for determining surface water pCO 2 from alkalinity and total inorganic carbon[END_REF]Millero et al., 2002). At sea studies have been complemented by further work using laboratory set-ups [START_REF] Lee | The reliability of the thermodynamic constants for the dissociation of carbonic acid in seawater[END_REF][START_REF] Lueker | Ocean pCO 2 calculated from dissolved inorganic carbon, alkalinity, and equations for K 1 and K 2 : validation based on laboratory measurements of CO 2 in gas and seawater at equilibrium[END_REF] and more recently investigations into inconsistencies in culture experiments [START_REF] Hoppe | Implications of observed inconsistencies in carbonate chemistry measurements for ocean acidification studies[END_REF][START_REF] Koeve | Potential impact of DOC accumulation on fCO 2 and carbonate ion computations in ocean acidification experiments[END_REF].

The conclusions from these studies is that the first and second dissociation constants of carbonic acid (K' 1 and K' 2 ) calculated by [START_REF] Mehrbach | Measurement of the apparent dissociation constants of carbonic acid in seawater at atmospheric pressure[END_REF], (refit by [START_REF] Dickson | A comparison of the equilibrium constants for the dissociation of carbonic acid in seawater media[END_REF], are the most reliable when pCO 2 is being used in the calculations or is being calculated [START_REF] Wanninkhof | The optimal carbonate dissociation constants for determining surface water pCO 2 from alkalinity and total inorganic carbon[END_REF][START_REF] Lee | The reliability of the thermodynamic constants for the dissociation of carbonic acid in seawater[END_REF]. These constants were the only ones determined in NSW and the difference in accuracy has been attributed to the lack of boric acid in the ASW (Mojica Prieto and Millero, 2002). In solution the majority of dissolved boron (BT = total boron) is present in two forms; boric acid and borate, following the equilibrium relationship below:

B(OH) 3 + H 2 O K B  B(OH) 4 -+ H + (equation 3) where K B =[B(OH) 4 - ][H + ]/[B(OH) 3 ] (equation 4)
In NSW boric acid interacts with HCO 3 -and CO 3 2-to lower the activities, and concentrations of HCO 3 -, and increase the activities and concentrations of CO 3 2-, thus lowering pK 1 and increasing pK 2 . Calculations of pCO 2 depend directly on the ratio of K 2 /K 1 , thus errors in this ratio will be evident in a lack of internal consistency of pCO 2 . In contrast, the lack of use of this ratio in calculations with pH, A T or DIC results in other dissociation constants being more appropriate, notably those of [START_REF] Roy | Determination of the ionization constants of carbonic acid in seawater[END_REF] or [START_REF] Goyet | New determination of carbonic acid dissociation constants in seawater as a function of temperature and salinity[END_REF] [START_REF] Lee | The carbon dioxide system in the Atlantic Ocean[END_REF][START_REF] Clayton | The role of pH measurements in modern oceanic CO 2 -system characterizations: precision and thermodynamic consistency[END_REF].

Despite this, in an attempt to homogenize the way carbon chemists carry out these calculations, it is now the recommended standard procedure to use the constants of [START_REF] Mehrbach | Measurement of the apparent dissociation constants of carbonic acid in seawater at atmospheric pressure[END_REF], refit by [START_REF] Dickson | A comparison of the equilibrium constants for the dissociation of carbonic acid in seawater media[END_REF], on the Total pH scale [START_REF] Dickson | Guide to best practices for ocean CO 2 measurements[END_REF].

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

Beyond the scope of dissociation constants, additional errors have been associated with systematical differences in equipment [START_REF] Körtzinger | The international at-sea intercomparison of fCO 2 systems during the R/V Meteor Cruise 36/1 in the North Atlantic Ocean[END_REF][START_REF] Ribas-Ribas | Intercomparison of carbonate chemistry measurements on a cruise in northwestern European shelf seas[END_REF], high DIC/A T ratios (>1.10) [START_REF] Lee | The reliability of the thermodynamic constants for the dissociation of carbonic acid in seawater[END_REF], the interference of fatty acids [START_REF] Gripenberg | On the Alkalinity of Baltic Waters[END_REF], and contribution from organic acids and bases to A T [START_REF] Hernandez-Ayon | Estimating the contribution of organic bases from microalgae to the titration alkalinity in coastal seawaters[END_REF][START_REF] Koeve | Potential impact of DOC accumulation on fCO 2 and carbonate ion computations in ocean acidification experiments[END_REF][START_REF] Ulfsbo | Modelling organic alkalinity in the Baltic Sea using a Humic-Pitzer approach[END_REF]. The latter errors are only thought to be a problem where dissolved organic carbon (DOC) reaches particularly high concentrations [START_REF] Koeve | Potential impact of DOC accumulation on fCO 2 and carbonate ion computations in ocean acidification experiments[END_REF], thus this is not an issue for open ocean measurements, although presents a major problem for calculations in freshwaters [START_REF] Abril | Technical Note: Large overestimation of pCO 2 calculated from pH and alkalinity in acidic, organic-rich freshwaters[END_REF]. Similarly, the contributions of fatty acids and high DIC/A T ratios in the open ocean are sufficiently low that internal consistency remains within the same order of magnitude as the measurement precision. In culture experiments, the differences between measured and calculated carbonate parameters have been much larger, and thus coastal areas where high primary production takes place, could also be susceptible to large discrepancies in calculations of the carbonate system.

Furthermore, in the original work of the recommended constants of [START_REF] Mehrbach | Measurement of the apparent dissociation constants of carbonic acid in seawater at atmospheric pressure[END_REF], the constants were only determined to salinities above 26. The salinity range was extended in further work by [START_REF] Millero | Dissociation constants of carbonic acid in seawater as a function of salinity and temperature[END_REF], however, the results have not yet been applied to further internal consistency studies. Coastal regions can also exhibit non-oceanic-like, chemical relationships with salinity, for example, the atypical borate-salinity relationship in the Baltic Sea [START_REF] Kremling | Untersuchungen ueber die chemische Zusammensetzung des Meerwassers aus der Ostsee[END_REF], which is not taken into account when performing carbonate system calculations and can lead to further error. Here we present the first internal consistency study performed in coastal waters with seasonal resolution over a range of salinities and biogeochemical regimes.

Methods

Hydrography

The hydrographic properties in the surface waters of the North Sea demonstrate the presence of 3 different water masses, acknowledged in the literature as North Atlantic water, Baltic water, and German Bight water, with the resulting composite water mass known as Central North Sea Water [START_REF] Kempe | Sinks and sources of CO 2 in coastal seas: the North Sea[END_REF][START_REF] Bozec | Assessment of the processes controlling seasonal variations of dissolved inorganic carbon in the North Sea[END_REF]; Salt et al., 2013). In the northern North Sea (nNS), North Atlantic water flows over the Shetland Shelf bringing characteristic warm and high saline waters into the North Sea.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
The inflowing water follows an anti-clockwise trajectory around the North Sea, with additional North Atlantic inflow occurring through the English Channel [START_REF] Otto | Oceanic CO 2 System: An evaluation of ten methods of investigation[END_REF]. In the southeast the Skagerrak water mass, a mixture of Central North Sea water and Baltic water, is easily identifiable through much lower salinities and temperatures than the central North Sea. The German Bight water comes from the coastal region bordering northwest continental Europe in the southeastern North Sea,

where the riverine influence is greatest in the North Sea.

The North Sea can be divided into two biogeochemical regimes, within which behavior is determined by depth [START_REF] Thomas | Controls of the surface water partial pressure of CO 2 in the North Sea[END_REF][START_REF] Prowe | Mechanisms controlling the air-sea CO 2 flux in the North Sea[END_REF][START_REF] Omar | Spatiotemporal variations of fCO 2 in the North Sea[END_REF]. The nNS (north of 56°N) is deeper than the south (50-400 m) and experiences more ocean-like behavior, e.g. higher salinity and seasonal stratification. The southern North Sea (sNS) is much shallower (<50 m) and remains well-mixed throughout the year. Both regions receive nutrient and organic carbon inputs from terrestrial sources, facilitating high primary production that leads to overall annual autotrophy in the nNS and near balanced trophic status in the sNS [START_REF] Bozec | The continental shelf pump for CO 2 in the North Sea -evidence from summer observation[END_REF][START_REF] Schiettecatte | High temporal coverage of carbon dioxide measurements in the Southern Bight of the North Sea[END_REF]).

Measurements and comparisons

The data used in this study comes from 6 separate re-occupations of the North Sea, which took place in August 2001, November 2001, February 2002, May 2002, August 2005and August 2008. All four carbonate parameters were measured on two of these cruises; February 2002, and August 2005, with a precision of ±2-3 mol A T kg -1 for A T , ±1-2 mol DIC kg -1 for DIC, ±0.001 for pH, and ±1 atm for pCO 2 , respectively.

The remaining cruises measured A T , DIC, and pCO 2 with the same level of precision.

Carbonate system measurements

For all of the aforementioned cruises, A T was measured by potentiometric titration to the second end-point (pK 2 ) of carbonic acid, using 0.1M hydrochloric acid as the values are fitted to a Gran titration plot using least squares [START_REF] Gran | Determination of the equivalence point in potentiometric titrations of seawater with hydrochloric acid[END_REF]. An uncertainty of ±2-3 mol A T kg -1 was found as the combined precision for all cruises.

Samples for DIC were not filtered in order to prevent gas exchange with the atmosphere. Any sediments or particles collected in the samples would likely settle to the bottom of the sample bottle before analysis and thus interference in the measurements should be minimal. All of the DIC data collected in the North Sea was analyzed by the coulometric method of [START_REF] Johnson | Coulometric total carbon dioxide analysis for marine studies: maximizing the performance of an automated gas extraction system and coulometric detector[END_REF], with an overall uncertainty of ±1-2 mol DIC kg -1 for all cruises. For quality control and calibration of both DIC and A T , Certified Reference Material (CRM), of known DIC and A T , acquired from Prof. Andrew Dickson (Scripps Institute of Oceanography, USA) were used.

In 2002 and 2005, discrete samples were taken from the CTD and analyzed for pH potentiometrically on the total hydrogen ion scale, calibrated using TRIS (2-amino-2hydroxymethyl-1,3-propanediol) and AMP (2-aminopyridine) buffers prepared at a salinity of 35 according to [START_REF] Dickson | pH buffers for sea water media based on the total hydrogen ion concentration scale[END_REF]. The accuracy of this method is strongly dependent on the sensitivity of the pH meter and how accurately the two buffer solutions are made up. For the measurements used in this work a precision of ±0.001 pH units was attained. For all North Sea cruises pCO 2 was measured continuously, underway, from surface waters, providing the greatest coverage of data. The instrument used continuous equilibration with a headspace and infrared detection with an accuracy of ±1 µatm [START_REF] Körtzinger | At-sea intercomparison of two newly designed underway pCO 2 systemsencouraging results[END_REF]. The calibration of pCO 2 measurements was carried out regularly (every 2-3 hours) using National Oceanic and Atmospheric Administration (NOAA, Global Monitoring Division, Carbon Cycle Greenhouse Gases Group) calibrated gas standards. For all cruises a minimum of 3 calibration gases were used, spanning a concentration range from 0 ppm to 750 ppm.

The concentrations of silicate and phosphate are also required to account for their

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

contribution to A T , which can be of the order of up to 1 mol A T kg -1 in the North Sea. The inorganic nutrients nitrate (NO 3 ), ammonium (NH 4 ), dissolved silicate (H 4 SiO 4 ), and phosphate (PO 4 ), were analyzed following the methods of Grasshof et al. (1983). Reference nutrient samples (RMNS Kanso, Japan) were analyzed in every run as a control, giving a precision of ±0.2, 0.1, 0.2 and 0.01 for the four nutrients, respectively, for each cruise. The salinity and temperature values were recorded by CTD (Conductivity-Temperature-Depth) and salinity values were calibrated postcruise to salinity samples taken onboard.

Calculations

The internal consistency was assessed using the measured in-situ salinity, temperature, silicate, phosphate and pressure values. The KSO 4 of [START_REF] Dickson | Standard potential of the reaction: AgCl(s) + 1/2H 2 (g) = Ag(s) + HCl(aq), and the standard acidity constant of the ion HSO 4 in synthetic seawater from 273.15 to 318[END_REF] is used and pH values are given throughout on the total hydrogen ion scale. Not all station positions were occupied and measured in all four seasonal cruises that are compared here. In order to overcome this, the percentage of data points considered internally consistent is often used, and given in tables 3 and 4, to allow direct comparison between the seasons without the loss of data. The calculations were performed using the CO2SYS program [START_REF] Lewis | Program developed for CO 2 system calculations[END_REF] [START_REF] Hansson | A new set of acidity constants for carbonic acid and boric acid in sea water[END_REF], refit by [START_REF] Dickson | A comparison of the equilibrium constants for the dissociation of carbonic acid in seawater media[END_REF](H), [START_REF] Mehrbach | Measurement of the apparent dissociation constants of carbonic acid in seawater at atmospheric pressure[END_REF], refit by [START_REF] Dickson | A comparison of the equilibrium constants for the dissociation of carbonic acid in seawater media[END_REF](M), Hansson & Mehrbach (1973), refit by [START_REF] Dickson | A comparison of the equilibrium constants for the dissociation of carbonic acid in seawater media[END_REF](HM), Mojico-Prieto & Millero (2002)(MM) and Millero et al.

(2006)(ML). In the remaining text, for simplicity, the constants will be referred to by the abbreviations (R, GP, H, M, HM, MM and ML). The parameterizations of these constants, along with the medium and salinity ranges for which they were determined, are shown in Table 1.

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

Assessment

The internal consistency of the carbonate system can be considered in a number of different ways. Here, we assess how accurate and precise it is by comparing measured parameters to calculated parameters. Comparing the measured and calculated values, an accurate but imprecise system would be indicated by an average difference (∆)

close to 0, but with a high level of variance, denoted by the standard deviation (). A precise but inaccurate system, however, is indicated by a small , but with a bias, indicated by an average ∆ varying from 0. Thus, an inaccurate and imprecise system would have a high standard deviation and a large average ∆, whilst an accurate and precise system would show the inverse. The measured and calculated values of each parameter, for each set of apparent dissociation constants, was fitted to the form Y = a.X + b, by means of total least squares using Matlab 2012®, where X is the calculated value and Y is the observed value. By fitting the data using total least squares, we take into account the errors in both the x and y variables.

Error Propagation

When calculating the carbonate parameters in CO2SYS, the uncertainties associated with the observations and calculations are carried through (propagated). Errors are introduced via the measurement techniques of the input parameters, the measurement techniques of the output parameters against which we compare our calculations, and the determinations of the K 1 and K 2 values of carbonic acid. The measurement errors carried through in the calculations depend on the input parameter combinations. These margins vary insignificantly over the range of conditions (e.g. temperature, salinity) observed in the North Sea. The error propagation for each calculated parameter, with all variations of input parameter, was calculated following equation 5.

(z) 2 = {z/x} 2 .(x) 2 + {z/y} 2 .(y) 2 (equation 5)

The partial derivatives were determined seasonally, using average temperature and salinity values, and determined over the range of values observed. The determined errors for each set of input parameters for each season are given in brackets in the

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
standard deviation section of table 2a and 2b. The error calculation performed here does not take into account errors in the acid dissociation constants, as we later compare the performance of these constants. When the term 'internally consistent' is used, it refers to whether the observed values and associated measurement errors fit within the limits of the calculated propagated error.

Total Borate Calculations

To independently calculate the total borate contribution from each water mass, the individual water mass fractions were first calculated, following [START_REF] Friis | The salinity normalization of marine inorganic carbon chemistry data[END_REF]. To calculate the Elbe contribution we used the salinity and A T data from the stations closest to the German Bight (at latitudes >54.4°N and longitudes >8.0°E) to form a regression line with the Central North Sea end member. We then extrapolated this to a salinity of 0 to get a representative freshwater end member. This assumes that all the fresh water present comes from the same source, however, due to other rivers in the region having a similar catchment lithology we assume that the chemical signal of the water will be similar. This resulted in a 0-salinity end member of 2616 µmol A T kg -1 in February 2002, and 1678 µmol A T kg -1 in August 2005. These values comparing well with recent seasonal A T data from the Elbe [START_REF] Amann | Carbon dynamics in the freshwater part of the Elbe estuary, Germany: Implications of improving water quality[END_REF], which reports A T values of ~2500 µmol A T kg -1 in winter and ~1750 µmol A T kg -1 in summer at salinities <1, thus we assume our estimation is sufficient for the purpose of the calculation.

The borate alkalinity (BAlk) was then calculated for the water mass end members from total boron concentrations. The total boron concentrations (BT) were calculated from (a) the equations of Barth (1998) (equations 6 and 7) for the Elbe (German Bight) salinity-boron relationship, (b) [START_REF] Gripenberg | On the Alkalinity of Baltic Waters[END_REF] (equation 8) for the Baltic salinity-boron ratio, (c) [START_REF] Uppström | The boron/chlorinity ratio of deep-sea water from the Pacific Ocean[END_REF] for the North Atlantic ratio and also the revised North Atlantic ratio of Lee et al. (2010). At a salinity of 35 the two seasonal linear relationships of Barth (1998), for summer (September, 1994) and winter (February, 1995), yield an approximate difference of 6 μmol kg -1 in the calculated The water mass fractions were then multiplied by the theoretical BAlk contribution and summed to obtain the actual BAlk of the sample. All of the following equations (6 -9) are given for total borate concentration (BT) in mol kg -1 seawater.

BT = (0.000127/10.811) . S + 1.66e -5 (for September, 1994) (equation 6) BT = (0.000136/10.811) . S + 1.57e -5 (for February, 1995) (equation 7)

The Baltic Sea has demonstrated a number of chemical relationships varying from standard ocean behaviour. The chlorinity-boron relationship has been documented by [START_REF] Gripenberg | On the Alkalinity of Baltic Waters[END_REF]), Dyrssen and Uppstrom (1974), and Kremling, (1970), who all noted a positive anomaly in the boron/chlorinity ratio. There is little variation in their findings, with the ratio varying from 0.22 (Kremling, 1970) to 0.244 [START_REF] Gripenberg | On the Alkalinity of Baltic Waters[END_REF]. Here, we use the following equation from the data presented by Gripenberg 2a and2b, respectively. For average conditions in the North Sea, the HCO 3 -and CO 3 2-species were by far the greatest contributors to A T supplying ~2212 mol A T kg -1 , followed by borate (~68 mol A T kg -1 ), hydroxide (~3 mol A T kg -1 ), phosphate (~0.07 mol A T kg -1 ) and silicate (~0.04 mol A T kg -1 ). The largest changes between the summer and winter situation were in carbonate alkalinity (HCO 3 -+ 2*CO 3 2-

), which was significantly lower in August, when photosynthesis was ongoing. The contribution of silicate and phosphate were also much smaller then, due to lower concentrations resulting from uptake by phytoplankton. The borate contribution remains the second largest contributor to A T for both seasons, but was slightly higher in August. Isolating the various contributors to this change, we find that the increases in temperature, salinity and pH from February to August all affect an individual increase of ~3 mol kg -1 in B(OH) 4 -(total increase ~9 mol kg -1 ).

The change in average pCO 2 , from February 2002 to August 2005, was just 8 atm, however, this is not representative of the entire North Sea, as the pCO 2 range was two times larger in August (287 -455 µatm) compared to February (315 -396 µatm). The range reflected the two biogeochemical regimes within the North Sea: the stratified

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
nNS, where stratification enables prolonged phytoplankton growth causing a large drawdown of CO 2 , shows very low pCO 2 values in August, whereas in the sNS, very low pCO 2 values are observed only in spring, and a combination of remineralization and higher temperatures lead to much higher pCO 2 in August [START_REF] Bozec | The continental shelf pump for CO 2 in the North Sea -evidence from summer observation[END_REF][START_REF] Thomas | Enhanced Open Ocean Storage of CO 2 from Shelf Sea Pumping[END_REF][START_REF] Thomas | Controls of the surface water partial pressure of CO 2 in the North Sea[END_REF][START_REF] Schiettecatte | High temporal coverage of carbon dioxide measurements in the Southern Bight of the North Sea[END_REF]). In February, the water column is well-mixed throughout, and biological activity is lower thus leading to a more homogenous surface distribution of pCO 2 .

The results from all internal consistency calculations using the 4 carbonate parameters, with various dissociation constants, in the February 2002 and August 2005 data sets are shown in Tables 2a and2b, respectively. A regression slope of >1

indicates that the measured value is larger than the calculated values. A regression slope of 1.000 is achieved with several combinations in both seasons. Overall, in February, when the three parameters involved are DIC, A T and pCO 2 , (e.g.

DIC(A T ,pCO 2 ), A T (DIC,pCO 2 ) or pCO 2 (DIC,A T )) we find that the constants of M are the most suitable, whereas when DIC, A T and pH are involved the constants of GP outperform the others. The largest deviations are consistently observed when the input parameters are pH and pCO 2 , with MM producing the most accurate and precise reproduction of measurements. We also observe that the observed standard deviations of differences are largely comparable with the propagation error, with the exception of pH, in which we observe notably larger discrepancies than can be accounted for from the propagation error. In August the constants of ML are the best choice when the parameter combination is A T , DIC and pCO 2 and, additionally, when A T , DIC and pH are all involved. However, for calculations using both pH and pCO 2 the constants of GP and H are the most accurate.

In February and August, the variability in the standard deviations between the dissociation constants (e.g. comparing  A T , R with  A T , ML ) were relatively minor compared to the differences obtained with various input parameter combinations (e.g.

comparing  A T , (pH,pCO2) with  A T , (DIC,pH) etc.). As such, beyond assessing the parameter combinations with the largest degree of noise, the standard deviations are

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

not further discussed in this section. The results for both February and August are consistent, showing that the best regression coefficients for A T were attained with the inputs (DIC,pCO 2 ), for DIC it was (A T ,pCO 2 ) and for pCO 2 it was (DIC,A T ). For A T and DIC, this is in line with smaller error propagations, as similarly found by Dickson and Riley (1974), however, the much improved result for calculating pCO 2 using (DIC,A T ) is indicative of the improvements in measurement accuracy obtained with these parameters with the use of CRMs. When calculating pH, a regression coefficient of 1.000 was attained with all combination parameters. Similarly, consistency was found in the parameter input combinations resulting in the greatest precision, with A T being best represented by (DIC,pH), DIC by (A T ,pH), pH by (A T ,pCO 2 ) and pCO 2 by (DIC,pH). The combination pair resulting in the lowest residuals varied between the two seasons, possibly indicating a bias in some of the measurements. The prevalence of better results being obtained when using pH as an input parameter, indicate that these measurements are reliable, however, the difference in the resulting residuals from e.g. DIC(A T ,pH) and DIC(A T ,pCO 2 ), are still directly comparable.

The results from February show a slight overall decrease in precision, represented by a higher standard deviation, compared to that obtained in August, however, the average differences show a similar level of accuracy for DIC and pCO 2 and a mildly lower bias in February for A T and pH. The standard deviations of the differences between the calculated and measured parameters () for February 2002 and August 2005 (shown in brackets, respectively) are approximately double to those found in the open ocean by [START_REF] Lee | The carbon dioxide system in the Atlantic Ocean[END_REF], of ±3 mol A T kg -1 in A T (±7, 6 mol kg -1 ), ±3 mol DIC kg -1 in DIC (±7, 6 mol kg -1 ), and ±1.3% in pCO 2 (±10, 8 atm; ±3%, 2%).

The effects of salinity on dissociation constants

Different salinity and temperature parameterizations used in the determination of pK 1 and pK 2 (Table 1) are largely responsible for the deviations in the calculations between dissociation constants. Plotting the difference in ∆A T between results from different dissociation constants (e.g. the differences between the residuals using the constants of R and the residuals using the constants of M: ∆A T R -∆A T M ; or ∆DIC, ∆pH

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

and ∆pCO 2 ) shows an exponential-like increase, or decrease, around a salinity of ~34.5 (Fig. 3), indicating that at higher salinities the choice of dissociation constants becomes more significant. We compared the differences in residuals of all dissociation constants to those of M, and found the same relationship to that shown in Fig. 3 for the constants of R, GP, H and HM. The differences in ∆A T between M and MM showed a gradual linear increase with salinity, with a slope of 0.3 (R 2 =0.9, p-value<0.01), indicating that across the range of 5 salinity units, there is an increase of 1.5 mol A T kg -1 difference in ∆A T . The residuals resulting from the use of constants . This ratio is used directly in the calculation of pH and A T from DIC and pCO 2 . Due to the different parameterizations of K 1 and K 2 for the different constants, the difference in this ratio between constants may become larger or smaller at lower salinities. This will lead to larger or smaller differences between the measured and calculated parameters. As such, at lower salinities, one would expect, overall, lower variation between the calculations from the different dissociation constants, as is observed here.

To investigate the performance of the six dissociation constants across a range of different salinities, the North Sea data from August ('05) and February ('02) were split into two groups; low salinity waters, with salinities less than 34.5, and high salinity waters, with salinities greater than or equal to 34.5. The 1 standard deviation (σ) of ∆A T , ∆DIC, ∆pH and ∆pCO 2 is larger, for all parameters and in both seasons, for low salinity waters compared to high salinity waters. In February, in high salinity waters, the standard deviation of ∆A T is 2 μmol A T kg -1 lower than that of the total dataset; ∆DIC is 1 μmol DIC kg -1 lower, ∆pH is 0.002 lower and ∆pCO 2 is 2 μatm lower. In the less saline waters there is an increase in the standard deviation of ∆A T , ∆DIC, ∆pH and ∆pCO 2 by 2 μmol A T kg -1 , 1 μmol DIC kg -1 , 0.003, and 2 μatm, respectively. The same pattern is observed in August, however the differences are marginally smaller.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

Furthermore, in February, there is a clear dichotomy between the different dissociation constants in low and high salinity waters (Fig. 4). As previously mentioned, in winter, when DIC, A T , and pH are involved in calculations (that is to say, the calculation of pH(DIC,A T ) or DIC(A T ,pH) or A T (DIC,pH)), the constants of GP produce the least deviation between measured and calculated parameters.

However, we find that in high salinity (S>34.5) waters the constants of M become more appropriate and in the lower salinity waters (S<34.5) the constants of R outperform the others. Despite these changes, the constants of MM and ML show little variation in performance between the two salinity environments. For August, there was no difference between the most suitable constants for the entire dataset compared to the different salinity regimes regarding measurements with DIC, A T and pH, with the constants of GP still most closely representing the measurements. With regard to measurements of pH, pCO 2 and DIC, however, in the lower salinity waters the constants of H outperformed those of GP. Similarly to February, the constants of MM and ML also perform well and show the best consistency across the entire salinity range, with the constants of ML generally displaying a greater level of accuracy. The salinity range does not change between the two seasons, however, in February, 2002 a 3.5 fold higher Elbe river discharge was recorded compared to August, 2005 (GKSS, 2012, http://coast.gkss.de/staff/kappenberg/elbe/abfluss/elbe.abfluss). A larger proportion of freshwater in the North Sea could explain why a greater difference is observed between the two salinity environments during this time.

From this application of the various dissociation constants of carbonic acid, we show that the range of salinity directly affects the robustness of carbonate system calculations with different constants. Specifically, the constants R, H, GP and HM show much larger errors in mediums with a larger freshwater constituent compared to those of MM and ML. The larger range of salinities for which these latter constants were determined is likely to contribute to more accurate determinations of K 1 and K 2 , thus also leading to more accurate K 1 /K 2 ratios. Based on these results, for datasets that cover a large range of salinities, for consistency, one can achieve the best overall results using the constants of ML in carbonate system calculations. Despite this, and bearing in mind that the current recommended dissociation constants for carbonate

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

system measurements are those of M, all further calculations and comparisons are performed using this set of constants.

The Role of Borate Alkalinity

The ratio of K 1 /K 2 varies between different parameterizations of carbonic acid dissociation constants, partly depending on whether the constants were determined in artificial seawater (ASW) or natural seawater (NSW). In NSW interactions between borate and carbonate cause an increase in K 1 and a decrease in K 2 (Mojico-Prieto and Millero, 2002), which does not occur in ASW due to the lack of boric acid. As such, we suggest that finding the constants of R, GP and H to perform comparably with those of M in less saline waters may be indicative of the lower borate concentrations in freshwater and thus the errors in pK 2, arising from the lack of boric acid in their ASW medium, are minimized. Borate has a long residence time in the oceans leading to a conservative relationship with salinity, thus in carbonate system calculations the total boron contribution is calculated using the linear relationship to chlorinity as found by [START_REF] Uppström | The boron/chlorinity ratio of deep-sea water from the Pacific Ocean[END_REF], and the borate then calculated using a thermodynamically derived K B . The data used by [START_REF] Uppström | The boron/chlorinity ratio of deep-sea water from the Pacific Ocean[END_REF] is based on data obtained from the Pacific Ocean. However, recent findings suggest that this relationship is not consistent in all the worlds oceans (Lee et al., 2010). Lee et al. (2010) recently found that this relationship is not as ubiquitous in the oceans as originally thought and has made adjustments which, on provisional tests, led to an increase in the internal consistency of carbonate system calculations in the North Atlantic. The [START_REF] Uppström | The boron/chlorinity ratio of deep-sea water from the Pacific Ocean[END_REF] In the North Sea, the Baltic input through the Skagerrak, and the Elbe input in the German Bight, have also revealed non-conventional salinity-borate relationships [START_REF] Kremling | Untersuchungen ueber die chemische Zusammensetzung des Meerwassers aus der Ostsee[END_REF]Barth, 1998). To evaluate the potential error introduced into calculations of the carbonate system in the North Sea by unaccounted-for borate, the borate concentrations in the Elbe-and Baltic-influenced areas of the North Sea were calculated using the water mass fractions and borate-salinity equations from the literature (see section 2.5). This 'pre-calculated' borate contribution to alkalinity (BAlk pre ) was then implemented into CO2SYS, replacing the normal salinity-borate calculation [START_REF] Uppström | The boron/chlorinity ratio of deep-sea water from the Pacific Ocean[END_REF]. For comparative purposes BAlk was also calculated using the traditional [START_REF] Uppström | The boron/chlorinity ratio of deep-sea water from the Pacific Ocean[END_REF] equation but using various [H + ] determinations, including that from direct pH measurements, using DIC and A T as input parameters and using A T and pCO 2 as input parameters (DIC/pCO 2 produce the same results as A T /pCO 2 , hence have not been included). The BAlk pre calculations showed little change from the BAlk concentrations, the difference showing a normal distribution with 0 as the mean of the absolute values and  of 3 μmol kg -1 . The equations of Barth (1998) are very similar to the oceanic relationships, thus no significant difference was observed in the German Bight between the two calculations. The greatest difference occurs in the Skagerrak where the [START_REF] Uppström | The boron/chlorinity ratio of deep-sea water from the Pacific Ocean[END_REF] ratio underestimates the BAlk compared to the equation derived from

Gripenberg's data. The differences in calculation resulted in a maximum difference in BAlk of 8 μmol kg -1 in August 2005. Despite these differences, using the new calculated BAlk pre values in CO2SYS did not improve the internal consistency; on the contrary the changes caused greater deviations between the measured and calculated parameters, reducing the number of internally consistent stations for all input combinations. Thus we conclude that the non-conservative salinity-borate relationships in the North Sea are not limiting factors for the internal consistency.

When the [H + ] from pH measurements is used to calculate BAlk directly, and this

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

value is implemented into the thermodynamic calculations, the overall internal consistency is again unimproved. On the other hand, when the BAlk calculated using

Seasonal comparison (A T , DIC, and pCO 2 )

Seasonal Internal Consistency

The results from the internal consistency study from all 6 data sets, using just three carbonate parameters (A T , DIC and pCO 2 ) show a clear dominance of the constants of M and ML outperforming the others. Both of these sets of constants show comparable results. The results of the seasonal calculations with the use of M constants are given in Table 3. The statistics indicate that the internal consistency of the North Sea follows a seasonal pattern, with Maythe most biologically productive time periodproducing less consistent results and the rest of the year displaying a good level of consistency. These findings are in agreement with those of [START_REF] Koeve | Potential impact of DOC accumulation on fCO 2 and carbonate ion computations in ocean acidification experiments[END_REF], who similarly noted significant inconsistencies in the carbonate system in waters of high

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT biological activity.
Additionally, the largest fraction of Baltic water found in the Skagerrak is present in May (17%). Further analysis shows that during this season there is a significant (p=0.008, R 2 =0.55), positive correlation between the Baltic fraction and the deviation between calculated and measured parameters (Figure 5(b)). Despite this, in May, the internal consistency could also not be improved by changing the BAlk calculation as done above (see section 3.1.2). However, [START_REF] Gripenberg | On the Alkalinity of Baltic Waters[END_REF] concluded, from back titrations of alkalinity, that whilst borate did not directly affect the accuracy of the A T titrations, the interactions between boric acid and organic material could lead to underestimated alkalinity values. Recent work by [START_REF] Ulfsbo | Modelling organic alkalinity in the Baltic Sea using a Humic-Pitzer approach[END_REF] similarly find that in the Baltic, dissolved organic matter can contribute up to 30 mol A T kg -1 to alkalinity, thus rendering the A T equation in CO2SYS invalid. It thus seems likely that the lack of internal consistency observed in Baltic-influenced waters in May stems from additional organic material in the water at this time of year, with possibly further errors introduced from the interaction of this organic material with boric acid.

The North Sea receives a high input of allochthonous dissolved organic matter, for which analysis of colored dissolved organic matter (CDOM) indicates that it contains a high proportion of terrigenous sources (Stedmon et al., 2009). Thus may contain terrestrial "titratable" compounds, which, are also not accounted for in our calculations. It has been suggested that fatty acids and humic substances may contribute to measured A T and in the Baltic Sea there is a higher than average concentration of both. The contribution of fatty acids to A T has been fairly extensively examined by [START_REF] Osterroht | Extraction of dissolved fatty acids from sea water[END_REF] and maximum contributions of <1 mol A T kg -1 to alkalinity were found, which is insufficient to cause the inconsistencies observed here.

In contrast, more recently, the presence of humic substances has been shown to contribute to large excess A T in the Baltic Sea [START_REF] Ulfsbo | Modelling organic alkalinity in the Baltic Sea using a Humic-Pitzer approach[END_REF]. The contribution of organic matter to A T has been observed elsewhere, as well. [START_REF] Muller | Estimating the organic acid contribution to coastal seawater alkalinity by potentiometric titrations in a closed cell[END_REF], using the technique of Hernandez-Ayon et al. ( 2005) on fjord surface waters, compared the A T estimation of DIC to the measured value and found that organic

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
alkalinity can contribute from 2-22 μmol A T kg -1 in regions of high organic matter.

These errors are sufficiently large to account for the noted decrease in internal consistency observed in the North Sea in May, compared to other seasons.

As the North Sea consists of two different biogeochemical regimes, the nNS and the sNS, we investigated whether these two environments contributed towards a difference in the internal consistency of the carbonate system. The North Sea was subsequently divided along the 56°N latitude and the internal consistency reexamined in both areas, applying just the constants of M (Table 4). In May there is a significant increase in the standard deviation of ∆DIC and ∆pCO 2 in the nNS compared to the sNS, however, ∆A T shows no difference. As previously mentioned there is greater freshwater influence in the sNS. However, the difference observed between the nNS and sNS is larger than the difference previously mentioned between oceanic and estuarine waters. In the nNS larger blooms of phytoplankton are maintained for long periods due to seasonal, thermal stratification, which is in contrast to the sNS where rapid remineralization takes place in the permanently well-mixed water column. Errors have previously been associated with intense biological activity [START_REF] Kim | Significant contribution of dissolved organic matter to seawater alkalinity[END_REF][START_REF] Hoppe | Implications of observed inconsistencies in carbonate chemistry measurements for ocean acidification studies[END_REF]. Thus we are inclined to suggest that this may have a similar effect here. However, it is interesting that the errors in A T are not affected on a north-south latitudinal gradient, like pCO 2 and DIC.

The different regimes in the nNS and sNS also influence the dominant phytoplankton species present in each location. Calcifying organisms and the presence of coccolithophores are well documented in the nNS, where the development of coccolithophore blooms is controlled by the North Atlantic influence and the stratification of the water column [START_REF] Houghton | Coccolith sedimentation and transport in the North Sea[END_REF][START_REF] Burkill | Dimethyl sulphide biogeochemistry within a coccolithophore bloom (DISCO): an overview[END_REF]. The presence of coccolithophores in waters sampled for DIC could conceivably lead to a degree of CaCO 3 dissolution during the DIC titration. This would introduce a degree of noise to the DIC values in the NNS, and possibly account for some of the larger errors observed in this region. Furthermore, [START_REF] Fraga | On the variation of alkalinity during phytoplankton photosynthesis[END_REF] quantified the change in A T during photosynthesis and found that the differences are

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

caused predominantly by CaCO 3 formation. Calcification influences the A T /DIC (X) ratio of the waters as it removes one mole of DIC per 2 moles of A T [START_REF] Wolf-Gladrow | Total alkalinity: The explicit conservative expression and its application to biogeochemical processes[END_REF]. [START_REF] Lee | The reliability of the thermodynamic constants for the dissociation of carbonic acid in seawater[END_REF][START_REF] Lee | The carbon dioxide system in the Atlantic Ocean[END_REF] found that high X ratios (>1.10) can impact the sensitivity of pK 2 , thus affecting the consistency between measured and calculated carbonate parameters. The range of X (A T /DIC) found in the North Sea varies from 1.06 -1.16 throughout the year, with the smallest range (1.06 -1.09) occurring in February and the largest range (1.07 -1.16) found in May. During May, X also shows a positive linear correlation with ∆DIC, ∆A T and ∆pCO 2 with R 2 values of 0.6, 0.6 and 0.7, respectively (Figure 6). The size of the error introduced by varying X can be up to 5 mol A T kg -1 , as found by [START_REF] Lee | The reliability of the thermodynamic constants for the dissociation of carbonic acid in seawater[END_REF], which would contribute a large proportion of the errors we observe in the North Sea.

Conclusions

On a basin-wide scale in the North Sea we find that the carbonic acid dissociation constants of Mehrbach et al., (1973, refit by Dickson and[START_REF] Dickson | A comparison of the equilibrium constants for the dissociation of carbonic acid in seawater media[END_REF], and the constants of [START_REF] Millero | Dissociation constants of carbonic acid in seawater as a function of salinity and temperature[END_REF] provide the best internal consistency. In the case of the former constants, this is in agreement with previous findings [START_REF] Lee | The reliability of the thermodynamic constants for the dissociation of carbonic acid in seawater[END_REF][START_REF] Wanninkhof | The optimal carbonate dissociation constants for determining surface water pCO 2 from alkalinity and total inorganic carbon[END_REF]. In regions of lower than average ocean salinity, there is markedly less deviation between the internal consistency of the different constants and the constants of ML outperform those of M. At high salinities these two sets of constants are comparable in performance, thus over datasets with a wide variety of salinities, we recommend the use of ML constants for more consistent results.

In the North Sea in general, we find that the internal consistency of the carbonate system is much less than that found in the open ocean, by approximately 2-fold, which is unsurprising considering the additional input of terrestrial compounds in coastal areas, and generally higher productivity. Application of the new boratesalinity relationship noted by Lee et al. (2010), and further accounting for the nonconventional salinity-borate relationships in the water masses contributing to the North Sea, did not result in any improvement in the internal consistency. The internal consistency is relatively constant throughout the year, with the exception of May,

A C C E P T E D M A N U S C R I P T
where there is a marked decrease in our ability to accurately calculate the carbonate system. Waters in which the Baltic Sea constituent is larger than 2% show lower internal consistency in May driven by the presence of excess organic alkalinity unaccounted for in carbonate system calculations. In the entire North Sea there is a significant correlation between the internal consistency and X (A T /DIC), similar to that found in previous studies. Higher X values in the nNS compared to the sNS are likely to contribute to the larger errors observed in this area through the increased sensitivity on pK 2 .

A [START_REF] Hansson | A new set of acidity constants for carbonic acid and boric acid in sea water[END_REF], refit by [START_REF] Dickson | A comparison of the equilibrium constants for the dissociation of carbonic acid in seawater media[END_REF](H), [START_REF] Mehrbach | Measurement of the apparent dissociation constants of carbonic acid in seawater at atmospheric pressure[END_REF], refit by [START_REF] Dickson | A comparison of the equilibrium constants for the dissociation of carbonic acid in seawater media[END_REF](M), Hansson & Mehrbach (1973), refit by [START_REF] Dickson | A comparison of the equilibrium constants for the dissociation of carbonic acid in seawater media[END_REF] Table 3. Derived linear regressions, using total least squares, and standard deviations of differences between measured and calculated parameters (A T , DIC, and pCO 2 ) with the percentage of internally consistent stations for all six cruises. The calculations were performed using the acid dissociation constants of [START_REF] Mehrbach | Measurement of the apparent dissociation constants of carbonic acid in seawater at atmospheric pressure[END_REF], refit by [START_REF] Dickson | A comparison of the equilibrium constants for the dissociation of carbonic acid in seawater media[END_REF]. Figure 3. The difference in ∆A T (µmol A T kg -1 ) across the measured salinity range.

The ∆A T R and ∆A T M are calculated using (DIC,pCO 2 ) as input parameters.

A C C E P T E D

M A N U S C R I P T A C C E P T E D M A N U S C R I P T -The consistency of the carbonate system follows a seasonal pattern, with the largest errors in May.

ACCEPTED MANUSCRIPT

ACCEPTED MANUSCRIPT

-The carbonic acid dissociation constants of [START_REF] Millero | Dissociation constants of carbonic acid in seawater as a function of salinity and temperature[END_REF] are most consistent across a large salinity range (0-36).

  the cruises taking place in August 2001, November 2001, February 2002, May 2002 and August 2005 the samples for A T were filtered on GF/F shortly after sampling, prior to analysis. The resultant acid volumes and corresponding e.m.f

[

  such, for our calculations, equation 6 was used for the calculations of A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT August 2005 (summer), and equation 7 for the calculation of February 2002 (winter).

  ML showed no discernible difference to those of M over the salinity range investigated. At lower salinities the ratio of K 1 /K 2 is larger[START_REF] Mosley | Modelling of pH and inorganic carbon speciation in estuaries using the composition of the river and seawater end members[END_REF] Zeebe and Wolf-Gladrow, 2001) leading to relatively greater HCO 3 -concentrations compared to CO 2 and CO 3 2-

  Figure 2a and 2b clearly show that borate is the second largest non-carbon contributor to A T , constituting on average ~4% (Zeebe and Wolf-Gladrow, 2001), thus deviations from the true borate concentration can potentially cause substantial errors in calculations involving A T . Additionally, the K B of the boric acid system is very close to ocean pH, thus any miscalculation of [H + ] can also lead to significant errors in the calculation of [B(OH) 3 ] and [B(

  borate-salinity relationship was substituted in CO2SYS for the new Lee et al. (2010) relationship, to see if this improved the internal consistency. The new calculations led to a maximum difference A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT in borate alkalinity of 4 µmol A T kg -1 , during both seasons. Despite this change, the influence on internal consistency was almost imperceptible, and resulted in no significant improvement.

  concentration that contributes significantly to a lack of internal consistency for (A T , pCO 2 ) and (DIC, pCO 2 ). This can be illustrated by Figure5(a), where the change in A T between the measured and calculated values are plotted against the change in [H + ] between that calculated from (DIC, A T ) and (DIC, pCO 2 ). When the [H + ] from (DIC, A T ) is used to calculate BAlk the A T values are shifted towards the x-axis (0), indicative of the improvement in internal consistency. The reason for such a bias is unknown, and may be caused by an undetected bias in one of the measured parameters.

  ., 1998. Evaluation and comparison of the global carbon cycle in the coastal zone and in the open ocean. In: Brink, K. H., Robinson, A. R. (Eds.), The Sea. John Wiley & Sons, New York, pp. 213-252. Wootton, J. T., Pfister, C. A. and Forester, J. D., 2008. Dynamic patterns and ecological impacts of declining ocean pH in a high-resolution multi-year dataset. P.Natl Acad. Sci. 105, 18,848-18,853. Zeebe, R. E. and Wolf-Gladrow, D., 2001. CO 2 in Seawater: Equilibrium, Kinetics, Isotopes. Elsevier, Amsterdam.

  (HM), Mojico-Prieto & Millero (2002)(MM) and Millero et al. (2006)(ML).

Figure Captions Figure 1 .

 Captions1 Figure Captions

Figure 2 .

 2 Figure 2. The distribution of total alkalinity in the North Sea and respective contributions to total alkalinity by (i) carbonate species, (ii) borate, (iii) silicate, (iv) phosphate and (v) ammonia in February 2002 (panel a) and in August 2005 (panel b).

Figure 4 .

 4 Figure 4. Radar plot showing (a) the regression coefficients for A T (DIC,pCO 2 ) and (b) pCO 2 (DIC,A T ) in different salinity environments, and their respective standard deviation (c), (d) of the fit.

Figure 5 .

 5 Figure 5. (a) The change in A T caused by the difference in BAlk calculated using the [H + ] from (DIC, A T )(red) and (DIC, pCO 2 )(black). (b) Regression of DIC with the calculated water mass fraction of Baltic water, selectively displaying those with a fraction of 2% or higher.

Figure 6 .

 6 Figure 6. The change in DIC and A T (both in µmol kg -1 ) and pCO 2 (in µatm) observed with a varying ratio of A T /DIC (X) during May 2002.

Figure 1 A-

 1 Figure 1

  The sampling program of the North Sea is shown in Figure1, with the different distinct, biogeochemical regions labelled. The two seasons in which data was collected of all four measurable carbonate parameters were winter (February) and summer (August). These two data sets allow us to compare the internal consistency in two very different conditions, in terms of values of the carbonate variables, and water temperature. As the North Sea is a temperate coastal sea, it is characterized by high primary production during summer and very low, or no, production in winter. The seasonal cycle is fully assessed in the works of[START_REF] Bozec | Assessment of the processes controlling seasonal variations of dissolved inorganic carbon in the North Sea[END_REF] andThomas et al.

		Results and Discussion
	3.1	A seasonal comparison
	(1960):
		BT = (0.000244/10.811) . (S/1.80655)	(equation
	8)
	For the North Atlantic fraction the traditional seawater relationship of Uppström
	(1974) was used:
		BT = (0.000232/10.811) . (S/1.80655)	(equation 9)
	Furthermore, Lee et al. (2010) recently found boron-salinity relationships
	significantly different to those reported by Uppström (1974), which were based on
	measurements solely from the Pacific Ocean. Here, we also apply this new boron-
	salinity equation (equation 10), to test whether this improves our internal consistency
	in the North Sea.
		BT= (0.0001336/10.811) . S	(equation 10)

(2009)

, thus here the discussion is kept brief. The average DIC in February 2002 was 64 mol DIC kg -1 higher than that in August. The average A T changed by just 2 mol A T kg -1 , however, the various contributions of chemical species to A T showed larger variability. The changes in A T contributions in the North Sea from February to August are shown in Figures

  [H + ] from (DIC, A T ) is implemented, the internal consistency improves notably if DIC or A T are being calculated, for both seasons. The number of internally consistent DIC, A T ) is only used to calculated the BAlk, whereas the [H + ] used in the other CO2SYS calculations is calculated from the stated input parameters. As such, we can conclude that it is the effect of miscalculation of [H + ] on the BAlk

	(IC) stations in February and August, for A T (DIC,pH) and DIC(A T ,pH), increased by
	11% (7 stations) and 36% (29 stations). The internal consistency of A T (DIC,pCO 2 )
	increased by 8% and 19% for February and August, respectively, and DIC(A T ,pCO 2 )
	increased by 6% and 16% stations, respectively. It should be noted that the [H + ]
	calculated from (

Table Captions Table 1 .

 Captions1 Table of parameterizations of the carbonic acid dissociation constants, with temperature and salinity range for each set, and the medium in which measurements were performed.

	Table 2. Calculated statistics for all calculated values of parameters compared to the
	measured values during February 2002 (a), and for August 2005 (b). Both datasets
	were tested using the dissociation constants of Roy et al. (1993)(R), Goyet & Poisson
	(1989)(GP),

Table 4 .

 4 Calculated internal consistency statistics for all data sets divided into

	northern North Sea stations (nNS)(latitude>56N) and southern North Sea stations
	(sNS)(latitude<56N).

Table 1 .

 1 Table of parameterizations of the carbonic acid dissociation constants, with temperature and salinity range for each set, and the medium in which measurements were performed.

	GP	pK 1 = 812.27/T + 3.356 -0.00171.S.ln(T) + 0.000091.S 2 pK 2 = 1450.87/T + 4.604 -0.00385.S.ln(T) + 0.000182.S 2	-1 -40	10 -50	Artificial SW
		pK 1 = 851.4/T + 3.237 -0.0106.S + 0.000105.S 2
	H	pK 2 = -3885.4/T + 125.844 -	2 -35	20 -40	Artificial SW
		18.141.ln(T) -0.0192.S +
		0.000132.S 2
		pK 1 = 3670.7/T -62.008 +
		9.7944.ln(T) -0.0118.S +
	M	0.000116.S 2	2 -35	26 -43	Natural SW
		pK 2 = 1394.7/T + 4.777 -0.0184.S + 0.000118.S 2
	HM	pK 1 = 845/T + 3.248 -0.0098.S + 0.000087.S 2 pK 2 = 1377.3/T + 4.824 -0.0185.S + 0.000122.S 2	2 -35	20 -40	Artificial SW
		pK 1 = -43.6977 -0.0129037.S + 1.364e -4 .S 2 + 2885.378/T +
		7.045159.ln(T)
	MM	pK 2 = -452.0940 + 13.142162.S -8.101e -4 .S 2 + 21263.61/T +	0 -45	5 -42	Natural SW
		68.483143.ln(T) + (-581.4428.S
		+ 0.259601.S 2 /T -
		1.967035.S.ln(T)

Acknowledgements

We thank the captains and crews of the Research Vessel Pelagia and an anonymous reviewer for constructive comments on a previous version of the manuscript. This program was supported by the Research Council for Earth and Life Sciences (ALW) of the Netherlands Organization for Scientific Research (NWO), project number 817.01.004 and contributes to CARBO-OCEAN, an integrated project of the European Union (contract no. 511176-2). This work contributes to IGBP/IHDP LOICZ. AVB is a senior research associate at the FRS-FNRS.

A C C E P T E D

M A N U S C R I P T

ACCEPTED MANUSCRIPT

(atm) A T ,pH (±3) 9 (11) 14 ( 11) 16 ( 11) 6 ( 10) 12 ( 11) 4 ( 10) 10 (10)


Olivier Kermorgant
email: kermorgant@unistra.fr

A dynamic simulator for underwater vehicle-manipulators

In this paper we present a dynamic simulator for intervention autonomous underwater vehicles. Prototyping and testing of such robots is often tedious and costly, and realistic simulation can greatly help validating several aspects of the project. In order to benefit from existing software, the presented system is integrated with ROS, through the Gazebo dynamic simulator, and the underwater image rendering UWSim. The whole approach allows realistic rendering of dynamic multirobot simulation, with contact physics, buoyancy, hydrodynamic damping and low-level PID control. This paper details the modeling choices that are done and exposes how to build its own AUV model. Integration with other ROS programs is exposed, and a simulation shows an example of behavior during a black box recovery mission.

Introduction

A strong trend in underwater robotics is the use of autonomous underwater vehicles instead of the classical Remotely Operated Vehicles and maned submersibles. Risk and cost are highly reduced, as it is never easy to deploy a team on a surface vessel or operators in a submersible. However, experimentation with AUV's is very difficult, because of the environment and the nature of the vehicle. Small experiments can be carried in water tanks, for example low-level control and basic prototyping, but this already requires space and resources. Higher-level experiments like navigation, waypoint following, seabed mapping or sensor-based control are designed to be carried in open environments, which involves high costs, human resources and are highly time consuming. Furthermore, in most autonomous underwater experiment the researchers do not have full knowledge of what is happening underwater, which is another difficulty during early development stages.

For these reasons, simulators have been developed in order to help prototyping AUV control and design. They allow having a full real-time access to all data during an experiment, and thus make it possible to greatly improve the AUV before going into the real experiments. A survey of AUV simulator has been carried in 2008 [START_REF] Matsebe | A review of virtual simulators for autonomous underwater vehicles (auvs)[END_REF]. Since then, the robotic community got used to the ROS framework [START_REF] Cousins | Welcome to ros topics[END_REF] which acquired its own AUV simulator, called UWSim [START_REF] Prats | An open source tool for simulation and supervision of underwater intervention missions[END_REF]. This simulator has been used extensively in the Trident project, and renders realistic images through OpenSceneGraph1 (OSG) and osgOcean2 . OSG is an open source 3D graphics application, while osgOcean was developped to render realistic underwater images in OSG. Fig. 1 gives an example of comparison between classical and underwater-oriented rendering. Embedded cameras, basic sonar and other AUV sensors are supported, and multiple underwater or surface vehicles can be present in the same simulation, which allows carrying experiments that would be very difficult in the real environment.

A drawback still present in the current state of the art is the absence of dynamic AUV simulation. Indeed, UWSim is only a kinematic simulator, with an external dynamic module that has to be coded in Matlab and only handles single-body vehicles. It is thus limited to kinematic control in the case of intervention AUV's (I-AUV), that carry a robotic arm. Several works have been using UWSim to show advanced whole-body control schemes [START_REF] Casalino | Agility for underwater floating manipulation: Task & subsystem priority based control strategy[END_REF][START_REF] Kermorgant | A global control scheme for free-floating vehicle-manipulators[END_REF], but the simulations lack realism. On the other hand, the ROS community is used to the Gazebo simulator [START_REF] Koenig | Design and use paradigms for gazebo, an open-source multi-robot simulator[END_REF], that was designed mostly for ground robots. This simulator handles dynamics, contact physics and is very versatile through its plugin-based design. A good example of Gazebo extension is the recent Hector Quadrotor package [START_REF] Meyer | Comprehensive simulation of quadrotor uavs using ros and gazebo[END_REF] that proposes quadrotor UAV simulation. However, as seen in Fig. 1a the rendering is of course far less satisfactory for underwater environments.

The purpose of this paper is to present an integration between Gazebo and UWSim, in order to achieve both dynamic and visually realistic I-AUV simulation. In this work we focus on the low-level dynamics, as higher-level control schemes or navigation tasks can already be implemented through the numerous ROS packages. The proposed simulator, called free-floating Gazebo 3 , aims to be flexible, letting each user specifying the characteristics of his own (I-)AUV. A typical scenario would be to use CFD software or water tank experiments to estimate the hydrodynamics coefficients of the vehicle. This parameters can then be used in the proposed simulator in order to allow real-time simulation of navigation or sensor-based missions, which leads to a good compromise between computation cost and realism compared to pure kinematics simulators.

The paper is organized as follows. In Section 2 the simulated dynamics are described together with their implementation in the robot model. Section 3 presents the low-level PID control that is included in the simulator, allowing to control the I-AUV in position, velocity or effort depending on the controller design. The integration with other programs such as high-level control is then exposed in Section 4. Finally, an experiment is detailed in Section 5, showing the recovery of a black box.

Description and simulation of a submerged object

The simulation is carried through two separate Gazebo plugins:

-A world plugin that simulates the overall buoyancy direction and water surface, together with the water current. The current may be varying but is assumed to be the same in the whole environment. -A model plugin that handles low-level control through thruster effort (vehicle body) or joint torque (embedded arm).

We first describe the overall simulated dynamics before going into implementation details.

Dynamic model

Gazebo performs dynamic simulation and takes as inputs the forces and torques that acts on a body. Through plugins, it is possible to add any kind of effort and thus to consider an hydrodynamic model. In this section we describe the dynamic model of a robot link and the approximations that are carried in the implementation.

We denote p the pose of a robot link expressed in the world frame, and v the velocity screw expressed in the local frame. In this case, the general form of the state-space dynamic model of the link yields [START_REF] Fossen | Guidance and control of ocean vehicles[END_REF]:

M v = G(v)v + d(v) + τ g + τ u + τ e (1) ṗ = Jv (2)
where the parameters are defined as:

-M is the inertia matrix. It is the sum of the rigid body inertia matrix M d and the added inertia matrix M a due to the fluid. Besides of the mass that can be measured directly, the values of M d and M a are usually determined empirically. The error on M a can be quite large and we will not consider this matrix in the model.

-G(v) represents the Coriolis and centrifugal forces. It is also the combination of the corresponding dynamic G d and hydrodynamic G a matrices, that are computed from the inertia matrices. For the same reason as M a , the matrix G a is considered null. -d(v) is the hydrodynamic damping force, or parasitic drag. Its values are also usually poorly known, but this effect cannot be ignored in the underwater case. -τ g is the combination of the gravity and the buoyancy forces and induced torques. In most vehicles, the center of gravity and the center of buoyancy are vertically aligned in order to produce a keel effect that stabilizes the roll of the vehicle. -τ u represents the forces and torques that are added by the user. In our case, we consider additional forces through several thrusters that can be placed on the AUV body. No other kind of motors (for example fins) are considered. For the links related to the robotic arm, the applied torque and forces correspond to the joint efforts.

-τ e regroups the other forces and torques applied to the robot link. This corresponds to the interactions with other robot links or through contact points with other objects, which are modeled by the Gazebo simulator. -J is the transformation matrix between the world frame and the link frame.

Because of the usual modeling uncertainties and the low velocity of the vehicle, the added inertia and Coriolis forces are not taken into account. All additional efforts related to the hydrodynamic effects are considered into the damping force d(v), which is a quadratic function of the velocity between the robot link and the water velocity:

d(v) = -k (v -J -1 v c ) (3)
where k expresses the damping coefficients and v c = (v c , 0) is the water velocity expressed in the world frame and has null angular values. The same modeling is used for non-actuated objects that may be present in the simulation. A free-floating buoy or an object lying on the seabed are also subject to buoyancy and hydrodynamic drag, and a robot that interacts with such objects will have to cope with the corresponding force. We now detail the implementation of the given model.

Implementation

World plugin for hydrodynamics The water velocity v c and the buoyancy direction are taken into account in a world plugin, The Gazebo simulator subscribes to a ROS topic that gives the intensity and direction of the water current. It is thus possible to have the current vary if needed. The buoyancy direction is of course opposed to the gravity, and becomes null as the vehicle goes above the surface. A Gazebo world file is given and handles this plugin.

We use additional tags in the URDF 4 to handle buoyancy and damping related to a given link or object. This allows giving the position of the center of buoyancy with regards to the link frame, together with the compensation coefficient that indicates the percentage of gravity that is compensated. In most I-AUV this coefficient is slightly greater than 1, meaning that without external effort the vehicle tends to come to the surface. Here is an example of these additional tags: <buoyancy> <compensation>1.01</compensation> <origin xyz= ".5 0 .2"/> <limit radius=".5"/> <damping xyz="40 80 80" rpy="20 50 50"/> </buoyancy>

The <limit> tag indicates where the buoyancy force begins to decrease when approaching the water surface. In I-AUV simulation we are mostly interested in controlling the vehicle when it is entirely underwater, therefore no complex computation are carried to handle the surface approach: the buoyancy force applied to a given link simply decreases to zero when the link is near the surface.

Model plugin for thruster control Thrusters are also defined in the URDF. For each thruster the maximum effort is given, together with the mapping between the thrust direction and the body frame: <link>base_link</link> <bodyCommandTopic>body_command</bodyCommandTopic> <jointCommandTopic>joint_command</jointCommandTopic> <!--for each thruster, give the map to XYZ+RPY and the maximum effort --> <thruster> <map>-1 0 0 0 0 ${-body_width/2}</map> <effort>30</effort> </thruster> It is thus possible to simulate thruster saturation and to handle it in the control law. For now all thrusters are assumed to be fixed on the same link, which is defined as base link in the given example. This tag is read by the model plugin to allow low-level control. In the given example, we see that the Gazebo simulator will listen to thruster forces on the body command topic, and to joint efforts on the joint command topic. These efforts can be generated from a high-level force controller, or more usually from the low-level controller that is included in the simulator. We now detail the implementation of this controller. 3 Low-level PID control

In this section we describe the PID controller of the AUV simulator. Because of the numerous uncertainties on the dynamic model, it is very tedious to design a AUV advanced force controller. The usual approach is to have a high-level controller that generates position or velocity setpoints for the body and the arm, and a low-level PID controller that ensures the given setpoints are followed.

General design A independent ROS node handles the global PID controller, which is implemented as a set of parallel or cascaded controllers depending on the user choice. The two configurations are detailed in Fig. 2 (cascaded) and Fig. 3 (parallel). p * and v * correspond to the desired position and velocity of a joint value or a body linear or angular direction. p and v represents the estimation of the position and of the velocity, coming from the sensors or from an estimation algorithm. The velocity PID is the same in both configuration, but the position PID outputs the velocity setpoint in cascaded mode, and the effort in parallel mode. Dynamic Reconfigure5 allows real-time gain tuning, and a configuration file is used to store the tuned gains. A ROS service is used to switch between position and velocity control, which are defined independently for the body and the arm.

Anti-windup and body thrusters mapping The maximum velocities and efforts that are given in the URDF are used for anti-windup inner loops in each PID. In particular, the maximum thruster efforts are mapped to the body directions in order to get the actual controllable directions and their associated maximum effort. No PID are instanciated for non-controllable directions (typically roll, but sometimes also pitch and sway). The position setpoint of the AUV is given in the world frame, while the velocity setpoint are given in the vehicle frame. The several PID outputs represent the desired wrench τ * u to apply to the AUV body. This desired wrench is then mapped to the thruster efforts τ * t in order to apply the saturation, and mapped back to the actual wrench τ u to be applied in Gazebo:

τ u = T.τ t = T.sat(τ * t ) = T.sat(T + τ * u ) ( 4
)
where T is the mapping matrix between the body directions and the thrusters and T + is its Moore-Penrose pseudo-inverse. Another possible strategy in case of thruster saturation is to scale all thruster efforts instead of only applying the saturation:

τ u = T. nt min i=1 τ max i τ * i τ * t ( 5
)
where n t is the number of thrusters and τ max i is the maximum effort of thruster i. The real-time load of the thrusters is available for logging or high-level control, as we will see in Section 5. We now detail the integration of the simulator with Gazebo, UWSim and higher-level control nodes.

Integration

In this section we detail the interactions between the I-AUV simulation in Gazebo, the PID controller, the UWSim visualization and other potential higherlevel controllers. First the file integration is mentioned as it is necessary to synchronize the Gazebo and UWSim description files. We then expose and detail an example of graph of nodes.

File synchronization

Gazebo and UWSim both have their own approach to describe the robot model, additional moving objects and the world terrain. As UWSim is only a kinematic simulator, the URDF used to describe the robots are usually simpler that the Gazebo ones, that contains all the inertial and collision-related data. On the opposite, the scene file, used to describe the whole world setup in UWSim, contains all information about the considered simulation while the same data is separated in several files when using Gazebo. In this section we briefly expose the synchronization script that helps building lesser-information files from higher ones.

The starting point is the UWSim launchfile that contains the UWSim scene file, that describes the whole setup. This allows retrieving the list of robots URDF files and moving objects, their starting position and the terrain 3D mesh that is used. The URDF files may already exist or need to be generated from existing xacro files, that are used in Gazebo. All starting positions are found in the UWSim scene file, which allows generating a Gazebo model spawner: this ensures that the same objects will be present at the same position in Gazebo and UWSim. Finally, for now UWSim and Gazebo do not use the same convention for file pathes. The pathes are thus adapted on the fly when generating or updating missing files. With such synchronization, an experiment can first be carried only in UWSim with kinematic control for simplicity and debugging, then with the dynamic simulator and low-level PID loops to have a realistic behavior. We now detail the interactions between the several ROS nodes in a typical simulation.

Node interactions

This section details the interactions between the simulation node, integrated through plugins in Gazebo, and the other nodes of a typical I-AUV simulation. The graph of nodes is represented in Fig. 4. We can see all moving objects have their own namespace. Here /g500arm5e is the robot namespace and regroups all nodes and topics related to its low-level control, while /blackbox is related to a moving black box. Moving objects are considered as non-actuated robots, thus have their own namespaces. The gazebo node subscribes to the /gazebo/current topic, which allows defining and modifying the water current.

The main controller is called /main control. It reads body and joint positions from the corresponding topics. In this simulation these values are the true ones and come directly from Gazebo. Water current, as well as body and joint position setpoints are generated. These setpoints are handled by the pid control node, which generates the joint effort joint command and the body wrench body command. These topics are read by Gazebo in order to apply the corresponding forces and torques to the vehicle and its arm.

As we can see on the graph, UWSim only receives data in order to visualize the current pose of the I-AUV and of the black box. Still, an embedded camera is simulated in UWSim and the generated images could be used in the main controller to perform visual servoing as shown in [START_REF] Kermorgant | A global control scheme for free-floating vehicle-manipulators[END_REF].

We now expose the simulation results.

Experiments

In this section we show an example of behavior of a black box recovering mission.

(a) End of approach. The hook is lying on the seabed and the arm cannot follow its setpoint.

(b) Beginning of lifting. One of the vertical thruster is saturating.

(c) Recovery complete. Experimental setup The considered I-AUV is the Girona 500 [START_REF] Ribas | The girona 500, a multipurpose autonomous underwater vehicle[END_REF], equipped with a 3-degrees of freedom robot arm mounted with a hook on the end-effector.

The AUV body has 5 thrusters, allowing to control all degrees of freedom except for the roll. A keel effect is induced by the relative position of the center of gravity and the center of buoyancy such that the roll is stable. The AUV is slightly lighter than water, thus the two vertical thrusters are always used to stabilize the depth.

Simulation scenario As the goal is to show the dynamic behavior, position control is in open loop and the PID gains are not well tuned. The scenario includes several steps that illustrate difficulties of underwater interventions:

1. At first only joint control is performed to set the arm in an idle position. 2. At iteration 70, body position control is used to approach the black box. At the end of the approach, around iteration 350, the hook will touch the seabed and thus the arm will not be able anymore to follow its setpoint. 3. At iteration 520, the joint setpoint is changed to prepare the recovery 4. At iteration 820, the body position setpoint is changed to go up and recover the box, which is entirely carried by the AUV after iteration 840. Because of the weight of the box, some joints and some thrusters will be saturated at this time.

Fig. 5 depicts the I-AUV in various steps. At the end of the approach (Fig. 5a), the hook slightly hits the seabed, which will appear on the joint measurements. Fig. 5b shows the beginning of the recovery, which corresponds to some oscillations and to the saturation of a vertical thruster. In Fig. 5c the I Body position control The overall behavior of the body control is shown in Fig. 6. The vertical dotted lines correspond to the key-iterations that are mentioned above. The measurement rate is 30 Hz. The low-level PID controllers ensure the body position and orientation setpoints are followed after iteration 70 (beginning of body position control). Interaction with the black box is clearly visible from iteration 840. As soon as the box is entirely recovered, small oscillations appear on the orientation error (Fig. 6b). This is due to the fact that the box behaves like a pendulum when carried by the recovering hook, and induces some perturbations on the vehicle position control.

Thruster load The described body behaviors also appear in Fig. 7, that represents the thruster load in percent. When gived a new position setpoint at iteration 70, we see one of the vertical thruster (cyan) is saturated. This explains the slow reaching of the heave setpoint (red) in Fig. 6a. From iteration 350 the hook is lying on the seabed. The load is thus reduced on the vertical thruster. It is saturated again after iteration 820, when the new position setpoint makes the I-AUV go up. The oscillations of the orientation error in Fig. 6b are also visible in Fig. 7 after iteration 820. The thruster load oscillations are more important than the orientation error, which means the PID controllers of the thrusters manage to greatly reduce the perturbation induced by the carrying of the black box.

Joint behavior

The behavior of the 3-degrees of freedom arm is represented in Fig. 8. First the setpoints are followed but when the hook comes in contact with the seabed at iteration 350, the arm is stuck between the ground and the I-AUV body and it becomes impossible to follow the position setpoint. The setpoint 8. Joint positions and setpoints. Slew (blue), shoulder (green) and elbow (red). The hook touching the seabed is clearly visible at iteration 350, as the arm cannot ensure position setpoint. Oscillations appear after iteration 820 because of effort limits.

q 2 q 3 q * 1 q * 2 q * 3 Fig.
changes at iteration 520, making it possible to be followed again. Once the arm begins to lift the black box at iteration 820, the joint effort limits clearly appear as the shoulder and elbow cannot cope with the weight of the black box.

Discussion

The whole experiment demonstrates the dynamic simulation, together with contact physics between the hook and both the seabed and the black box. Dynamic coupling is clearly visible when lifting the box, as its pendulum-like oscillation are passed on the whole I-AUV. The PID controllers of the thrusters manage stabilizing the body position, but the arm still oscillates and cannot follow its position setpoint due to the limits on joint efforts. The overall simulation shows that dynamic simulation of a I-AUV raises realistic difficulties that are possible to try and solve in the simulated environment, before having to perform real experiments.

Conclusion

We have presented a dynamic simulator for intervention autonomous underwater vehicles. Through the ROS framework, it was possible to use both the dynamic simulation capabilities of Gazebo and the realistic underwater rendering of UWSim. The considered modeling is focused on the main effects due to hydrodynamic forces, that are drag and buoyancy. A possible improvement is to take into account the added inertia and Coriolis, but these effects are difficult to quantify precisely. That is why all uncertainties are here considered in the damping coefficients d(v). Additional URDF tags allows defining the buoyancy and damping parameters of each robot link or other moving object, such as the black box that is considered in the example. The presented experiment was designed to show that dynamic simulation can be used to reproduce realistic difficulties such as thruster or joint effort limit, collision and uncertainties on setpoint following. As in real I-AUV, this can be due to non-optimal tuning of the low-level controllers, or on the uncertainties on the hydrodynamic parameters. The proposed simulation framework can thus be used to validate new algorithms for sensor-based control or state estimation. In future works we will focus on underwater sensors in order to perform realistic sensor-based missions.

Fig. 1 .

 1 Fig. 1. Rendering comparison between Gazebo (a) and UWSim (b). osgOcean allows UWSim to render automatically water color, visibility, floating particles and other underwater characteristics.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Cascaded PID, here in position control. The position PID outputs a velocity setpoint that is then followed with the velocity PID.+

Fig. 4 .

 4 Fig. 4. Auto-generated graph of nodes (round-shaped) and topics (rectangular).

Fig. 5 .

 5 Fig. 5. UWSim view of the simulation.

  Body orientation error in degree. Roll (magenta), pitch (cyan) and yaw (yellow).

Fig. 6 .

 6 Fig.6. 3D behavior of the I-AUV body. The heave (red) follows its setpoint after iteration 70, while the orientation error (b) is null even if small oscillations appear when lifting the black box after iteration 840.

2 Fig. 7 .

 27 Fig.7. Thruster load in percent. Surge/yaw (blue and green), sway (red) and heave/pitch (cyan and magenta).

  -AUV has recovered the black box and is going up. Because of the weight of the black box, some joints are saturated.

	5
						t x
	0					t y t z
						t * x
	-5					t * y
						t * z
	-10
	0	200	400	600 iterations	800	1000
	(a) Body position (plain) and setpoint
	(dotted). Surge (blue), sway (green) and
	heave (red).

http://www.openscenegraph.org

http://code.google.com/p/osgocean

https://github.com/freefloating-gazebo

Universal Robot Description File

http://wiki.ros.org/dynamic reconfigure