Application of a Lagrangian transport model to organo-mineral aggregates within the Nazaré canyon

Introduction Conclusions References


Introduction
Understanding the exchange of energy and matter between the shelf and the open ocean has been the focus of several European research programmes such as OMEX (Wollast and Chou, 2001), EUROSTRATAFORM (Weaver et al., 2006) and HERMES (Weaver and Gunn, 2009).Most recently the HERMES and HERMIONE programmes have addressed the distribution of organic matter, carbon flow and biodiversity in European continental margins (e.g.García et al., 2007;Ingels et al., 2009;Van Oevelen et al., 2011).In these studies submarine canyons are identified as important transport systems of sedimentary organic matter from the continental shelf to the deep ocean (Monaco et al., 1999;Schmidt et al., 2001;Canals et al., 2006), as important depocentres of sediments and organic matter of often higher quality (Epping et al., 2002;Van Weering et al., 2002;De Stigter et al., 2007;García et al., 2008García et al., , 2010) ) as well as hotspots of biodiversity (Ingels et al., 2009;Tyler et al., 2009;Cunha et al., 2011).Consequently, the transport of organic particles in submarine canyons is relevant in terms of global carbon budgets (Thomsen et al., 2002;Accornero et al., 2003;Masson et al., 2010).
Most of the present understanding on the transport of organic particles within submarine canyons has been derived from field observations which have subsequently been summarized in conceptual models of canyon dynamics.The downward transport and the redistribution of sediments and organic particles is controlled by hydrodynamic processes interacting with the bottom topography, such as internal tide circulation, internal waves, the formation of nepheloid layers, down and along slope bottom currents, intermittent gravity flows or the cascading of dense water (e.g.Van Weering et al., 2002;Canals et al., 2006;De Stigter et al., 2007).Hence, submarine canyons which are dominated by the formation of nepheloid layers and internal tides circulation, for example, should mostly concentrate organic material close to the canyon walls; while canyons dominated by down canyon circulation or cascading should mostly transfer organic particles to greater water depths.
The characteristics of marine aggregates collected from depth have been determined for the Iberian continental margin (Thomsen and Gust, 2000;de Jesus Mendes and Thomsen, 2007) but have not been used to date for numerical modelling.The application of lagrangian transport models linked to hydrodynamic models has a high potential to predict various environmental scenarios.At the western Iberian margin, lagrangian transport models have been applied to the Galician coast (Carracedo et al., 2006), Ria de Vigo (Huhn et al., 2012;Abascal et al., 2012), Rio Lima estuary (Vale and Dias, 2011), Ria de Aveiro lagoon (Dias et al., 2001), Óbidos lagoon (Malhadas et al., 2009) and Tagus estuary (Braunschweig et al., 2003).The operational model MOHID-PCOMS (MOdelac ¸ão HIDrodinâmica Portuguese Coast Operational Modelling System) (Mateus et al., 2012) runs in full operational mode for the western Iberian coast with daily hydrodynamic and ecological results.The model adequately represents the hydrodynamic features of the region and the seasonal variations in the dynamical processes.In this study, the MOHID model simulated the dispersion of OMAs within the Nazaré canyon by coupling the hydrodynamic model to a lagrangian transport model.The simulations in this current study cover several months of spring, a period during which surface phytodetritus production and subsequent transport to benthic communities is of ecological interest.The numerical model was assessed to determine whether it agreed with the flux passing through the upper and middle part of the canyon as described by the current conceptual model of organic matter transport.Lastly, the hypothesis that the Nazaré canyon acts as a conduit for OMAs, and therefore an enhanced carbon flux, was tested.

Study area
The western Iberian shelf and slope are intersected by several submarine canyons.The Nazaré canyon is the largest of these extending ∼ 210 km offshore, from the Nazaré beach running down to 5000 m depth (Tyler et al., 2009).According to Lastras et al. (2009), the canyon can be divided into three sections based on the hydrography and its physical characteristics.The upper section embraces a V-shaped valley incised into the shelf starting at the canyon's head and extending to a depth of 2700 m and is branched by a short side-valley called Vitória tributary.The middle section is characterized by a broad meandering valley with terrace slopes descending from 2700 m to 4000 m depth and the lower section, a flat floored valley which descends to a depth of 5000 m.The canyon cuts the entire Portuguese continental shelf and slope.The hydrodynamic processes are intensified by the rugged topography because the internal waves are preferentially formed in the canyon (Quaresma et al., 2007) and trapped as internal tidal energy.This mechanism is responsible for sediment resuspension and transport at the shelf (Quaresma et al., 2007) and in the upper section of the canyon (De Stigter et al., 2007).Martín et al. (2011) analysed the near bottom particle dynamics for the upper and middle Nazaré canyon and determined two contrasting dynamic environments.In the upper section (1600 m depth) high current speeds with spring tides up to 80 cm s −1 were registered and also high mass fluxes of particulate matter (mean 65 g m −2 d −1 ; maximum 265 g m −2 d −1 ), while at the deepest station (3300 m) the mass fluxes were below 10 g m −2 d −1 .The authors also concluded that storms can trigger sediment transport at the middle Nazaré canyon.

Organo-mineral aggregate data
The dispersion patterns, residence time estimation and travel trajectories of organic particles of different sizes under spring hydrodynamic conditions were studied.The OMAs of three different size classes (i.e.429 µm, 2000 µm and 4000 µm) were sampled during OMEX I, OMEX II, EU-ROSTRATAFORM and HERMES cruises to the northeastern Atlantic continental margin (RV Pelagia 95; RV Pelagia 1998;RV Meteor 1998/1999;RV Pelagia 2004) (Thomsen et al., 2002;de Jesus Mendes and Thomsen, 2007).The 429 µm aggregates belong to a dominant class of aggregates with the same median aggregate parameter size observed at the western Barents Sea, the northeast Greenland Sea, the Celtic Sea, and the Nazaré and Setúbal canyons (Thomsen and Graf, 1995;Thomsen and Ritzrau, 1996;Thomsen and Van Weering, 1998;de Jesus Mendes and Thomsen, 2007).Frequently these aggregate sizes were found at the shelf and at depths > 2500 m, while aggregates with larger dimensions (> 900 µm) were found at 3400 m depth at the northwest Iberian continental margin (Thomsen et al., 2002).The median aggregate sizes (429 µm) were constituted of organic matter (≤ 80 % wt) and lithogenic material (≥ 20 %) while the aggregates with larger dimensions (2000 µm, 4000 µm), also known as fluffy phytodetrital aggregates, were constituted of small amounts of lithogenic material and were highly transparent (> 80 % organic matter).
Critical shear velocities (U * cr ), critical deposition velocities (U * d ) and particle settling velocities (W s ) were determined for the three different aggregate sizes (Thomsen et al., 2002) (Table 1).These velocities were mandatory for the lagrangian model, and their units were converted into the model requirement units.

Hydrodynamic module
A high-resolution hydrodynamic model was used to simulate the evolution of the 3-D physical structure of the Iberian coast, and its influence on OMA transport to and within the Nazaré canyon.The model is an open source software under continuous development, named MOHID Water (http: //www.mohid.com),and a component of the MOHID Water Modelling System (MWMS), an integrated water modelling software that simulates water dynamics in water bodies, porous media and watersheds (Mateus, 2012).The MWMS is able to simulate broad processes and scales in marine systems ranging from coastal areas to the open ocean (Coelho et al., 2002;Santos et al., 2002;Santos et al., 2005;Mateus et al., 2012).The hydrodynamic model solves 3-D incompressible primitive equations considering hydrostatic equilibrium and the Boussinesq approximation, and its description can be found in Martins et al. (2001).The turbulent vertical mixing coefficient is determined using the General Ocean Turbulence Model (GOTM).

Lagrangian transport module
The lagrangian transport module of MOHID was used to simulate particle transport following the methodologies proposed in previous works (Braunschweig et al., 2003;Saraiva et al., 2007;Malhadas et al., 2009).The lagrangian module simulates the movement of aggregates located at specific water depths using the current fields calculated by the hydrodynamic module, thus solving the equation of transport independent of momentum balance equations.The lagrangian module derives the hydrodynamic information (current fields) from the system and updates the calculations without having the need to solve all variables at the same time.It uses the concept of passive tracers, characterized by their spatial coordinates, area and others properties such as settling velocities, and critical and depositional bottom shear stress (Table 1).(Thomsen et al., 2002;de Jesus Mendes and Thomsen, 2007).In our study, the model simulates the OMA trajectories using the concept of settling velocity, and each particle is assigned a time to perform random movement.These particles are placed at origins which emit the tracers at a specific depth and at one instant in time.The dispersion and distribution field of the particles is monitored using monitoring boxes to compute their residence time.For this project we use the term "residence time" for the temporal interval required by the OMAs to leave each monitor box.This is a new and alternative approach to the previous concept proposed by Braunschweig et al. (2003).The lagrangian results to characterize the OMA behaviour showed the average distance, displacement and velocity of OMAs of different sizes for each box .The distance was related to the total length that the OMAs travelled (km); the displacement was the difference between the initial and final position of the OMAs (km) and the velocity of the OMAs (km yr −1 ).

Model set-up for the Nazaré canyon
The domain configuration of the Nazaré canyon includes three levels of nested models using a one-way coupling (Fig. 1).This nesting methodology is described in Leitão et al. (2005).The first level covers the west coast of Iberia between 5.5-12.6 • W and 34.4-45.0• N with a resolution of 5.6 km.The boundary conditions of this level are provided by the 3-D operational model MOHID-PCOMS (Mateus et al., 2012).The operational model is forced by data from PSY2v2 Mercator Ocean solution for the North Atlantic and by MM5 atmospheric forecast model with 9 km resolution operated at IST (http://meteo.ist.utl.pt).Tide is imposed from 2-D barotropic model forced by the FES2004 global solution.
The second level covers the stretch from Figueira da Foz to Ericeira between 8.86-10.38• W and 39.02-40.08• N with a constant grid spacing of 2 km.The third grid has a resolution of 400 m for the Nazaré canyon area between 9-10.22 • W and 39.3-39.8• N. The vertical resolution of the three different levels adopted in this one-way nested modelling scheme is with 50 vertical layers, 43 Cartesian coordinates on the bottom and 7 sigma coordinates on the upper 10 m.The bathymetric data for the levels construction were provided by the National Oceanography Centre, The lagrangian module was run with tracers originating from 10 boxes of same dimensions distributed along the Nazaré canyon area at water depths between 59 and 3189 m (Fig. 2) (Table 2).Each box corresponds to a geographic domain of 3 × 3 cells of the model grid leading to a total of 1.44 km 2 .The boxes were filled with aggregates at a height of 0.5 m from the sea floor.While applying the module, properties such as the monitor box area and spatial coordinates, OMA settling velocities, and critical and depositional bottom shear stress were taken into consideration.For the OMAs of three different size classes, the settling velocities (W s ) increased with increasing aggregate size, while the critical shear velocities (U * cr ) decreased over the same aggregate size spectrum.However, the depositional bottom shear stress (T d ) is highest for the medium-sized aggregates (2000 µm), and has lower values for the 429 µm and 4000 µm aggregates (Table 1).The boxes were located in the upper (canyon head -2700 m) and middle (2700-4000 m) part of the Nazaré canyon according to Lastras et al. (2009).The first box was located at the canyon's head (59 m) and box 2 at the shelf    and 5), whereas the 429 µm OMAs continuously decreased with time inside the box (Fig. 3).The OMAs in boxes 3, 4, 5, and 6 had high residence times, indicating a reduced transport of aggregates in this part of the canyon.Box 7 at 1498 m showed a decrease in the residence times particularly for the 429 µm and 4000 µm (Figs. 3 and 5).The model predictions for boxes 8, 9, and 10 located offshore showed a very active transport for the OMAs of different size classes.After 74 days, there was a sudden decrease in OMAs escaping from box 8, and this loss was more pronounced for the 429 µm (Fig. 3) and 4000 µm (Fig. 5) than for the 2000 µm OMAs (Fig. 4).The residence times of the 4000 µm showed a significant depletion in box 9 (Fig. 5).On the 46th day of simulation, there was an abrupt decrease of aggregate fraction followed by another significant escape on the 74th day.
The 429 and 2000 µm OMAs however showed a gradual and less pronounced depletion with time (Figs. 3 and 4).Box 10 at 3189 m depth showed a significant depletion in the residence time of the 4000 µm OMAs on the 74th day (Fig. 5), whereas 429 and 2000 µm OMAs presented a gradual and less pronounced depletion as was the case of box 9 (Figs. 3  and 4).

OMA dispersion patterns
A higher percentage of OMAs escaped from the shelf break box 2 and from the offshore boxes 8, 9, and 10 for size classes 429 and 4000 µm when compared to the 2000 µm size class (Table 2).Very few 2000 µm OMAs escaped from the boxes along the canyon axis depth gradient, with box 2 and 10 showing a slightly higher escape percentage.When com-paring the 429 µm and 4000 µm OMA size classes, a higher percentage of 429 µm OMAs escaped from box 1 to 6, while from box 7 to 10 the 4000 µm OMAs showed higher percentages of escape (Table 2).Figures 6 and 7 represent the dispersion patterns for the 429 µm and 4000 µm OMAs in each box predicted by the model for an initial period of 22 days, the half-life of fresh phytodetritus (Sun et al., 1991).The figures show the aggregate trajectories along the depth gradient during that period.The 429 µm OMAs from box 2 at the shelf break were dispersed and transported in different directions (Fig. 6).OMAs travelled southward along the coast with the Portugal current, up-canyon in the direction of the coast and downcanyon (Fig. 6).The 4000 µm size class OMAs were only dispersed down-canyon (Fig. 7) for the first 22 days simulated.Within the lower canyon region, the 429 µm OMAs from boxes 8 and 10 were mainly dispersed up-canyon after 22 days, with those from box 9 showing a symmetric dispersion on the up-down canyon direction (Fig. 6).The dispersion of the 4000 µm OMAs from the same boxes was not appreciable when compared with the 429 µm OMAs (Fig. 7).Boxes 1, 3, 4, 5, 6 and 7 for both 429 and 4000 µm OMAs did not show considerable dispersion.

OMA behaviour
The 429 µm OMAs at the shelf break (box 2) and in the lower region of the canyon (boxes 8, 9, and 10) travelled longer distances (Fig. 8) and at higher velocities (Fig. 10) than the 2000 µm and 4000 µm OMAs.The highest distance value for the 429 µm was in box 2, while for the two classes of phytodetrital aggregates it was in box 9 (Fig. 8).The displacement was higher in box 2 for the 429 and 2000 µm size classes and in box 2 and 9 for the 4000 µm size class (Fig. 9).The velocities of phytodetrital aggregates were higher in box 9, while for the 429 µm in box 2 (Fig. 10).The 2000 µm OMAs travelled the shortest distance and at the lowest velocities.On average the 429 µm OMAs travelled 2.5 times farther away and with a speed 8 times higher than the 2000 µm and 2.2 times farther away and 7 times faster than the 4000 µm OMAs.In terms of displacement, the 2000 µm travelled a net distance 0.34 km and 0.47 km less than the 4000 µm and 429 µm OMAs respectively.OMAs at the remaining boxes generally showed short travelling distances, displacements and velocities.

Discussion
The conceptual model of OMA transport drawn from the model results mostly agrees with what other authors have described for the Nazaré canyon.This holds especially true for the 429 µm sized particles.In comparison to the transport of small lithogenic particles (De Stigter et al., 2007), large aggregates do not travel over long distances due to their different transport behaviour.The different sections of the canyon show different patterns of resuspension, transport and deposition of OMAs.From the upper to middle canyon regions, tidal currents are an important mechanism of resuspension and transport of sedimentary particles (De Stigter et al., 2007), and the residence time of the OMAs showed a sinusoidal pattern for boxes 1 to 5 at the upper canyon (Figs.3-5), also indicating a close match with the semidiurnal peaks of the tides (Vitorino et al., 2002).
The canyon head was characterized by active transport of OMAs, particularly the 429 µm size class.Larger amounts of OMAs escaped from box 2 (Table 2) and travelled up to 168 km (Fig. 8), and at maximum velocities of 568 km yr −1 (Fig. 10).Here, longest displacements (Fig. 9) and dispersion were observed within the canyon, up and down the canyon, as well as southwards along the coast (Fig. 6).A large percentage of the 4000 µm size class OMAs also escaped from box 2, and showed long displacement (Fig. 9) and some dispersion down canyon (Fig. 7).However, these OMAs travelled at much lower velocities with maximum distances of 9.8 km yr −1 and for much shorter distances of maximum values of 2.9 km.At the canyon head, the 429 µm OMAs exhibited the highest lateral carbon flux with the 4000 µm class being the next.This active lateral transport could be associated with the formation of nepheloid layers at these depths (Van Weering et al., 2002;Oliveira et al., 2002;De Stigter et al., 2007).
In the middle of the upper canyon (from box 3 to 6), OMA transport slowed down as indicated by the small percentages of the three OMA classes escaping from the boxes (Table 2), the very high residence times in the canyon (Figs.3-5), the lack of dispersion (Fig. 6 and 7), and no appreciable travel distances (Fig. 8) and displacements (Fig. 9), which occur at the slowest velocities (Fig. 10).Hence, the large amounts of OMAs remaining in the boxes and the lack of lateral transport indicated that the transport of OMAs in this region is dominated by short travel times followed by rapid deposition due to their increased settling velocities.This supports the idea that this section of the canyon is a deposition area of sedimentary organic matter (Schmidt et al., 2001;Van Weering et al., 2002).The OMAs with the highest residence times were the 2000 and 4000 µm sizes (Figs.3-5), which barely moved as indicated by their extremely short distances, displacements and velocities .Hence, these large phytodetrital aggregates are the major contributors in terms of carbon flux to the sediments at this region of the canyon, which may fuel the benthic communities with a food source.Although faunal abundances and biomass generally show a decreasing trend with increasing water depth, higher amounts of fresher phytodetritus and labile organic matter characterize this region of the canyon (García and Thomsen, 2008;Pusceddu et al., 2010), where the higher faunal abundances and biomasses have been found (García et al., 2007, Koho et al., 2007).Farther down, also within the middle upper canyon, the model simulations show a slight increase in the lateral carbon fluxes at box 7 at 1498 m depth.This box showed a slight increase in velocities (Fig. 10), displacements (Fig. 9) and distances (Fig. 8) of particularly the 429 µm and 4000 µm OMA size classes, and a slight increase of the percentages of particles escaping from the box (Table 2).We barely identify dispersion of OMAs though (Figs. 6 and 7), and the aggregates with 2000 µm systematically show low travelling velocities, displacements, distances and box escape percentages.We therefore conclude that this region acts as a transitional zone and is mostly characterized by a depositional regime, but where a certain amount of lateral transport occurs.Indeed, favourable conditions for sediment resuspension have been described for this region of the canyon (De Stigter et al., 2007;Oliveira et al., 2007;Martín et al., 2011).High current speeds have been observed at ∼ 1600 m depth in combination with high mass fluxes of particulate matter (Martín et al., 2011), which may explain the slight increase of lateral transport in our results.The model simulations were only carried out for a spring period and do not consider the possible role of intermittent sediment gravity flows in the transport of material.Thus the OMA transport predictions could be underestimations.If the highly energetic winter conditions were taken into account, enhanced resuspension and transport of OMAs through this part of the canyon might have been more conspicuous.This could be further evaluated in future model simulations.
The offshore region of the canyon was characterized by resuspension of OMAs and acceleration of the transport.At boxes 8, 9 and 10, the OMAs residence times were low (Figs.3-5), accompanied by high escape percentages from the boxes, particularly of the 4000 µm size class (Table 2).There was an increase in OMA travelling distances (Fig. 8), displacements (Fig. 9) and velocities (Fig. 10) that reached similar values to the ones observed at the canyon head (box 2).The 429 µm size class was again the driver of the particle flux reaching maximum velocities of 230 km yr −1 and distances of 68 km in box 9 (Figs. 10 and 8).The phytodetrital aggregates, particularly the 4000 µm, showed also an active transport but not as pronounced as the 429 µm.Box 10, located in the middle canyon (3189 m), showed a slight decrease in the carbon flux with average velocities ranging from 2.6 km yr −1 to 0.6 km yr −1 (Fig. 10) and average distances ranging from 0.8 km to 0.2 km (Fig. 8) for the three different OMA classes.These boxes were located between 2077 and 3189 m water depth in a steep section of the canyon under the influence of high bottom currents and internal waves (De Stigter et al., 2007;Martín et al., 2011).High current speeds and variable seasonal fluxes were observed in springsummer at ∼ 3300 m (Martín et al., 2011), which would explain the more active nature of this part of the canyon in terms of sediment resuspension and horizontal carbon flux.

Conclusions
Exploring the potential of the operational modelling, the MOHID-PCOMS was used to give the necessary boundary conditions to apply a hydrodynamic model in the Nazaré canyon.A lagrangian transport model was successfully coupled to the MOHID model, giving an overview of the OMA transport patterns along the Nazaré canyon bottom depth.The model simulations were performed during the spring season when phytodetritus production is high.With respect to our original hypothesis, the model results show that, during the specific time of investigation (spring 2009), the canyon did not function as a conduit of organo-mineral aggregates to the deep sea.Rather, it acts as a temporary depocentre of sedimentary organic matter during spring conditions.Previous studies (De Stigter et al., 2007;Martín et al., 2011) indicate that this may not always be the case, and the canyon is an active conduit of sediment transport to the deep sea.The model results show that the carbon flux in the canyon is not constant and unidirectional within areas of resuspension, transport and deposition.For instance, large, carbon-rich aggregates with their specific transport behaviour are not exported but rather remain in a given area for long periods of time.These aggregates, however, are frequently resuspended into the BBL and therefore allow mineralization to occur under turbulent conditions of the BBL.This is in agreement with other studies carried out within the canyon.The differences between transport patterns of the median OMAs and phytodetrital aggregates were also predicted by the model, and the lateral transport of the larger OMAs is less pronounced than for the median OMAs resulting in the carbon deposition.On the other hand, the model results did not include the possible role of sediment gravity flows, and therefore may underestimate the rates of OMA transport.Nevertheless, the model could also be applied to evaluate the transport patterns of other substances in the canyon such as pollutants.Further studies are required to analyse the differences in the carbon flux transport in an autumn-winter season and the impact of the river discharges on the increasing carbon fluxes in Nazaré canyon.

Fig. 3 .
Fig. 3.The residence time of the 429 µm OMAs in each monitor box over the simulated period.The monitor boxes are set at increasing depths from box 1 to box 10: 59 m, 262 m, 357 m, 575 m, 331 m, 945 m, 1498 m, 2077 m, 2657 m, and 3189 m.

Fig. 4 .
Fig. 4. The residence time of the 2000 µm OMAs in each monitor box over the simulated period.The monitor boxes are set at increasing depths from box 1 to box 10: 59 m, 262 m, 357 m, 575 m, 331 m, 945 m, 1498 m, 2077 m, 2657 m, and 3189 m.

Fig. 5 .
Fig. 5.The residence time of the 4000 µm OMAs in each monitor box over the simulated period.The monitor boxes are set at increasing depths from box 1 to box 10: 59 m, 262 m, 357 m, 575 m, 331 m, 945 m, 1498 m, 2077 m, 2657 m, and 3189 m.

Fig. 8 .
Fig. 8. Distance predicted by the model for the three classes of OMAs.

Fig. 9 .
Fig. 9. Displacement predicted by the model for the three classes of OMAs.

Fig. 10 .
Fig. 10.Velocity predicted by the model for the three classes of OMAs.

Table 1 .
OMA characteristic data used in the model simulations: aggregate size d (µm); settling velocity W s (cm s −1 ); critical and depositional velocities U * cr and U * d (cm s −1 ); critical and depositional bottom shear stresses T cr and T d (N m −2 )

Table 2 .
Percentage of OMAs escaping from the monitor boxes predicted by the model.