**Application of Photo-Identification and Lengthened Deployment Periods to Baited Remote Underwater Video Stations (BRUVS) Abundance Estimates of Coral Reef Sharks**

## **Mauvis Gore 1,2, Rupert Ormond 1,2,3,\*, Chris Clarke 4, Johanna Kohler 1,2,5, Catriona Millar 1,2 and Edward Brooks <sup>6</sup>**


Received: 18 June 2020; Accepted: 23 October 2020; Published: 3 November 2020

**Abstract:** Baited Remote Underwater Video Stations (BRUVS) are widely used for monitoring relative abundances of fishes, especially sharks, but only the maximum number of individuals seen at any one time (MaxN) is usually recorded. In both the Cayman Islands and the Amirante Islands, Seychelles, we used photo-ID to recognise individual sharks recorded on BRUVS videos. This revealed that for most species the actual numbers of separate individuals (IndN) visiting the BRUVS were significantly higher than MaxN, with, for example, ratios of IndN to MaxN being 1.17 and 1.24 for Caribbean reef, *Carcharhinus perezi*, and nurse, *Ginglymostoma cirratum*, sharks in the Cayman Islands, and 2.46 and 1.37 for blacktip reef, *C. melanopterus*, and grey reef, *C. amblyrhynchos*, sharks, respectively, in the Amirantes. Further, for most species, increasing the BRUVS deployment period beyond the 60 min normally used increased the observed IndN, with more than twice as many individuals in the Cayman Islands and >1.4 times as many individuals in the Amirantes being recorded after 120 min as after 60 min. For most species, MaxN and IndN rose exponentially with time, so data from different deployment periods cannot reliably be compared using catch-per-unit-effort (CPUE) calculated as catch-per-unit-time. In both study areas, the time of first arrival of individuals varied with species from <1 min to >2 h. Individually identifiable sharks were re-sighted after up to 429 days over 10 km away in the Cayman Islands and 814 days over 23 km away in the Amirantes, demonstrating that many individuals range over considerable distances. Analysis of Cayman re-sightings data yielded mean population estimates of 76 ± 23 (SE) and 199 ± 42 (SE) for *C. perezi* and *G. cirratum*, respectively. The results demonstrate that, for sharks, the application of both photo-identification and longer deployment periods to BRUVS can improve the precision of abundance estimates and provide knowledge of population size and ranging behaviour.

**Keywords:** capture–recapture; Cayman Islands; Seychelles; monitoring; endangered species; maximum number of individuals; photo-identification; mark-recapture; movement ecology

## **1. Introduction**

Globally, many species of shark have suffered a dramatic decline in abundance over the past 30 years [1–3], principally as a result of by-catch and finning in response to the demand for shark fin by the Far-East restaurant trade [4,5]. As a consequence, 11 shark species are now considered Critically Endangered, 15 Endangered and 48 Vulnerable [6]. Most large-bodied sharks are threatened or endangered [2,7], while even once common reef species have become scarce in many areas with, for example, the long-established shark fishery in the Seychelles recording a significant decrease in catch per effort from 16 to 6 kg day−<sup>1</sup> between 1992 and 2016 [8].

As a result, monitoring the population abundance of sharks is much needed [9]. Baited Remote Underwater Video Stations (BRUVS), underwater video cameras set to record the fishes attracted to bait, is now the method most frequently used for surveying and monitoring sharks in particular [10–12]. Stereo-BRUVS, in which two cameras are focused towards the bait, can provide further data by enabling both the distance to a fish and its length to be determined [13], although the associated equipment costs are necessarily higher if similar numbers of units are to be deployed. Cappo et al. [14] have reviewed BRUVS techniques for camera, bait, and orientation of bait and camera, as well as potential auditory, olfactory and behavioural cues, while Harvey et al. [15] concluded that BRUVS provide a statistically robust and cost-effective method of monitoring diverse assemblages of fishes in a number of habitat types.

BRUVS are not the only visually-based method regularly employed for monitoring fishes. Several forms of Underwater Visual Census (UVC) are widely used by divers [16–18]. However, remote video has the advantage of being practical in conditions where divers cannot easily operate, as in very deep water or rough conditions, and of being able to detect species, such as large sharks, which are often too scarce or too wary to be regularly sighted by divers [13]. Video methods are only fully effective in reasonably clear water, and baited video will only detect species that either are resident close to the BRUVS, or passing through, such as roaming schools of carangids, or like sharks have a well-developed olfactory sense enabling them to home-in on the bait [19].

To date, the most commonly used measure of the numbers of fishes coming to the BRUVS has been the maximum number of fishes recorded in a single frame during the viewing period (referred to as MaxN or Nmax). While simple for an observer to record, this measure suffers a number of limitations, most obviously that the observed MaxN may be much lower than the total number of separate individual fish attracted to the bait. Harvey et al. [15] discussed alternative approaches, including the use of cumulative MaxN (recorded each time MaxN is updated for each species). Campbell et al. [20] compared MaxN with MeanCount (the mean number of fishes observed in a series of frames or stills). They concluded that these two measures were comparable for abundance estimation, although MeanCount resulted in less precision for all species analysed. Similarly, Schobernd et al. [21] compared the advantages of using MaxN and MeanCount in both simulations and laboratory experiments, as well as in empirical data. They concluded that an advantage of MeanCount is that it is linearly related to actual abundance, whereas MaxN is not.

Understanding the accuracy and precision of monitoring and surveying methods is important for interpreting any apparent change or lack of change in species number or relative abundance. A weakness of the above counting methods is that they do not record whether different individuals are visiting a BRUVS at different times, as opposed to a single individual making repeat visits. Nor do they detect when the same individual visits two or more different BRUVS. This limitation is most obvious with less abundant but not uncommon species, such as many reef sharks and rays, in which different individuals may visit the same BRUVS, but at different times.

To address this issue, studies on other taxa have experimented with photo-identification as a tool for distinguishing different individuals. Stobart et al. [22] used BRUVS to record spiny lobster, *Palinurus elephas*, and compared three measures of abundance: MaxN, mean MaxN (using 5 min sampling periods) and the true number of individuals per recording based on individual identification (referred to in the present study as IndN). The authors concluded that while all three measures could distinguish areas of contrasting population density, MaxN and MeanN were not appropriate for

documenting changes in population abundance because they were liable to sample saturation. That is, an increase in population abundance may not result in any increase in the MaxN. More recently, Irigoyen et al. [23] used three relative abundance indices in a survey of broadnose sevengill shark, *Notorynchus cepedianus*: Nmax, NmaxIND (cumulative number of different sharks and similar to IndN) and Nocc (total number of occurrences in a BRUVS session). The index NmaxIND appeared to show a greater abundance of sharks than Nmax. Sherman et al. [24] used photo-identification to determine the true numbers of two species of ray visiting BRUVS and found that the actual abundances were for *Neotrygon orientalis* 2.4, and for *Taeniura lymma* 1.1, times greater than the recorded MaxN values. However, to date this approach has not been generally applied to BRUVS-based monitoring of sharks, in part because the earlier cameras used were of lower resolution, so that it was not feasible to distinguish between most individuals of the same species.

Photo-identification has been found to be practical in some other studies of sharks. Population estimates of basking shark, *Cetorhinus maximus* [25], and great white shark, *Carcharhinus carcharius* [26], have been made based on individual recognition through quality images of dorsal fins. Photo-identification based on images captured by divers or fishers has also been used in ecological studies of blacktip reef shark, *C. melanopterus* [27], Greenland shark, *Somniosus microcephalus* [28], nurse shark, *Ginglymostoma cirratum* [29], sicklefin sharks, *Negaprion acutidens* [30], broadnose sevengill shark, *Notorynchus cepedianus* [23], and whale shark [31,32], as well as flapper skate, *Dipturus intermedius* [33]. In addition, BRUVS have been used to monitor the presence and size [34], as well as site fidelity [11], of individually identified sharks. Here, the feasibility of applying photo-identification methodology to the video footage recorded by BRUVS is assessed to monitor more effectively the abundance of coral reef shark species and test whether the data can also be used to generate additional information on population size or ranging behaviour.

The deployment period used to record individuals visiting a BRUVS (i.e., the effective length of deployment) has also varied between studies. The majority of previous studies have used deployment periods of 60 min (for example [24,35,36]), often because this was the maximum recording time that could be guaranteed by the batteries in older cameras. Some studies used shorter times, such as 30 min ([19] fish species), while Brooks et al. ([10] shark species) used 90 min, Harasti et al. ([34] white shark) 5 h and Stobart et al. ([22] lobster) 7 h. After using BRUVS to record coral reef fishes in the Hawaiian Islands, Asher [37] concluded that while a short deployment of 20 min was sufficient to capture fast-reacting species and residents, longer periods of 60 min were required to capture macropiscivores, including sharks. Our experience suggested that additional individual and species of shark could arrive at a BRUVS even after this time. A study by Torres et al. [38] exceeded other deployment times in using a duration of 24 h for sharks in the western Mediterranean. As a result of the wide range of deployment times used in different studies, the effect of deployment period on the numbers of individuals and species of shark recorded at a site was also tested in the present study.

In the present study, the BRUVS video recordings obtained during two shark conservation projects, one in the Cayman Islands (western Caribbean) and the other in the Amirante Islands (south-western Seychelles) were reviewed. Photo-identification methodology was applied to distinguish the separate individuals visiting BRUVS over different periods. The aims of the study were to assess whether (a) the use of photo-identification can provide better estimates of the numbers of sharks being recorded, (b) the approach can be used to generate estimates of population size by applying mark–recapture analyses to the data, and (c) the method can also be used to provide information about the movements of individuals. In addition, (d) it was investigated whether for reef sharks the numbers of species and individuals recorded varied significantly with BRUVS deployment period. We found that photo-identification could enhance the value of data obtained in relation to all three questions. However, the present study is not intended to compare abundances between areas or habitats, or across years (hence locations of individual stations are not included); those analyses are being reported on separately.

## **2. Materials and Methods**

### *2.1. Study Areas*

The study was based on data collected from two widely separated island groups, one in the western Atlantic and the other in the Indian Ocean. The Cayman Islands are located in the North-West Caribbean Sea (Figure 1a) and include three islands, Grand Cayman, Little Cayman and Cayman Brac. These islands are emergent sections of the Cayman Ridge which runs adjacent to the 7 km deep Cayman Trench. Sampling was conducted around the largest island, Grand Cayman (in a study area bounded by Latitudes 19.25◦ N and 19.42◦ N and Longitudes 81.05◦ W and 81.44◦ W), and also around Little Cayman (in a study area bounded by Latitudes 19.64 and 19.73◦ N and Longitudes 79.95 and 80.13◦ W), 130 km to the northeast. Both islands have well-developed fringing coral reefs beyond which a narrow coastal shelf drops steeply to very deep water. A series of marine parks and conservation/replenishment areas, occupying about 25% of the 160 km long coastline, were designated over 30 years ago. Shark abundance, although a little higher than on many other Caribbean Islands, was still considered low in relation to relatively pristine reef areas [39]. In consequence sharks were given full protection throughout Cayman waters in 2015 [40].

**Figure 1.** Locations of study areas: (**a**) the three islands constituting the Cayman Islands in the Caribbean Sea, and (**b**) D'Arros Island and St. Joseph's Atoll in the Amirante Islands, Seychelles, in the Indian Ocean. Data analysed for this study resulted from deployments of BRUVS at varied intervals of 0.5–2 km around the whole of both Grand Cayman and Little Cayman in the Cayman Islands, and D'Arros Island and St. Joseph's Atoll in the Seychelles. The faint red lines around each island indicate the extent of territorial waters [41].

*Oceans* **2020**, *1*

The Amirante Islands in the south-western part of the Seychelles, in the western Indian Ocean (Figure 1b), are a series of widely separated small islands and coral reefs sitting on the 180 km by 40 km Amirante Bank, beyond which the plateau gives way to very deep water. The islands include D'Arros Island (53.29◦ E, 5.41◦ S) and the adjacent St. Joseph Atoll (53.33◦ E, 5.43◦ S), around both of which the data used here were collected (the combined study area being bounded by Latitudes 5.27◦ S and 5.47◦ S and Longitudes 53.17◦ E and 53.38◦ E). D'Arros and the outer faces of St. Joseph Atoll are surrounded by well-developed fringing reefs that slope down from a reef crest to the seabed at 15 to 20 m, while the atoll itself encloses a 3.5 km wide sandy lagoon, 1–6 m deep. Adjacent bank areas (5–35 m deep) that were also surveyed have a mainly sandy seabed, which in its shallower parts are colonised by dense seagrass beds (mostly *Thalassodendron ciliatum*) and occasional coral patches. D'Arros and St. Joseph Atoll are more remote from significant human populations than the Cayman Islands, and the adjacent reefs have been protected from fishing for almost three decades; in contrast, however, more distant bank areas are regularly fished, including by targeted long-lining for sharks. Nevertheless, sharks are more abundant in the Amirante Islands than in the Cayman Islands, and a comparison between the two areas provided an indication of how different levels of abundance influence the value of applying photo-identification to BRUVS recordings. The precise locations of sampling points are not relevant to the analysis, and for both ethical and legal reasons are not included in this report.

## *2.2. BRUVS*

In both study areas the BRUVS body was built from an open weave plastic crate, upturned and weighted down by a series of dive weights attached close to the rim (now the lower side). The bait arm was a 1.5 m long PVC tube fixed to the top of the crate, with the bait enclosed in a double-layered plastic mesh bag placed at the far end of the bait arm (Figure 2). The amount of bait used was standardised at 300 g of mackerel, *Scomber scombrus,* in the Cayman Islands and bonito, *Euthynnus a*ffi*nis*, in the Amirante Islands. A camera was fixed on the crate next to the base of the bait arm and pointed towards the bait bag. The cameras used were models 3 and 4 Silver GoPros, which recorded in high definition (HD) and could be fitted with double battery packs, allowing recording periods of 2+ h. The availability of HD cameras was key to the present work, as the high quality of images extracted from video permitted the detection of fine detail on the individual sharks.

**Figure 2.** Image of a BRUVS as used in the present study, showing the open weave plastic body, a high-definition camera (GoPro Hero 3 or 4) under a roll cage, a bait arm with outer bait bag (lower left) and lines (upper right) leading to a surface marker buoy. (Photo with permission by James Lea).

BRUVS were normally deployed in both study areas on the fore-reef at depths of 10–18 m, but at 2–4 m deep in the lagoon of St. Joseph Atoll and in the Cayman Islands' Sounds (equivalent to lagoons). Deployment by boat took place between 07:35 and 15:34 h in the Amirante Islands and 08:26 and 15:04 in the Cayman Islands. BRUVS were positioned on the seabed in open spaces among corals on reefs, or in gaps within seagrass beds. Appropriate placement of the BRUVS was achieved with the aid of either a viewing scope or of a snorkeler who signalled when the BRUVS was over a suitable area. In both locations BRUVS were deployed in sets of four, one at each of four stations that were located at variable intervals of 500 m to 2 km in a line along each section of coast (i.e., coasts orientated in a particular direction); the arrangement thus approximated a stratified random sampling design. In repeat deployments over successive years BRUVS were placed at the same stations at the same sites. Survey work was undertaken when the Douglas sea scale was ≤ 4 and wave or swell height ≤ 1.5 m. Following retrieval of the BRUVS, memory cards were removed from the cameras and the files copied to computers on which the videos could be viewed. Visual analysis was carried out for the full length of each video, with a minimum of 25% of each video being checked by a second researcher to ensure that comparable standards were being observed. For each shark that was detected, the time(s) of arrival in view and loss from view were recorded on each occasion, along with species, sex, maturity and size, as estimated by comparison with the size of the bait bag and graduated marks on the bait arm [42].

In the Cayman Islands BRUVS were deployed once or twice a year, in or around June and December, and for present purposes data from 2015 to 2017 have been analysed. In the Amirante Islands BRUVS were deployed once per year, in or around February or March, and for this study data from 2015 to 2017 were similarly selected.

## *2.3. Photo-Identification*

To facilitate individual recognition, so far as possible detailed information was recorded on the appearance of each newly observed shark. The characteristics recorded included the shapes of fins and other body parts (if these appeared in anyway atypical), and any variation in natural patterning or colouration, as well as the shape and location of any naturally acquired markings, such as wounds, injuries or scars (Figures 3 and 4). Details were recorded of distinguishing features on both sides of the body, so that with few exceptions individuals could be recognised on viewing from either side. Acquired marks could be temporary, as with minor abrasions or injuries, but most used in individual identification were more permanent, as when part of a fin was missing or torn. Evidence from longer-term studies suggested that while many non-fatal injuries heal, scars and serious fin injuries in sharks can remain evident for many years [25,26].

**Figure 3.** *Cont.*

**Figure 3.** Examples of variations in natural patterning and injuries on dorsal fins of sharks recorded as screen grabs from BRUVS videos. (**a**) Differences in natural patterning on right and left sides of first dorsal fins of two different blacktip reef sharks, *Carcharinus melanopterus*, (first shark i and ii, second shark iii and iv), and (**b**) dorsal fins of four different grey reef sharks, *C. amblyrhynchos.* The yellow circles highlight the location of distinctive features; the yellow lines are added to allow variation in the shape of the black fin tip to be seen more clearly. Note that the different sides of an individual blacktip reef shark's dorsal fin tend to be similar, but are not identical to each other.

**Figure 4.** Examples of variation in natural markings (circled) on different shark species taken from BRUVS videos of two different Caribbean reef sharks, *Carcharinus perezi* (top), a nurse shark, *Ginglymostoma cirratum*, (centre), and a set of four images (bottom) showing distinctive features on different individual whitetip reef sharks, *Triaenodon obesus*.

Features observed on individual sharks were marked on a standardised cartoon (sketch-like drawing) for each shark (see [25] for further discussion of the method and for examples). Characteristics recorded on each cartoon were then entered into a computer database, along with the details of each deployment and selected images showing the individual's distinctive features. The cartoons and database could be interrogated for individual characteristics to facilitate matching with any other recorded individual. Potential matches were then either confirmed or rejected by comparison of the images or video clips of each shark. Potential matches were independently verified by a second researcher.

## *2.4. IndN and MaxN*

Data on the numbers of sharks were used to compare the results of using two different methods for recording the numbers of sharks of each species attracted to each BRUVS. First, MaxN, the greatest number of sharks of a given species observed in a single frame during a set period, was determined. Second, photo-identification was used to determine so far as possible IndN, the total number of different individuals of each species visible at any time in the field of view. In practice, not all the sharks detected on screen approached the BRUVS close enough for the individual to be identified. Individuals that lacked distinguishing features enabling them to be recognised in the long term, but that could still be separated from other sharks seen on the same recording (for example, if the shark was a different size), were added to the tally. Nevertheless, on some occasions two or more similar individuals may not have been seen clearly enough to distinguish, and on others not every shark attracted may have entered the field of view. Thus, in some instances, the IndN recorded will have been less than the actual number of individuals visiting.

Catch per unit effort (CPUE) has been used to compare the numbers of sharks visiting different BRUVS [10], but this metric can be misleading because, as evident from our results, the numbers of new sharks arriving at a BRUVS tends to decline with time. That is, the plot of the accumulative number of individuals against time is asymptotic. Therefore, longer deployments may generate lower CPUEs, even if more sharks are recorded. For this reason, the values of MaxN and IndN were compared with respect to three different deployment periods commencing from the time when the BRUVS arrived at the seabed; these were 60, 120 and 150 min. This comparison was undertaken to investigate the effects of using longer deployment periods by determining if these result in the detection of more individuals and/or more species. Deployment period was the term used to describe the recording time on the seabed.

## *2.5. Statistical Analysis and Mark–Recapture Population Estimates*

A series of tests were applied to the data to determine statistical significance as detailed in the relevant parts of the Results section; these include Friedman ANOVA, Wilcoxon matched pairs test, ANOVA, paired t and post-hoc LSD tests. The statistical package used in the study for parametric and non-parametric tests (described in Results) was Statistica 64 V. 12.

Although large numbers of BRUVS videos were analysed (>500) and a large number of individually identifiable sharks recorded, the numbers of individuals re-sighted on successive campaigns were too low for open population mark–recapture models to be employed, such as the Schnabel and Jolly–Seber methods which we have applied to basking shark in Scotland [25]. However, the numbers of re-sightings in the Cayman Islands for six surveys over three years were high enough for successive population estimates of two species to be made using the Chapman derivative of the basic Lincoln–Peterson estimator [43,44], which is more robust for small sample sizes [45]. For each "re-capture" occasion *t*, the population size *Nt* was calculated where *nt* = the number of sharks sighted on sample occasion *t*, *rt* = the number of marked sharks re-sighted in sample *nt*, and *mt* = the number of marked sharks at the beginning of sample occasion *t*:

$$N\_l = \left[ (m\_l + 1) \* (n\_l + 1) / (r\_l + 1) \right] - 1 \tag{1}$$

Assuming that *k1* number of sharks were initially marked (i.e., identified) at sample occasion 1 and that *k2* = *n2* − *r2* was the number of unmarked (previously unknown) sharks sighted by the end of sample occasion 2, then *mt*, the number of effectively marked sharks at the next sample occasion *t*, equals *kt*. The total population size (*N*) was estimated as the mean across sample occasions *s* [41,42] according to:

$$
\overline{N} = \sum \text{Nt/(s-1)}\tag{2}
$$

#### **3. Results**

A total of 661 BRUVS were deployed, 451 in the Cayman Islands (at 92 locations) and 210 BRUVS in the Amirante Islands (at 60 locations). These BRUVS recorded a total of six species in the Cayman Islands and eight in the Amirante Islands (Table 1), representing two orders (Carcharhiniformes and Orectoloboformes) and three families (Carcharhinidae, Ginglymostomatidae and Sphyrnidae) of sharks.

**Table 1.** The numbers of individually identifiable sharks of different species (i.e., those carrying distinguishing features enabling them to be recognised again as individuals) and the corresponding percentages these represent of all records of that species observed on video recordings from Baited Remote Underwater Video Stations (BRUVS) (a) in the Cayman Islands (451 BRUVS) and (b) in the Amirante Islands (210 BRUVS).


## *3.1. Recognisability of Individual Sharks*

In most species the majority of sharks observed could be recognised as individuals (Table 1), although the proportions of sharks observed for which sufficient features were evident to permit this varied from 100% for tiger shark, *Galeocerdo cuvier,* to 9.1% for bull shark, *Carcharhinus leucas*. In some species, individuals were relatively easy to distinguish. For example, the dorsal fins of *C. melanopterus* (85.3% individually identifiable) are capped by an area of black pigmentation, the exact shape of which varies markedly both between individuals and between opposite sides of the same individual (Figure 3). Similarly, whitetip reef sharks, *Triaenodon obesus*, (63.9% individually identifiable) have spot patterns on the body that are unique to individuals (Figure 4). Other species, such as *C. amblyrhynchos*, (77.9% individually identifiable) and *G. cirratum* (82.8% individually identifiable) (Figure 4), required more careful examination of the video footage to discern individually distinctive markings. *C. leucas* and silvertip, *C. albimarginatus*, sharks were the most difficult to identify as individuals on recordings since they tended not to approach the camera, making detection of individual features often impracticable. The mean percentage of distinguishable sharks across species was 60% and the overall percentage of recognisable individuals (irrespective of species) was 66%.

## *3.2. Comparison of IndN with MaxN*

For all shark species combined in both the Cayman Islands and the Amirante Islands, the mean IndN (1.88, 2.71, 2.67) was significantly greater than mean MaxN (1.20, 1.57, 1.51) for each of the three deployment periods (60, 120 and 150 min, respectively) (Table A1). Considering the species separately, in the Cayman Islands for both *C. perezi* and *G. cirratum*, mean IndN was significantly greater than mean MaxN for each deployment period (60, 120 and 150 min) (Figure 5, Table A1), the overall ratio of IndN to MaxN being 1.17 for *C. perezi*, and 1.24 for *G. cirratum*. Similarly, in the Amirante Islands, mean IndN was significantly greater than mean MaxN for all three deployment periods for *C. amblyrhynchos*, *C. melanopterus* and *N. acutidens,* but only for 120 min for tawny nurse shark, *Nebrius ferrugineus*, (Figure 6, Table A1). The overall ratios of IndN to MaxN were 2.46 for *C. melanopterus*, and 1.37 for *C. amblyrhynchos*. For other species, the sample sizes were too low for significant differences to be detected (Figures 5 and 6). For all species in which there were significant differences between IndN and MaxN, the ratio of IndN to MaxN increased with time, as it did also for *C. limbatus* and *T. obesus*, in which differences were not significant (Table A2).

**Figure 5.** Cayman Islands: The mean number (bar) and 95% confidence interval (error bar) of individually identified sharks observed for deployment periods of 60, 120 and 150 min for IndN compared with the same periods for MaxN. Significant differences found between pairs of IndN and MaxN using paired *t* tests are noted by for *p* ≤ 0.05, - for *p* ≤ 0.001 and --for *p* ≤ 0.0001 (see Table A1).

**Figure 6.** Amirante Islands: The mean number (bar) and 95% confidence interval (error bar) of individually identified sharks observed for deployment periods of 60, 120 and 150 min for IndN compared with the same periods for MaxN. Significant differences found between pairs of IndN and MaxN using paired *t* tests are noted by for *p* ≤ 0.05, - for *p* ≤ 0.001 and -- for *p* ≤ 0.0001 (see Table A1).

The number of BRUVS deployments in which the application of individual identification revealed the presence of more sharks than would otherwise have been evident (i.e., where IndN > MaxN) is shown in Table 2. While the use of photo-identification revealed the presence of additional sharks in the more abundant species (three in the Cayman Island and six in the Amirante Islands), this was not the case for the less common species.

**Table 2.** Proportion of BRUVS deployments in which IndN was greater than MaxN in (a) the Cayman Islands (*n* = 451) and (b) the Amirante Islands (*n* = 210). The proportion as a percentage of BRUVS is the percentage of BRUVS with sharks in which IndN was greater than MaxN.


## *3.3. E*ff*ect of Deployment Period for All Species Combined*

The mean IndN and MaxN of all species combined after different recording times in each of the two study areas are significantly different in each case (Table 3). In both the Cayman Islands and the Amirante Islands there was a significant increase in IndN and in MaxN both from 60 to 120 min, and from 120 to 150 min in the Amirante Islands (Table A3).

**Table 3.** Comparison of IndN and MaxN for all shark species combined ± standard deviation (S.D.) recorded on BRUVS in (a) Cayman Islands (BRUVS *n* = 109) and (b) Amirante Islands (BRUVS *n* = 55) for 60, 120 and 150 min deployment periods. Only cases where BRUVS ran for 150 min were used to compare across the three periods with Friedman ANOVA χ2, *p* = probability value, df = 2.


Overall, using IndN there were 2.1 times as many sharks recorded after 120 min as after 60 min in the Cayman Islands, and 1.4 times as many recorded in the Amirante Islands after 120 min than after 60 min. Similarly, there were 1.6 and 1.4 as many individual sharks recorded after 150 min as after 120 min in the two areas, respectively.

## *3.4. E*ff*ect of Deployment Period for Separate Shark Species*

The mean IndN (i.e., the mean number of separate individuals visiting a BRUVS) was significantly different over the periods 60, 120 and 150 min for *C. limbatus*, *C. perezi*, *G. cirratum* and lemon shark, *Negaprion brevirostris*, in the Cayman Islands, and for *C. amblyrhynchos*, *C. melanopterus*, *N. ferrugineus*, *N. acutidens* and *T. obesus* in the Amirante Islands (Table 4). Among these species, *C. perezi*, *G. cirratum*, *C. amblyrhynchos, C. melanopterus* and *N. acutidens* showed a significant increase in mean number from 60 to 120 min, while *C. perezi* and *G. cirratum* also showed an increase from 120 to 150 min (Table A4).

**Table 4.** Comparison of mean ± standard deviation (S.D.) values of IndN and MaxN for each shark species for BRUVS deployment periods of 60, 120 and 150 min for (a) the Cayman Islands and (b) the Amirante Islands. Shown too are the results of a Friedman ANOVA (χ2, *p* = probability value, *n* = number of BRUVS in sample; df = 2) for the differences in the means across the three deployment periods. Post hoc results are available in Table A4.


The mean MaxN was significantly different across different deployment periods for *C. perezi*, *G. cirratum* and *N. brevirostris* in the Cayman Islands and for *C. amblyrhynchos*, *C. melanopterus*, *N. ferrugineus*, *N. acutidens* and *T. obesus* in the Amirante Islands (Table 4). Post-hoc tests indicated that in the Amirante Islands, *C. amblyrhynchos*, *C. melanopterus* and *N. acutidens* all showed a significant increase in mean MaxN from 60 to 120 min, as well as *G. cirratum* in the Cayman Islands, while *C. perezi* showed significant increases in mean MaxN both from 60 to 120 and from 120 to 150 min (Table A4).

## *3.5. Arrival Times at BRUVS*

As a further means of assessing effective BRUVS deployment periods, the mean times of first arrival (at those BRUVS on which the species were recorded) were determined for each species in the Cayman Islands and the Amirante Islands (Figure 7). For two species in the Cayman Islands, the mean first arrival time was close to 60 min, while the remaining four species took longer (mean of 71–76 min). In the Amirante Islands, the mean first arrival time was less than 60 min for four species and was about an hour for two other species, but in two species, *C. albimarginatus* and great hammerhead (*Sphyrna mokarran*) shark, the mean arrival time was greater than an hour. It was notable that some individual *C. amblyrhynchos* and *C. melanopterus* (in the Amirante Islands), and *G. cirratum* and *C. perezi* (in the Cayman Islands), took more than 2 h to arrive.

**Figure 7.** Mean and 95% confidence interval of time taken (hh:mm) by shark species to arrive at the BRUVS camera. Top: Cayman Islands, Bottom: Amirante Islands. The red line indicating a mean time taken to arrive of 1 h (60 min) is included for convenience.

In the Amirante Islands, the variation in first arrival time across species was statistically significant (ANOVA: F(7,533) = 5.05, *p* = 0.00001), with *C. leucas* when present arriving significantly earlier than *C. albimarginatus*, *C. amblyrhynchos*, *C. melanopterus*, *N. ferrugineus, N. acutidens* and *T. obesus* (Post hoc LSD tests: *p* < 0.006, respectively). *C. melanopterus's* first arrival time was significantly shorter than those for *C. amblyrhynchos*, *N. ferrugineus* and *N. acutidens* (Post hoc LSD tests: *p* <0.04, respectively). There was no significant variation in first arrival time across species in the Cayman Islands (ANOVA: F(5,411) = 0.19, *p* = 0.9).

## *3.6. Re-Sighting Rates*

The application of photo-identification methodology allowed many individuals to be recognised again if they were re-sighted. The proportions of the more commonly recorded species that were re-sighted during the study periods on a different BRUVS deployment ranged from 0% for *C. limbatus* to 10.1% for *G. cirratum* in the Cayman Islands, and from 4.5% for *N. ferrugineus* to 9.8% for *C. melanopterus* in the Amirante Islands (Table 5). In the Cayman Islands, a third of individual *G. cuvier* and *S. mokarran* were re-sighted, but only three individuals of either species were recorded in total. The mean of the re-sighting rates for each species across both locations was 12%; the combined re-sighting rate (from pooling all individuals irrespective of species) was 10.9%.

**Table 5.** The re-sightings rates on BRUVS for shark species (a) in the Cayman Islands and (b) in the Amirante Islands, both within years and between years. Only re-sights recorded on BRUVS videos are included here (as opposed to other re-sightings that were recorded during other related activities). Note that the total number of individuals re-sighted over the whole period is not necessarily the sum of the re-sights in each year, since some re-sights may be of individuals that have already been re-sighted.


## *3.7. Movements of Individual Sharks*

Re-sightings of individuals on a BRUVS different from those on which they were originally recorded provided information about shark movement patterns. In the Cayman Islands, individuals were re-sighted at distances of up to 10 km for *C. perezi* and 15 km for *G. cirratum* from the location of the original sighting (Table A5). In the Amirante Islands, individuals were re-sighted at distances of up to about 25 km for *C. melanopterus* shark, 33 km for *C. amblyrhynchos* shark, and about 27 km for *N. acutidens*, *N. ferrugineus* and *T. obesus* (Table A5), although both these last two species were also often re-sighted on BRUVS close to the stations at which they were originally recorded. While a majority of re-sights occurred within the same season, in the Amirante Islands some re-sights of *C. amblyrhynchos*, *C. melanopterus* and *N. acutidens* occurred after longer periods of time, the longest being of an individual *C. melanopterus* re-sighted after 814 days. In the Cayman Islands, some individuals were also detected after long periods of time, the longest being 429 days for *C. perezi* and 184 days for *G. cirratum*.

## *3.8. Population Estimates*

For most species, the numbers of re-sightings were too low to permit the reliable estimation of population size, but using the Chapman derivative of the basic Lincoln–Peterson estimator successive estimates were possible for the numbers of *C. perezi* and *G. cirratum* present in the Cayman Islands study area. The mean estimated number of individually recognisable *C. perezi* over 2015–2017 was 76 ± 23 (S.E.), and of individually recognisable *G. cirratum* over the same period was 199 ± 42 (S.E.). Comparable estimates for *C. melanopterus* and *C. amblyrhynchos* in the Amirante Islands study area are in hand.

## **4. Discussion**

Comparison of the present data with those obtained by other researchers in the same ocean regions is possible (Table 6 [10–12,42,46–48]), but confounded by differences in method, such as the time and depth ranges at which the BRUVS were deployed. In particular, while most researchers report MaxN for a deployment period of 60 min, relevant studies have used deployment periods ranging from 54 to 90 min, with one using 6 h [48]. Some authors (Brooks et al. [10]) have corrected for differences in deployment period by reporting numbers as catch per unit effort (CPUE) of one hour, on the assumption that MaxN increases linearly with time. However, in the present study, the values of both MaxN and IndN rise exponentially with time, hence longer deployment periods will appear to yield lower values of CPUE if presented this way. Nevertheless, some comment is possible. In the Cayman Islands, the abundance of two species of sharks was greater and of four species less than that reported in Eleuthera, the Bahamas. The abundances of *N. acutidens* appear markedly greater in the Amirante Islands than in the eastern Indian Ocean. Abundances of *C. melanopterus* and *N. ferrugineus* appear broadly comparable. More detailed analysis of the spatial and temporal patterns of abundance for the different species recorded in our study areas will be reported separately.

It is also important to note that when fish visible on BRUVS videos are counted (e.g., to determine MaxN), normally the results provide information only on relative abundance (that is, differences in abundance between locations or between occasions), but not on overall population size or density. Even though such population indices are widely used for wildlife monitoring and are generally assumed to be directly proportional to population density, in many field studies this assumption has been found to be invalid [49]. It is a weakness of BRUVS methodology as generally employed that the most commonly used measure, MaxN, does not increase linearly with actual abundance, but tends to saturate [22]. While this would be less of an issue for schooling or non-territorial species, it becomes more of a concern in species which are territorial or keep their distance from each other.

**Table 6.** A comparison of the values of MaxN reported by authors for BRUVS deployed (A) in the Caribbean Sea and tropical western mid-Atlantic Ocean, and (B) inthe Indian Ocean. In (A), the area in the Bahamas was Eleuthera (1: Brooks et al. [10]) and in Belize it was Glover's Reef (6: Bond et al. [11]). (B) included an area inthe BMR in the BIOT, Indian Ocean (3: Tickler et al. [12]), MOU74 and Rowley Shoals, E. Australia (4: Meekan et al. [42]), Raja Ampat in Indonesia (5: Jaiteh et al.[46]; 6: Beer [47]), and Houtman Abrolhos Islands, W. Australia (7: Santana et al. [48]). \* denotes unfished or protected areas. Catch per unit effort (CPUE) values forthe present study were calculated as the shark MaxN per hour on BRUVS.


In the present study, we sought to address this issue by identifying sharks as individuals and determining the actual number of sharks visiting each BRUVS. As this study confirmed, most individual sharks carry distinctive marks or features allowing them to be identified as individuals. Some marks are natural while others are acquired, with older individuals tending to have more characteristics through injuries and scars than neonates and juveniles. In 12 of the 14 species recorded, at least two-thirds of the observed sharks could be recognised as individuals through such features. In most cases where it was not possible to discern any distinguishing features, this was because the shark's appearance on camera was too brief or too distant. Some species, notably *C. albimarginatus* and *C. leucas*, in which the percentages of distinguishable individuals were the lowest (9.1% and 16.7%, respectively), never closely approached the bait even though they may have been attracted to the vicinity of the BRUVS by olfactory cues. With other species the proportion of individually identified sharks was lower than it might otherwise have been because of high turbidity and poor visibility at some stations. This was especially the case in the sandy lagoon in St. Joseph Atoll, where only a minority of individuals could be identified. In particular, this affected the percentage of individually identifiable *C. melanopterus* (85.3%), which otherwise are readily identified as individuals because of the very variable shape of the black patches on their dorsal fins (see also [50]). In summary, while many sharks can be recognised as individuals with certainty, a proportion of them viewed on BRUVS cannot be distinguished. Thus, in many cases the difference between the true number of visiting sharks and MaxN will be even greater than between MaxN and IndN.

The species for which for most BRUVS IndN was greater thanMaxN are those that are largely coastal by habitat. In contrast, more open-water species *C. albimarginatus*, *C. limbatus, G. cuvier* and *S. mokarran*, appeared rarely and only as single individuals. Mean IndN was significantly greater than MaxN for each deployment period overall (60, 120, 150 min). Further, in most species, including *C. limbatus, C. perezi* and *G. cirratum* in the Cayman Islands and *C. amblyrhynchos*, *C. melanopterus*, *N. ferrugineus* and *T. obesus* in the Amirante Islands, the mean ratio of IndN/MaxN increased with recording period. These findings indicate that in at least these two study areas, for most species an increased deployment period resulted in visits by significantly more individuals, not simply in return visits by the same individuals. Thus, IndN may at times be a more useful index than MaxN, especially for assessing ongoing population trends in the least abundant species, such as *S. mokarran*.

A comparison of the data obtained after different BRUVS deployment periods also showed that while most shark species present might be recorded during the course of survey work using 60 min deployments, nevertheless both the number of species recorded and the numbers of individuals of most species increased with longer BRUVS deployment periods. Interestingly, in a very recent study by Torres et al. [38], a deployment time of 24 h showed that while no additional species were observed after 60 min deployment, additional individuals continued to arrive for a long time. In the present study, the 120 and 150 min deployment periods showed a greater number of both individuals and species, with some individuals of four species taking over 2 h to arrive. An example of a non-shark elasmobranch arriving late during a recording session was a bowmouth guitarfish, *Rhina ancylostoma*, which arrived at a BRUVS only 2 h 17 min 7 s after deployment. Thus we concluded that extended recording periods may be advantageous, provided the fieldwork schedule permits this.

Besides deployment period, the number of sharks attracted to a BRUVS, and hence the precision of monitoring, is also influenced by other factors. Our experience over 10 years suggests that the quality, as well as the quantity, of bait can be important. The oil in bait appears to be a component particularly attractive to sharks, but the oil content can vary greatly not only between bait species, but with the season in which it is caught (relative to fish maturity and reproductive state) for fat content, where they are fished, and with catch storage conditions. For example, the fat content of fished tuna has been reported to vary from 0.4 to 13.5%, and of sardine from 1.2 to 18% [51–54]. The quantity of bait may in our experience be less critical. We used less bait than many researchers, but in a pilot study we found that using over five times the amount of bait used here (1550 g as opposed to 300 g) resulted in only a small, non-significant increase in the numbers of sharks arriving at the BRUVS.

The method of containment of bait on the BRUVS is also important. When, during earlier work, bait was enclosed only in a large mesh-size bag, most or all of it could be quickly removed by sharks, but also by scavenging fishes such as species of snapper (*Lutjanus*), emperor (*Lethrinus*), durgon (*Melichthys*) and octopus. In such circumstances, all of the bait could be consumed in well under an hour, in which case extended deployment would not result in further sharks being recorded. In both studies described here, double bags of fine mesh inside larger mesh were employed, so that the scent of the bait could still disperse readily, while the bait itself was difficult for fishes to extract. Other workers have used perforated cans as bait containers [23,32]. As in any monitoring program, the data obtained from BRUVS-based shark monitoring programs will be more reliable if such considerations are taken into account and all aspects of the method standardised as far as possible.

As well as providing more precise estimates of abundance, the application of photo-identification to BRUVS can provide information on the behaviour of individuals. Re-sightings showed that some individuals of most species moved distances of up to 33 km across the study area. *C. melanopterus*, *G. cirratum* and *N. ferrugineus*, however, were also often re-sighted close to the stations where they were originally recorded. Data from acoustic tagging work undertaken in both study areas and elsewhere have provided evidence that some of the species studied can move over even greater distances than evident from the BRUVS data. For example, one *C. perezi* was recorded moving 125 km from Grand Cayman to Little Cayman and the same distance back again to Grand Cayman [39]. Similarly, in the present study individual *C. amblyrhynchos* in the Amirante Islands were observed on different BRUVS up to 20 km apart. Both longand short-range movements of this species have also been documented in New Caledonian waters by Bonnin et al. [55], who suggested that the long-range movement of adults was motivated by mating opportunity. In contrast, only a few of the >50 tagged *C. melanopterus* in the Amirante Islands have been recorded crossing the <1 km wide channel between D'Arros and St. Joseph Atoll [56], even though in Polynesia this species has been recorded crossing open water between islands up to 50 km apart [57].

If the low rates of re-sighting that were observed in the present study are typical of such photo-identification work, then the movement data generated will be more limited than can be obtained using acoustic tagging, which frequently records the locations of individuals if they are present within a study area, provided that a dense network of receivers is installed over the area of interest. Information derived from photo-identification work may nevertheless be of value where the considerable funding required to mount an effective acoustic tagging study is unavailable. Furthermore, where BRUVS are deployed on a regular basis within a specific area, photo-identification data can reveal whether individual sharks make frequent use of an area or are resident within it.

It was also anticipated that re-sightings data might allow the application of mark–recapture modelling to estimate the population of different species within the study areas. In practice, the numbers of individuals of even the more frequently recorded species re-sighted on a different BRUVS during the same sampling campaign were low (ranging from 2 to 9%), and the numbers re-sighted in the following year even lower (0.4–2.7%). Encouragingly, however, low rates of re-sighting suggest that known individuals must be mixing with a larger population of other individuals. Thus the model used indicated for the Cayman Islands a mean effective population sizes of 76 ± 23 (S.E.) *C. perezi* and 199 ± 42 (S.E.) *G. cirratum*. These estimates are of interest in light of shark conservation efforts in Cayman, where all shark species have been protected since April 2015 [39,40]. It will be important to determine whether estimated population sizes of these species increase with time. Local opinion surveys indicated that the shark protection legislation had overwhelming public support [39], being opposed by only a small minority of boat-based fishers who may be tempted to ignore it.

There are also two factors which confound the interpretation of the population estimates. Firstly, the calculated population sizes are only of individually identifiable sharks. Following Gore et al. [25], these numbers need to be corrected upwards to allow for the proportions of sharks (of the species concerned) that lacked distinguishing marks. However, a simple proportional correction based on the proportions of individually identifiable sharks may result in an over-estimate, because some of the sightings of non-identifiable sharks will also have been re-sightings. Furthermore, while some sharks

could not be identified as individuals because they lacked distinguishing features, others could not be identified as individuals since they did not approach the BRUVS closely enough or because the water was too turbid, and so could in fact have been known individuals.

A second uncertainty relates to the size of area to which these population estimates apply. The three Cayman Islands have a combined terrestrial area of only 264 km2 and an estimated length of coast of 160 km, most of which is fringed by coral reef. These reefs drop steeply into water >500 m deep within 2 km of the shore. However, *C. perezi*, for example, have been recorded at depths of 150 m or more [58–60], and in the Cayman Islands were recorded moving between Grand Cayman and Little Cayman, islands 125 km apart separated by water > 1000 m deep [39]. Further, it seems from the results of conventional and acoustic tagging studies of sharks on reefs that, while some individuals may be resident or semi-resident in an area, many are only occasionally present in the area or may never be recorded again [39,58,61,62]. Thus, some of the individuals observed in the present study may be foraging and mixing over a much larger region than the 160–320 km2 shelf area immediately surrounding the individual Cayman Islands. Comparable analysis is underway of the Amirante Islands data, but is confounded by the facts that the study area covers only a third of the archipelago and that the area did not achieve marine protected area status until 2019.

## **5. Conclusions**

The results of the present study demonstrate that the photo-identification of individual sharks can be applied to BRUVS videos and that in most species this approach can result in a greater ability to distinguish between different abundances. This was most evident for species in which there were usually too few sharks for more than one individual to appear on screen at a time. Extended recording periods enhanced this benefit, increasing both the number of shark species and the numbers of individual sharks detected per BRUVS, at least up to 150 min. The application of photo-identification to BRUVS can also provide information about the movement ecology of a species; this is helpful since it indicates for management the appropriate size of the protected area or the foraging range of specific individuals of interest or concern. Importantly, the method clarifies whether within a study area the sharks recorded on BRUVS are frequently the same individuals or usually different ones. Additionally, the method can generate indicative estimates of population size, which can be of value in provoking either concern or reassurance (depending on the numbers) in relation to local shark conservation programs. However, the process of matching images is very time consuming (in the present study the videos from 661 BRUVS were analysed, an undertaking which required about 1500 h of desk time), and the reliable detection of distinguishing features on BRUVS videos requires patience and experience. Sherman et al. [24] similarly noted that identifying individuals does provide better accuracy and additional information, but the determination of actual numbers takes time. MaxN is a simpler and quicker method, and is particularly appropriate for large-scale surveys or routine monitoring. On the other hand, the more sensitive estimates attained by determining IndN are valuable if attempting to monitor population trends of a scarce species or in a smaller area, such as an island or an MPA. IndN is sensitive to differences in abundance that MaxN would not detect. The balance of costs and benefits in applying this approach will depend on the aims and circumstances of any particular project.

## **6. Research Ethics**

The work on endangered shark species was carried out in the Cayman Islands with permission from and on behalf of the Cayman Islands Department of Environment, and, in the Amirantes, under a letter of approval from the Republic of Seychelles Ministry of Environment, Energy and Climate Change.

**Author Contributions:** Conceptualisation, M.G. and R.O.; Data curation, M.G., J.K. and C.M.; Formal analysis, M.G., R.O., J.K. and C.M.; Funding acquisition, M.G., R.O. and C.C.; Investigation, M.G., R.O., C.C., J.K., C.M. and E.B.; Methodology, M.G., R.O. and E.B.; Project administration, M.G.; Resources, C.C.; Supervision, M.G. and R.O.; Validation, R.O.; Writing—original draft, M.G., R.O. and J.K.; Writing—review and editing, M.G., R.O., C.C., J.K., C.M. and E.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** The work in the Cayman Islands was supported by a Darwin Initiative grant (DPLUS036) from the UK Department of Environment, Food & Rural Affairs, by the Department of Environment, Cayman Islands, and by the CayBrew Whitetip Shark Fund. Work in the Seychelles was supported by the Marine Research Facility, North Obhur, Jeddah, Saudi Arabia. RO was part-supported through a research chair at King Abdulaziz University, Jeddah, funded by HRH Prince Khalid bin Sultan.

**Acknowledgments:** We thank the Cayman Islands Department of Environment and all the staff on D'Arros Island, Seychelles, for their support and assistance. We also thank Peter Davies, Martin Eaton, Frida Lara and Robert Nestler for assistance with fieldwork and Tim Austin and James Lea for advice.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **Appendix A**

**Table A1.** Comparison between mean IndN and mean MaxN for shark species in (a) the Cayman Islands and (b) the Amirante Islands for 60, 120 and 150 min deployment periods, with differences tested using post hoc Wilcoxon matched pairs test (statistic Z), *n* = number of cases and *p* = probability values for Figures 5 and 6.


**Table A2.** The ratio of IndN to MaxN for shark species for deployment periods of increasing length (60, 120 and 150 min) in the Cayman Islands and Amirante Islands. *n* = number of BRUVS included in analyses. IndN = number of individuals identified. MaxN = maximum number of individuals observed on screen at any one time. *n* decreases with increasing time periods, because not all recordings ran for 120 min, and a minority for 150 min.


**Table A3.** Statistical comparison of pairs of deployment periods of IndN and MaxN for all shark species combined in (a) Cayman Islands and (b) Amirante Islands, using post hoc Wilcoxon matched pairs tests (statistic Z), *n* = number of cases and *p* = probability values for Table 3.



**Table A4.** Statistical comparison for pairs of BRUVS deployment periods (60, 120, 150 min) of IndN and of MaxN for each shark species separately in (a) Cayman Islands and (b) Amirante Islands, using post hoc Wilcoxon matched pairs tests (statistic Z), *n* = number of cases and *p* = probability values for Table 4.

**Table A5.** The range of distances and the respective time periods taken for each shark species.


Data availability: 10.6084/m9.figshare.13061546.

## **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Use of Polyphosphates and Soluble Pyrophosphatase Activity in the Seaweed** *Ulva pseudorotundata*

## **Juan J. Vergara \*, Patricia Herrera-Pérez, Fernando G. Brun and José Lucas Pérez-Lloréns**

Departamento de Biología, Facultad de Ciencias del Mar y Ambientales, Universidad de Cádiz, E-11510 Puerto Real, Cádiz, Spain; patricia.herrera@uca.es (P.H.-P.); fernando.brun@uca.es (F.G.B.); joselucas.perez@uca.es (J.L.P.-L.)

**\*** Correspondence: juanjose.vergara@uca.es

Received: 1 October 2020; Accepted: 14 December 2020; Published: 16 December 2020

**Abstract:** The hydrolytic activity of different types of polyphosphates, and the induction of soluble pyrophosphatase (sPPase; EC 3.6.1.1) activity have been assessed in cell extracts of nutrient limited green seaweed *Ulva pseudorotundata* Cormaci, Furnari & Alongi subjected to different phosphorus regimes. Following a long period of nutrient limitation, the addition of different types of (poly)phosphates to artificial seawater enhanced growth rates on fresh weight and area, but not on dry weight bases. Chlorophyll and internal P content were affected by P supply. In contrast, internal soluble reactive P was kept low and was little affected by P additions. Soluble protein content increased in all treatments, as ammonium was added to prevent N limitation. The C:N:P atomic ratio revealed great changes depending on the nutrient regime along the experiment. Cell extracts of *U. pseudorotundata* were capable of hydrolyzing polyphosphates of different chain lengths (pyro, tripoly, trimeta, and polyphosphates) at high rates. The sPPase activity was kept very low in P limited plants. Following N and different kind of P additions, sPPase activity was kept low in the control, but slightly stimulated after 3 days when expressed on a protein basis. The highest activities were found at the end of the experiment under pyro and polyphosphate additions (7 days). The importance of alternative P sources to phosphate and the potential role of internal soluble pyrophosphatases in macroalgae are discussed.

**Keywords:** phosphorus; polyphosphates; pyrophosphate; pyrophosphatase activity; seaweed; *Ulva*

## **1. Introduction**

Although nitrogen (N) is considered to be the main limiting nutrient in the ocean [1], phosphorus (P) can also limit primary production of algae in some coastal areas [2]. The major P form for algae is inorganic ortho-phosphate (Pi), whose uptake in seaweeds occurs through specific transport systems [3–5]. Besides the main Pi source (i.e., from the water column and sediments), the hydrolysis of organic P monoesters by the periplasmic enzyme alkaline phosphatase also renders Pi [6]. However, although Pi is the main form of P used by algae, their growth can also be supported by more complex inorganic P compounds [7]. Accordingly, pyrophosphates (PPi) and polyphosphates (polyP) could be important sources of P [8]. However, these compounds are not molybdenum sensitive, not being quantified in standard analytical measurements of soluble reactive phosphorus (SRP) [9,10]. These Pi polymers, which can contribute to the P pool in the coastal ocean, are mostly derived from anthropogenic activities [9]. However, since P is taken up as Pi, enzyme cleavage of Pi polymers at the cell surface is needed. Data on the role of these alternative sources of P in seaweeds are scarce.

In contrast to P uptake kinetics and transmembrane transport mechanisms, little is known about P internal cellular metabolism. Marine algae, as opposed to most terrestrial plants, are able to store Pi polymers as polyP [2,8]. Studies of 31P-NMR on the chlorophyte *Ulva lactuca*, identified polyP as the main intracellular storage for P [11]. This polyP pool was accumulated under external Pi availability, while it was mobilized under P limitation. The PolyP polymers have also been reported in other *Ulva* spp., in the Rhodophyta *Ceramium* sp. [11] and in *Porphyra purpurea* [12], among other species.

Soluble pyrophosphatase (sPPase) (EC 3.6.1.1) is a ubiquitous enzyme essential for cell anabolism [13]. It is located in cellular organelles (plastids and mitochondria) of a variety of photosynthetic autotrophs [14], where many anabolic and biosynthetic reactions occur [15]. This enzyme is responsible for PPi cleavage and Pi regeneration in vivo, thus allowing biosynthetic reactions to proceed [16]. Its function is also essential to replenish Pi for phosphorylation [17]. Despite the importance of this enzyme, the few data available are mostly from photosynthetic prokaryotes and some microalgal groups [14], whereas studies in seaweeds are much scarcer, e.g., [7].

Under this framework, the aim of this work was double: (1) to assess the capacity of cell extracts of the green macroalgae *U. pseudorotundata* to hydrolyze Pi polymers, and (2) to study the effect of different P sources (i.e., Pi, PPi and PolyP) on the response of nutrient-limited algae (growth, internal C:N:P composition, chlorophyll, protein content and internal sPPase activity). Results indicate the capacity of *U. pseudorotundata* to exploit alternative sources of P for growth, as well as the induction of internal sPPase activity in cell extracts of nutrient-limited specimens following P additions.

## **2. Materials and Methods**

## *2.1. Plant Material and Preculture Conditions*

Seaweed specimens, identified as *Ulva pseudorotundata* Cormaci, Furnari & Alongi, according to [18], were collected from earthern pools at an intensive fish farm (Acuinova S. L., San Fernando, Cádiz, southern Spain). Once in the laboratory, thalli were rinsed with seawater and cleaned of mud and epiphytes. Some of these field-collected specimens were used both, to optimize sPPase assays and to determine the hydrolytic activities of different kinds of polyphosphates in cell extracts (see below). Pretreatment: the remaining biomass was grown in aquaria with constant aeration and seawater pumped from a natural underground well (pH 7.3, salinity 39.6, 35 μM NH4 <sup>+</sup>, 0.8 μM PO4 3−, and negligible NO3 − levels, mean values) for 70 days before starting the experiments (see below). Lighting (145 μmol photons m−<sup>2</sup> s<sup>−</sup>1, which is saturating for this species [19]) was provided with two fluorescent lamps (Phillips cool white) in a 12 h:12 h light-dark photoperiod. Mean pH and temperature were 8.4 and 25.1 ◦C, respectively, along this pre-treatment phase. Culture was kept at a high biomass density (16 g FW L<sup>−</sup>1) in a total volume of 20 L, and seawater was renewed weekly. The aquarium was cleaned at every seawater change with distilled deionised water plus sodium hypochlorite. Germanium dioxide (2 mg L−<sup>1</sup> , final concentration) was added to inhibit diatom growth [20]. The objective of this long preculture period at a high biomass density was to obtain nutrient-depleted thalli (i.e., with very low N and P quota), despite there not being a complete absence of PO4 <sup>3</sup><sup>−</sup> and NH4 <sup>+</sup> in seawater. In fact, the growth rate was very low along this period (data not shown) and N and P quotas approached the subsistence ones for this species (see initial conditions in results), with a N:P atomic ratio of about 15, close to the Redfield ratio, indicating an incipient N and P co-limitation for *Ulva* species [21].

## *2.2. Experimental Design*

The N and P-limited thalli were grown under four different P treatments: (1) without P (control), (2) with orthophosphate (Pi, 5 μM final concentration), (3) with pyrophosphate (PPi, 5 μM final concentration) and (4) with polyphosphate (PolyP, 5 μM final concentration). Inorganic N (as NH4 +, 5 μM final concentration) was also added to all treatments (including the control) to avoid a parallel N co-limitation. All chemical reagents for P analytics were from Sigma, and polyphosphate used is referred to as "phosphate glass", with a mean length chain of 13–18 Pi molecules. Thallus discs (2 cm diameter) were punched with a cork borer and grown in 1 L flasks (see below) within an incubation chamber (Koxka model EC−540-F) at 20 ◦C and at a photon flux of 200 <sup>μ</sup>mol m−<sup>2</sup> s−<sup>1</sup> in a 12 h:12 h light-dark photoperiod provided by fluorescent tubes (cool white, Phillips). Each of the four treatments (*n* = 3 per treatment) were run for seven days in 1 L-aerated flasks (12 in total) at a biomass density about 1 g FW L−<sup>1</sup> (c.a. 20 discs per flask) with synthetic seawater [22], which allowed us to control nutrient availability. During the experiment, seawater and *Ulva* samples were taken at days 0, 2, 4 and 7. Seawater was renewed at days 2 and 4.

## *2.3. Analyses*

Three discs were randomly sampled in each of the 12 flasks, measured and weighted (fresh and dry, after drying in an oven at 60 ◦C for 48 h) to estimate growth rates assuming an exponential growth model. Fresh weight/area (g FW m−2) and DW/FW ratios were also estimated. Total chlorophyll (*a* + *b*) concentration was determined according to [23], after extracting liposoluble pigments in N, N, dimethyl formamide for 24 h at 4 ◦C in darkness. Internal C and N concentrations were determined from dried biomass at the beginning and at the end of the experiment in a Perkin-Elmer 240-C elemental analyser. The time course for the C:N ratio was not monitored since there was not enough algae material to be collected. The SRP was determined according to [24]. Total internal P was also determined according to [24], after acid digestion in an autoclave (121 ◦C, one hour, Raypa Sterilmatic AE-150). In both cases, absorbance was measured at 885 nm (Hitachi, U-1100 spectrophotometer, Tokyo, Japan).

Cell extracts were prepared by grinding *Ulva* discs with liquid N2 with a mortar and a pestle. Cell extract was resuspended in 50 mM Tris buffer pH 7.5 at an approximate ratio of 100 mg FW in 1.5 mL buffer. Tris buffer was supplemented with 1 mM DTT and 1 mM PMSF to protect active enzyme groups from oxidation and to inhibit proteases. Cell extracts were centrifuged at 15,000× *g* for 15 min at 4 ◦C. Cell debris and membrane fractions were discarded and soluble components were recovered, frozen in liquid N2, and stored at −80 ◦C for enzyme assays. An aliquot of cell extracts was taken for soluble protein determination, to scale enzyme activities. Soluble protein content of cell extracts was determined according to [25].

The sPPase activity was measured in cell extracts according to [26], monitoring the Pi production at 25 ◦C for 10 min. Assay medium contained Tris buffer 50 mM pH 8.0, 20 mM of substrate (either PPi, tripolyP, trimetaP or polyP depending on the treatment), 20 mM Mg2<sup>+</sup>, and an appropriate quantity of cell extract. The ion Mg2<sup>+</sup> is an essential cofactor for sPPase activity for most of the assayed organisms [14]. In preliminary tests, a 200 μL extract was assayed in a total volume of 1 mL, and 500 μL of the assay medium was assayed for Pi production. Given the high activity of the enzyme, 50 μL cell extract was used instead of 200 μL, and 100 μL of the assay medium was diluted with 400 μL Tris buffer pH 8.0 before measuring Pi, to reduce the sensibility of the assay, yielding a dilution factor of 20× in comparison to preliminary assays. Immediately, an aliquot was taken and mixed with ammonium molybdate and ascorbic acid to measure Pi production.

A unit of sPPase activity was defined as the quantity of enzyme activity that yields 1 μmol Pi per minute in the assay conditions stated above. This activity was scaled either to proteins or DW. In all enzyme assays two blanks were run in parallel: (1) in absence of the substrate (PPi), which indicates the amount of internal SRP in the sample, and (2) in the absence of cell extract, to determine spontaneous PPi hydrolysis and/or SRP contamination of the substrate. These assays were done with 200 μL extract, and no further dilutions were performed for the measurements. Internal SRP values were quite low compared to the rate of Pi production by the enzyme activity, and the spontaneous hydrolysis was almost undetectable (data not shown).

#### *2.4. Statistics*

The effect of P source and time on thallus area, fresh weight:area, dry weight:fresh weight ratio, internal SRP, internal total phosphorus, soluble protein and sPPase activity were analyzed by a 2-way ANOVA (Table A1). Prior to any statistical analysis, data were checked for normality (Shapiro–Wilk normality test) and homocedasticity (Bartlett test of homogeneity of variances test). The effects of different P sources on enzymatic activity and different P sources on C:N:P tissue composition were analyzed by means of 1 way ANOVA (significant Tukey's test differences indicated in Tables 1 and 2). Post hoc comparisons among means were tested by the Tukey test [27]. In all cases, the null hypothesis was rejected at the 5% significance level and data were presented as mean ± SE.

**Table 1.** *Ulva pseudorotundata*. Enzyme soluble hydrolytic activities for different polyphosphate substrates in cell extracts of field collected samples. Data are expressed as mean (*n* = 3–4) ± SE. Different letters indicate significant differences among treatments at α = 0.05 following Tukey's test and 1 way ANOVA.


**Table 2.** *Ulva pseudorotundata*. Internal C:N:P composition at the beginning of the experiment and after 7 days under different phosphorus treatments. Data are expressed as mean (*n* = 3) ± SE. Different letters indicate significant differences among treatments at α = 0.05 following Tukey's test and 1 way ANOVA.


### **3. Results**

The hydrolytic capacity of several polyphosphate compounds assayed with soluble cell extracts from *Ulva pseudorotundata*, as a first approach, is shown in Table 1. Hydrolytic activities were recorded for all the substrates tested. Both sPPase and tripolyphosphatase activities showed similar values, whereas the cyclic trimetaphosphate was hydrolysed at a lower rate. In contrast, polyphosphatase activity was 3–4 times higher than that of the sPPase.

*Ulva pseudorotundata* was maintained for a long period under nutrient limitation (preculture) and, subsequently, grown under different P sources, but with NH4 <sup>+</sup> addition to avoid N-limitation. Seaweed growth rates were affected by P treatments, but their values depended largely on the plant attributes used in the calculations (Figure 1a). Thus, normalization by fresh weight (FW) resulted in higher values, with those obtained under different P additions being higher than the controls. In contrast, growth rates were minimal when normalized by dry weight (DW), being unaffected by P treatments. On area basis, the growth rates were also lower than those estimated on a FW basis, but were affected by P additions. The time-course of *U. pseudorotundata* disc area grown under different P treatments is shown in Figure 1b. A positive thallus expansion was observed regardless of P treatment, with the highest values under PPi and polyP enrichment, and the lowest in the control. As a consequence of the different water gain and area expansion, the FW/area ratio increased along time in all treatments (Figure 2a). In contrast, the DW/FW ratio showed the opposite trend, decreasing along time in all the treatments (Figure 2b).

Total chlorophyll content was also affected by P sources (Figure 3), with the lowest values recorded in the controls. Chlorophylls *a* and *b* displayed the same trend, as the chlorophyll *a*/*b* ratio was constant along the experiment (data not shown).

**Figure 1.** *Ulva pseudorotundata*. (**A**) Specific growth rates on fresh weight (FW), dry weight (DW) and area basis of *seaweed* discs cultured under different P regimes (control without P, Pi, PPi, and PolyP) for 7 days. Data are means ± SE (*n* = 3 independent cultures). (**B**) Disc area of *seaweed* discs along the experiment as a function of the different P regimes. Data are means ± SE (*n* = 6–9 discs).

**Figure 2.** *Ulva pseudorotundata*. (**A**) Time course of fresh weight/area ratio and (**B**) dry weight/fresh weight ratio, as a function of different P regimes. Data are means ± SE (*n* = 6–9 discs).

**Figure 3.** *Ulva pseudorotundata*. Time course of total (*a* + *b*) chlorophyll content as a function of different P regimes. Data are means ± SE (*n* = 3).

Whereas total internal P content was unaffected in the controls along time, it doubled the initial values under different P sources, especially for PPi and PolyP (Figure 4a). In contrast, internal SRP content was kept at quite low levels, (1–4 μmol Pi g−<sup>1</sup> DW). Although a transient increase after 2 days was recorded for the PPi treatment, all values were fairly constant along the experiment in all treatments (Figure 4b).

**Figure 4.** *Ulva pseudorotundata*. (**A**) Time course of total internal P and (**B**) internal soluble reactive P (SRP), as a function of different P regimes. Data are means ± SE (*n* = 3).

Seawater SRP was also monitored along the experiment (Figure 5). At the beginning, it was undetectable in both the control and the PPi treatment, but values of 7.5 μM and 0.5 μM were recorded for Pi treatment and PolyP treatments, respectively, indicating, in the latter case, the possibility of a certain degree of Pi contamination in the added substrate (added PolyP was about 5 μM). Significant Pi levels were detected in PolyP-enriched cultures, probably as a consequence of the partial hydrolysis of long chain Pi molecules. In contrast, Pi was not detected in either Pi or PPi treatments, except for a high value in the PPi treatment on the fourth day. To check for the spontaneous hydrolysis of PPi and PolyP in seawater, 1 L aerated flasks were maintained without algae for 3 days, and SRP from PPi and polyP hydrolysis was measured. Values ranged from 0.06 to 0.16 μM Pi (almost undetectable) for PPi and from 0.50 μM to 0.66 μM for polyP, indicating the lack of a significant polyP hydrolysis and suggesting Pi contamination in the polyP substrate used for enrichment.

**Figure 5.** *Ulva pseudorotundata*. Time course of SRP in seawater, as a function of different P regimes. Data are means ± SE (*n* = 3).

In contrast to total internal P content, the soluble protein content increased along the experiment regardless of the P treatment as expected, since all cultures were supplemented with N (Figure 6). This enhancement was especially noticeable after 4 days.

**Figure 6.** *Ulva pseudorotundata*. Time course of soluble protein content as a function of different P regimes. Data are means ± SE (*n* = 3).

Tissue C, N and P composition, as well as the C:N, C:P and N:P atomic ratios are shown in Table 2. These data showed great variations caused by different P supplies. In parallel with thallus expansion and growth, internal C content increased at the end of the incubations. Tissue N content increased by 3–4 times at the end of the experiment, with slightly higher values in the P enriched treatments than in the control. Tissue P content almost doubled the initial values in all P treatments, whereas no changes were observed in the control. As a consequence of the noticeable changes in N and P composition, and, to a lesser extent in C, C:N:P atomic ratios displayed great variations along the experiment. Thus, the C:N atomic ratio dropped from 22.8 down to about 8 in the control, and even more in the P enriched treatments. The C:P values were driven by P changes, being higher in the control than in P enriched treatments. As expected, the largest increase in the N:P atomic ratios was recorded for the control, since no P source was added. However, the N:P ratio also showed higher values following N and P additions than initial thalli, denoting a greater N enrichment than P in the cell composition.

The sPPase activity in cell extracts was assayed along the experiment and scaled on both a protein and DW basis (Figure 7). Remarkable low activities were recorded in the control at the beginning of the experiment (almost undetectable on a DW basis, Figure 7B). There was a transient increase of activity in the control after 2 days when expressed on a protein basis, then this activity was maintained at low levels. The maximum level for Pi treatment was recorded after 4 days, and the highest levels were recorded for PPi and specially polyP enrichment at the end of the experiment (7 days). These values for sPPase activity were at the same order of magnitude as those recorded in field-collected plants (Table 1).

**Figure 7.** *Ulva pseudorotundata*. (**A**) Time course of sPPase activity on a protein basis and (**B**) on a dry weight, DW basis, as a function of different P regimes. Data are means ± SE (*n* = 3).

#### **4. Discussion**

#### *4.1. Growth and Phosphorus Use*

*Ulva pseudorotundata* grew and accumulated P regardless of the P source supplied in the incubations. The long preculture period resulted in nutrient limited thalli, with internal N and P cell quotas close to subsistence values for *Ulva* species [28] and a N:P ratio close to the Redfield ratio. Nutrients (N, P) resupply at the onset of the experimental phase (culture), resulted in thallus growth and cell C and N gain regardless of the P form supplied (treatment) even in the control without P supply, as N was also added. The highest growth rates were achieved under PPi and PolyP additions, which could be explained by the higher content of Pi per mol of substrate of both substrates (i.e., PPi and PolyP, 2 and 13–18 Pi molecules, respectively). The metabolic stimulation driven by N and/or P additions was also reflected by the increase of the FW/area and the drop of the DW/FW ratios. This could be a result of the increase in cell turgor (increased water content) and cell division. Therefore, in the short term, there was a stimulation of growth (on fresh weight and area basis), which can be indicative of cell division in this bi-stromatic sheet like species [19]. Initial tissue C content was lower than mean values for *Ulva* species (about 27 %DW, [29]), but increased in all treatments as a result of the metabolic stimulation caused by N and/or P supply. The growth rates recorded in this study were much lower than maximum ones measured in this species both in the laboratory and in the field (0.2 to 0.3 d<sup>−</sup>1) [19,29]. These low values were probably because thalli were nutrient limited, and the immediate response was, once nutrients were supplied, to replenish cell nutrient quotas before achieving maximum growth rates [30]. In this way, cell N and P increased in all treatments (excepting P for the control, as expected).

Growth supported by different P sources has been previously reported in macroalgae [7]. Phosphorus is taken up through specific Pi transporters. Pi usually comes from SRP or from dissolved organic P molecules broken by alkaline phosphatase enzymes. But if the available P sources are pyro or polyposphates, an external hydrolytic activity is needed. This activity can be mediated through plasma membrane specific pyro or polyphosphatases, and/or from non-specific alkaline phosphatases located in the periplasmic space [31]. These non-specific phosphatases can even be released to culture medium, as shown by [7] in the green seaweed *Cladophora glomerata*. Although we used synthetic seawater, epiphytic bacteria occurring on *Ulva* thalli surfaces, can account for some of the PPi and polyP hydrolysis, also leaving some free Pi available for the macroalgae.

#### *4.2. Soluble Pyrophosphatase Activity*

Soluble cell extracts of field-collected *U. pseudorotundata* showed hydrolytic capacity for a variety of polyphosphates of different chain lengths. The activities measured in this study correspond to the enzymes located in the soluble fraction, as cell extracts were centrifuged and, thus, cell wall and membrane components were eliminated in the pellet. The high catalytic activity values, typical for this enzyme, were on the same order of magnitude or even higher than those recorded in several photosynthetic microorganisms [14]. Hydrolysis of PPi and tripolyposphate yielded similar rates. In contrast, lower values were obtained for trimetaphosphate, probably because of its cyclic structure. Polyphosphatase activity was higher than sPPase as a result of the breakdown of a polymeric substrate, causing a cascade of hydrolytic reactions, as it was observed in the cyanobacterium *Synechocystis* sp. strain PCC6803 [14].

As far as we are aware, there are few data on sPPase activity in marine macroalgae [7]. This sPPase activity has been successfully measured, displaying large changes depending on nutrient status and growth state. Thus, sPPase activity was almost undetectable in long-term nutrient-limited macroalgae. The transient increase in the control, at least on a protein basis, indicated that the response is not only dependent on P supply, but a more general response to metabolic stimulation caused by N addition. The sPPase induction has been observed in *Synechocystis* sp. strain PCC6803. This cyanobacterium showed a biphasic curve: a first sPPase maximum during the exponential growth period, and a second one, of lower magnitude, in the stationary phase [32]. As stated in the introduction, this is a key enzyme to maintain the biosynthetic reactions, and therefore, it was stimulated in response to nutrient supply in *U. pseudorotundata*. In contrast, quite low levels of sPPase activities were recorded in nutrient limited algae (initial values). It should be taken into account that preculture lasted 70 days, leaving sPPase activities at basal levels. After this long period of nutrient limitation, the maximum levels of sPPase activity registered for the Pi, PPi treatments after 4 days and polyP after 7 days showed the stimulation of sPPase activities in parallel to the stimulation of protein synthesis and growth after several days of supply of N and different P sources.

#### *4.3. Ecological Implications*

Some studies have revealed the importance of polyP molecules in contributing to the total available P pool of aquatic sediments [9,10]. However, the measurement of SRP misses these important Pi polymers than are non-molybdenum sensitive. The green macroalga *Ulva pseudorotundata* occurs in dense mats in the intertidal mudflats of the Palmones river estuary (Southern Spain) [29], reducing not only underneath light availability [33], but also (more than 40 times) the net phosphorus flux from sediment to the water column [34]. As cited above, this sediment may contain not only Pi, but also inorganic polymeric forms, that, as evidenced in our study, can support *U. pseudorotundata* growth. In fact, *U. pseudorotundata* eventually exhibits huge biomass accumulations, with canopies that can reach a thallus area index of 17 layers [29]. This is a great barrier for the flux of P between sediment and water. Therefore, it can be concluded that this overlooked P pool must be analysed and considered when studying biogeochemical cycles in coastal ecosystems and nutrient-mass balances, especially when green tides of macroalgae develop. Methods such as 31P-NMR and/or PPi-dependent enzymatic assays

to measure Pi polymers should be adopted in current protocols to understand the phosphorus cycle, especially in human impacted coastal areas.

**Author Contributions:** Conceptualization, J.J.V., F.G.B. and J.L.P.-L.; methodology, J.J.V., F.G.B. and P.H.-P.; lab research, J.J.V., F.G.B. and P.H.-P.; data curation, P.H.-P., J.J.V. and F.G.B.; writing—original draft preparation, J.J.V., F.G.B. and J.L.P.-L.; writing—review and editing J.J.V.; supervision, J.J.V.; funding acquisition, J.L.P.-L. and F.G.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research has been supported by the project CTM2017-85365-R from the Spanish Ministry of Economy, Industry and Competitiveness.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **Appendix A**

**Table A1.** *Ulva pseudorotundata*. Statistic results of 2-way ANOVAs analyzing the effects of Pi treatments and time on the different variables analyzed. \* = *p* < 0.05; \*\* = *p* < 0.01; \*\*\* = *p* < 0.001. n.s. Non-significant. df. Degrees of freedom.


#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Structure and Dynamics of the Ras al Hadd Oceanic Dipole in the Arabian Sea**

**Adam Ayouche 1, Charly De Marez 1, Mathieu Morvan 1, Pierre L'Hegaret 1, Xavier Carton 1,\*, Briac Le Vu <sup>2</sup> and Alexandre Stegner <sup>2</sup>**


**Abstract:** The Ras al Hadd oceanic dipole is a recurrent association of a cyclone (to the northeast) and of an anticyclone (to the southwest), which forms in summer and breaks up at the end of autumn. It lies near the Ras al Hadd cape, southeast of the Arabian peninsula. Its size is on the order of 100 km. Along the axis of this dipole flows an intense jet, the Ras al Had jet. Using altimetric data and an eddy detection and tracking algorithm (AMEDA: Angular Momentum Eddy Detection and tracking Algorithm), we describe the life cycle of this oceanic dipole over a year (2014–2015). We also use the results of a numerical model (HYCOM, the HYbrid Coordinate Ocean Model) simulation, and hydrological data from ARGO profilers, to characterize the vertical structure of the two eddies composing the dipole, and their variability over a 15 year period. We show that (1) before the dipole is formed, the two eddies that will compose it, come from different locations to join near Ras al Hadd, (2) the dipole remains near Ras al Hadd during summer and fall while the wind stress (due to the summer monsoon wind) intensifies the cyclone, (3) both the anticyclone and the cyclone reach the depth of the Persian Gulf Water outflow, and (4) their horizontal radial velocity profile is often close to Gaussian but it can vary as the dipole interacts with neighboring eddies. As a conclusion, further work with a process model is recommended to quantify the interaction of this dipole with surrounding eddies and with the atmosphere.

**Keywords:** ras al hadd oceanic dipole; arabian sea; cyclonic and anticyclonic eddies; altimetric data; angular momentum eddy detection and tracking algorithm (AMEDA); HYCOM model; ARGO floats

## **1. Introduction**

Over the last three decades, there has been growing interest in the Arabian Sea and in its marginal gulfs and seas, for geopolitical, economic, and scientific reasons. These seas form a complex region in terms of oceanographic variability as they are influenced by the North Equatorial Current and by the Southwest Monsoon Current. The monsoon winds, blowing over this region, are highly variable with the seasons, coming strong and northeastward from May to September, weaker and southwestward from December to February, and quite weak and variable in direction during the two intermonsoon periods (March–April and October–November). These winds drive coastal and regional ocean currents, upwellings, and downwellings near the coast. These winds also play a major role in the vertical mixing of the upper ocean water masses [1,2]. The alongshore currents are therefore seasonally variable too. During the summer (southwest) monsoon, the Oman Coastal Current (OCC) flows northeastward, southeast of the Arabian Peninsula, and the West India Coastal Current (WICC) flows southward. In winter, the WICC flows northward (see Figure 1a).

**Citation:** Ayouche, A.; De Marez, C.; Morvan, M.; L'Hegaret, P.; Carton, X.; Le Vu, B.; Stegner, A. Structure and Dynamics of the Ras al Hadd Oceanic Dipole in the Arabian Sea. *Oceans* **2021**, *2*, 105–125. https://doi.org/10.3390/ oceans2010007

Academic Editor: Michael W. Lomas Received: 1 December 2020 Accepted: 27 January 2021 Published: 4 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Satellite observations and high resolution, primitive equation model simulations reveal that a strong turbulent activity takes place near the ocean surface in the Arabian Sea and in its marginal gulfs. The mesoscale (the mesoscale in the ocean characterizes motions of horizontal size from 10 to 250 km and time scales of 2 to 30 days) surface eddies (eddies are oceanic vortices) can be generated by the instability of alongshore currents [3,4] and are influenced by the wind stress curl and by Rossby waves [5,6]. These Rossby waves are generated at the western coast of India [7] and propagate westward towards the south Arabian or Somali coasts. There, they can break and intensify pre-existing eddies. Since both these Rossby waves, and the coastal currents, are surface intensified, the eddies thus formed are also concentrated near the ocean surface. Some of these eddies have already been studied, some of them have received less attention.

**Figure 1.** (**a**) Regional map of the Arabian Sea—in summer—displaying the main oceanographic phenomena: Currents, outflows, upwellings, and eddies and (**b**) snapshot of the Ras al Hadd dipole in summer 2020 with the cyclone to the northeast and the anticyclone to the southwest (from "Ocean Virtual laboratory web portal (ESA, OceanDataLab, https://ovl.oceandatalab.com (accessed on 1 December 2020)")).

Among the most often studied eddies in the region is the Great Whirl, a large surface eddy near the coast of Somalia. It is generated, and its intensity is modulated, by the summer monsoon wind stress curl, and by the incoming Rossby waves [8–10]. In the central Arabian Sea, Fisher et al. [11] observed several surface eddies via a moored currentmeter array. These eddies were generated during the summer monsoon wind stress. In the Gulf of Aden, the formation of eddies and their westward propagation are often due

to wind stress during the summer monsoon [12] but current instabilities play an essential role in eddy dynamics during the fall intermonsoon and Rossby waves play a role during the winter monsoon. In November, eddies in the eastern Gulf of Aden are formed by a northward inflow of the Somali Current through the Socotra Passage [13]. In the Gulf of Aden, Bower and Furey [14] discovered two eddies, the Summer Eddy and the Lee Eddy, formed by the splitting of the Gulf of Aden Eddy at the end of the Spring intermonsoon. The Summer Eddy propagates westward into the Gulf of Aden, while the Lee Eddy remains fairly stationary north of Socotra Island and these two eddies are influenced by the wind stress curl.

The importance of these mesoscale eddies lies both in their interaction with the regional dynamics, or with the atmosphere, and in their effect on the marginal sea outflows. Indeed, the Red Sea and the Persian Gulf produce salty waters which are exported into the Arabian Sea via the straits of Bab el Mandeb and of Hormuz. These two water masses are the RSOW (Red Sea Outflow Water) also called RSW (Red Sea Water) and the PGW (Persian Gulf Water). After exiting the marginal seas, the RSW and PGW outflows mix with the surrounding water masses in the Gulf of Aden and in the Gulf of Oman, respectively, then, the RSW plume stabilizes at 600–1000 m and the PGW plume at a 250–350 m depth [15–17]. Due to the influence of mesoscale eddies or bottom topography [4,18–20], the PGW or RSW outflows are destabilized and they often shed smaller-scale eddies (of radii ranging between 5 and 20 km). These eddies, formed at a depth near the coast, drift away into the open ocean, carrying their warm and salty waters (see [3,21,22] for the small RSW eddies, and [23–25] for the PGW eddies). These small eddies contribute to the regional budgets of heat, salt, and momentum.

A less studied surface mesoscale eddy, or more precisely a pair of eddies, is the Ras al Hadd oceanic dipole. This association of a cyclone and of an anticyclone, exists from May to October, near Ras al Hadd, the southeastern-most cape of the Arabian Peninsula (see Figure 1b). An intense offshore jet, the Ras al Hadd jet, lies between the two eddies. It was studied by Böhm et al. [26] using high resolution radiometric data and Acoustic Doppler Current Profiler observations. Though Böhm et al. mentionned the presence of the two eddies near the jet, and noted their size of about 150 km, they mostly concentrated on the jet which advects cold water offshore during the summer monsoon. Indeed, this dipole marks the eastern edge of the Oman upwelling system which contains cold water. Apart from the study by Böhm et al., the Ras al Hadd dipole has not been studied in detail. Taking advantage of a new, high-resolution altimetric dataset, analyzed via vortex identification and tracking software (AMEDA: Angular Momentum Eddy Detection and tracking Algorithm), of ARGO profiling float data in the region (from the ARGO program), of a regional primitive equation model (HYCOM), the structure and the evolution of the Ras al Hadd dipole are examined here.

This paper is organized as follows: First, we present satellite and profiling float data, the AMEDA algorithm, and the HYCOM model simulations (Section 2). Then we analyze the Ras al Hadd dipole by studying its life cycle over a year (April 2014 to April 2015), and its horizontal and vertical structures using the altimetric and the float data, and the model results (Section 3), evaluating its inter-annual variability over a 15 year period. These results are discussed to present the impact of this dipole on its environment (Section 4). Finally, conclusions are drawn.

## **2. Materials and Methods**

Time series of sea surface height data were obtained via a combination of up-to-date measurements by altimetric satellites (Topex/Poseidon, ERS-2, GFO, Jason-1, Envisat, Jason-2, Cryosat-2, Altika et HY2A). Using all available missions, an ADT (Absolute Dynamic Topography) product for the Arabian Sea was processed by CLS-Argos on a <sup>1</sup> 8 ◦ Mercator grid, with time intervals of 24 h, for the 2000–2015 period. Those products were processed by SSALTO/DUACS and distributed by AVISO+, https://www.aviso.altimetry. fr (accessed on 1 December 2020)) with support from CNES, Centre National d'Etudes

Spatiales, Toulouse, France The geostrophic velocity fields were then derived from absolute dynamic topography. The MDT used to calculate the ADT from SLA is described on the AVISO/CNES website https://www.aviso.altimetry.fr/en/data/products/auxiliaryproducts/mdt/mdt-description.html (accessed on 1 December 2020) and results from an average of sea surface height above the geoid over the 1993–2012 period. It takes into account the last geoid model GOCO05S (based on the complete GOCE mission and on 10.5 years of the GRACE mission) and 25 years of altimetric and in situ data (hydrological profiles, drifting buoys). The spatial resolution of this disseminated regional data set resolves the internal deformation radius, *Rd* ≈ 50–60 km, in the Arabian Sea. The spatiotemporal interpolation performed on the multiple altimetric tracks in order to build gridded products provides this <sup>1</sup> 8 ◦ resolution while the signal resolution of altimetry is about 100 km [27,28]. The sensitivity test performed on the AMEDA algorithm [29] shows that only eddies with a characteristic radius larger than 25 km could correctly be detected with the regional AVISO/CMEMS data set. In addition, a recent work [30] compared eddies identified in a <sup>1</sup> 60 ◦ resolution model with those found in the <sup>1</sup> 8 ◦ AVISO product used here and concluded that eddies with radii smaller than 25 km could not reliably be detected. The Ras al Hadd dipole eddies have radii larger than 45 km. Thus these two eddies can be detected with confidence in this AVISO/CMEMS product.

To relate vortex generation or intensification with the wind stress curl, the ECMWF (European Center of Medium-range Weather Forecast) meteorological model re-analyses (ERA-Interim) were used [31]. These re-analyses provided the regional wind every 6 h, with a <sup>1</sup> 4 ◦ resolution.

To support our analysis of the Ras al Hadd dipole, we also used HYCOM (Hybrid Coordinate Ocean Model) model simulations. These model simulations were performed by SHOM (Service Hydrographique and Oceanographique de la Marine, France). The model has a 5 km horizontal resolution, 40 vertical levels (the upper ones being sigma levels—or terrain following levels, and the lower ones, isopycnic levels—that is, isodensity levels). The atmospheric forcing for the ocean model was extracted from the Meteo-France atmospheric model results over the region with a <sup>1</sup> 4 ◦ horizontal resolution and a 6 h period for the year 2011. This year was chosen as typical of the meteorological variability in this region. The HYCOM model was initialized using a situation obtained from a global, lower resolution ocean model (MERCATOR PSY3, with <sup>1</sup> 4 ◦ resolution). The results of MERCATOR PSY3 were also used for the external boundary conditions of HYCOM (temperature, salinity, and currents). The geographical domain of the HYCOM model covers the Arabian Sea, the gulfs of Oman and of Aden, the Red Sea, and the Persian Gulf, north of 10◦ S, and west of 80◦ E. The tidal forcing of HYCOM was provided by TOPEX and was included in the model at the open boundaries. Vertical mixing was parameterized with the KPP (K-profile parameterization) scheme [32]. This numerical model has been described and validated in [18].

To identify vortices, and in particular the Ras al Hadd dipole, in the altimetric data, we used the AMEDA algorithm. The AMEDA (Angular Momentum Eddy Detection and tracking Algorithm) method is an algorithm that detects and tracks oceanic vortices in the horizontal (geostrophic) velocity fields derived from satellite measurements of the sea surface [29]. The AMEDA method is available and decribed on the website of the DYNED ATLAS project https://www.lmd.polytechnique.fr/dyned/methods (accessed on 1 December 2020) with its strengths and weaknesses, and in particular, a comparison of the geostrophic and cyclogeostrophic velocities. This comparison is also performed in [33]. The vortex centers are identified as extreme in local angular momentum, and each vortex boundary is determined via the connexity of the vortices thus found (the boundary must be a closed streamline). Tracking vortices in time relies on the spatial closeness and on the shape similarity between two vortices successively identified in time. Nevertheless, short gaps in the time series are allowed provided that the temporal continuity can be ensured (the distance between two successive positions of a given vortex must not exceed a threshold). The model also identifies merging (connection) and splitting of vortices, via the

connection or separation of contours. This algorithm contains a minimal number of tunable parameters. It has been successfully tested on AVISO altimetry datasets, on ROMS model output, and on velocity fields issued from laboratory experiments. Here the AMEDA algorithm was applied to the high resolution altimetric maps.

Finally, ARGO profiling float data were extracted from the Coriolis database with only validated data with a high quality index that corresponded to the period of study, used. The ARGO profilers are deep reaching oceanic floats, equipped with temperature, conductivity (salinity) and possibly oxygen concentration sensors. These floats first dive to a given depth (the parking depth, usually between 500 and 1000 m), then they drift at this depth under the influence of the oceanic currents, for a prescribed duration (their cycle period, usually 5 or 10 days). Once every cycle, the ARGO floats move up to the sea surface, thus performing a vertical profile of temperature, salinity, and sometimes oxygen concentration, of the ocean. At the surface, the floats are positioned via GPS and they transmit their data back to the scientific institutes, via a satellite link. Then, they dive back to their parking depth for a new cycle. Therefore the vertical profile data are positioned but the exact trajectory of the ARGO floats at depth is unknown. The ARGO data used here (temperature and salinity only) were collected by floats with cycling periods of 5 or 10 days, and a parking depth of 1000 m. The Ras al Hadd AMEDA outputs for the studied period with their respective Argo profiles given in https://github.com/aayouche/RASALHADD (accessed on 1 December 2020).

## **3. Results**

#### *3.1. Seasonal Variability of the Ras al Hadd Dipole from Mid-April 2014 to Mid-April 2015*

3.1.1. Origin and Evolution of the Two Eddies Composing the Dipole, Obtained from Surface Data Analysis

The AMEDA algorithm applied to the altimetric data (the satellite measurements of the sea surface elevation) allows the tracking in time of the eddies detected. This tracking can thus be reversed in time to determine the origin of the two eddies (the cyclone and the anticyclone) which form the Ras al Hadd dipole. These two eddies are detected in the altimetric data as early as April 2014. At that time, the two eddies lie far from each other, in the Arabian Sea. The cyclone is then at the mouth of the Gulf of Oman, near the Iranian/Pakistani coast. The anticyclone is then located east of Ras al Hadd. Altimetric maps show that the two eddies slowly migrate during the 2014 summer monsoon (see Figure 2). The cyclone (CE) drifts southeastward (early to mid-June), then it remains at the same location and elongates until mid-August. The anticyclone (AE) drifts southeastward early to mid June and then remains steady (at the same location). It intensifies in July due to the interaction with a few neighboring eddies and then weakens in mid-August.

During their migrations, the cyclone and anticyclone are altered by their mutual interaction, and by interactions with other eddies. For instance, in June 2014, the anticyclone becomes strongly elongated under the influence of a cyclone (see Figure 3a). The aspect ratio of the anticyclone (obtained from a best elliptical fit) increases beyond a value of 3, in early June 2014 (Figure 3d). Hydrodynamical theory states that an elliptical vortex with an aspect ratio larger than 3 is prone to breaking into two separate vortices [34]. And indeed, a temporary breaking of the anticyclone is observed in June.

In April 2015, the cyclone merges with another cyclone, which is then located to its east (see Figure 3b). Hydrodynamical theory states that two eddies can merge when the distance between them becomes smaller than 3 times their radii. In April 2015, the two cyclones are closer than this distance (see Figure 3c).

AMEDA also provides the time evolution of the cyclone and anticyclone radii, for a year (note that these two eddies may be associated in the dipole or separate during that period). Notable time variations of these radii reflect in part the successive interactions of the dipole with neighboring eddies (see Figure 4). In particular, we can see the growth

of the cyclone by merger with a cyclone nearby, in August 2014 and in addition, what is noticeable is the growth of the cyclone radius from September to November 2014. A possible mechanism explaining this intensification is presented below (the effect of the wind). Between mid-January and early March 2015, the cyclone is elongated and loses mass due to the shear exerted by other eddies. In April 2015, it merges again with a cyclone nearby.

**Figure 2.** Snapshots of the Absolute Dynamic Topography (ADT) anomaly for the 2014 summer monsoon; the blue and red contours are for anticyclones and cyclones respectively; the dashed contours circle the eddy core; and solid lines indicate the outermost closed contours. The black dots represent the anticyclone and cyclone instantaneous center positions.

**Figure 3.** *Cont*.

**Figure 3.** (**a**) Sea surface height map showing the elongation of the Ras al Hadd anticyclone (green contour, eddy 227) under the influence of cyclone 35, on 5 June 2014; (**b**) sea surface height map showing the merging of the Ras al Hadd cyclone with eddy 35 on 7 April 2015 (green contour); (**c**) ratio of mutual distance to the mean radius of the two cyclones for the merging event; and (**d**) aspect ratio of the anticyclone for the splitting event.

**Figure 4.** Time series of the Ras al Hadd cyclone and anticyclone radii from April 2014 to 2015.

The anticyclone of the dipole evolves as follows: In early October 2014, it grows in a favorable environment (with neighboring anticyclonic eddies) then, it splits and decreases in size. It intensifies again at the end of December 2014, early January 2015, in a favorable anticyclonic environment (via mass transfer between eddies). It decreases thereafter by elongation due to the shear of its cyclonic partner. Finally, it intensifies again via a merger with a neighboring anticyclone (in February–March 2015). These evolution of the two eddies of the dipole show the importance of vortex-vortex interactions in this region.

To explain another part of the eddy radii and intensity variations with time, we investigate the effect of the wind on the cyclone and on the anticyclone. We correlate the absolute dynamic topography anomaly (the Absolute Dynamic Topography (ADT) anomaly is the total ADT minus its instantaneous geographical average in a subdomain of 5 × 5 degrees around the eddy) and the wind stress curl for the Ras al Hadd dipole. In summer, the wind stress curl correlation was often positive with the cyclone and negative with the anticyclone (Figure 5). Indeed, during this period, the summer monsoon winds blow from Africa to Asia, and veer northward around Ras al Hadd. At this location, the wind stress curl is therefore positive in summer and early fall. During the winter monsoon and spring inter-monsoon, the wind direction is more variable and therefore the wind correlation with the eddies fluctuates more (see again Figure 4).

Now, we study the internal structure of each eddy. The internal velocity-pressure balance in a nearly circular eddy is the cyclogeostrophic balance:

$$V\_{\theta}^{2}/r + f\_{0}V\_{\theta} = \begin{array}{c} (1/\rho) \, dp/dr \end{array}$$

where *V<sup>θ</sup>* is the azimuthal (tangential) velocity, *r* is the local radius of the point considered in the vortex, *f*<sup>0</sup> is the Coriolis parameter, *ρ* is seawater density, and *dp*/*dr* is the radial gradient of pressure. The first term in this equation is the centrifugal acceleration, the second is the Coriolis acceleration, and the right hand side in the pressure gradient. For large vortices, the balance lies essentially between the Coriolis acceleration and the pressure gradient and the equilibrium is geostrophic. For smaller and faster spinning eddies, the first term, the cyclostrophic correction to geostrophy, is not negligible. We calculate the ratio of the centrifugal to the Coriolis accelerations for the two eddies, when they are well developed (when their radius is close to 100 km):

$$\text{Če/}\text{Čo} = V\_{\theta}/f\_0R$$

to estimate how nonlinear this balance is. This ratio is shown on Figure 6 with respect to the eddy radius, for each of them. The cyclostrophic correction reaches 20% in the eddy core (10 km from the eddy center) and decreases to about 11% for cyclones and 7% for anticyclones beyond the radius of maximal velocity. We note that it never reaches a large value (e.g., 50%). Therefore, the geostrophic balance is a good approximation for the internal dynamics of these two eddies.

**Figure 5.** (**a**) Cross-correlation between the wind stress curl and the cyclone of the Ras al Hadd dipole with days starting in April 2014; (**b**) same as (**a**) for the correlation of the wind with the anticyclone of the dipole.

**Figure 6.** Ratio of the centrifugal to the Coriolis accelerations *Vθ*/ *f*0*R* for the cyclone and the anticyclone of the Ras al Hadd dipole.

## 3.1.2. Vertical Structure of the Dipole from the HYCOM Model Results

To complement this surface analysis of the dipole, we present its vertical sections of temperature, salinity, density, and perpendicular current issued from the HYCOM model (see Figure 7). The HYCOM model reproduces well the location and intensity of the Ras al Hadd dipole during the summer monsoon period (see Figure A1). This numerical model has been validated with in-situ observations in previous studies and showed that the main water masses (e.g., the PGW) are well reproduced in the Gulf of Oman and in the northern Arabian Sea [24]. The sections presented in Figure 7a are snapshots extracted from the model. Firstly, the upwards or downwards curvature of the main pycnocline is clearly seen for each vortex, in particular, near a 200 m depth though a few small scale density features are observed near the surface (Figure 7a,b). The 200 m depth corresponds to a marked thermocline and halocline (associated with the presence of Persian Gulf Water), above which the Ras al Hadd eddy velocities are intensified and can reach 0.5 m/s. This magnitude compares well with that measured by ADCP [26]. The velocity field is slightly asymmetric with respect to each eddy vertical axis, but this is due to the deformation effect induced by the surrounding flow (Figure 7c). The vortex flow has a baroclinic component but no flow reversal occurs, at least down to 1500 m depth (not shown). The flow coherence above a 250 m depth indicates that the surface information (eddy drift, interaction or deformation) correctly characterizes the deeper structure of each eddy (Figure 7d).

**Figure 7.** (**a**) Snapshot of the sea surface height (SSH) anomaly (HYCOM). The black dashed line represents the location of the vertical section. (**b**–**d**) Represent the temperature, salinity, and the cross-section velocity for the cyclone and anticyclone of the Ras al Hadd dipole, respectively. The black solid lines represent the isopycnals.

## *3.2. Interannual Evolution of the Dipole Over the 2000–2015 Period*

#### Variations of the Surface Features of the Dipole

The AMEDA algorithm is now applied to altimetric data over the 2000–2015 period, to study the long-term variations of the Ras al Hadd dipole structure. During this period, years 2005 and 2009 are not presented because of the dipole position (southwest of the Ras al Hadd cape; see Figures A2 and A3). This dipole is long-living and similar in intensity (radius) to the the Ras al Hadd dipole observed during the other studied years. In this paper, we study the dipole only when it is located at the Ras al Hadd cape at the entrance of the Gulf of Oman (GO), which is not the case for these years. Another dipole is observed at the entrance of GO but is short lived (less than 3 months), therefore for these reasons the dipole is not presented for years 2005 and 2009. For both eddies composing the dipole, we

show the time variation of: (1) The radius *Rmax* of (2) the maximal azimuthal velocity *Vmax* (defining each eddy core), (3) the steepness *α* of each eddy velocity profile fitted on:

$$V\_{\theta}(r) = \frac{V\_{\text{max}}r}{R\_{\text{max}}} \exp(\frac{1}{\alpha}) \exp(-(r/R\_{\text{max}})^{\alpha}/\alpha).$$

and (4) the dynamical Rossby number *Ro* = *<sup>ζ</sup><sup>c</sup>* <sup>2</sup> *<sup>f</sup> <sup>e</sup>*<sup>−</sup> <sup>1</sup> *<sup>α</sup>* where *ζ<sup>c</sup>* is the maximal relative vorticity of each eddy and *f* is the Coriolis parameter. The calculation of the steepness *α* is accurate only for large eddies; this is the case of the two eddies composing the dipole. Figure 8 shows that, over 14 years, the anticyclone of the Ras al Hadd is slightly wider and spins often faster than the cyclone. A direct correlation with the interannual variability of the summer wind stress averaged over the Northwestern Arabian Sea (see Figure 7a of [35]) was not found for the dipole. This may be due to the fact that the dipole is both a local feature and that its intensity also varies by interaction with the surrounding eddies. For the anticyclone, the velocity profile is most often Gaussian (*α* ∼ 2), but for 2011 when this anticyclone was deformed, it created locally intense (*α* ≥ 3) and locally weaker (*α* ∼ 1) vorticity gradients. For the cyclone, the velocity profile was also often Gaussian. In 2004 and 2012, it was less steep (*α* ∼ 1), either because the cyclone was weaker (2012) or because it was not interacting much with neighboring eddies. The nonlinearity of oceanic eddies is also measured by the dynamical Rossby umber: During the summer monsoons of these 14 years, the values of *Ro* were below 0.3; note that the velocity deduced from altimetry is weaker than that measured at sea because of the low spatial resolution of altimetric data. Nevertheless, we can state that the geostrophic balance holds to a good accuracy for both eddies.

**Figure 8.** *Cont*.

**Figure 8.** *Cont*.

**Figure 8.** 2D (surface) statistics for both eddies of the Ras al Hadd dipole with error bars (vertical segments at the end of the colored lines): (**a**) Maximal azimuthal velocity *Vmax*, (**b**) steepness of the radial profile of azimuthal velocity, (**c**) dynamical Rossby number *Ro*, and (**d**) radius of maximal velocity *Rmax*. The horizontal gray line in the steepness plot indicates a Gaussian profile of velocity and the horizontal gray line in the *Rmax* plot indicates the value of the internal Rossby radius of deformation.

#### *3.3. Variations of the Vertical Structure of the Dipole*

The vertical structure of the Ras al Hadd dipole is now evaluated using ARGO float co-localization in the two eddies. The co-localization is performed during the summer monsoon for 6 years (from 2010 to 2015). In the inner core of each eddy, the colocalized profiles give the structure of potential density and a few background profiles are obtained outside of the eddy (at 150 km from the eddy outermost closed contour). At each time step (when ARGO data are available), the inner core density anomaly is computed (*ρeddy* − *ρbackground*). The background state is calculated out of the outermost eddy closed contour (using an average over 5 points). Then a time average is calculated (over the summer monsoon) and over all the detected profiles. The number of ARGO floats co-located with the dipole varies from year to year. It is indicated on Figure 9b. The eddy depth is defined here as the depth of the maximum density anomaly. On Figure 9b, we can see that the anticyclone is globally deeper than the cyclone. The anticyclone reaches a 280 m depth for most of the 6 years while the cyclone depth is closer to 225 m. The year 2010 is of interest for the study; the number of co-localized profiles reaches 30 (Figure 9b) for both the cyclone and the anticyclone. For this year, the vertical distribution of the density anomaly (Figure 9a) shows a strong intensification near the surface (in the upper 50 m); these are mixed layer variations of the eddy density anomalies. They are negative for the anticyclone and positive for the cyclone. The density anomaly remains noticeable for the anticyclone and cyclone even below a 250 m depth. The thickness of each eddy can be computed from the vertical profiles of density anomalies. Choosing a 0.1 kg/m<sup>3</sup> value (an e-folding scale from the maximum) as a threshold for each eddy thickness, we obtain a thickness of 250 m for the anticyclone and of 300 m for the cyclone.

**Figure 9.** (**a**) Mean density anomaly (error bars for standard deviation (STD)) in the cyclone and anticyclone and (**b**) dipole depth as defined by the depth of the maximal density anomaly (the number of co-localized ARGO float profiles is indicated in white bold characters).

#### **4. Discussion**

The Ras al Hadd dipole is a recurrent feature composed of two eddies (a cyclone to the northeast and an anticyclone to the southwest) near the southeasternmost cape of the Arabian Peninsula, Ras al Hadd. It appears in the summer monsoon and breaks up in winter. Its position (at the cape or west of the cape) can vary from year to year. Its radius and azimuthal velocity also vary from year to year. We have found that, seasonally, the local wind stress curl intensifies the cyclone and weakens the anticyclone. It is likely that the anticyclone is strengthened both by the circulation at the upwelling front and by an upwelling Rossby wave (like the Great Whirl); a formal proof of this should be searched in forthcoming studies. Interannually, a correlation of the dipole intensity with the regionally averaged wind stress has not been evidenced. However, as Piontkovski et al. [35] indicate, the wind stress is spatially heterogeneous so that a regional average may not represent the local effect. The recurrence of this eddy dipole is also related to the Ras al Hadd front which separates the warm waters of the Gulf of Oman, from the colder water of the upwelling south of Oman. By inspecting the dipole life cycle, it was shown that each eddy evolved via strong interactions with neighboring eddies. The two eddies are fairly wide (100 km) and their Rossby number is small, so that the cyclostrophic correction to geostrophy is small. The radial profile (the velocity steepness profile) also evolves in time with the eddy interactions. It becomes steeper as the eddy periphery is eroded under the influence of neighboring eddies. It becomes smoother as the eddy undergoes few or no interaction and as dissipation slowly acts. On average over the years, the anticyclone is slightly steeper than the cyclone. In most cases, the two eddy velocity profiles are close to Gaussian (with a steepness in the range [2.25, 2.4]). The Burger Number of each eddy (*Bu* = ( *Rd Rmax* )2) is of order unity, and the Rossby number is small, which shows that quasi-geostrophic dynamics should apply to model the evolution of these eddies, which are stable [34].

An order of magnitude of the kinetic energy and heat content stored in this dipole can be obtained. Assuming that the azimuthal velocity also decays vertically as a Gaussian, we obtain, from our eddy model, and with the maximal velocities observed, a kinetic energy on the order of:

$$K = \frac{1}{2} \rho \,\,\pi \int\_0^\infty V\_0^2 \, (r/R)^2 \, \exp[-2(r/R)^2] \, r \, dr \int\_{-\infty}^0 \exp[-2(z/H)^2] dz$$

with *<sup>R</sup>* <sup>=</sup> <sup>√</sup>2*Rmax*, *<sup>V</sup>*<sup>0</sup> <sup>=</sup> <sup>√</sup>2*eVmax*, *<sup>H</sup>* <sup>=</sup> <sup>250</sup>*m*. For each eddy, this leads to:

$$K = \frac{\sqrt{\pi}}{8} K\_{0\prime} \text{ with } K\_0 = \pi R^2 H \rho V\_0^2 / 2.$$

Taking *<sup>ρ</sup>* = 1025 kg/m3, *<sup>R</sup>* = 5 × 104 m, *<sup>H</sup>* = 250 m, and *<sup>V</sup>*<sup>0</sup> = 1.15 m/s, we have *<sup>K</sup>* ≈ <sup>10</sup><sup>14</sup> <sup>J</sup> for each eddy.

Via geostrophy and hydrostatics and assuming a vertical dependence of the azimuthal velocity as exp(−(*z*/*H*)2), we obtain:

$$\rho'(r,z) = -\text{Ro}/\text{Bu}\,(z/H)\,\exp(-\left(r/R\right)^2)\,\exp(-\left(z/H\right)^2).$$

where *ρ* is the density anomaly. Such a vertical dependence can indeed fit the vertical variation of the density anomaly shown in Figure 9a. Via a linear equation of state for the upper ocean, assuming that temperature variations follow density variations, we have:

$$T'(r, z) = T\_0 \left( z / H \right) \exp(-\left( r / R \right)^2) \exp(-\left( z / H \right)^2),$$

where *T* is the temperature anomaly associated with each eddy, and *T*<sup>0</sup> = 2 ◦C its maximum. The total heat (anomaly) content of each eddy is therefore:

$$Q = \rho C \int\_0^\infty r \, dr \int\_{-\infty}^0 dz \, T'$$

or *<sup>Q</sup>* <sup>=</sup> <sup>√</sup>*πMeCT*<sup>0</sup> where *Me* is the eddy mass and *<sup>C</sup>* is the heat capacity of seawater. This yields *<sup>Q</sup>* = 2.85 × <sup>10</sup><sup>19</sup> J for each eddy.

Both values of eddy kinetic energy and of heat content are one order of magnitude smaller than those of the Gulf Stream rings. However the lifetime of the Ras al Hadd dipole (about 6 months) renders its contribution to the regional budgets of eddy kinetic energy and heat non negligible. Its heat fluxes towards the atmosphere should be evaluated in further work.

Furthermore, the Ras al Hadd eddies have a dynamical influence on deeper water masses, though the azimuthal velocity of these eddies is intensified in the upper 250 m of the water column. This surface intensification is due to at least two mechanisms: The stronger energy of the wind forced ocean above the thermocline, and the presence of a heterogeneous water mass at depth (PGW), which creates a strong pycnocline, attenuating the vertical penetration of the energy. The values of the azimuthal velocity for both eddies found in the present study are comparable with those measured at sea (with ADCP's) by [26]. The anticyclone spins slightly faster than the cyclone below a 250 m depth. However both eddies have a noticeable velocity at depth (about 0.25–0.3 m/s), which can influence the subsurface path of the PGW outflow in the region [25]. As noted by [25] and by [1], the strong surface-intensified eddies of the Gulf of Oman and of the Arabian Sea pull PGW fragments away from the coast. This has been modeled hydrodynamically by [19,20]. The surface eddies distort the alongshore outflow and also creates an unstable boundary layer vorticity sheet, which detaches offshore and breaks into subsurface submesoscale eddies containing PGW. These small eddies play a role in the transport of heat and salt from the Gulf of Oman to the Arabian Sea.

#### **5. Conclusions**

The objective of this paper was the description of the seasonal life cycle and of the inter-annual variability of the two eddies composing the Ras al Hadd oceanic dipole. Using altimetric data and eddy identification and tracking software (AMEDA), we were able to determine the characteristics of the two eddies composing it, for a year (2014–2015). We determined their origin and drift from their initial position to their final junction near Ras al Hadd. We showed the importance of eddy-eddy interactions during the lifetime of the dipole. The evolution of the dipole was also driven by the local wind stress curl, by the upwelling circulation, and by the incoming Rossby waves.

On a longer period (about 15 years), we noted that the anticyclone was larger than the cyclone of the Ras al Hadd dipole. Though the maximal internal velocities could reach 50 cm/s, in particular for the anticyclone, the cyclostrophic correction could be neglected as the dynamical Rossby number did not exceed 0.3. This is due to the large horizontal size of the eddies. On average, the two eddies have a Gaussian radial profile of velocity, but the steepness of this profile varies with the interactions undergone by each eddy. The vertical structure of both eddies is intensified above the layer of the Persian Gulf Water but influences it dynamically.

This leads us to propose process studies to deepen our understanding of this oceanic dipole. The quasi-geostrophic approximation holds for both eddies and this could be used first to perform a stability analysis of this dipole, with respect to velocity profile steepness, to the Burger number, assuming that the cyclone and anticyclone of the RAH dipole are isolated. Secondly, the interaction of this dipole with neighboring eddies and with the Rossby waves originating from the eastern boundary of the Arabian Sea deserves a specific analysis. Finally, the dipole interaction with the atmosphere should be assessed in detail, in particular for the heat and freshwater exchanges.

**Author Contributions:** Conceptualization, X.C., C.D.M., A.S. and A.A.; methodology, B.L.V., P.L. and A.S.; software and validation, B.L.V. and M.M.; analysis and investigation, A.A., C.D.M., M.M., P.L. and X.C.; writing—original draft preparation, A.A.; writing—review and editing, X.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by SHOM contract HYCOM to UBO.

**Data Availability Statement:** The data sources are specified in the text.

**Acknowledgments:** This work is a contribution to the PHYSINDIEN research program.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **Appendix A. Representation of the Ras al Hadd Dipole in the HYCOM Model**

To validate the ability of the HYCOM model to represent the Ras al Hadd dipole, we present here snapshots of the sea surface height anomaly in the model compared with the AVISO MADT anomaly. The dipole location (near Ras al Hadd) was well represented in HYCOM model. A small difference (∼5–10 cm) in the amplitudes of each eddy of this dipole is observed. The HYCOM model represents a more intense dipole compared to AVISO data. Differences may be due to the different horizontal resolution of HYCOM and AVISO, to the HYCOM model atmospheric forcings, or to the parameterization of momentum or heat diffusion in the model. Still, HYCOM reproduces well the known location and radius of the cyclone and anticyclone of the Ras al Hadd dipole. It also provides a proper positioning and extent of the upwelling. Finally, the surrounding eddies in the model and in the data have similar sizes and intensities. At deeper levels (in particular for the Persian Gulf outflow), the HYCOM model has been validated by [24]. Thus we can properly use the HYCOM model results in this region.

**Figure A1.** Snapshots of the ADT anomaly (AVISO) (**top row**). Black red and blue solid contours represent the external shape (closed contour) of the cyclone and anticyclone, respectively. The dashed contours characterize the inner core of these eddies. The black dot represents their centers. Snapshots of the SSH anomaly from HYCOM model (**bottom row**).

#### **Appendix B. Details for Specific Years 2005 and 2009**

In this appendix, the Ras al Hadd dipole locations and intensities are represented for the two specific years 2005 and 2009.

In our study, this dipole is considered to be long-living, intensified in summer monsoon and located near Ras al Hadd. Another important criterion for this dipole is the existence of the jet of Ras al Hadd between the cyclone (north) and anticyclone (south). This jet is also considered intense in summer monsoon season.

**Figure A2.** (**a**) Absolute Dynamic Topography anomaly for the 2005 year (snapshot); blue and red contours are for anticyclones and cyclones respectively, dashed contours are for eddy core, and solid contours for the last closed contour. The black dot represents the anticyclone and cyclone instantaneous center position. The solid red and blue trajectories are for cyclone and anticyclone, respectively during their lifetime. The black triangle and square represent the initial and final position of each eddy, respectively; (**b**) time series of the cyclone and anticyclone radii for 2005.

As can be seen in the following two figures, in 2005 and 2009, only a weak (and intermittent) dipole is observed near Ras al Hadd. The largest and most stable dipole lies southwest of Ras al Hadd. Thus, these two years are considered as "anomalous".

**Figure A3.** (**a**) Absolute Dynamic Topography anomaly for the 2009 year (snapshot); blue and red contours are for anticyclones and cyclones respectively, dashed contours are for eddy core, and solid contours for the last closed contour. The black dot represents the anticyclone and cyclone instantaneous center position. The solid red and blue trajectories are for cyclone and anticyclone, respectively during their lifetime. The black triangle and square represent the initial and final position of each eddy, respectively and (**b**) time series of the cyclone and anticyclone radii for year 2009.

We can suggest several reasons for these anomalous years which would have to be confirmed: A weaker summer monsoon wind would lead to a smaller extent of the upwelling region. Then, the eastern edge of this upwelling, which corresponds to the position of the anticyclone, would be west of its usual position near Ras al Hadd. Another possibility is that the eddy field in this region has interacted with the cyclone and anticyclone to the west of their usual positions and progressively intensified them there. Finally it is possible that the Oman Coastal Current is weaker these years and thus is unable to advect the dipole towards the cape.

#### **References**


## *Article* **Ocean Circulation Drives the Variability of the Carbon System in the Eastern Tropical Atlantic**

**Nathalie Lefèvre 1,\*, Carlos Mejia 1, Dmitry Khvorostyanov 1, Laurence Beaumont <sup>2</sup> and Urbain Koffi <sup>3</sup>**


**Abstract:** The carbon system in the eastern tropical Atlantic remains poorly known. The variability and drivers of the carbon system are assessed using surface dissolved inorganic carbon (DIC), alkalinity (TA) and fugacity of CO2 (fCO2) measured in the 12◦ N–12◦ S, 12◦ W–12◦ E region from 2005 to 2019. A relationship linking DIC to temperature, salinity and year has been determined, with salinity being the strongest predictor. The seasonal variations of DIC, ranging from 80 to 120 μmol kg<sup>−</sup>1, are more important than the year-to-year variability that is less than 50 μmol kg−<sup>1</sup> over the 2010–2019 period. DIC and TA concentrations are lower in the northern part of the basin where surface waters are fresher and warmer. Carbon supply dominates over biological carbon uptake during the productive upwelling period from July to September. The lowest DIC and TA are located in the Congo plume. The influence of the Congo is still observed at the mooring at 6◦ S, 8◦ E as shown by large salinity and chlorophyll variations. Nevertheless, this site is a source of CO2 emissions into the atmosphere.

**Keywords:** carbon cycle; tropical Atlantic; dissolved inorganic carbon; alkalinity

## **1. Introduction**

Tropical regions are strong sources of CO2 to the atmosphere due to the equatorial upwellings, bringing CO2-rich waters to the surface, and to high sea surface temperatures. Among all tropical regions, the Pacific Ocean is the most studied as El Niño events develop there and they are the dominant process governing the interannual variability of the air– sea CO2 flux (e.g., [1]). The Atlantic Ocean has been less studied with relatively fewer observations and modelling analyses [2]. The Western Tropical Atlantic (WTA) and the Eastern Tropical Atlantic (ETA) receive the discharge of the two largest rivers of the world with the Amazon, near the equator in the west, and the Congo near 6◦ S in the east. Numerous studies have highlighted the impact of the Amazon outflow on productivity and carbon uptake in the WTA [3–10]. The dynamics of the Congo plume and its influence on the coastal ocean of the ETA are less documented than the Amazon plume [11]. The Congo plume, like large river plumes, supplies inorganic nutrients, which stimulates biological production. There is evidence of high chlorophyll concentrations (>10 mg m<sup>−</sup>3) near its delta [12]. A few studies reported carbon uptake in the outer plume of the Congo river [13,14]. However, the extent of its influence on carbon properties in the ETA is poorly known.

Here, we focus on the 12◦ N–12◦ S, 12◦ W–12◦ E region defined as the ETA. This region is dominated by the seasonal variation of equatorial and coastal upwellings. The

**Citation:** Lefèvre, N.; Mejia, C.; Khvorostyanov, D.; Beaumont, L.; Koffi, U. Ocean Circulation Drives the Variability of the Carbon System in the Eastern Tropical Atlantic. *Oceans* **2021**, *2*, 126–148. https://doi.org/ 10.3390/oceans2010008

Academic Editor: Michael W. Lomas

Received: 14 December 2020 Accepted: 3 February 2021 Published: 8 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Atlantic cold tongue, which spreads between the African coast and 20◦ W south of the equator, usually forms in June and lasts nearly 5 months [15]. Its formation coincides with the intensification of the southeasterly trade winds. Its spatial extent is maximum in July–September as well as its cooling [15]. The coastal upwelling along the coast between Cape Palmas (Ivory Coast) and Cotonou (Benin) is also seasonal with the main cooling period occurring at the same time, from July to September and a minor upwelling occurring from January to February [16,17].

Earlier studies showed that the tropical Atlantic Ocean has a low productivity except in areas where nutrients are brought to the surface such as the equatorial region and the coastal upwelling [18]. In the Gulf of Guinea, along the northern African coast, the highest chlorophyll concentrations are observed in July–September and coincide with the presence of the upwelling, which leads to large fish catches [19,20]. This productive region is part of the Guinea Current Large Marine Ecosystem extending from the Guinea-Bissau to the Cape Lopez in Gabon [21].

The equatorial upwelling is also a productive region with a main phytoplankton bloom in July–August, when the Atlantic cold tongue is most developed, and a secondary bloom in December as observed by chlorophyll concentrations from ocean color satellite [22]. The surface productivity is explained by the vertical supply of nutrients that fuel the phytoplankton development [23]. When the thermocline is deeper, the vertical supply is more limited and the surface productivity remains low [24]. Chlorophyll and nitrate are the biochemical variables for which extensive data exist and they were used for modelling studies in this region [25].

Carbon observations collected during oceanographic cruises in the 1980s enabled estimates of the source of CO2 in this region (e.g., [26,27]). The seasonality of the flux and the CO2 outgassing have been further documented by hourly monitoring of the fugacity of CO2 (fCO2) at the time-series station at 6◦ S, 10◦ W, e.g., [28]. More recently, a coupled physical-biogeochemical model was used to compare the air–sea CO2 flux between the equatorial Atlantic and the equatorial Pacific [2]. One of the main result was that the variability of the sea surface partial pressure of CO2 (*p*CO2) appears to be driven by sea surface temperature (SST) in the equatorial Atlantic whereas it is driven by dissolved inorganic carbon (DIC) in the equatorial Pacific. However, the authors pointed out that the driving parameter in the equatorial Atlantic was subject to high uncertainty due to the lack of observations. They recommended regional analyses of key variables that regulate oceanic *p*CO2 such as DIC.

The objective of this paper is to analyze the seasonal and interannual variability of the carbon system in the Eastern Tropical Atlantic and the impact of the Congo outflow, using surface carbon measurements from 2005 to 2019 and satellite data. We determine the driving factors of DIC and alkalinity (TA) spatial variations and provide empirical relationships for the ETA.

## **2. Materials and Methods**

As part of the US–France–Brazil prediction and research moored array in the tropical Atlantic (PIRATA) program, five moorings have been deployed to monitor temperature, salinity, wind speed and precipitation in the eastern Tropical Atlantic [29]. The moorings are serviced and replaced every year during the French cruises. Two of these moorings have been equipped with a CO2 CARIOCA sensor to monitor seawater fCO2 at about 1.5 m depth. The fCO2 has been measured hourly by spectrophotometry at 6◦ S, 10◦ W since 2006 [30] and at 6◦ S, 8◦ E from 2017 to 2019. The accuracy of the CARIOCA sensor is estimated at ±3 μatm. The fCO2 data are archived by the Surface Ocean CO2 Atlas project (SOCAT, www.socat.info (accessed on 31 October 2020)). The buoy at 6◦ S, 8◦ E drifted in 2018 and again in 2019; therefore, this site has been abandoned. During the cruises for servicing the PIRATA moorings, seawater samples were taken for inorganic carbon and alkalinity (National Center for Environmental Information (NCEI), Ocean Carbon Data System (OCADS), https://www.ncei.noaa.gov/access/ (accessed on 31 October 2020), accession numbers 0108086–0108091, 0171189–0171191, 0171193–0171197). The samples were analyzed by potentiometric titration with a closed-cell using the method of Edmond et al. [31]. The system was calibrated by Certified Reference Materials (CRMs) provided by Prof. A. Dickson (Scripps Institution of Oceanography, San Diego, CA, USA). The accuracy of TA and DIC is estimated at ±<sup>3</sup> <sup>μ</sup>mol kg−<sup>1</sup> and ±<sup>5</sup> <sup>μ</sup>mol kg<sup>−</sup>1, respectively. In June 2006 and March 2019, seawater and atmospheric fCO2 were measured underway during the PIRATA cruises (EGEE 3 and PIRATA FR-29), using a CO2 system including an infrared LiCor 7000 analyzer [14]. The accuracy of fCO2 is estimated at ±2 μatm. The PIRATA cruises have provided a good coverage of carbon measurements from ship and moorings in the 10◦ S–6◦ N, 12◦ W–12◦ E region (Figure 1).

**Figure 1.** Distribution of DIC observations (2005–2019) and location of the five PIRATA moorings (×). The main zonal surface currents are the eastward Guinea Current (GC) and the westward South Equatorial Current (SEC). The location of the Angola Dome (AD) is indicated.

The sampling covers the two main surface currents: the eastward Guinea Current (GC) in the north and the westward South Equatorial Current (SEC) in the south. The boundary between the two currents is located approximately along 3◦ N [32]. The carbon system includes four parameters: fCO2, DIC, TA and pH. The knowledge of two parameters allows the calculation of the other two using the equations of the carbon system. At the deployment of the moorings and during the EGEE 3 and PIRATA FR-29 cruises, the three parameters fCO2, DIC and TA were measured so that the carbon system was overdetermined. Two parameters are used to calculate the remaining one that is subsequently compared with the measurement. The calculation is performed with the CO2SYS software for Matlab [33] using the dissociation constants of Mehrbach et al. [34] refitted by Dickson and Millero [35].

A number of 119 concomitant measurements of DIC, TA and fCO2 are used to check the consistency of the carbon system (Table 1).

**Table 1.** Statistics of the calculation of a carbon parameter from a given pair of parameters given by the Root Mean Squared Error (RMSE) and the correlation coefficient (*r*) between the calculated and the 119 measured values.


From the three possible pairs of parameters, the fCO2–TA pair gives a relatively low uncertainty to calculate DIC. The pair DIC–TA gives a high uncertainty to calculate fCO2, which is explained by the strong co-variation between DIC and TA (*r* = 0.96). This is in agreement with previous works estimating combined uncertainties in calculated parameters from different measurement pairs and reporting DIC–TA as the worst pair of parameters to calculate fCO2 [36,37].

The air–sea flux of CO2, FCO2, is calculated at the PIRATA moorings as follows:

$$\text{FCO}\_2 = k \ast sol \ast \Delta \text{fCO}\_2 \tag{1}$$

where *k* is the gas transfer velocity depending on the wind speed [38], *sol* is the solubility of seawater CO2 [39] and ΔfCO2 is the difference of fCO2 between the ocean and the atmosphere. A positive flux means a source of CO2 to the atmosphere. The PIRATA sites are equipped with an anemometer but, due to the failure of the sensor, there are some data gaps. Therefore, we use the wind speeds from the Cross-Calibrated Multiplatform (CCMP) product version 2 [40]. They are usually in good agreement with the PIRATA winds.

Empirical relationships are a useful means to estimate data from more accessible parameters, such as temperature and salinity, when no observation is available. In order to calculate DIC and pH from fCO2 at the moorings, we use the following empirical TA–sea surface salinity (SSS) relationship determined by Koffi et al. [41] with 190 observations collected from 2005 to 2007:

$$\text{TA } (\pm 7.2) = 65.52 \, (\pm 0.77) \* \text{ SSS} + 2.50 \, (\pm 27.22), r^2 = 0.97,\tag{2}$$

Using 638 observations of TA from 2005 to 2019, the RMSE is 8.5 μmol kg−<sup>1</sup> and *r* is over 0.99, which confirms the robustness of this TA–SSS relationship for the eastern tropical Atlantic. The salinity-normalized TA given by TA∗35/SSS gives an average of 2298.1 ± 8.1 <sup>μ</sup>mol kg−<sup>1</sup> and is close to the value of 2291 ± <sup>4</sup> <sup>μ</sup>mol kg−<sup>1</sup> given by [42] for the Atlantic between 30◦ S and 30◦ N for SST > 20 ◦C.

Relationships between DIC and physical parameters are investigated through various techniques such as MLR (multiple linear regressions), decision trees, random forest and feed forward neural networks (Appendix A).

Statistical tests are used to compare different groups of data. When the parameters do not follow a normal distribution, we use the Wilcoxon rank sum test, equivalent to the Mann–Whitney U test. In particular, we tested whether the data located north of the equator were significantly different from the data located south of it.

We examine the link between the year-to-year variability and the tropical southern Atlantic (TSA) index. The TSA index is an indicator of the SST in the eastern tropical South Atlantic. It is the average of SST in the box 30◦ W–10◦ E, 20◦ S–0◦ [43].

Satellite images of SST, SSS, chlorophyll and precipitation provide the spatial environmental conditions in the ETA, and their seasonal variability in the region. The chlorophyll concentrations are monthly composites from MODIS (Moderate-Resolution Imaging Spectroradiometer) Aqua on a 4 km grid. The monthly composites of SST are on a 9 km grid. MODIS data are available from July 2002 to January 2020. The sea surface salinity fields are 18-day Gaussian means from SMOS (Soil Moisture and Ocean Salinity Satellite) at a resolution of 25 km (CATDS CEC LOCEAN debias V4.0, doi:10.17882/52804). The data from SMOS cover the period from January 2010 to September 2019. A monthly climatology of SSS is calculated over this period. The precipitations are from the Tropical Rainfall Measuring Mission (TRMM) and are available monthly, on a 0.25◦ grid. The monthly climatology of rain rate is calculated for the 2005–2019 period.

#### **3. Results**

#### *3.1. Environmental Setting*

In the 12◦ S–12◦ N, 12◦ W–12◦ E region, the MODIS climatological SST varies from 24.55 ◦C in July to 27.01 ◦C in March. During the warm season, warm water spreads over most of the basin as shown for the month of March (Figure 2a). From June to September, the surface water is cooler (<25.5 ◦C) as the equatorial and coastal upwellings are well

developed. The months of March and August are chosen to illustrate the contrasts between the two seasons (Figure 2a,b). During the cold season, a central West African upwelling occurs in the northern part of the Gulf of Guinea between approximately 8◦ W and 4◦ E, north of 2◦ N, from June to October [17]. The cooling is visible in the SST map of August as a small band of cooler SST along the coast (Figure 2b). Along the coast of Gabon, around 10◦ E near the equator, the coastal upwelling merges with the equatorial upwelling forming the Atlantic cold tongue (ACT). In the southernmost part of the basin, austral winter cooling leads to the coldest SST.

**Figure 2.** Distribution of the climatological sea surface temperature (SST) for (**a**) March and (**b**) August and chlorophyll a (Chla) for (**c**) March and (**d**) August from Moderate-Resolution Imaging Spectroradiometer (MODIS). The five PIRATA moorings are indicated by the black stars (\*).

During the warm season, the productivity is low and the chlorophyll concentrations remain close to 0.2 mg m−<sup>3</sup> except along the coastline, and in the Congo plume near 6◦ S, 10◦ E (Figure 2c). During the cold season, the chlorophyll concentrations increase along the coast, especially in the coastal upwelling of the Gulf of Guinea, and in the equatorial upwelling (Figure 2d). In the south, near 10◦ S, 9◦ E, the Angola Dome is maintained by the dynamic uplift of the thermocline, and it is also a productive region [44]. Signorini et al. [12] explain the surface distribution of the chlorophyll in the tropical Atlantic by the upwelling, defined as large vertical excursions of the thermocline, and by the river plumes. Both processes supply nutrients that lead to biological activity. Productive waters are advected by surface currents and arrive at the PIRATA moorings with the exception of the mooring at 10◦ S, 10◦ W that remains in very low chlorophyll waters throughout the year (Figure 2c,d). From 2010 to 2019, the MODIS chlorophyll concentrations during the July–September months are significantly higher than during the remaining months of the period (*p*-value < 0.0001), in agreement with previous studies [22]. This confirms that the cold season corresponds to the productive season.

The region is also affected by freshwater input from precipitation and river runoff. The maximum precipitation is an indicator of the position of the Intertropical Convergence Zone (ITCZ). The ITCZ migrates seasonally and is located at its northernmost position in July and at its southernmost position in March. In the ETA, the ITCZ is over the ocean during the warm season from approximately December to June. In March, it is located slightly north of the equator (Figure 3a) whereas in August it is over the land (Figure 3b). The precipitation feeds the numerous rivers that discharge into the Atlantic, which affects the distribution of SSS in the basin. The northern part of the basin is fresher than the southern part because of the presence of the ITCZ, and of the runoff of the rivers (e.g., Niger, Volta) located in the eastern Gulf of Guinea. North of about 2◦ N, the eastward Guinea Current, an extension of the North Equatorial Counter Current (NECC), is relatively fresh as it receives large precipitation due to the presence of the Intertropical Convergence Zone (ITCZ).

**Figure 3.** Distribution of the Tropical Rainfall Measuring Mission (TRMM) rain rate for (**a**) March and (**b**) August and Soil Moisture and Ocean Salinity Satellite (SMOS) SSS for (**c**) March and (**d**) August DIC. The PIRATA moorings are indicated by the black stars (∗).

South of the equator, the areas of freshwater are located along the coastline, associated to high precipitation, and near the Congo mouth at 6◦ S (Figure 3). The Congo River has a bimodal hydrological cycle with its maximum river discharge in December and a secondary peak of discharge in May, whereas the minimum discharge is in August and March [45]. In March, the SSS distribution follows the precipitations with low SSS (<35) associated to high precipitation in the northern and southeastern parts of the basin (Figure 3a,c). In August, the salinity is higher and low salinity (<34) is encountered near the rivers mouths of the Volta and Niger in the north and of the Congo in the south (Figure 3d). During this month of minimum discharge, the extension of low SSS water near the Congo mouth is limited.

The physical conditions in this region show a north–south difference between warmer and fresher waters in the north and colder and saltier waters in the south because of the precipitations, river outflows and the coastal and equatorial upwellings. In the northern part of the basin, the climatological SST ranges between 25.26 ◦C in July and 28.82 ◦C in February and the SSS is below 35.2, reached in August. In the southern part, the SST varies from 22.41 ◦C in July to 27.41 ◦C in February, and the SSS varies from 35.2 in April to 35.88 in July. The north–south gradient is visible on the satellite SST throughout the year. Even during the warm season, as illustrated in March (Figure 2a), the surface water is warmer north of the equator than south of it.

## *3.2. Variability and Drivers of the Carbon System*

The DIC concentrations of the samples collected from 2005 to 2019, west of 10◦ E, range from 1618 to 2133 μmol kg−<sup>1</sup> with the lowest concentrations (<1800 μmol kg<sup>−</sup>1) close to the Congo outflow and the highest (>2100 μmol kg−1) in the southwestern part of the region. The distribution of the DIC samples highlights the lower concentrations north of the equator compared to the south (Figure 1). As TA is strongly correlated to SSS, TA exhibits also a north–south difference. In the north, TA is lower with values ranging from 2114 to 2369 <sup>μ</sup>mol kg−<sup>1</sup> and a mean of 2278 ± <sup>38</sup> <sup>μ</sup>mol kg<sup>−</sup>1. In the south, TA has a mean value of 2324 ± <sup>66</sup> <sup>μ</sup>mol kg−<sup>1</sup> with a maximum of 2433 <sup>μ</sup>mol kg−<sup>1</sup> reached south of 8◦ S, near 4◦ W where evaporation dominates. The lowest TA value of 1832 μmol kg−<sup>1</sup> is found in the Congo plume. Warm and fresh waters are located in the GC system and exhibit lower TA and DIC concentrations whereas colder and saltier water mass transported by the SEC have higher alkalinity and DIC content. The distributions of SST, SSS, DIC and TA north and south of the equator are significantly different (*p* < 0.001).

DIC concentrations vary according to the main system of currents and depend on the water masses present in the region, which are characterized by their SST and SSS. This led Koffi et al. [41] to propose an empirical relationship using data over 2005–2007. However, due to the atmospheric CO2 increase and air–sea exchange, CO2 increases in the ocean over time and the relationship presents a bias. The difference between the observations and the regression increases over time, which confirms that the relationship is no longer valid for recent years. According to the atmospheric CO2 measurements during EGEE 3 and PIRATA FR-29, the atmospheric fCO2 has increased from 367.6 ± 1.7 μatm in June 2006 to 392.2 ± 2.4 <sup>μ</sup>atm in March 2019, which corresponds to a rate of 1.9 <sup>μ</sup>atm yr<sup>−</sup>1.

In order to take into account the CO2 increase over time, we develop a new multiple linear regression (MLR) by adding a time variable:

$$\text{DIC} = -4585.8 \left( \pm 371.16 \right) - 12.46 \left( \pm 0.53 \right) \ast \text{SST} + 54.36 \left( \pm 0.93 \right) \ast \text{SSS} + 2.49 \left( \pm 0.19 \right) \ast \text{Year} \tag{3}$$

This relationship is based on 637 samples collected in the 12◦ W–10◦ E, 10◦ S–10◦ N region from 2005 to 2019. The RMSE is 15.3 μmol kg−<sup>1</sup> and *r*<sup>2</sup> = 0.93. The negative sign of the SST coefficient in Equation (3) reflects the contribution of the upwelling with cold SST leading to higher concentrations of DIC. The same regression is performed after normalization of the data (by subtracting the mean and dividing by the standard deviation) in order to determine the weight of each predictor. The coefficients for SST, SSS and Year are −0.42, 0.75 and 0.20, respectively. The stronger predictor is SSS, which is the dominant factor driving the DIC variations. However, without taking into account the SST, the variance explained would drop from 93% to 87% and the RMSE would increase to 21 μmol kg<sup>−</sup>1.

Equation (3) is determined using discrete samples of DIC and needs to be validated with independent data. Surface fCO2 has been measured underway during two cruises (EGEE 3 and PIRATA FR-29) and at two moorings (6◦ S, 10◦ W and 6◦ S, 8◦ E) on an hourly basis. In order to check the robustness of Equation (3), DIC is calculated from the alkalinity–SSS relationship and the fCO2 measurements. The calculated DIC is then compared with the predicted DIC (Table 2). DIC given by Equation (3) is in good agreement with DIC calculated from underway fCO2. The correlation coefficient is always over 0.96 and the RMSE is within the uncertainty given by Equation (3). The highest RMSE is found at 6◦ S, 8◦ E.

Using monthly MODIS SST and SMOS SSS fields on a 0.25◦ grid, the monthly regional distribution of DIC can be estimated using Equation (3) to map its spatial variations (Figure 4). The monthly climatology of DIC is built for the January 2010 to September 2019 period. Overall, the monthly maps of DIC are similar to the SSS maps with lower values in the northern part of the basin and close to the river mouths.

**Table 2.** Comparison between the DIC predicted by Equation (3) and DIC calculated from measured fCO2 and alkalinity from Equation (2). The statistics include the Root Mean Squared Error (RMSE), the correlation coefficient (*r*) and the number of data (N). The data of the moorings correspond to daily means.


\* RMSE using the coefficients of Equation (3) with all decimal places are calculated in Appendix A.

**Figure 4.** Climatology of DIC (in μmol kg<sup>−</sup>1) obtained from Equation (3) with the monthly fields of MODIS SST and SMOS SSS. The DIC climatology is calculated over the January 2010 to September 2019 period.

Lower DIC concentrations are associated with the northern water mass and range from a mean value of 1924 μmol kg−<sup>1</sup> in April to 2027 μmol kg−<sup>1</sup> in August. Higher DIC concentrations are found in the southern water mass with a mean value of 2001 μmol kg−<sup>1</sup> in April to 2096 μmol kg−<sup>1</sup> in August. The monthly north–south DIC concentrations are significantly different throughout the year (*p*-value < 0.0001). During the upwelling period (July–September), the DIC concentrations are significantly higher than during the warm season (*p*-value < 0.0001). Thus, for the upwelling period, the mean DIC is <sup>2070</sup> ± <sup>70</sup> <sup>μ</sup>mol kg<sup>−</sup>1, whereas it is 2007 ± <sup>28</sup> <sup>μ</sup>mol kg−<sup>1</sup> during the warm season over the 2010–2019 period. In the upwelling season, the north–south contrast in DIC corresponds to a north–south SST difference. For example, in August, warm surface waters (>26 ◦C) are separated from cold waters (<24 ◦C) with the frontal zone located around the equator (Figure 2b). This front is also observed on the DIC maps from June to September (Figure 4).

### *3.3. Impact of the Congo Plume*

The Congo has by far the strongest river discharge in this region. Near the Congo mouth, the carbon system is affected by the Congo discharge and is examined from nine observations taken east of 10◦ E, around 6◦ S, over the period 2013–2018 from March to June. The lowest salinity (27.3 psu) was observed in March 2018, at the time of the secondary minimum of the Congo discharge. There was no sampling during the period of high Congo discharge, in December and May. Using the nine samples near the Congo outflow, the mixing between the river and oceanic waters gives the following equation:

$$\text{DIC} = \\$0.60 \left( \pm \ 1.97 \right) \* \text{SSS} + 231.73 \left( \pm \ 62.19 \right) \tag{4}$$

with a RMSE of 16.0 μmol kg−<sup>1</sup> and *r*<sup>2</sup> = 0.99. The relationship covers the region close to the Congo mouth from 6◦ S to 5◦ S, east of 10◦ E. Further west, carbon data are available at the mooring at 6◦ S, 8◦ E. Hourly fCO2 has been monitored at 6◦ S, 8◦ E from 2017 to 2019 with a CARIOCA sensor. This dataset includes observations from March to August. In addition, DIC samples have been taken at this site. These DIC samples and the DIC calculated from daily averages of fCO2 and TA–SSS tend to depart from the mixing line (Figure 5). The RMSE is 32 μmol kg−<sup>1</sup> (N = 358) and 24 μmol kg−<sup>1</sup> (N = 11) for the calculated DIC and for the DIC samples at 6◦ S, 8◦ E, respectively.

**Figure 5.** DIC–SSS relationship obtained with the samples near the Congo outflow (dashed line) with DIC samples used for the regression (red stars), DIC samples at 6◦ S, 8◦ E (blue squares) and DIC calculated (cyan dots) from daily fCO2 measured at 6◦ S, 8◦ E and TA calculated from SSS at 6◦ S, 8◦ E.

The DIC data at 6◦ S, 8◦ E are better estimated by Equation (3), giving a RMSE of 14 μmol kg−<sup>1</sup> for both the 11 samples taken at the site and the DIC calculated from the daily fCO2. This suggests that the conservative mixing applies near the River mouth but it is no longer valid further west, at 8◦ E. Nevertheless, the influence of the Congo at the mooring is noticeable. Using satellite based data products, Hopkins et al. [46] show that the Congo plume extends northwest along the coastline or west-southwest into the Atlantic transported by the SEC. From MODIS chlorophyll images (Figure 2c,d), high chlorophyll concentrations are observed at the mooring at 6◦ S, 8◦ E in both seasons and originate from the coast near the Congo mouth.

Using all the available data measured at 6◦ S, 8◦ E over the 2006–2019 period, the SSS are grouped by month (Figure 6a). A data gap occurred from 2007 to 2013 as no mooring was deployed. This site exhibits large salinity variations. The lowest salinity of 30.22 psu,

occurs in May 2018 during the secondary maximum discharge of the Congo. The numerous outliers usually point towards lower salinity than the median. In December, during the maximum Congo discharge, the salinity of the outliers ranges from 32 to 29 psu. The large range of salinity of the outliers (>1 psu) is usually observed around the maximum discharge periods, December and May. The median salinity of December is over 35 psu, which indicates that the plume has not reached the mooring yet, except for some years as in 2006 when a salinity close to 29 psu was observed. The lowest SSS median is in January, which suggests that the Congo plume, transported by the westward SEC, reaches the mooring. Based on satellite images, Hopkins et al. [46] found that the maximum offshore extent of the Congo plume occurs between December and April with a January peak. The freshwater supplied by the secondary peak of discharge, in May, tend to move towards the Gulf of Guinea [47]. The minimum Congo discharge occurs in August, which corresponds to the highest salinity at 6◦ S, 8◦ E. The measured values of SSS and DIC, calculated from fCO2 at the mooring in 2017, 2018 and 2019, and averaged for each month, are also indicated (Figure 6).

**Figure 6.** (**a**) Distribution of SSS at 6◦ S, 8◦ E grouped by month using SSS observations from 2006 to 2019 (the data at 6◦ S, 8◦ E are for the periods 28 June 2006 to 10 June 2007 and 5 July 2013 to 8 August 2019). The central mark of the box corresponds to the median with the bottom and top edges of the box indicating the 25th and 75th percentiles. The whiskers extend to the extreme data points and the crosses (×) are considered as outliers. (**b**) Same for the distribution of DIC at 6◦ S, 8◦ E calculated with Equation (3) using the SSS and SST recorded at the PIRATA mooring. The monthly measurements of SSS and DIC calculated from fCO2 at the PIRATA mooring are indicated for 2017 (black stars), 2018 (green stars) and 2019 (red stars).

The DIC variations follow closely the SSS variations with the highest salinity corresponding to the highest DIC in August, the largest variations in January and the lowest DIC concentrations occurring during the January–May period (Figure 6b). The year-to-year variability is very pronounced as shown by the samples taken at this site. In April 2017, the salinity is lower than 1 psu and DIC is higher than 50 μmol kg−<sup>1</sup> compared to April 2018 and 2019. Likewise, in June 2018 a low salinity (34 psu) corresponds to a low DIC (1968 μmol kg−1) whereas in June 2017 and 2019 the salinity is over 35 psu and the DIC concentrations over 2029 μmol kg<sup>−</sup>1.

The distribution of the chlorophyll from MODIS is examined at the 6◦ S, 8◦ E mooring (Figure 7). The chlorophyll concentrations at this site are relatively high, which is in contrast with the other four PIRATA moorings where concentrations are always lower than 1.5 mg m−<sup>3</sup> at 6◦ S, 10◦ W and lower than 3 mg m−<sup>3</sup> at the two equatorial moorings. The lowest chlorophyll concentrations (<0.25 mg m−3) are encountered at 10◦ S, 10◦ W as this site is outside the area of influence of both the Atlantic cold tongue and the Congo plume.

**Figure 7.** Distribution of MODIS chlorophyll at 6◦ S, 8◦ E grouped by month using data from 2002 to 2019. The central mark of the box corresponds to the median with the bottom and top edges of the box indicating the 25th and 75th percentiles. The whiskers extend to the extreme data points and the crosses are considered as outliers.

At 6◦ S, 8◦ E, the highest concentrations of MODIS chlorophyll *a* occur in July– September and, to a lesser extent, in January–February, which is in agreement with the results of Signorini et al. (1999) using SeaWiFS (Sea-viewing Wide Field-of-view Sensor) chlorophyll *a* from 1997 to 1998. The concentrations of the years 2017, 2018 and 2019 are indicated (Figure 7) to compare with the DIC variations (Figure 6b). The chlorophyll *a* from MODIS is highly variable at 6◦ S, 8◦ E with the highest concentrations in August, as shown by the whisker reaching over 7 mg m−<sup>3</sup> (Figure 7), when the concentrations of DIC are the highest. There is no clear link between DIC and chlorophyll concentrations as the highest chlorophyll are not associated with low DIC that would occur if strong carbon uptake were taking place. In July–September, both chlorophyll and DIC concentrations are high. This suggests that biological activity is not the dominant driver of carbon variations at this site.

#### *3.4. Year-to-Year Variability of the Carbon Parameters*

In the ETA, the longest time-series of carbon parameters is the monitoring of hourly fCO2 at 6◦ S, 10◦ W that started in 2006. At 6◦ S, 8◦ E, the monitoring of fCO2 is rather limited with observations from 2017 to 2019 and a large data gap in 2018–2019 due to the drift of the buoy. The daily DIC at these two sites, calculated from daily fCO2, is compared with the DIC from the MLR, and with the samples taken during the PIRATA cruises at both locations (Figure 8a,b). The MLR reproduces rather well the variations of DIC at 6◦ S, 10◦ W but overestimates the DIC concentrations in some years, such as in July–September 2011 where the difference is over 35 μmol kg−1. The timing of the decreases and increases in DIC is in good agreement between the DIC calculated by the MLR and the DIC calculated from underway fCO2. At 6◦ S, 10◦ W, DIC values are higher in July–November and lower in April–June each year. At 6◦ S, 8◦ E the daily variations of DIC tend to show the lowest concentrations in April and the highest in August but with large variations. For example, DIC concentrations vary by more than 150 μmol kg−<sup>1</sup> within days in April 2017 (Figure 8b). In April 2019, the lowest DIC is much higher than in April 2017 with a difference reaching 125 μmol kg<sup>−</sup>1. Nevertheless, the MLR captures well the high frequency variations of DIC.

**Figure 8.** (**a**) Daily DIC determined from the multiple linear regressions (MLR) using SST and SSS recorded at the 6◦ S, 10◦ W mooring since 2006 (in black) and calculated from measured seawater fCO2 and TA–SSS (in cyan) at 6◦ S, 10◦ W The DIC samples are in red. (**b**) as in (**a**) for the mooring at 6◦ S, 8◦ E from 2017 to 2019. (**c**) Monthly DIC calculated with the MLR using SST and SSS recorded at the mooring (in black) and using satellite SST and SSS (in blue) collocated at 6◦ S, 10◦ W. (**d**) same as in (**c**) for the mooring at 6◦ S, 8◦ E from 2017 to 2019.

As expected, given the discrepancy observed in the daily variations, monthly DIC derived from satellite data tend to give higher concentrations than DIC calculated with the PIRATA data but reproduce the seasonal variations of DIC (Figure 8c,d). The monthly DIC smooth significantly the variations. It is especially visible at 6◦ S, 8◦ E on the MLR with a reduced range of variations of DIC, decreasing from 460 to 320 μmol kg−1. The seasonal cycle of DIC at 6◦ S, 8◦ E is more distinct than on the daily variations. Given that the satellite pixels correspond to a much larger area than the observations at the mooring, there is expected discrepancy between the two products but overall, the MLR reproduces well the seasonal cycle and the year-to-year variations of DIC.

After verifying the DIC variations over time at the two sites where observations are available, we examine the monthly DIC variations given by the MLR using satellite data at regional scale. The DIC concentrations are lower in the north than in the south with the maximum values occurring from July to September each year in both regions (Figure 9). The monthly means of DIC over the entire region are closer to the DIC variations observed in the south in Figure 9 because the region north of the equator includes a large section of land and, hence, the weight of this region is much smaller in the global average. The lowest concentrations of DIC occur in March–April in the south, whereas, north of the equator, low DIC concentrations (<1950 μmol kg<sup>−</sup>1) occur over a much longer period, usually from November to May. The difference between the highest and the lowest DIC concentrations range from about 80 to 125 μmol kg−<sup>1</sup> in both regions. This range is much higher than the year-to-year variability of DIC that is less than 50 μmol kg−<sup>1</sup> over the 2010–2019 period. On annual average, over the whole region, the DIC variations are less than 35 μmol kg−1. In both the northern and southern regions, the seasonal variability is responsible for the large variations of DIC (Figure 9).

**Figure 9.** Monthly DIC determined from the MLR using satellite SST and SSS for the region north and south of the equator between January 2010 and September 2019.

The year 2010 exhibits the lowest annual average of DIC in the north (*p*-value = 0.006) and in the south but, in the southern region, it is not significantly different from the other years. The salinity data retrieved from SMOS are of lower quality from January to March 2010, at the beginning of the series. If we remove these 3 months in the time series, the annual mean of DIC in 2010 is no longer significantly different from the other years in both regions. The highest DIC annual average occurs in 2015 north of the equator and in 2017 south of it, but these years are not significantly different from the other years. Overall, the year-to-year variability of DIC does not highlight any anomalous year over the 2010–2019 period.

#### **4. Discussion**

## *4.1. Main Features of the Carbon Parameters*

Both DIC and TA are strongly correlated to physical parameters. The TA–SSS relationship is robust and does not vary much with new observations. The DIC can be estimated by taking into account the water masses, characterized by their temperature and salinity, and an increase in oceanic surface CO2 over time due to the atmospheric CO2 increase. In stratified regions, like the tropics, the oceanic CO2 tends to follow the rate of the atmospheric increase. The MLR used here shows that SSS, SST and year describe well the DIC variations in this region. Other methods, taking into account non-linear processes, do not perform significantly better than the MLR (Appendix A), which implies that the DIC dependency on the SSS, SST, and the year is essentially linear. However, feed forward neural networks show slightly better results with values of the first and last quantiles of DIC, suggesting that a non-linear dependency exists between DIC and input variables (Appendix A). Low concentrations of DIC are mostly encountered in coastal area in river plumes. As most samples are in the open ocean region, the MLR is adequate here. However, a non-linear method would be required when more samples in the coastal zone become available.

The DIC distribution in the 12◦ W–10◦ E, 12◦ S–12◦ N region (Figure 4) shows high DIC concentrations in most part of the ETA during the upwelling season. Overall, the DIC values range from 1945 to 2090 <sup>μ</sup>mol kg−<sup>1</sup> with an average of 2024 ± <sup>37</sup> <sup>μ</sup>mol kg−<sup>1</sup> from 2010 to 2019. The highest DIC values occur during the cold and productive season, in July–September each year (Figure 9). In this season, the MODIS chlorophyll is significantly higher. The vertical supply of subsurface waters is responsible for bringing CO2-rich waters to the surface that spread over the basin. It is the dominant process and the biological activity is not strong enough to counterbalance the carbon supply.

The DIC concentrations in the ETA are similar to DIC observations collected further west. In the Western Tropical Atlantic, in the SEC, DIC concentrations average <sup>2045</sup> ± <sup>40</sup> <sup>μ</sup>mol kg−<sup>1</sup> in April–September and 2058 ± <sup>22</sup> <sup>μ</sup>mol kg−<sup>1</sup> in October–March from analyses made over the period 1989–2014 [48]. In the ETA, south of the equator, the DIC concentrations, between 2010 and 2019, average 2059 ± <sup>39</sup> <sup>μ</sup>mol kg−<sup>1</sup> in April–September and 2036 ± <sup>32</sup> <sup>μ</sup>mol kg−<sup>1</sup> in October–March.

The lowest DIC concentrations in the ETA are observed east of 10◦ E near the Congo mouth. In this area, DIC is better estimated with the conservative mixing Equation (4) than with the MLR. The slope of Equation (4) is different from the ones given by published relationships (Table 3). As previously noticed [49], the value of the slope is sensitive to the salinity range. Using all the data or data with salinity higher than 33 changes the slope from 46.5 μmol kg−1/psu to 50.6 μmol kg−1/psu [49]. It also changes the end-member (DIC at SSS = 0) for the Congo from 355 ± <sup>48</sup> <sup>μ</sup>mol kg−<sup>1</sup> to 221 ± <sup>97</sup> <sup>μ</sup>mol kg−1. As no sample was taken during the highest Congo discharge, our lowest salinity is above 27. Nevertheless, the end-member for the Congo River of 231.7 ± 62.2 <sup>μ</sup>mol kg−<sup>1</sup> is in good agreement with the average riverine DIC concentrations of 258 ± <sup>29</sup> <sup>μ</sup>mol kg−<sup>1</sup> measured at Brazzaville–Kinshasa by [45].

**Table 3.** Mixing equations, including Equation (4), between the Congo River and oceanic waters.


TA values are also lower near the Congo mouth and are better estimated by the following equation than by Equation (2):

$$\text{TA} = 61.89 \ (\pm 0.70) \* \text{SSS} + 137.51 \ (\pm 22.08) \tag{5}$$

which gives a TA end-member of 137.51 ± 22.08 <sup>μ</sup>mol kg<sup>−</sup>1, in agreement with the riverine TA values ranging from 85 to 235 μmol kg−1, the variations being strongly dependent on the Congo discharge unlike the DIC end-member that is relatively constant [45].

The chlorophyll images show the high biological activity near the Congo mouth, which decreases DIC concentrations and leads to carbon uptake as previously proved by the estimates of the CO2 flux in this region [13,14]. The mouth of the Congo is located at 6◦ S, 12.3◦ E and its plume can be traced 400 to 1000 km from the mouth. The influence of the Congo plume is evidenced at 6◦ S, 8◦ E by the large salinity (Figure 6a) and chlorophyll variations (Figure 7), as well as the spatial distributions of chlorophyll (Figure 2c,d) and salinity (Figure 3c,d). In order to relate the DIC distribution (Figure 4) to the chlorophyll climatology, the chlorophyll fields were re-gridded on a 0.25-degree map and the correlation between DIC and chlorophyll concentrations was calculated at each grid point from January to December (Figure 10).

The blue regions correspond to negative correlations and coincide with high chlorophyll concentrations and low DIC, which are caused by carbon uptake. The positive correlations (in red) correspond to regions where chlorophyll and DIC follow the same seasonal cycle. In the equatorial upwelling region, high DIC and high chlorophyll concentrations are observed during the upwelling season (July–September) and lower concentrations occur during the rest of the year, which explains the positive correlation between DIC and Chla. Outside the upwelling region, low chlorophyll concentrations are observed but follow the seasonal cycle and the DIC–Chla correlation is also positive. Near river outflows, high chlorophyll concentrations are usually associated with low SSS and low DIC concentrations leading to negative DIC–Chla correlations. The spatial extension of

the negative correlations indicates the areas of influence of the river discharge. As SSS is a strong driver of DIC variations in the ETA, the negative DIC–Chla correlations correspond to areas with low SSS and high chlorophyll concentrations.

**Figure 10.** Map of correlations between chlorophyll and DIC concentrations calculated over 12 climatological months.

At 6◦ S, 8◦ E, low DIC and TA are observed especially in March–April during the period of maximum offshore extension of the Congo plume [46]. Nevertheless, the DIC observations at this site are better reproduced by the MLR, which suggests rapid mixing of the plume with oceanic waters during its westward propagation.

Observations made at 6◦ S, 8◦ E are compared with those at 6◦ S, 10◦ W. Unfortunately, measurements of fCO2 are available at both sites for 6 months only, from March to August 2017. The CO2 flux is much lower at 6◦ S, 8◦ E than at 6◦ S, 10◦ W (Table 4). However, the difference of fCO2 between the ocean and the atmosphere, ΔfCO2, is similar at both sites. DIC and SSS, and hence TA, are significantly lower at 6◦ S, 8◦ E. An increase in DIC would increase fCO2 but the increase in TA would decrease fCO2. The increase in DIC from 6◦ S, 8◦ E to 6◦ S, 10◦ W is compensated by the increase in alkalinity, so that overall, the mean fCO2 (and ΔfCO2) remains similar to the mean fCO2 at 6◦ S, 8◦ E.

**Table 4.** Mean and standard deviation of the carbon parameters, salinity and temperature at 6◦ S, 8◦ E and 6◦ S, 10◦ W from March to August 2017 from observations at the PIRATA moorings. The winds are from the Cross-Calibrated Multiplatform (CCMP).


The flux is significantly lower at 6◦ S, 8◦ E because the wind speed is much weaker close to the coast than offshore. However, the CO2 flux remains positive with only slight negative values, or close to zero, that may occur in March–April in some years (Figure 10). The flux and the ΔfCO2 are strongly correlated at this site (*r* = 0.95). During this time, the offshore spatial extension of the Congo plume is maximum [46] and the SSS at 6◦ S, 8◦ E is low (Figure 6a). The surface waters coming at the mooring are from the Congo plume and have low CO2 content because of the biological activity near the Congo mouth and the low CO2 content of the Congo waters. Although DIC–TA is the worst pair to calculate fCO2, the flux has been reconstructed at 6◦ S, 8◦ E to examine its annual variations (Figure 11, MLR sat). Some differences can reach more than 2 mmol m−<sup>2</sup> d−<sup>1</sup> as in June 2017 but the reconstruction suggests that the site is a source of CO2 with some possible CO2 uptake occurring during the first months of the year.

Δ

**Figure 11.** CO2 flux from daily fCO2 observations at 6◦ S, 8◦ E (in cyan), averaged monthly (in blue) and using the MLR with satellite SST and SSS (in red).

### *4.2. Year-to-Year Variability*

Using the MLR and the SST and SSS satellite fields from 2010 to 2019 at regional scale, the anomalies of DIC are calculated and compared with the anomalies of the physical parameters (Figure 12). As salinity is a strong predictor of DIC, the monthly variations of DIC anomalies from 2010 to 2019 are correlated with SSSA (*r* = 0.87). The anomalies of DIC are negatively correlated with SSTA and the correlation is much weaker (*r* = −0.38). The anomalies of SST in the ETA are related to the TSA index that cover a much larger area of the South Atlantic (30◦ W–10◦ E, 20◦ S–0◦), with a correlation coefficient between SSTA and TSA of 0.77.

In the tropical region, the year 2010 is characterized by much warmer temperatures whereas 2012 is a cold year [50]. The positive SSTA in 2010 are observed in the northern Tropical Atlantic and affect the air–sea CO2 flux [51]. In the ETA, SSTA are close to zero suggesting that the anomalous conditions of 2010 are limited to the Northern Tropical Atlantic. However, the ETA presents negative SSS anomalies throughout the year 2010 that result in negative DIC anomalies (Figure 12a). There was no anomaly of precipitation in the ETA in 2010 but, during the first months of 2010, the position of the ITCZ in the tropical Atlantic remained north of the equator instead of migrating south of it. As the ITCZ is associated with high precipitations, lower-than-usual salinities advected by the NECC might explain the negative SSSA in the ETA. When comparing the SSS anomalies in 2010 north and south of the equator, the SSSA are significantly lower north of the equator (mean of −0.36) than south of it (mean of −0.16, *p*-value = 0.0024). Nevertheless, the impact of the SSSA anomalies on DIC is moderate as it does not exceed 40 μmol kg<sup>−</sup>1, which is much less than the seasonal DIC variations in this region, as shown, for example, in Figure 9.

**Figure 12.** (**a**) Monthly DIC and SSS anomalies (DICA and SSSA), (**b**) SST anomalies (SSTA) and Tropical Southern Atlantic (TSA) index from January 2010 to September 2019.

The other strong anomaly in the tropical Atlantic is the cold event of 2012 mainly observed in the northeastern Atlantic and the southern hemisphere [50]. From November 2011 to early 2012, a strong cooling anomaly (SSTA < 0) occurs in the region and is associated with a negative TSA index (Figure 12b). The distribution of fCO2 at 6◦ S, 10◦ W was affected by the cooling with a decrease in fCO2 associated to a decrease in SST during January– March 2012 compared to other years [52]. However, the SST variations have a moderate effect on DIC and this cooling anomaly does not affect much the DIC concentrations (Figure 12a).

According to [2], interannual fCO2 is mainly driven by SST in the equatorial Atlantic extending from 50◦ W to 5◦ E and by DIC in the equatorial Pacific. The fCO2 monitoring is still sparse in the tropical Atlantic to conclude. Nevertheless, the interannual events of 2010 and 2012 affected fCO2 and suggest an important role of SST anomalies whereas in the ETA, DIC is mainly driven by SSS variations. The use of DIC and TA to calculate fCO2 does not reproduce well the fCO2 variations observed at 6◦ S, 10◦ W, which might be explained by the strong co-variation of fCO2 and SST. Using pH that has a strong dependency with SST with either DIC or TA would better reproduce observed fCO2. Monitoring fCO2 and calculating TA, with the robust TA–SSS relationship, allow the calculation of pH. This would be a substitute for direct pH measurements. Given the strong seasonality of carbon parameters in this region, long term monitoring is necessary to be able to quantify the acidification rate in this region.

#### **5. Conclusions**

The DIC distribution has been examined in the ETA (12◦ W–12◦ E, 12◦ S–12◦ N) as well as the factors controlling its distribution. DIC strongly depends on the water masses and increases over time due to the atmospheric CO2 increase. In the northern part of the basin, relatively fresh and warm waters are associated with lower DIC concentrations compared to the southern part where colder and saltier waters are enriched in CO2. A multiple linear regression between DIC and SST, SSS and year has been determined with SSS as the main predictor of DIC. Other regression techniques, such as decision tree, random forest and feed forward neural networks, do not significantly improve the DIC estimation except for low DIC concentrations corresponding to coastal data. However, a neural network would be required when more coastal data become available as non-linearity occurs at low values.

East of 10◦ E, the strong influence of the Congo plume is evidenced with a different relationship which corresponds to conservative mixing between the river and oceanic waters. The end-members of 231.7 ± 62.2 <sup>μ</sup>mol kg−<sup>1</sup> for DIC and 137.5 ± 22.1 <sup>μ</sup>mol kg−<sup>1</sup> for TA are in good agreement with published riverine measurements. The Congo plume reaches its more offshore spatial extension during December–April and affects the distribution of DIC at the 6◦ S, 8◦ E mooring. Nevertheless, the conservative mixing is no longer valid as the plume rapidly mixes with oceanic waters. Consequently, the carbon uptake is restricted to the region close to the Congo mouth and the carbon supplied by the upwelling season is the dominant process in the ETA. The seasonal cycle dominates the variations of DIC in the region, whereas the year-to-year variability of DIC is less than 40 μmol kg−<sup>1</sup> over the 2010–2019 period. SSS is the main driver of DIC in this region and the cold event of 2012 does not lead to significant DIC anomalies.

Although no pH measurements are made, the observations of fCO2 and the calculation of TA from SSS allow the calculation of pH. Pursuing the monitoring of fCO2 at time-series stations will contribute to the quantification of the rate of ocean acidification in this region.

**Author Contributions:** Conceptualization, N.L.; methodology and validation, N.L., C.M. and D.K.; writing—original draft preparation, N.L.; writing—review and editing, N.L., C.M., D.K., L.B. and U.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Institut de Recherche pour le Développement (IRD), the European integrated project AtlantOS (grant agreement 633211), and the French Integrated Carbon Observation System (ICOS-France).

**Data Availability Statement:** The data presented in this study are openly available and the references are given in the paper.

**Acknowledgments:** We are very grateful to the volunteers who collected the samples during the PI-RATA cruises and to the US IMAGO from the IRD, especially F. Roubaud, J. Grelet and P. Rousselot for the installation and deployments of the CO2 sensors. DIC and TA have been analyzed at the service national d'analyses des paramètres océaniques du carbone in LOCEAN. The L3\_DEBIAS\_LOCEAN\_v4 Sea Surface Salinity maps have been produced by LOCEAN/Institut Pierre Simon Laplace (UMR CNRS/UPMC/IRD/MNHN) laboratory and ACRI-st company that participate in the Ocean Salinity Expertise Center (CECOS) of Centre Aval de Traitement des Donnees SMOS (CATDS). This product is distributed by the Ocean Salinity Expertise Center (CECOS) of the CNES-IFREMER Centre Aval de Traitement des Donnees SMOS (CATDS), at IFREMER, Plouzane (France). We thank two anonymous reviewers for their constructive comments.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Regression Methods and Results**

Four data-driven regression methods (Multiple Linear Regression, Decision Tree, Random Forest and Feed Forward Neural Networks) have been applied to estimate DIC in the Eastern Tropical Atlantic. Figure A1 presents the dataset with the relationships between DIC and the predictors (SST, SSS and Year) and the distribution divided by DIC quantiles (in color).

Decision Trees (DT) are a non-parametric supervised learning method used for classification and regression. The goal is to create a model that predicts the value of a target variable by learning simple decision rules inferred from the data features. During the tree construction, the dataset is iteratively split by the parameters maximizing the variance of the predicted variable.

Random Forest (RF) algorithms are ensemble methods using so called perturb-andcombine techniques [53] designed for decision trees. This means a diverse set of regressors is created by introducing randomness in the model construction. The prediction of the ensemble is given as the averaged prediction of the individual regressors. Each tree in the ensemble is built from a sample drawn with replacement (i.e., a bootstrap sample) from the training set. In addition, when splitting each node during the construction of a tree, the best split is found either from all input features or a random subset of a predefined size.

**Figure A1.** Scatter plots between variables and density plots for each variable (on the main diagonal). Data are colored according to the intervals of the quantiles of DIC (25%, 50%, 75%).

Feed forward neutral networks (NN) are part of artificial neural networks. NN are non-linear functions able to uniformly approximate a wide class of functions encountered in physics [54], which include continuous and derivable functions with respect to each variable. The minimization used allows a global optimization of the set of parameters. Feed forward neural networks are a broadly applied technique in classification and regression tasks. Table A1 summarizes the characteristics of the feed forward Neural model. We use keras API [55] in Python using TensorFlow [56] as back-end.

**Table A1.** Characteristics of the neutral networks (NN).


The relationship between the DIC and its potential predictors, such as the sea surface salinity (SSS), sea surface temperature (SST), the interannual trend (expressed by the year variable) have been explored using four regression models and four independent observation datasets.

For the NN, the intra-annual variations were taken into account by adding two variables: sinDoY and cosDoY, representing the sine and cosine of the day of the year, respectively, normalized between 0 and 2π. Having a validation dataset is fundamental to avoid overfitting on trained data. From the data available for training, i.e., the DIC (2005–2019), two datasets were created for training and validation processes. We applied the following method to get a homogeneous representation of data in the validation dataset: the 637 measurements of DIC (2005–2019) were separated into 4 groups according to the 25%, 50%, 75% quantiles of DIC (this corresponds to the 1941.1 μmol kg−1, 1992.2 μmol kg−<sup>1</sup> and 2044.6 μmol kg−<sup>1</sup> DIC values). Finally, data in each group were sorted by SST and a sample was chosen for validation, leaving the rest for training. Best results were obtained by selecting one data out of six. The sizes of these datasets are given in Table A2.


**Table A2.** Data distribution for the NN showing the total number of data (N), the number of data for training (NTrain) and for validation (NVal) according to the 25%, 50%, 75% quantiles (Q25, Q50, Q75) of the DIC. The corresponding values of DIC are indicated in μmol kg<sup>−</sup>1.

The NN is then compared with the MLR using the same groups of data (Table A3). The two methods give similar RMSE but the NN tends to perform better than the MLR for low DIC values.

**Table A3.** RMSE for each DIC interval for both MLR and NN methods tested on the two time-series stations and the two cruises.


The different methods, MLR, NN, DT and RF are compared without separating the data into groups (Table A4). The NN usually performs slightly better than the other methods but, as it does not lead to significant improvement, we use the MLR, which is the simplest method. This also means that the dependencies of the DIC on its main predictors, SSS, SST, and year, are mostly linear. Among these features, the SSS has the largest impact. For random forests and decision trees, the relative importance of SSS is more than 90% (up to 95% for 6◦ S, 10◦ W), while that of the SST is about 5% (calculated as total variance reduction brought by the corresponding variable).

Although not very different, the DT is slightly less efficient than the other methods. The DT tend to group the predicted values into clusters, due to relatively small tree depths and consequently small numbers of leaves (i.e., discrete cases with observations attributed to them). DTs typically exhibit high variance and tend to lead to overfitting. RFs decrease the variance of the ensemble estimator. By taking an average of tree predictions, with injected randomness, some errors can cancel out. In our examples RF yield slightly better results that DT for all stations, in most cases being close to the MLR model. Nevertheless, all methods used here give similar results.


**Table A4.** Same as Table 2, but including a comparison between the multiple linear regression (MLR), the decision tree (DT), the random forest (RF) and the feed forward neural networks (NN) methods.

#### **References**


## *Article* **Millennial-Scale Environmental Variability in Late Quaternary Deep-Sea Sediments from the Demerara Rise, NE Coast of South America**

**Steve Lund \* and Ellen Platzman**

Department of Earth Sciences, University of Southern California, Los Angeles, CA 90007, USA; platzman@usc.edu **\*** Correspondence: slund@usc.edu

**Abstract:** We carried out a rock magnetic study of two deep-sea gravity cores from the Demerara Rise, NE South America. Our previous studies provided radiocarbon and paleomagnetic chronologies for these cores. This study presents detailed rock magnetic measurements on these cores in order to characterize the rock magnetic mineralogy and grain size as indicators of the overall clastic fraction. We measured the magnetic susceptibility, anhysteretic remanence, and isothermal remanence and demagnetized the remanences at several alternating field demagnetization levels. The magnetic intensities estimate the magnetic material concentration (and indirectly the overall clastic fraction) in the cores. Ratios of rock magnetic parameters indicate the relative grain size of the magnetic material (and indirectly the overall clastic grain size). Rock magnetic intensity parameters and rock magnetic ratios both vary systematically and synchronously over the last 30,000 years in both cores. There is a multi-millennial-scale cyclicity, with intervals of high magnetic intensity (high magnetic and clastic content) with low magnetic ratios (coarser magnetic and clastic grain size), alternating in sequence with intervals of low magnetic intensity with high magnetic ratios (finer grain size). There is also a higher-frequency millennial-scale variability in intensity superposed on the multi-millennial-scale variability. There are nine (A–I) multi-millennial-scale intervals in the cores. Intervals A, C, E, G, and I have high magnetic and clastic content with coarser overall magnetic and clastic grain size and are likely intervals of enhanced rainfall and runoff from the NE South American margin to the coastal ocean. In contrast, intervals B, D, F, and H have lower clastic flux with finer overall grain size, probably indicating lower continental rainfall and runoff. During the Holocene, high rainfall and runoff intervals can be related to cooler times and low rainfall and runoff to warmer times. The opposite pattern existed during the Pleistocene, with higher rainfall and runoff during interstadial conditions and lower rainfall and runoff during stadial conditions. We noted a similar pattern of Pleistocene multi-millennial-scale variability in a transect of deep-sea sediment cores along the NE Brazilian margin, from the Cariaco Basin (~10 N) to the NE Brazilian margin (~1◦ N–4◦ S). However the NW part of this transect (Cariaco Basin, Demerara Rise, Amazon Fan) has an out-of-phase relationship with the SE part of the transect (NE Brazilian margin) between warm–cold and wet–dry conditions. One possible cause of the high–low rainfall and runoff patterns might be oscillation of the Intertropical Convergence Zone (ITCZ), with higher rainfall and runoff associated with a more southerly average position of the ITCZ and lower rainfall and runoff associated with a more northerly average position of the ITCZ.

**Keywords:** Demerara Rise; millennial-scale variability; Late Quaternary paleoclimate

## **1. Introduction**

The circulation pattern within the western Equatorial Atlantic Ocean and the large sediment supply from the Amazon River create a complex sedimentation system off the coast of NE South America. The Amazon continental shelf accumulates about five hundred million tons of sediment annually (Kuehl et al., 1986) [1], carried by about 20% of the world's

Millennial-Scale Environmental Variability in Late Quaternary Deep-Sea Sediments from the Demerara Rise, NE Coast of South America. *Oceans* **2021**, *2*, 246–265. https://doi.org/10.3390/ oceans2010015

**Citation:** Lund, S.; Platzman, E.

Academic Editor: Michael W. Lomas

Received: 30 December 2020 Accepted: 22 February 2021 Published: 5 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

freshwater entering into the ocean from the Amazon River (Franzinelli and Potter, 1983) [2]. Over 80% of the sediment discharged by the Amazon River plume originates from the Andes (Milliman and Meade, 1993; Meade, 1994) [3,4]. This large flux of sediment creates the Amazon Fan, the third largest modern deep-sea fan, an expanded continental shelf found off the coast of NE South America, as well as a distal sediment blanket on the adjacent Demerara Rise and Abyssal Plain. These sediments and their stratigraphic variability provide a unique record of Equatorial Atlantic ocean circulation and Amazon River runoff and their climate and environmental variability (e.g., Damuth and Fairidge, 1970; Damuth and Kumar, 1975; Damuth and Flood, 1985; Flood et al., 1995; Maslin et al., 2000; Maslin and Burns, 2000) [5–10].

Cruise KNR197-3 (chief scientists D. Oppo and W. Curry) collected multicores, gravity cores, and long piston cores along a transect of the Demerara Rise perpendicular to the local coastline (Figure 1). The cores were collected in varying water depths from ~500 m to 3200 m. The sediment cores were taken to better understand the Late Quaternary sequestration of carbon and related physical sedimentation. This paper examines the Late Quaternary pattern of physical sedimentation in the uppermost portion of the Demerara Rise (core depths of less than 1100 m) based on rock magnetic measurements. This part of the Demerara Rise transect has the highest sediment accumulation rates throughout the Late Quaternary. We have previously published a paleomagnetic secular variation record for some of these cores (Lund et al., 2017) [11], along with radiocarbon dating (Huang et al., 2014) [12], which provide a chronostratigraphic framework for our study. We think this study of Late Quaternary physical sediment variability on the Demerara Rise can be used to estimate the environmental pattern of Amazon River plume sedimentation into the deep ocean adjacent to NE South America for the last ~30 ka.

**Figure 1.** *Cont*.

**Figure 1.** (**A**) Map of the northeastern continental margin of South America showing the Demerara Rise and general flow of the surface North Brazil Current. (**B**) Demerara Rise bathymetry and locations of our studied cores (stars).

## **2. Regional Ocean Circulation**

The ocean circulation pattern in the western Equatorial Atlantic Ocean is dominated by the Atlantic Meridional Overturning Circulation (AMOC) (Berger and Wefer, 1996; Schmitz and McCartney, 1993) [13,14]. This system brings warm water into the North Atlantic Ocean, where it is transformed into North Atlantic Deep Water (NADW), which then migrates to the South Atlantic Ocean (Metcalf and Stalcup, 1967; Kirchner et al., 2009) [15,16]. The principal surface ocean current along the northeastern margin of South America is the North Brazil Current (NBC) (Figure 1). The NBC forms off the easternmost margin of Brazil, where the westward-flowing South Equatorial Current (SEC) bifurcates to form the northwardflowing NBC and the southward-flowing Brazil Current (Stramma et al., 1995) [17]. The NBC crosses the equator as a warm-water current, where it connects to the North Equatorial Current (NEC), Caribbean Current, and Gulf Stream (Schmitz and McCartney, 1993) [14]. The NBC plays an important role in the cross-equatorial transfer of heat to the North Atlantic Ocean (Metcalf and Stalcup, 1967) [15]. The NBC also distributes the sediment plume associated with Amazon River discharge northward along the continental shelf and Demerara Rise–Abyssal Plain (Ruhlemann et al., 2001) [18].

Three deep-ocean currents flow along the NE margin of South America below the NBC (e.g., Huang et al., 2014) [12]. Antarctic Intermediate Water (AAIW) flows northward at depths of ~500–1500 m, directly below the NBC, and is an important source of nutrients and carbon to the North Atlantic Ocean. NADW flows southward at depths of ~2000–4000 m below AAIW. The deepest parts of the Demerara Abyssal Plain may see northward flow of Antarctic Bottom Water (AABW) at depths below 4000 m.

#### **3. Regional Deep-Sea Sedimentation**

All of the currents carry varying amounts of suspended sediment, depending partly on flow velocities and partly on advection of sediment at the surface by aeolian processes or coastal runoff. The most significant source of advected sediment is the Amazon River discharge. This large flux of sediment creates the Amazon Fan, the third largest modern deep-sea fan, and an expanded continental shelf. The Amazon delta and adjacent shelf currently trap perhaps two-thirds of the discharged sediment (Figueiredo and Nittrouer, 1995; Nittrouer et al., 1995) [19,20], however the finest sediment fraction is probably dispersed beyond the shelf by the NBC. This flux of Amazon River water and sediment beyond the shelf environment brings vast amounts of nutrients to the surface waters of the western Equatorial Atlantic Ocean, causing large blooms of phytoplankton (Smith and DeMaster, 1996) [21].

The western Equatorial Atlantic has two major modes of sedimentation caused by sealevel changes. Higher sea levels associated with interglacial periods foster the deposition of most sediment from the Amazon River onto the north Brazil continental shelf (Nittrouer and DeMaster, 1986; Nittrouer et al., 1995) [20,22]. A lesser portion of the sediment is transported as a suspended load either by longshore currents toward the northwest, where it is deposited along the continental shelf as far as French Guinea, or by the NBC, where it is deposited onto the Amazon Fan and adjacent Demerara Rise–Abyssal Plain (Lund et al., 2018) [23]. Lower sea levels associated with glacial periods cause sediment to bypass the exposed continental shelf and be deposited directly into the NBC and head of the Amazon Submarine Canyon (Damuth and Fairbridge, 1970) [5] and onto the fan.

The "on–off" supply of sediment to the Amazon Fan results in distinctive sediments for glacial versus interglacial intervals. Glacial sediments on the abyssal plain of the most recent glacial period consist of mottled and dark grey hemipelagic clays, a low concentration of carbonate and foraminifera, and occasional redeposited terrigenous sand and silt layers. Above the glacial sediments, Holocene pelagic sediments consist of tan foraminiferal marls or oozes with diverse foraminifera assemblages (Damuth, 1977; Damuth and Flood, 1985) [7,24]. Holocene and interglacial sediments in the western Equatorial Atlantic are characterized by lower concentrations of magnetic material (mostly magnetite) and by finer-grained magnetic material (Meinert and Bloemendal, 1989) [25]. In contrast, the glacial sediments have higher concentrations of magnetic material (mostly magnetite) that is coarser-grained.

Climate also plays a major role in the amount and distribution of sediment deposited on the Amazon Fan, adjacent continental shelf, and more distal Demerara Rise–Abyssal Plain at shorter timescales. The seasonal migration of the Intertropical Convergence Zone (ITCZ) creates distinctive wet and dry seasons and has a large influence on environmental patterns of the coastal region (e.g., Wainer and Soares, 1997) [26]. Boreal summer placement of the ITCZ near 10◦ N causes dryness in the Amazon Basin and strong precipitation further to the north, including the Cariaco Basin (Haug et al., 2001) [27]. During the boreal winter, the ITCZ moves over the Amazon Basin (0◦–10◦ S), creating a pronounced wet season (austral summer monsoon from December to February).

There is evidence at millennial timescales that the average annual placement of the ITCZ may be biased in its position. During the Early–Middle Holocene (approximately 8000–4000 years BP), there is widespread evidence that the climatic conditions in the tropical Andes were significantly drier than present (e.g., Mayle et al., 2004; Mayle and Power, 2008) [28,29]. More southerly areas experienced similar increased aridity in the Early Holocene. For example, the driest climatic conditions at Lake Junin (11◦ S), Lake Titicaca (14–17◦ S), and Sajama Mountain (18◦ S) were centered at ca. 10,000, 5500, and 4000 years BP, respectively [30–33] (Thompson et al. 1998; Seltzer et al., 2000; Abbott et al., 2003; Bush et al., 2005). In contrast, more northerly sites in Venezuela and Columbia (including Cariaco Basin) showed significant increases in rainfall in the Early Holocene, as the ITCZ migrated even more northerly during summer months (Haug et al., 2001) [27]. Overall, the sediments of the northern South American margin and their stratigraphic variability provide a correlative record of Amazon Basin climate and environmental variability (e.g., Maslin et al., 2000, Lund et al., 2018) [9,23].

#### **4. Core Descriptions**

Cruise KNR197-3 collected a sequence of cores at ~500 m to ~3000 m water depth in a transect outward from the South American margin (Figure 1) as part of a paleoceanographic study of Atlantic Meridional Overturning Circulation (AMOC) in the Atlantic Ocean. Typically, multicores (MC), gravity cores (GGC), and long piston cores (CDH) were collected to provide a complete, composite sedimentologic record along the transect. The focus of our rock magnetic study is 9GGC (405 cm) and 25GGC (352 cm) gravity cores with companion 10CMC (31 cm) and 25DMC (32 cm) multicores (Figure 1, Table 1). We also make use of a 46CDH piston core because it has a large number of radiocarbon dates and XRF measurements that can be transferred into the gravity cores based on magnetic susceptibility (CHI) profile correlations. The radiocarbon dates (Huang et al., 2014) [12] indicate that our paleomagnetic or rock magnetic study covers the last ~19 ka in core 25GGC and ~31 ka in core 9GGC.

**Table 1.** Site locations and descriptions of Demerara Rise cores discussed in this paper.


#### **5. Rock Magnetic Measurements**

U-channel samples were collected from both gravity cores for our rock magnetic studies. The multicores were sampled discretely with contiguous 2 × 2 × 2 cm cubes because the cores were too short (<50 cm) for good-quality u-channel measurements. Paleomagnetic measurements of these cores have been previously described by Lund et al. (2017) [11]. All cores were subjected to whole-core magnetic susceptibility measurements while onboard the ship. The u-channeled sediments were given a series of rock magnetic measurements after completion of the paleomagnetic studies. All samples were given an artificial anhysteretic remanent magnetization (ARM) sequence with a bias field of 0.05 milliTesla (mT) and an alternating field of 100 mT. After sample measurement, all samples were subjected to a sequence of ARM af demagnetization and measurements at 10 mT, 20 mT, 40 mT, and 60 mT. Finally, all samples were given another artificial isothermal remanent magnetization (IRM) sequence with a pulsed 1 T field. (This 1-T artificial remanence is usually referred to as a saturation IRM (SIRM), because all magnetite and titanomagnetite grains, which hold almost all of the natural remanence (NRM), are saturated.) After sample measurement, all samples were subjected to a sequence of SIRM af demagnetization and measurements at 10 mT, 20 mT, 40 mT, and 60 mT. One u-channel from each of the gravity cores also had their ARMs and SIRMs subjected to 80 mT and 100 mT af demagnetization. Lund et al. (2017, 2018) [11,23] summarized the AF coercivity patterns of the NRM, ARM, and SIRM. Those papers also considered variations in the ARM/Chi ratio, which showed no evidence of fine-grained (bacterial) magnetite. These studies both suggest that detrital magnetite is the dominant magnetic mineral in these sediments. The resulting CHI, ARM, and SIRM intensities for both gravity cores are shown in Figures 2 and 3. The magnetic intensity patterns in the two multicores and neighboring multicores (Lund et al., 2018) [2] are consistent with the pattern in the uppermost 50 cm of the gravity cores. This provides corroborating evidence that the gravity cores did not lose significant surface sediment in the coring process.

The magnetic intensity variations reflect changes in the percentage of detrital magnetic material (and presumably total clastic material) in the sediment. Magnetite is the most common magnetic mineral noted in several other studies of deep-sea sediments from the Equatorial Atlantic–Amazon Fan region (Meinert and Bloemendal, 1989; Maslin et al., 2000; Franke et al., 2007; Bendle et al., 2010) [9,25,34,35]. Our previous rock magnetic studies of deep-sea sediments from the neighboring Demerara Abyssal Plain also support this view (Lund et al., 2018) [23]. These studies also indicate that ferric oxides (mostly hematite and goethite) are present in varying small amounts in the sediments. The magnetic intensity variations may be due to varying clastic flux (presumably mostly from the Amazon drainage) or varying dilution by biogenic components due to variable upwelling and biological productivity.

Our ARM and Chi measurements and ratios of ARM or SIRM under different levels of demagnetization (SIRM20/SIRM0; SIRM40/SIRM20, ARM20/ARM0, ARM40/ARM20) estimate the grain size distribution of the magnetic mineral fraction (e.g., Stacy and Banerjee, 1974; Bloemendal et al., 1988) [36,37]. For all rock magnetic ratios, higher (lower) values indicate a finer-grained (coarser-grained) set of magnetic grains. Variations in magnetic grain size may be due simply to a finer-grained or coarser-grained source of overall clastic material; they might also reflect changes in the distribution of magnetic minerals if more than one is significant in concentration.

**Figure 2.** Rock magnetic intensity variability of gravity core 9GGC. Individual rock magnetic intensities are shown as Chi, anhysteretic remanence (ARM), and saturation isothermal remanence (SIRM) values. See text for further discussion.

**Figure 3.** Rock magnetic intensity variability of gravity core 25GGC. Individual rock magnetic intensities are shown as Chi, anhysteretic remanence (ARM), and saturation isothermal remanence (SIRM) values. See text for further discussion.

Selected samples from gravity core 9GGC were also given an extra series of magnetic hysteresis measurements. The samples were cycled through a MicroMag up to 1 Tesla (T). Values of remanent magnetization (Mr), saturation magnetization (Ms), coercivity field (Hc), coercivity of remanence field (Hcr), and field intensity (H) were calculated to estimate the grain size distribution and percentage of ferric magnetic minerals (goethite, hematite) versus titanomagnetite minerals (primarily magnetite). H represents the amount of magnetic material that acquires a remanence between 0.3 and 1 T (ferric magnetic minerals, primarily hematite). The primary purpose of these measurements was to assess the percentage of ferric magnetic minerals that primarily develop in Amazon Basin soils and reflect the relative importance of erosion from the Andes (primarily magnetite and titanomagnetite) versus the Amazon Basin itself (a more significant portion of goethite and hematite) on total Amazon drainage to the coastal ocean.

#### **6. Development of a Revised Magnetic Chronostratigraphy Technique**

Lund et al. (2017) [11] developed chronostratigraphies for cores 9GGC and 25GGC by first correlating the two cores based on environmental magnetic and paleomagnetic variability and then using that correlation to assess the radiocarbon dates in cores 9GGC and 25GGC. (Selected radiocarbon dates from 46CDH were also tied into the two gravity cores using magnetic susceptibility correlations.) The final chronostratigraphies for 9GGC and 25GGC were then derived from a selection of radiocarbon dates tied to each core individually and connected with linear interpolation. We then identified 52 paleomagnetic and rock magnetic features that both cores have in common. We estimated the ages of each feature based on the chronologies outlined by Lund et al. (2017) [11]. We found age differences between the same correlatable features in both cores of as much as 1000 years because of differences in linear interpolation between radiocarbon dates in the two cores.

We developed a revised magnetic chronostratigraphy technique for these two cores that provides more consistent ages for environmental and paleomagnetic features. The starting point was to date the 52 environmental magnetic or paleomagnetic correlation features in the two cores (summarized in Supplementary Table S1) based on the original chronos-

tratigraphies developed by Lund et al. (2017) [11]. The 52 feature ages in 9GGC and 25GGC were then averaged. Those feature ages were then used to rebuild chronologies for both cores. The results are shown in Figure 4. Figure 4 shows the original Lund et al. (2017) [11] chronologies based on calibrated radiocarbon dates (large open circles in Figure 4), with linear interpolation and the revised chronology based on the 52 feature ages summarized in Supplementary Table S1. This procedure makes PSV and environmental variability more synchronous between the two cores. It does not, however, improve the overall chronology of the sediments. It simply makes the two records chronologically consistent. There are two underlying assumptions in this procedure that we think are warranted. First, we presume that the sediment remanence acquisition process is comparable in both cores. Second, we presume that the sedimentation process near our study cores is not time transgressive. Both assumptions seem reasonable based on the strong (±2 cm) correlation between both PSV and rock magnetic variability among the two cores over their entire lengths.

The revised chronostratigraphies were used to estimate sedimentation rate variability in both cores. The results are shown in Figure 5. There is a notable overall pulse of high sedimentation around the Pleistocene–Holocene transition (~10–16 ka). Superposed on this are multi-millennial variations in sedimentation rates that we can tie to our rock magnetic variability described below. In particular, note that colder intervals of the Late Pleistocene (Younger Dryas (YD) and Heinrich events H1 and H2) have lower sedimentation rates than the neighboring warmer intervals (Bolling–Allerod (B/A) and Dansgard–Oeschger (D/O) interstadial events 2–4).

**Figure 4.** Time–depth plots for cores 9GGC and 25GGC. The open circles are calibrated radiocarbon dates used by Lund et al. (2017) [11] with linear interpolation to dates the cores. We built a revised chronology, as summarized in Supplementary Table S1, to create the more detailed chronologies shown here (dark lines).

**Figure 5.** Sedimentation rates for GGC9 and GGC25. There is a multi-millennial-scale alternation in sedimentation rates (faster–slower) that is indicated by grey versus white zones labeled A–I. These correspond to labeled zones in Figures 6–9.

**Figure 6.** Rock magnetic intensity variability for core 9GGC versus time. Alternating light and dark grey zones are multi-millennial-scale oscillations in magnetic intensity (light = higher intensity, dark = lower intensity) labeled A-I. Vertical lines estimate finer-scale millennial intensity peaks that are superposed on the longer-term variability.

**Figure 7.** Rock magnetic intensity variability for core 25GGC versus time. Alternating light and dark grey zones are multi-millennial-scale oscillations in magnetic intensity (light = higher intensity; dark = lower intensity) labeled A–I. Vertical lines estimate finer-scale millennial intensity peaks that are superposed on the longer-term variability.

**Figure 8.** Rock magnetic coercivity variability for core 9GGC versus time. Alternating light and dark grey zones are multimillennial-scale oscillations in magnetic coercivity (light = lower coercivity or coarser-grained; dark = higher coercivity or finer-grained) labeled A–I. Vertical lines indicate finer-scale millennial intensity peaks noted in Figures 6 and 7. There is no significant evidence for coercivity variations associated with those intensity variations.

**Figure 9.** Rock magnetic coercivity variability for core 25GGC versus time. Alternating light and dark grey zones are multi-millennial-scale oscillations in magnetic coercivity (light = lower coercivity or coarser-grained; dark = higher coercivity or finer-grained) labeled A–I. Vertical lines indicate finer-scale millennial intensity peaks in Figures 6 and 7. There is no significant evidence for coercivity variations associated with those intensity variations.

## **7. Environmental Variability**

Figure 6 (9GGC) and Figure 7 (25GGC) display the pattern of magnetic intensity variation as a function of time over the last 30 ka. The two gravity cores have the same pattern of intensity variation, with six multi-millennial-scale intervals in common (A–F in Figures 6 and 7). Core 25GGC also displays three more intervals (G-I) that are older than the data in 9GGC. Each of these intervals lasts ~1500–4000 years. The intervals oscillate between high and low intensity. The intensity variations could be associated with variations in the rate of clastic flux to the coastal ocean or with variations in biological productivity causing biogenic carbonate or silica dilution.

Figure 8 (9GGC) and Figure 9 (25GGC) display the pattern of magnetic coercivity variations over the last 30 ka. The coercivity varies between high (finer-grained clastics) and low (coarser-grained clastics) values and is the same pattern in both cores. The timing of coercivity variations is also consistent with the pattern of intensity variations. High intensities generally go with lower coercivities. This relationship between intensity and coercivity strongly suggests that clastic flux variations in amount and grain size are the primary cause of the observed multi-millennial-scale variability. Higher rates of clastic flux probably indicate stronger outflow from the Amazon delta; this should also yield overall coarser clastic flux. This pattern is consistent with the sedimentation rate variations noted in Figure 5. High intensity–coarser clastic flux intervals that suggest stronger Amazon outflow are also intervals of higher than average sedimentation rate. Additionally, the lower intensity–finer clastic intervals that suggest lower Amazon outflow are also intervals of lower than average sedimentation rate.

Several other deep-sea sediment studies discussed below use Ti/Ca or Fe/Ca from XRF studies to assess clastic flux variability. Delia Oppo (unpublished) [38] supplied us with an XRF and magnetic susceptibility record for core 46CDH (Figure 1, Table 1). The results are plotted in Figure 10. We correlated our multi-millennial-scale variability (intervals A–I) to the Ti/Ca results using the magnetic susceptibility records. The high (low) Ti/Ca intervals (relatively high clastic flux) are consistent with our high-intensity–coarser (low intensity–finer) rock magnetic intervals.

Intervals A–I can be easily associated with known regional to global climate variations over the last 30 ka. Intervals A–C can be associated with Late Holocene (A) and Middle Holocene drought (B) in the Amazon Basin (Maslin and Burns, 2000; Lund et al., 2018) [10,23] and final deglaciation (C) with high Amazon outflow (Maslin et al., 2000; Maslin and Burns, 2000) [9,10]. These three intervals represent Holocene variations in clastic flux, most likely due to variations in multi-millennial-scale rainfall or glacial melting. The low intensity (and modestly higher coercivity, finer-grained) of intervals D, F, and H can be associated with cold intervals (stadials) of the Younger Dryas (YD) and Heinrich events 1 (F) and 2 (H). This observation is consistent with the results of Maslin et al. (2011) [39] for the YD and Heinrich event 1. These cold intervals can also be interpreted as relatively dry intervals. The intervening higher flux (coarser -rained) intervals are associated with the warmer (interstadial) Bolling–Allerod (D/O) cycle 1 (E), D/O cycle 2 (G), and D/O cycles 3–4 (I). These intervals probably represent less cold and wetter intervals within the Amazon Basin. These are all labeled on Figures 6–9. There is one complication to this analysis. It is likely that the colder Pleistocene intervals (YD, Heinrich events) are slower North Brazil Current intervals. This change in flow speed might have contributed to the variations in sedimentation rate. However, a number of published studies (discussed below) that show this same multi-millennial-scale pattern of sediment variability all associate the variation predominantly with variable rainfall or runoff from the Amazon Basin to the NE South American margin.

There is also a higher-frequency scale of variability in these records that is more subdued. This pattern is indicated by the solid lines in Figures 6–9. The lines indicate local intensity peaks. The relationship with coercivity is more complicated and subdued at best or perhaps there is no relationship at all. These cycles occur in all of the longer-duration intervals noted above. They suggest a higher-frequency oscillation in Amazon outflow, with the black lines indicating times of higher flow (higher clastic flux). These more subtle oscillations in intensity are consistent between cores 9GGC and 25GGC. However, the lack of any notable coercivity variability makes their interpretation more questionable. Table 2 presents the intensity peak ages and the durations between successive peaks. The average duration is ~1750 ± 450 years.

We finally assess the evidence for variations in ferric (primarily hematite)/ferrous (magnetite and titanomagnetite) ratios in the clastic sediments. These ratios indicate the relative importance of ferrous versus ferric components to the overall sediment magnetization. Hysteresis measurements of selected samples from 9GGC are shown in Figure 11. The rock magnetic data suggest that the magnetization is primarily magnetite–titanomagnetite in the size range of ~5–15 microns. This distribution could have varying percentages of coarser (MD) or finer (SD) grains and could still be within the pseudo-single domain (PSD) range for magnetite (Figure 11A). The hysteresis data let us estimate the amount of magnetization that is acquired with increasing magnetic fields; magnetic contributions between 0.3 and 1.0 T are usually referred to as hard components (H), which is associated with ferric minerals (primarily hematite). The ratio of total magnetic remanence to H (H/SIRM; Figure 11B) estimates the relative portion of ferric magnetization in any sediment sample. The time interval 6–9 ka has the lowest ratio of the last 25 ka. Bendle et al. (2010) [35] and Lund et al. (2018) [23] have both documented low H/SIRM ratios in the Early–Middle Holocene (6–9 ka) of deep-sea sediments from the Demerara Abyssal Plain; our results corroborate this. They all indicate that the Amazon Basin itself was relatively dry during the Early–Middle Holocene, with lower than average erosion of Amazon Basin soils (relative to contributions from the Andes).

**Figure 10.** Proxy data for core 46CDH (Figure 1, Table 1), situated directly between gravity cores 9GGC and 25GGC. Top—magnetic susceptibility record; the variability can be directly correlated to data from 9GGC (Figure 2) and 25GGC (Figure 3). Below are Ti and Ca counts from unpublished XRF data on 46CDH (Delia Oppo, personal communication). The rock magnetic intensity zonation is consistent with the XRF data. High chi goes with high Ti (a clastic indicator; low chi goes with high Ca (a biologic indicator).


**Table 2.** Timing of millennial-scale intensity peaks and durations between successive intensity peaks.

**Figure 11.** Magnetic hysteresis measurements of selected sediment horizons in core 25GGC. (**A**) Plot of remanent magnetization (Mr)/saturation magnetization (Ms) versus coercivity of remanence (Hcr)/coercivity (Hc). All samples fall within the pseudo-single domain (PSD) range. (**B**) Saturation isothermal remanence (SIRM) and field (H)/SIRM plotted versus time in 25GGC. See text for further discussion.

## **8. Comparison with Other Records**

The first evidence for millennial-scale variability in the tropical Atlantic region was noted by Grimm et al. (1993) [40] in Lake Tulane (27.6◦ N) in Florida and Curry and Oppo (1997) [41] in the Ceara Rise (5◦ N), also on the NE South American continental margin. The Lake Tulane study documented enhanced rainfall during Pleistocene North Atlantic cold (Heinrich) events. The Ceara Rise study linked coherent multi-millennial-scale tropical sea-surface temperature and high-latitude North Atlantic temperature variations. Locally, we can compare our deep-sea sediment records with others from the NE margin of South America.

Lund et al. (2018) [23] studied rock magnetic variability of Holocene deep-sea sediments from the Amazon Fan–Demerara Abyssal Plain. Our results from this study are the same. Our Pleistocene results (~7◦ N) can be compared with three other studies along the NE South American margin: the Cariaco Basin (~10◦ N), Demerara Rise–Amazon Fan (~4◦ N–7◦ N), and outboard of the NE Brazilian margin (~1◦ N–3◦ S). All of these studies show strong evidence of multi-millennial-scale variability from 12–30 ka, dominated by colder phases associated with the Younger Dryas and Heinrich events 1 and 2 versus warmer phases associated with the Bolling–Allerod and DO cycles 1, 2, and 3.

Our results are consistent with results from the Cariaco Basin (Haug et al. 2001; Deplazes et al., 2013) [27,42] showing lower sedimentation rates and clastic flux during the Younger Dryas, H1, and H2 associated with dryer conditions along NE South America, while the intervening intervals have higher clastic flux associated with heavier rainfall and runoff. These variations are associated with shifts in the long-term position of the Intertropical Convergence Zone (ITCZ).

Haggi et al. (2017) [43] studied the isotopic composition of plant waxes in a deep-sea sediment core from the same region (7◦ N) as our cores. They suggested that the Amazon rainforest was affected by intrusions of savannah or open framework vegetation types during H1 and H2 associated with drier conditions. The intervening intervals indicated higher Amazon rainfall. This pattern is completely consistent with our Pleistocene multimillennial-scale variability. They attributed these long-term variations to changes in atmospheric circulation.

Crivellari et al. (2018) [44] studied the H1 interval in more detail in the same deep-sea core as Haggi et al. (2017) [43]. They identified an early part of the H1 interval (HS1a: 17–18 ka) that had higher influx of sediment to the site, however the core of HS1 (HS1b: 15–17 ka) showed a decrease in terrestrial flux. Our results do not show evidence of this two-part structure, however our study and the study by Haggi et al. (2017) [43] do show a broadly drier response to H1.

A group of four deep-sea sediment studies further to the south (1◦ N–4◦ S) along the NE Brazilian margin also show the same pattern of multi-millennial-scale environmental variability but with mostly an opposite phasing: colder intervals had higher clastic flux and warmer intervals had lower clastic flux. However, there is some disagreement among these four studies that could be due to chronology issues. Arz et al. (1998) [45] used Ti/Ca XRF ratios to document multi-millennial-scale pulses of significant clastic sedimentation into the Atlantic slope along the NE corner of Brazil (3.5◦ S). One pulse correlates to H1, but two pulses correlate to warm Bolling–Allerod and DO interstadial 2 phases; the intervening Younger Dryas and H2 cold intervals had low clastic flux. Nace et al. (2014) [46] studied a marine core along the NE Brazilian margin at 0.3◦ N. They also identified multi-millennialscale pulses of high clastic flux associated with Ti/Ca XRF data, which they associated with increased precipitation or weathering on the adjacent continent. They correlated these events to multiple (cold) Heinrich events. However, their data actually showed high clastic pulses associated with H1 (~16 ka) and (warm) D/O cycles 2–3 (~30 ka) and low clastic flux at ~26 ka (cold H2). This pattern is consistent with the data from the study by Arz et al. (1998) [45]. Zhang et al. (2015) [47] studied Fe/Ca XRF data from a deep-sea core at 1.6◦ S and noted higher clastic flux during all three late Pleistocene cold intervals, namely the Younger Dryas, H1, and H2 events. Mulitza et al. (2017) [48] studied a deepsea sediment core (2◦ S) in the same region and concluded that both the Younger Dryas and H1 cold intervals were intervals of increased rainfall and higher clastic flux based on Fe/Ca XRF data. These results are opposite to the Younger Dryas and H2 results of Arz et al. (1998) [45] from the same region. Overall, these records do generally support an opposite phase of high clastic flux in cold intervals than we noted for the Demerara Rise.

Speleothem records from the Amazon Basin display similar complexity in Late Pleistocene rainfall variability. Cheng et al. (2013) [49] suggested that broadly the eastern Amazon region had lower rainfall during the late Pleistocene versus the western Amazon region. This might support our observation that NW marine records (Cariaco Basin, Demerara Rise) had an out-of-phase relationship for rainfall and runoff with SE marine records (NE Brazilian margin). Cheng et al (2013) [49] also noted that both regions had multi-millennial-scale variability on top of longer-duration orbital-scale variability with wetter intervals in NE Brazil (Rio Grande do Norte record) during the H1 and H2 stadials. This is consistent with the NE Brazilian margin marine records noted above. Wang et al. (2017) [50] also noted a similar pattern of multi-millenial-scale rainfall variability in their Paraiso speleothem record (~3◦ S) of NE Brazil.

#### **9. Discussion**

Cores 9GGC and 25GGC are dominated by millennial- to multi-millennial-scale variability throughout their entire extent. It is worthwhile noting that the timescales of variability do not significantly change between the Holocene and the Pleistocene. The dominant variability seems to follow a simple oscillatory multi-millennial pattern of wet versus dry conditions in the Amazon Basin (stronger and faster Amazon outflow with coarser clastic grains during wetter conditions versus lower and weaker Amazon outflow with finer clastic grains during drier conditions) in NE South America. Intervals A, C, E, G, and I are associated with relatively wet conditions (and strong Amazon flow), while intervals B, D, F, and H are associated with relatively dry conditions (and weaker Amazon flow). The sediments reach the Demerara Rise from the Amazon outflow under control of the North Brazil Current after passing the local shelf. Coarser (finer) clastic grains are more (or less) likely to reach the North Brazil current under conditions of higher (or lower) Amazon flow. The deep-sea sediment records south of the Amazon Fan (NE Brazilian margin, 1◦ N–4◦ S) suggest an out-of-phase relationship with our Demerara Rise results and suggest that NE Brazil to the east of the Amazon River delta receives most rainfall and runoff when the Demerara Rise and Cariaco Basin are drier. This also suggests that when we consider Amazon River outflow versus NE Brazil outflow, we are talking about western versus eastern Amazon Basin patterns as a whole.

Multi-millennial-scale intervals A–C are of the Holocene–Glacial transition in age. This interval seems to associate low rainfall and Amazon flow to warmer conditions (Middle Holocene warmth—B) and higher rainfall and Amazon flow to cooler conditions (Late Holocene—A; glacial transition—C). On the other hand, the Pleistocene intervals (e.g., Stuiver and Grootes, 2000 [51]; Cronin, 2009 [52]) D–I have an opposite relationship between Amazon flow and rainfall and temperature. Intervals D (Younger Dryas), F (H1), and H (H2) are all times of low rainfall and Amazon flow with colder conditions, while intervals E (Bolling–Allerod), G (D/O interstadial 2), and I (D/O interstadials 3–4) are times of strong rainfall and Amazon flow and warmer conditions. On this basis, we argue that multi-millennial-scale cyclicity in the Holocene is out of phase with multi-millennial-scale variability in the Pleistocene.

The higher-frequency magnetic intensity variability indicated by vertical black lines in Figures 6–9 suggests another level of wet–dry oscillation associated with Amazon Basin hydrology and Amazon outflow. However, the data do not show a consistent associated variation in coercivity, and so the interpretation of this more subtle higherfrequency variability is more questionable. Even so, there is a great deal of evidence for millennial-scale environmental variability (e.g., Bond et al., 1997 [53]), most dominantly with ~1500 year cyclicity. The average duration of our successive millennial-scale high intensity peaks is 1794 ± 482 years (Table 2). All of our cycles that are common to both 9GGC and 25GGC have an average of ~1600 years. The pre 20 ka interval in 9GGC alone shows three intervals with ~2600 year cyclicity (underlined in Table 2). As noted previously, our millennial-scale intensity cycles are not as notable as our longer-term cycles. If we are simply missing three cycles in the pre 20ka interval of 9GGC then it seems that our millennial-scale intensity variability is consistent with previous estimates of ~1500 year cyclicity.

Several of the published studies discussed above have commented on possible mechanisms for tropical environmental variability. One process that fits the overall data is multimillennial-sale variation in the placement of the Intertropical Convergence Zone (ITCZ). More southerly average placement will enhance NE South America and Amazon Basin rainfall. More northerly placement will lead to drier conditions. Haug et al. (2001) [27] suggested this mechanism caused a wet early Holocene wet interval in the Cariaco Basin (10.5◦ N) further to the northwest along the NE South American margin. This occurred at the same time as our interval B (early Holocene dry conditions), which was also noted in other nearby cores of the Demerara Rise by Lund et al. (2018) [23]. This pattern is also consistent with Lake Tulane wet intervals during Pleistocene Heinrich (cold) events versus our dry events in cold intervals D, F, and H.

The dominance of millennial- to multi-millennial-scale environmental variability in the tropical North Atlantic Ocean is also interesting to note given the recent study by Obrochta et al. (2012) [54] on variability of a similar scale in the high-latitude North Atlantic Ocean (ODP Site 609). They focused on the Late Quaternary ice-rafting evidence for a "1500-year cycle" in climate variability over the last glacial cycle. They suggested that much of that cyclicity might be better interpreted as an admixture of ~1000 and ~2000 year cycles (Obrochta et al., 2012) [54]. It seems clear that the entire Atlantic region has been subject to continuing millennial- to multi-millennial-scale variability over the Late Quaternary, however the cause(s) of that regional climate pattern and the relationship between tropical and high-latitude North Atlantic variability are still not certain.

### **10. Conclusions**

We have carried out a rock magnetic study of two deep-sea gravity cores (9GGC and 25GGC) from the Demerara Rise, NE South America. Previous studies (Huang et al., 2014 [12]; Lund et al., 2017 [11]) provided radiocarbon and paleomagnetic chronologies for these cores. In this study, we performed detailed rock magnetic measurements on these cores to characterize the rock magnetic mineralogy, magnetic concentration, and magnetic grain size as an indicator of the overall clastic fraction of the cores. We measured the magnetic susceptibility (chi), anhysteretic remanence (ARM), and isothermal remanence (SIRM) and demagnetized the remanences at several levels of af demagnetization. The magnetic intensities estimate the proportion of magnetic material (and indirectly overall clastic fraction) in the cores. The ratios of rock magnetic parameters (ARM/chi, ARM10/ARM0, SIRM10/SIRM0) indicate the relative grain size of the magnetic material (and indirectly the overall clastic grain size). Selected hysteresis measurements also gave us a sense of the variation in ferric (goethite, hematite)/ferrous (magnetite, titanomagnetite) ratios for the magnetic minerals.

The rock magnetic intensity parameters (chi, ARM, SIRM) and the rock magnetic ratios (ARM/chi, ARM10/ARM0, SIRM10/SIRM0) both vary systematically and synchronously over the last 30,000 years in both cores. There is clear evidence of a multi-millennialscale pattern of cyclicity, with intervals of high magnetic intensity (high magnetic and clastic content) and low magnetic ratios (coarser magnetic and clastic grain size) alternating in sequence. There is also evidence of finer millennial-scale variability in intensity superposed on the multi-millennial-scale variability, with an average cycle duration of ~1800 ± 500 years. There are nine (A–I) multi-millennial-scale intervals in the cores. Intervals A, C, E, G, and I have high magnetic and clastic content with coarser overall magnetic and clastic grain size and are likely intervals of enhanced rainfall and runoff from the NE South American margin to the coastal ocean. Alternatively, intervals B, D, F, and H represent periods of lower clastic flux with finer overall grain size, probably indicating lower rainfall and runoff from the continental margin. During the Holocene, high rainfall and runoff intervals can be related to cooler times and low rainfall and runoff to warmer times. The opposite pattern is true during the Pleistocene, with higher rainfall and runoff occurring during interstadial conditions and lower rainfall and runoff occurring during stadial conditions. The highest ferric/ferrous ratios occur during Holocene interval B (low rainfall and runoff) during the Early Holocene warm interval ~5–9 ka. We associated that interval with drought conditions in the Amazon region, consistent with previous estimates

of Amazon flow in this interval (Maslin and Burns, 2000 [10]; Bendle et al., 2010 [35]; Lund et al., 2018 [23]).

We noted a similar pattern of Pleistocene multi-millennial-scale variability in published marine deep-sea core studies from Cariaco Basin (~10◦ N) to the NE Brazilian margin (~1◦ N–4◦ S). These studies documented clastic flux variability using XRF measurements of Ti and Fe. However, it seems likely that the NW part of this transect (Cariaco Basin, Demerara Rise, Amazon fan) has out-of-phase warm–cold wet–dry relationships that support the notion of the eastern Amazon Basin and western Amazon being alternately wetter or drier. One possible cause of the high–low rainfall–runoff patterns might be oscillation of the Intertropical Convergence Zone (ITCZ), with higher rainfall and runoff associated with a more southerly average position of the ITCZ and lower rainfall–runoff associated with a more northerly average position of the ITCZ.

**Supplementary Materials:** Supplementary materials can be found at https://www.mdpi.com/2673 -1924/2/1/15/s1, Supplementary Table S1: Correlations for revised chronostratigraphy.

**Author Contributions:** Writing—S.L., Data acquisition—S.L. and E.P., data analysis—S.L. and E.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** NSF grant EAR1547605.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


## *Article*
