Potential vorticity conserving flows and vortex-wave interaction : the role of vertical velocity and isopycnal diffusion on plankton heterogeneity

POTENTIAL VORTICITY CONSERVING FLOWS AND VORTEX-WAVE INTERACTION: THE ROLE OF VERTICAL VELOCITY AND ISOPYCNAL DIFFUSION ON PLANKTON HETEROGENEITY Barcelona, Febrero de 2012 Mariona Claret Ilustraci´on de la contraportada realizada por Marc Gasser Rubinat. Anexo I D/Dª...............................................................SECRETARIO/A DEL DEPARTAMENTO DE FISICA DE LA UNIVERSIDAD DE LAS PALMAS DE GRAN CANARIA, CERTIFICA, Que el Consejo de Doctores del Departamento en su sesión de fecha.............................tomó el acuerdo de dar el consentimiento para su tramitación, a la tesis doctoral titulada “Potential vorticity conserving flows and vortex-wave interaction: the role of vertical velocity and isopycnal diffusion on plankton heterogeneity” presentada por la doctoranda Dª Mariona Claret Cortés y dirigida por el Doctor Álvaro Viúdez Lomba y la Doctora Yvette H. Spitz. Y para que así conste, y a efectos de lo previsto en el Artº 73.2 del Reglamento de Estudios de Doctorado de esta Universidad, firmo la presente en Las Palmas de Gran Canaria, a…......de.........................................de dos mil doce. Potential vorticity conserving flows and vortex-wave interaction: the role of vertical velocity and isopycnal diffusion on plankton heterogeneity (Flujos que conservan la vorticidad potencial e interacción onda-vórtice: el efecto de la velocidad vertical y la difusión isopicna en la heterogeneidad planctónica) Tesis doctoral presentada por Mariona Claret Cortés dirigida por el Dr. Álvaro Viúdez Lomba y codirigida por la Dra. Yvette H. Spitz para obtener el grado de Doctora por la Universidad de Las Palmas de Gran Canaria, Departamento de Física, Programa en Oceanografía (Bienio 2006-2008). En Barcelona, a febrero de 2012 La Doctoranda El Director La Codirectora Institut de Ciències del Mar Per a la meva `avia, ´unica en su especie. Per tenir un bomb´o sempre a punt, per viure el dia a dia amb ironia, per la il·lusi´o dels seus ulls. THESIS ABSTRACT This thesis investigates physical-ecological and vortex-wave interactions through potential vorticity (PV) considering a stratified and oligotrophic ocean. To this end, a NPZ (Nutrients-Phytoplankton-Zooplankton) model is coupled to a physical one that conserves PV explicitly on isopycnals. The physical-ecological coupled model is initialized using stationary NPZ solutions numerically stable with the fluid at rest. These solutions are implemented homogeneous both on horizontal and isopycnals levels to quantify the effect of horizontal and vertical advection caused by mesoscale and submesoscale vortex structures, and isopycnal mixing. At the interior of the vortex separatrix, plankton and PV distributions translate in phase at vortex propagation speed. Within cyclones, isopycnal doming enhances plankton biomass at the vortex center in different trophic conditions. Furthermore, isopycnal mixing associated to small-scale motions maximizes the phytoplankton (P) biomass in cyclones through a resonant response between Pand diffusive timescales. This Pincrease is significant in mesotrophic conditions and occurs where the vertical displacement of isopycnals is maximum, and hence where vertical gradients of PV are large. At the separatrix outer, horizontal and vertical advection are of the same order of magnitude than the ecological forcing and enhance Pthrough different mechanisms. Firstly, vertical velocity wuplifts nutrients and Pto better lit levels. Presponds with some time lag to this perturbation and the associated increase in biomass occurs far from the upwelling location due to the action of horizontal advection. As a result, Pcorrelates with w, and thus with horizontal gradients of PV, only at initial times. In the particular case of translating cyclones, this mechanism explains the development of a Ptrail at their wake. And secondly, the horizontal advection of a surface ecosystem patch by subsurface vortices decreases Pself-shading at the patch front in benefit of Pgrowth. Finally, interactions between vortex structures and pure inertial and gravity large amplitude waves are investigated. The advection of PV by waves causes vortices to be unsteady and modifies the upper and lower bounds of the wave frequency band. The advection of waves by vortices Doppler shifts the local wave frequency. When inertial waves are involved, a near-inertial right-handed helical wave is developed due to a non-linear interaction. As a result, total wincreases one order of magnitude and correlates with horizontal gradients of PV. These results aim to shed further light on the ecological impact of long-lived coherent vortices in the open ocean. LIST OF SYMBOLS Physics ρmass density ρ′density anomaly ρ0constant density α0constant specific volume ppressure p′pressure anomaly zstratification constant d(x, t) depth in the reference density configuration of isopycnal located at (x, t) D(x, t) vertical displacement of isopycnals in the spatial description D(s, d) vertical displacement of isopycnals in the isopycnal description fCoriolis frequency feeffective Coriolis frequency Nbackground buoyancy frequency Ntotal buoyancy frequency cPrandtl ratio ωllocal wave frequency ωpparticle wave frequency u= (u, v, w) velocity vector uggeostrophic velocity wqQG vertical velocity Qg hgeostrophic Q-vector Uddipole speed ω≡ωh+ζkrelative vorticity ωggeostrophic vorticity ω′ageostrophic vorticity PV anomaly Πtotal PV ϕ= (ϕ, ψ, φ) vector potential FFroude number RRossby number κdimensional diffusion coefficient Knon-dimensional diffusion coefficient Zmin minimum dimensional vertical coordinate z Lspatial conversion factor Ttemporal conversion factor Ecology Pphytoplankton Zzooplankton Nnutrients NTtotal nitrogen GPphytoplankton production rate GZzooplankton production rate K0half-saturation for Puptake Lphotosynthetic rate Iavailable radiation I0surface available radiation Awlight attenuation by sea water Aplight attenuation by P Ψ0initial slope P-Icurve V0Pmaximum uptake rate MPphytoplankton mortality MZzooplankton mortality Ξ0Pspecific mortality rate Rzooplankton grazing rate R0zooplankton maximum Λ0Ivlev constant Γ0fraction of Zgrazing egested Ψ0Zexcretion/mortality rate Θ0Zspecific excretion/mortality rate Chapter 1 Introduction Marine ecosystems are highly sensitive to ocean physics. Pelagic organisms are embedded in the fluid and thus they are explicitly affected by its dynamics. Since most of these individuals are non-motile, they drift in the vast ocean. These wanderers owe their name to the Greek equivalent, plankton. Two microscopic groups are enclosed within plankton, autotrophs and heterotrophs. The formers are called phytoplankton, where phyto denotes its plant-like character due to the performance of photosynthesis. The latters are the animal-like plankton, the zooplankton, since they feed from phytoplankton or the smallest zooplankton. Phytoplankton is of great importance because is the basis of the oceanic food web. Its growth is affected by fluid motion because growth limiting factors, light and nutrients, are vertically segregated in the ocean. On the one hand, the well-lit layer is akin to the tip of an iceberg, it is the upper fraction of an ocean which is on average twenty times deeper. On the other hand, nutrients slowly and incessantly sediment into the dark deep ocean. As a result, vertical velocities in the ocean play an implicit role in ecosystems by fertilizing the well-lit zone. Many works have been conducted to characterize the ecological patterns associated to specific physical structures. In this regard, Haury et al. (1978) first drew an ecological counterpart of the Stommel diagram, which plots physical ocean variables in a spatio-temporal logarithmic frame, to represent zooplankton biomass variability. Nowadays, physical-biological interactions have been extensively reported from planetary to Kolmogorov scales (Steele, 1978;Mann and Lazier,1991;Denman and Gargett,1995). The combination of these wide range of spatial scales results in complex and deterministic plankton patterns, which are reflected in the sea surface chlorophyll distributions as measured from space satellites. Highly 1 2 CHAPTER 1 productive regions are localized at western boundaries, sub-polar latitudes, and the Equator (Yoder et al.,1993), where coastal upwelling, water mass subduction, and Trade Winds forcing, respectively, overcome the light-nutrient vertical segregation. These areas are called eutrophic, because communities are not limited by nutrients, as opposed to oligotrophic zones. The largest part of the open ocean is oligotrophic. Its sea surface chlorophyll concentration (CC) has a geometric mean about two orders of magnitude smaller than that of eutrophic waters, and a variability dominated by submonthly mesoscale variance (Doney et al.,2003). Since CC images carry information about the ocean surface turbulent flow (Nieves et al., 2007), plankton variability in the oligotrophic open ocean may be likely caused by mesoscale phenomena, and vortical structures may be particularly involved. The dynamics and threedimensional structure of ecosystems in the open ocean remain largely unknown because synoptic and high spatial resolution sampling there is more costly and challenging than in the coastal ocean. In this context, the aim of this thesis is to provide further insight on how ubiquitous mesoscale eddies alter the basis of the food web in oligotrophic environments. To this end, we analyze first which factors shape plankton patterns at mesoscales in order to choose an appropriate physical descriptor (section 1.1). We then review the mechanisms through which vortices generate ecological heterogenity (section 1.2) and identify those which are related to our purpose (section 1.3). Finally, the specific objectives of this thesis are outlined by chapters (section 1.4). 1.1 Plankton patchiness at mesoscales: Ecological footprint of potential vorticity First synoptic maps of sea surface chlorophyll were obtained in the early 80’s (Gordon et al., 1980;Gower et al.,1980). At first glance, these type of maps show high correlation with sea surface temperature at mesoscales, which has led to consider plankton as a passive tracer to a first approximation. However, discrepancies between physics and ecology arise when they are quantified (see review Martin,2003). Two main facts account for these differences, nutrient pumping at smaller scales and plankton time lag response to perturbations. Submesoscale structures have gained increasing attention the last decade because they enhance vertical velocity wone order of magnitude with respect to mesoscale structures and are responsible for injection of allochthonous nutrients to the photic zone (see review Klein and Lapeyre, 1.1. PLANKTON PATCHINESS AT MESOSCALES: ECOLOGICAL FOOTPRINT OF POTENTIAL VORTICITY 3 2009), where light irradiance is greater than one per cent of that arriving at the sea surface. Since physical and ecological timescales are similar at the submesoscale, plankton couples to this nutrient upwelling and wcollocates with primary production (L´evy et al.,2001) and some photosynthetic indexes (Falkowski,1983;Cullen and Lewis,1988). However, exact correlation with wdepends on the ecological parameter timescale. Those parameters with timescales greater than that of horizontal advection, such as phytoplankton biomass, show a spatial lag with w(Lima et al.,2002). The importance of inherent ecological timescales in plankton patterns was nicely illustrated by Abraham (1998) using a two-dimensional turbulence numerical model. He observed an increasing patchiness from physics to phytoplankton, and to zooplankton. Therefore mesoscale plankton distributions result from an interplay between horizontal advection, vertical advection, and intrinsic plankton timescales. Whether these distributions are caused by ecological or physical phenomena depends on their relative timescale (Mahadevan and Campbell,2002). In order to relate both timescales, we seek to describe the fully three-dimensional flow nature with a single timescale. A very useful physical quantity that relates horizontal and vertical motion is potential vorticity (PV). The concept of PV was introduced by Beltrami in 1871, applied to an adiabatic inviscid multi-layered ocean by Rossby in 1936, and extended to baroclinic flows by Ertel in 1942 (check Vi´udez,2001, for cites and relation between PV definitions). According to the latter, specific PV is defined as Π≡ω+fk ρ·∇Tθ,(1.1) where the sum of vorticity ω= (ξ, η, ζ) relative to a reference frame rotating with Earth’s angular velocity and the vertical component of the planetary vorticity f, which is the Coriolis parameter, is the absolute vorticity. Additionally, ρis the mass density, Tθthe potential temperature, though it could be any scalar fluid property materially conserved (Pedlosky,1987, chapter 2), and ∇the three-dimensional gradient operator. The geophysical significance of PV lies on its material invariance, that is, dΠ dt= 0 .(1.2) One way to derive (1.2) comes from the conservation of circulation for frictionless motion considering material conservation of Tθand mass conservation (Pedlosky,1987, chapter 2). 4 CHAPTER 1 Thus (1.2) can be interpreted as the conservation of circulation along material circuits and angular momentum for a given volume, which let us to relate horizontal and vertical motions. When the fluid is barotropic this relation is straightforward since (1.1) becomes Π≡ζ+f h,(1.3) where his the vertical separation between neighboring material isosurfaces (eq. 3.4.9. Pedlosky,1987). Thus changes in himply adjustments in ζin order to conserve PV, which explains the well-known ballerina and ice skater effects. Instead, when the fluid is baroclinic the tilting of material isosurfaces introduces nontrivial changes in ω. Furthermore, if some balance condition is established in the momentum equations (along with proper boundary conditions), then for a given distribution of PV at a fixed time the velocity vector, pressure, and density three-dimensional fields associated to the balanced flow (void of inertia–gravity waves IGWs) can be recovered in a process known as PV inversion (Hoskins et al.,1985; McIntyre and Norton,2000;Vi´udez,2008a). PV is relevant to ecological dynamics for several reasons. Firstly, it describes the balanced velocity vector field and the vertical displacement of isopycnals D. Secondly, it relates vertical upwelling with horizontal advection, which is crucial to plankton patchiness (Martin et al.,2002). The importance of PV in ecology was first pointed out by Woods (1987) and Strass and Woods (1987). They observed chlorophyll increase in regions of large PV isopycnal gradients, which were thought to be areas of high w. Their hypothesis agrees with experimental (Pall`as-Sanz and Vi´udez,2005) and numerical (Vi´udez and Dritschel,2003, 2004b) works that demonstrate that zones with large horizontal gradients of ζand PV are related to high values of w. Large PV horizontal gradients involve large horizontal gradients in ρand therefore in D. Consequently, when these Dgradients are advected, large local rates of Dand woccur. However, vertical advection does not account for the whole plankton big picture aforementioned. 1.2 How do mesoscale vortices alter ecosystems? Many in-situ observations reveal eddies impact on marine biology. At the basis of the food web, vortices qualitatively affect community structure (Thompson et al.,2007;Huang et al.,2010), physiological processes (Bibby et al.,2008), and ecosystems transport (Batten 1.2. HOW DO MESOSCALE VORTICES ALTER ECOSYSTEMS? 5 and Crawford,2005). Additionally, they quantitatively alter biogeochemical balances by enhancing new primary productivity (NPP) (Mor´an et al.,2001). For instance, they are involved in the North Atlantic carbon balance, though their specific contribution is under scientific debate ranging from 50% to less than 10% of NPP (see review Oschlies,2008). As a result, vortices also perturb higher trophic levels (Mackas et al.,2005;Atwood et al.,2010). There are different mechanisms by which ocean eddies introduce the above mentioned ecological variability. Next, we analyze their ecological spatial signature in order to relate it with PV when this is possible. •Horizontal advection. Mesoscale two-dimensional turbulence induces a conservative transfer from large to small scales (Abraham,1998). If we consider that phytoplankton behaves as a passive tracer, then its pattern is determined by vorticity and strain. Inside vortices, vorticity dominates over strain and particles are trapped tracing nearly orbital trajectories. This explains long-term transport of chlorophyll rich waters offshore the Algerian (Arnone and LaViolette,1986), Alaskan (Batten and Crawford, 2005), and northwest African (Pelegr´ı et al.,2005) coasts among many others. Outside vortices, strain overcomes vorticity, resulting in chaotic motion (Provenzale,1999). As a result, eddies stir plankton patches into spirals and filaments (Lehahn et al.,2007) and even wave-like structures (Menkes et al.,2002) at their edges. •Eddy pumping. This term was coined by Falkowski et al. (1991) to denote NPP enhancement due to upwelling of the nutricline into the photic layer caused by the passage of a translating surface cyclone or subsurface anticyclone (McGillicuddy et al.,1999). On the one side, isopycnal uplift is greatest at the eddy center, which explains high values of sea surface chlorophyll at eddy cores (Ar´ıstegui et al.,1997;Barton et al.,1998; Mizobata et al.,2002;Siegel et al.,2008;Tew-Kai and Francis,2009;Siegel et al.,2011). On other side, PV is correlated with the isopycnal vertical displacement (Vi´udez and Dritschel,2003). Consequently, we expect chlorophyll enrichment at maxima absolute PV values. •Submesoscale. Fully-developed vortices are rarely spherical in the isotropic quasigeostrophic (QG) space. The QG space is the vertically stretched space of dimensions (x, y, cz), where c≡N/f is the Prandtl ratio, Nthe background buoyancy frequency, and fthe inertial frequency. Instead, a variety of elliptical geometries are shaped 6 CHAPTER 1 during their life-time, which lead to vortex rotation, and ultimately to a quadrupolar pattern distribution of vertical velocity w(Vi´udez and Dritschel,2003). Vertical velocity maxima are reached at vortex edges enhancing chlorophyll concentrations at the eddy periphery (Mizobata et al.,2002;Ladd et al.,2005;Siegel et al.,2011). Thus chlorophyll would be better related to PV horizontal gradients than to PV. However, to unveil this wpattern, resolutions smaller than 10 km, which are close to submesoscale sampling, are often required (L´evy et al.,2001). •Eddy-eddy interactions. Vortices interact in a number of ways through bounding, merging, or conforming complex structures (Voropayev and Afanasyev,1994). Since these interactions modify the flow, they also perturb ecological distributions. For instance, the ubiquitous presence of mushroom-like shapes of sea surface chlorophyll in frontal coastal zones (Sur et al.,1996;Stapleton et al.,2002) are imprints the interaction between two vortices of opposite PV, that is, the vortex dipole. Another example is the merging of anticyclones in the Kuroshio Current, which results in an increase of wwhich triggers phytoplankton enhancement at vortex edges (Yoshimori and Kishi,1994). •Isopycnal mixing. As aforemetioned, Dmaxima occur at the monopole center, which creates isopycnal gradient of light irradiance. In oligotrophic regimes phytoplankton growth is mainly limited by light and nutrients. Since isopycnal doming modifies the amount of light at which isopycnal confined phytoplankton is exposed, it affects its growth, and indirectly nutrient consumption. Thus isopycnal diffusion may increase the phytoplankton response to light at cyclone cores and anticyclone edges through a down-gradient replenishment flux of nutrients. •Vortex-wave interaction. On the one hand, interactions between vortices and IGWs are ubiquitous in the ocean. Some of these interactions imply inertial to near-inertial wave frequency shift (Perkins,1976) and trapping of wavepackets within anticyclones (Kunze,1986). On the other hand, phytoplankton distributions are perturbed by both eddies, and high frequency oscillations (Franks,1995a;Granata et al.,1995;G´omez et al.,2001;Sangr`a et al.,2001). Thus, we expect vortex-wave interaction to affect phytoplankton dynamics. 1.3. RELEVANT PHYSICAL-ECOLOGICAL INTERACTIONS IN COHERENT VORTICES 7 Several other vortical mechanisms have been proposed to explain plankton patterns. For instance, wind forcing increases primary production either at the vortex edges through eddy-wind interaction (Martin et al.,2001;Mahadevan et al.,2008) or anticyclone cores by deepening of the mixed layer (Thompson et al.,2007). The contribution of these mechanisms as well as diapycnal mixing to phytoplankton patchiness is beyond the scope of this work. 1.3 Relevant physical-ecological interactions in coherent vortices Vortices life cycle involve three different stages with different dynamics: formation, maturity, and decay (Sangr`a et al.,2005). High turbulent areas, such as coastal fronts, are hot-spots of vortex formation and decay. In contrast, the stratified open ocean is mainly dominated by long-lived mature vortices. Fully-developed coherent vortices translate across the ocean conserving their PV. The outer PV isosurface defines the vortex separatrix, which marks vortex limits and acts as an impermeable barrier. Within the separatrix the fluid rotates approximately as a solid body. Outside the separatrix, the horizontal velocity decreases exponentially with radius. Thus vortices can trap waters in their interior and transport ecosystems long distances. If those ecosystems are only affected by horizontal advection, the transporting distance will depend on the sinking rate of nutrients to the aphotic layer. This leads to long-term chlorophyll depletion within vortices unless some nutrient upwelling occurs. In this regard, eddy pumping assumes entrainment of new nutrients into the vortex (Martin and Pondaven,2003), which implies mixing across vortex boundaries and thus restricts the mechanism to developing cyclones and decaying anticyclones (Franks et al., 1986b). This mechanism causes the initial seed of ecological heterogeneity though does not explain how it is maintained during vortex translation. In contrast, the vertical velocity of mesoscale and submesoscale monopoles or complex vortical structures may contribute to ecosystem subsistence in one year living vortices. Additionally, nutrient injection by isopycnal diffusion is likely to occur in the open ocean, where the flow is mainly two-dimensional. Though the contribution of isopycnal mixing in upwelling new nutrients is minimal compared to the mentioned mechanisms (Siegel et al.,1999;Ledwell et al.,2008), it may play some role in upwelling regenerated nutrients. Finally, the way vortex-wave interaction affects ecology remains largely unknown. 2.1. INTRODUCTION 15 2.1 Introduction The dynamics of oceanic planktonic ecosystems is often investigated using nutrient-phytoplankton-zooplankton (NPZ)-type numerical models (Wroblewski,1977;Franks et al.,1986a; Franks,2002;Newberger et al.,2003). These models are discrete versions of the continuous partial differential equations for the field variables NPZ. In the absence of flow and horizontal gradients this system of equations has a number of one-dimensional (1D) steady analytical solutions which are very useful both to characterize different planktonic regimes and to serve as initial conditions in coupled physical-ecosystem numerical modeling. It may happen, however, that some of these 1D steady solutions be only continuous, but not continuously differentiable functions of the vertical coordinate zalong the water column. This potential lack of differentiability makes these NPZ solutions inappropriate as initial conditions in three-dimensional (3D) coupled physical-ecosystem models. These 3D coupled models, if formulated as it is usual in the spatial description, require existence of vertical derivatives for the vertical advective terms, present in the material rate of change of the ecosystem quantities, make sense. In this paper we present a particular example of a motionless, 1D steady NPZ solutions which, being continuous but not continuously differentiable at several depths, do however admit numerical equivalents that are continuously differentiable in the numerical sense (here meaning convergence of the vertical derivative with respect to increasing vertical resolution). In the next section the NPZ system of equations is briefly introduced. Though more sophisticated ecological models exist a simple one is used here because we seek to keep the number of free parameters as small as possible, while retaining the essential behavior of the ecological fields, particularly the development of NPZ anomalies due to the vertical advection of nutrients into the euphotic zone. The NPZ-type models are ecological bulk algorithms subjected to significant errors in the mathematical parametrization of the different ecological processes (see Anderson,2005;Flynn,2005;Mitra et al.,2007). The motionless 1D continuous steady ecological solutions and their numerical continuously differentiable equivalents are obtained in section 2.3. Convergence of this differentiable solution is reached however at vertical resolutions of a few centimeters. As an application example, in section 2.4 the differentiable steady NPZ solutions are used as initial conditions in a 3D coupled physical-ecosystem model (with poor vertical resolution) to address the role 16 CHAPTER 2 of the vertical velocity in a case of oceanic baroclinic instability. Concluding remarks are given in section 2.5. 2.2 NPZ Model The dependent variables of the NPZ model are the dissolved inorganic nitrogen (N), the phytoplankton (P), and the zooplankton (Z) biomass. These variables are expressed in units of concentration of nitrogen (here always in mmol N m−3), and satisfy the system of equations (Wroblewski,1977;Newberger et al.,2003) dP dt=N K0+NL P |{z } GP −R0(1 −e−Λ0P)Z |{z } R −Ξ0P |{z} MP ,(2.1) dZ dt= (1 −Γ0)R |{z } GZ −Φ0Z |{z} MZ ,(2.2) dN dt=−GP+ Γ0R+MZ+MP.(2.3) Above, the material rate of change dχ/dt≡∂χ/∂t +u·∇χis the sum of the local and advective rates of change of χ,u= (u, v, w) is the 3D velocity, and ∇is the 3D gradient operator. The Pincreases due to its production rate (GP), which depends on both the inorganic dissolved nitrogen uptake and the photosynthesis, and decreases due to the herbivore grazing (R) and the phytoplankton mortality rate (MP). Constant K0is the half-saturation concentration for phytoplankton uptake of nutrients. The Zincreases due to the ingestion of phytoplankton (GZ) and decreases due to the zooplankton specific excretion and mortality rate (MZ). The Pproduction (GP) follows the Michaelis-Menten kinetics and depends on the following photosynthetic rate, L(x, t)≡V0Ψ0I(x, t) pV2 0+ Ψ2 0I2(x, t)(2.4) (Newberger et al.,2003), which depends on the photosynthetically available radiation I(x, t)≡I0exp ½Awz−ApZ0 z P(x, y, z′, t)dz′¾,(2.5) 2.3. STEADY SOLUTIONS 17 where z≤0. Above, V0is the phytoplankton maximum uptake rate, Ψ0is the initial slope of the P-Icurve, Awis the extinction coefficient of seawater in the absence of phytoplankton (following the Lambert-Beer law and light variation with day time is not considered), I0exp{Apκ(x, t)}is the light attenuation by phytoplankton self-shading, where constant I0 is the surface photosynthetically available radiation, and Apis the extinction coefficient per unit concentration of phytoplankton. Adding (2.1)+(2.2)+(2.3) the total nitrogen NT≡N+P+Zis materially conserved, dNT dt= 0 .(2.6) It is convenient therefore to define the system of independent equations as (2.1)-(2.2)-(??), instead that the original set (2.1)-(2.2)-(2.3). In the next section two particular solutions of these equations are obtained in the case of steady and horizontally homogeneous ecosystem distributions in the absence of flow. 2.3 Steady Solutions 2.3.1 Analytical Steady Solutions In order to obtain analytical steady solutions of horizontally homogeneous distributions we set dP/dt= dZ/dt = 0 in (2.1) and (2.2), so that spatial functions depend only on z. The water column is divided into three layers. In the upper layer P6= 0 and Z6= 0, in the mid layer P6= 0 and Z= 0, and in the lower layer P=Z= 0. In the upper layer (z∈(z1,0]) the steady P, obtained directly from (2.2) (e.g., Busenberg et al.,1990, eq. 8; Newberger et al.,2003, eq.18) is independent of z P(z) = P1=−1 Λ0 ln ·1−Φ0 R0(1 −Γ0)¸, z ∈(z1,0] .(2.7) With the commonly used parameters for upwelling conditions given in Table 2.1,P1= 8.468 mmol N m−3. The values and a sensitivity analysis of these parameters are given in (Newberger et al.,2003). The maximum depth z1at which this solution is feasible is determined below. From (2.1) and (??) the zooplankton Z(z) = Z1(z) in the upper layer 18 CHAPTER 2 description value units Awlight attenuation 0.067 m−1 Aplight attenuation by P0.0095 m2mmol N−1 Ψ0initial slope of P-Icurve 0.025 m2(W d)−1 V0Pmaximum uptake rate 1.5 d−1 I0surface available radiation 158 W m−2 K0half-saturation for Puptake 1 mmol N m−3 Ξ0Pspecific mortality rate 0.1 d−1 R0Zmaximum grazing rate 0.52 d−1 Λ0Ivlev constant 0.06 m3mmol N−1 Γ0fraction of Zgrazing egested 0.3 Φ0Zexcretion/mortality rate 0.145 d−1 Table 2.1: List of constants (Newberger et al.,2003). Pand Zstand for phytoplankton and zooplankton, respectively. z∈(z1,0] is obtained solving the quadratic equation A0Z2 1+BZ1+C= 0 ,(2.8) where the coefficients A0≡Φ0 (1 −Γ0)P1 ,(2.9) B(z)≡ −L(z)−A0(K0+NT−P1) + Ξ0,(2.10) C(z)≡[L(z)−Ξ0] (NT−P1)−Ξ0K0.(2.11) Thus, Z1(z) = −B(z)−pB2(z)−4A0C(z) 2A0 .(2.12) The negative root solution above ensures that Z(z)< NT. To further simplify the problem we consider NTas constant, independent of z.NTmust be such that the discriminant B2(z)−4A0C(z)≥0. The maximum depth z1of the upper layer is defined as the shallower depth at which Z(z1) = 0. Using (2.1) or (??) this condition implies C(z1) = 0, and therefore L1≡L(z1) = Ξ0µ1 + K0 NT−P1¶.(2.13) 2.3. STEADY SOLUTIONS 19 Inverting L(z1) using (2.4) the depth z1defining the lower boundary of the upper layer is z1=1 2(ApP1+Aw)ln ·V2 0L2 1 (V2 0−L2 1)Ψ2 0I2 0¸.(2.14) In the mid layer, z∈(z2, z1], Z(z) = Z2(z) = 0 and P(z) = P2(z)6= 0. In this layer, due to (2.3) or (2.1), P2satisfies the relation P2(z) = NT−Ξ0K0 L(z)−Ξ0 ,(2.15) where L(z) is given by (2.4) and (2.5). The solution P2(z) is found here solving (??) numerically. This solution is feasible as long as P2(z)≤NTwhich, from (??), implies that L(z)>Ξ0. Consequently, the maximum depth z2of this mid layer is defined as that at which L(z2) = Ξ0. Inversion of this equation implies that z2satisfies the relation z2−Ap AwZz1 z2 P2(z′)dz′= =1 Aw ln ÃV0Ξ0 Ψ0I0pV2 0−Ξ2 0!−Ap Aw P1z1.(2.16) For its later use it is convenient to define the maximum depth z3of this mid layer as that obtained by neglecting the phytoplankton self-shading only in this mid layer, z3≡1 Aw ln ÃV0Ξ0 Ψ0I0pV2 0−Ξ2 0!−Ap Aw P1z1.(2.17) Clearly, z3< z2. Finally, in the lower layer, z∈[zmin, z2], there is neither Pnor Z(P3= Z3= 0), so that the dissolved inorganic nitrogen N(z) = N3=NT. The choice of NTis particularly important in this NPZ solution since it must be such that the resulting distance between z1and z2be large enough to be properly discretisized using a finite grid size. Based on the behaviour of z1and z3as functions of NT(Fig. 2.1), NT must be close to P1= 8.468 mmol N m−3. We select NT= 8.8 mmol N m−3. With this choice the values of the transition depths are z1=−15.26 m, z2=−25.13 m, and z3=−36.52 m. The maximum Z, obtained from (??) at z= 0, becomes Z0=Z(0) = 0.25 mmol N m−3. The solutions NPZ at the three layers are shown in Fig. 2.2. 20 CHAPTER 2 6 8 10 12 14 16 −50 −40 −30 −20 −10 0 z* 1 z* 3 N* 0 z1 z2 z3 P1N0 z (m) Figure 2.1: Depths z∗ 1(??) and z∗ 3(2.6) as functions of the variable total nitrogen N∗ 0. Units are mmol N m−3. Note that z1=z∗ 1(NT) and z3=z∗ 3(NT). In the upper layer the amount of Zdecreases with depth (from Z0to 0) which is compensated (since both Pand NTare constant) by an equally small increase of N. In the mid layer, Pdecreases from P1to 0, and is compensated (since Z= 0) by an equally large increase of N. In this steady solution, and in the upper layer, the large z-dependent GP(z) is mainly balanced by the large constant MP (Fig. 2.3). The remanent, smaller part of GP(z) is balanced by the small z-dependent grazing R(z). In the Zbalance, the small growth GZ(z) of zooplankton is balanced by MZ(z). In the mid layer R= 0, so that GP(z) and MP(z) exactly balance. These steady NPZ solutions are continuous functions of z, but they have been obtained without any requirement on differentiability conditions. It is clear, at least visually from Fig. 2.2, that these functions are not vertically differentiable at z=z1or z=z2. As a simple proof consider differentiability of P(z) at z2. The vertical derivative of (??) is ∂P2 ∂z =Ξ0K0 L(z)−Ξ0 ∂L ∂z =Ξ0K0 L(z)−Ξ0 V3 0Ψ0 p(V2 0+ Ψ2 0I2)3 ∂I ∂z =Ξ0K0 L(z)−Ξ0 V3 0Ψ0[Aw+ApP2(z)] p(V2 0+ Ψ2 0I2)3I(z).(2.18) As z→z2, we have P2(z)→0 and L(z)→Ξ0, and therefore ∂P2/∂z → ∞ in the mid layer, as observed in Fig. 2.2a. This limit does not match with the vertical derivative of P as z→z2in the lower layer, where ∂P3/∂z = 0. This lack of differentiability (the functions are continuous but not continuously differentiable) implies that these steady solutions are questionable as initial conditions in many coupled physical-ecosystem models. Most of these models are formulated in the spatial (Eulerian) description and require that solutions NPZ be continuously differentiable for the vertical gradients present in the advective derivative on the right hand side of (2.1)-(2.2) make sense. In the next section NPZ solutions, similar to the ones described above but continuously differentiable, are numerically obtained. 2.3. STEADY SOLUTIONS 21 −2 0 2 4 6 8 10 P(z) −60 −40 −20 0 z (m) z1 z2 (a) N0 P1 upper layer mid layer lower layer 0.0 0.1 0.2 0.3 Z(z) −60 −40 −20 0 Z0 (b) −2 0 2 4 6 8 10 N(z) −60 −40 −20 0(c) N0 Figure 2.2: Steady vertical profiles of (a) P(z), (b) Z(z), and (c) N(z) in the three layers. The constant concentrations P1,NT, and Z0, as well as the transition depths z1 and z2are indicated. Units are mmol N m−3. −0.2 0.0 0.2 0.4 0.6 0.8 1.0 −40 −30 −20 −10 0 z (m) GP RGZ (a) z1 z2 −0.2 0.0 0.2 0.4 0.6 0.8 1.0 −40 −30 −20 −10 0 PM ZM (b) Figure 2.3: Vertical profiles of (a) phytoplankton production GP(z), grazing R(z), and zooplankton growth GZ(z); (b) phytoplankton mortality MP(z) and zooplankton mortality MZ(z). Units are mmol N m−3d−1. 22 CHAPTER 2 2.3.2 Numerical Steady Solutions In order to obtain continuously differentiable steady P(z) and Z(z) solutions, the prognostic equations (2.1)-(2.2) are numerically integrated in time, using as initial conditions smooth profiles P(0)(z) and Z(0)(z), until a steady state is reached. The initial profiles P(0)(z) and Z(0)(z) are identical to the solutions obtained in the previous section (Fig. 2.2a,b) except that in the non homogeneous layers, that is, the mid layer in the case of Pand the upper layer in the case of Z, the z-dependent profiles are replaced by transition cosine functions. Since ∂cos(z)/∂z =−sin(z) = 0 at z={0, π},P(0) and Z(0) have continuous (zero) derivatives at z1and z3, and at 0 and z1, respectively. Specifically, we define the initial profiles P(0)(z) =          P1, z ∈(z1,0] 1 2P1h1 + cos (z−z1)π z1−z3i, z ∈[z3, z1] 0, z < z3, (2.19) and Z(0)(z) =    1 2Z1³1 + cos zπ z1´, z ∈[z1,0] 0, z < z1. (2.20) These initial profiles are shown in Fig. 2.4. Next, P(0)(z) and Z(0)(z) are vertically discretisized and integrated forward using (2.1)-(2.2) with u= 0. Nine different numerical resolutions, ranging from constant grid-size δz = 4 m (i= 1) to δz = 1.5625 cm (i= 9), are implemented (Table 2.2). 0 2 4 6 8 P(z) −40 −35 −30 −25 −20 −15 −10 z (m) (a) P(0) P(1) z3 z2 z1 −0.1 0.0 0.1 0.2 0.3 Z(z) −20 −15 −10 −5 0(b) Z(0) Z(1) z1 Figure 2.4: (a) P(z), and (b) Z(z). The initial profiles are P(0) and Z(0). The steady solutions for the different numerical resolutions are P(i)and Z(i)(i={1,...,9}). Close up views of P(i)at z∼z2, and of Z(i)at z∼z1are shown in Figs. 2.5a and 2.6a, respectively. Units are mmol N m−3. 2.3. STEADY SOLUTIONS 23 label grid points resolution (m) i ni= 25 ∗i+ 1 |Zmin|/(ni−1) 1 26 4 2 51 2 3 101 1 4 201 0.5 5 401 0.25 6 801 0.125 7 1601 0.0625 8 3201 0.03125 9 6401 0.015625 Table 2.2: List of numerical resolutions. Here Zmin =−100 m. Time integration is carried out using an explicit leap-frog scheme, together with a RobertAsselin time filter to avoid the computational mode (see e.g., Durran,1998, p. 62). The time integration was 104days, at the end of which the maximum forcing term in the local rate of change of Por Zwas of the order of 10−7mmol N m−3d−1. During the integration time Pand Zmonotonically converged to the steady solutions P(i)and Z(i)(i= 1,...,9) shown in Fig. 2.4. After a first look these solutions seem to be very similar to the non continuously differentiable solutions obtained in the previous section. However, a closer view around the layer boundary depths (zooms on z=z2and z=z3in Figs. 2.5a and 2.6a, respectively) reveals that for large resolutions (i≥7) the solutions P(i)and Z(i)become continuously differentiable functions of z. From the numerical point of view this means that the vertical gradients ∂P(i)/∂z and ∂Z(i)/∂z, here computed using a simple second order centered scheme, have converged to finite values and no longer depend on the numerical resolution (Figs. 2.5b and 2.6b). −0.2 0.0 0.2 0.4 0.6 0.8 1.0 P(z) (mmol N m −3) −25.4 −25.3 −25.2 −25.1 −25.0 (a) z (m) z2 3 4 5 6 7, 8, 9 −10123456 dP/dz (mmol N m −4) −28 −27 −26 −25 −24 −23 −22 (b) 123 4 56 7, 8, 9 Figure 2.5: A close up view at z∼z2of (a) P(i)(z), and (b) ∂P(i)/∂z. 30 CHAPTER 2 −1 −0.5 0 0.5 1 x (km) −1 −0.5 0 0.5 1 y (km) −5.0 −4.2 −3.3 −2.5 −1.7 −0.8 0.0 0.8 1.7 2.5 3.3 4.2 5.0 (a) P’ (mmol N m −3) −1 −0.5 0 0.5 1 x (km) −1 −0.5 0 0.5 1 −11.3 −9.4 −7.5 −5.7 −3.8 −1.9 0.0 1.9 3.8 5.7 7.5 9.4 11.3 (c) Z’ (mmol N m −3) (x10 2) −1 −0.5 0 0.5 1 y (km) −50 −25 0 z (m) (b) −1 −0.5 0 0.5 1 y (km) −50 −25 0 (d) Figure 2.10: (a) P′(x, y) at zb≃ −23.4 m (P′∈[−3.3,3.1] mmol N m−3). PV contours =±0.05 at z= 0 are included for reference. (b) P′(y, z) on vertical section x= 0 (P′∈ [−3.3,4.8] mmol N m−3). (c) Z′(x, y) at za=−12.5 m (Z′∈[−9.0,8.3] ×10−2mmol N m−3). (d) Z′(y, z) on vertical section x= 0 (Z′∈[−8.9,11.2] ×10−2mmol N m−3). Time t= 8 Tip. the total Pand Zchanges (the rhs of (2.1) and (2.2), respectively) forced by wand are independent of the effect of mere advection. Ascending fluid particles experience an increase of their Pcontent while descending particles experience a decrease of P. Note particularly the large Pdecrease at the northern side of the domain (Fig. 2.11a,b), where w < 0 (Fig. 2.8c,d). Positive Pbudgets occur at depths a bit shallower than negative Pbudgets, which explains why dP/dtis mostly negative at the depth shown in Fig. 2.11a. The material rate of change of Z(Fig. 2.11c,d) and ware also clearly correlated. However dZ/dtusually displays a minimum and a maximum along the water column, which is consistent with the two maxima in σ{dZ/dt}(Fig. 2.9d). The analysis of dP/dtand dZ/dtinto their local and advective changes (Fig. 2.12) shows that, as inferred from their standard deviations (Fig. 2.9c,d), there is a large cancelation between the local change and the horizontal advection of Pand Z. The vertical advection is smaller. Consistently also with the P′and Z′distributions (Fig. 2.10) the local change and horizontal advection of Zpresent patterns more elongated than those of P. This is a 2.4. COUPLED PHYSICAL-ECOSYSTEM NUMERICAL SIMULATIONS 31 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 y (km) −5.1 −4.3 −3.4 −2.6 −1.7 −0.9 0.0 0.9 1.7 2.6 3.4 4.3 5.1 (a) z=−18.7 m (mmol N m −3 d −1) (x10) −1 −0.5 0 0.5 1 x (km) −1 −0.5 0 0.5 1 y (km) (b) z=−23.4 m −1 −0.5 0 0.5 1 y (km) −50 −25 0 z (m) (c) Figure 2.11: Pproduction anomaly G′ P(x, y) on horizontal planes (a) z=−18.75 m (iz= 53, G′ P∈[−0.12,0.33] mmol N m−3d−1). (b) z=zb≃ −23.4 m (iz= 50, G′ P∈ [−0.51,0.035] mmol N m−3d−1). (c) G′ P(y, z) on vertical section x= 0 (G′ P∈ [−0.51,0.34] mmol N m−3d−1). Time t= 8Tip. consequence of the better material conservation of Zin comparison con P. Large local rates occur in the frontal areas, where both horizontal velocity and horizontal gradients of Pand Zare large. The vertical advection of Pand Zhave however similar patterns. This is so because ∂P/∂z ≃∂Ps/∂z > 0 at z=zband ∂Z/∂z ≃∂Zs/∂z > 0 at z=za, so that the vertical advection patterns (Fig. 2.12c,f) resemble the wpattern (Fig. 2.8c). The time evolution of σ{dP/dt},σ{dZ/ dt}, and σ{w}(Fig. 2.13) show that the ecosystem time response to wmaxima is about 5Tip. The second wmaximum at t≃37 Tip is related to the flow enhancement due to the fusion of two anticyclones. The vertical resolution used in this simulation (δz = 65 cm) is not good enough to fully resolve the large vertical gradients of Pand Zat transition depths. Based on Figs. 2.5b and 2.6b, vertical gradients are underestimated by a 50%. Larger vertical resolutions would correctly resolve the vertical advection of Pand Zwhich would cause an important increase in P′and Z′distributions. However, as another consequence of an increased vertical resolution, these larger anomalies would be restricted to thinner ocean layers, so that only quantitative changes are expected in the ecosystem variables. 32 CHAPTER 2 These numerical results correspond to a non diffusive NPZ ecosystem model coupled to an adiabatic inviscid physical model. These results will not apply when vertical mixing is added to the NPZ model (see Edwards et al.,2000) since in such a case the large vertical NPZ gradients found here would turn Psand Zsinto unsteady solutions. We note that the mere existence of vertical eddy diffusion in a numerical model already requires vertical differentiability. Vertical diffusion is not included here because the NPZ model is kept as simple as possible in order to analyze the vertical velocity forcing of NPZ anomalies. Including vertical diffusion will add new free parameters (the vertical diffusivity coefficients) to the already large list of NPZ parameters on Table 2.1. Furthermore, the absence of NPZ diffusion This is also consistent with the inviscid nature of the PV-conserving dynamical model (only a very small amount of numerical diffusivity is included to avoid grid-size noise). −1 −0.5 0 0.5 1 x (km) −1 −0.5 0 0.5 1 y (km) −8.2 −6.8 −5.5 −4.1 −2.7 −1.4 0.0 1.4 2.7 4.1 5.5 6.8 8.2 (a) dP dt (mmol N m −3 d −1) (x10) −1 −0.5 0 0.5 1 x (km) −1 −0.5 0 0.5 1 −2.0 −1.7 −1.3 −1.0 −0.7 −0.3 0.0 0.3 0.7 1.0 1.3 1.7 2.0 (c) dZ dt (mmol N m −3 d −1) (x10 4) −1 −0.5 0 0.5 1 y (km) −50 −25 0 z (m) (b) −1 −0.5 0 0.5 1 y (km) −50 −25 0 (d) Figure 2.12: (a) dP/dtat zb≃ −23.4 m (iz= 50,dP/dt∈[−0.82,0.081]). (b) dP/dtat x= 0 (dP/dt∈[−0.82,0.42]). (c) dZ/dtat z=−12.5 m (iz= 57,dZ/dt∈[−19.8,5.3] ×10−5). (d) dZ/dtat x= 0 (dZ/dt∈[−19.6,10.5] ×10−5). Time t= 8 Tip. 2.4. COUPLED PHYSICAL-ECOSYSTEM NUMERICAL SIMULATIONS 33 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 y (km) −9.5 −7.9 −6.3 −4.8 −3.2 −1.6 0.0 1.6 3.2 4.8 6.3 7.9 9.5 (a) (mmol N m −3 d −1) ∂P ∂t −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 −5.2 −4.3 −3.4 −2.6 −1.7 −0.9 0.0 0.9 1.7 2.6 3.4 4.3 5.2 (d) (mmol N m −3 d −1) (x10) ∂Z ∂t −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 y (km) (b) uh·∇hP −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 (e) uh·∇hZ −1 −0.5 0 0.5 1 x (km) −1 −0.5 0 0.5 1 y (km) −4.2 −3.5 −2.8 −2.1 −1.4 −0.7 0.0 0.7 1.4 2.1 2.8 3.5 4.2 (c) (mmol N m −3 d −1) w∂P ∂z −1 −0.5 0 0.5 1 x (km) −1 −0.5 0 0.5 1 −14.0 −11.6 −9.3 −7.0 −4.7 −2.3 0.0 2.3 4.7 7.0 9.3 11.6 14.0 (f) (mmol N m −3 d −1) (x10 2) w∂Z ∂z Figure 2.13: (a) ∂P/∂t at z≃zb=−23.4 m (∂P/∂t ∈[−8.3,7.4]). (b) uh·∇ hP(∈ [−9.4,8.3]). (c) w∂P/∂z (∈[−1.8,4.2]). (d) ∂Z/∂t at za=−12.5 m (∂Z/∂t ∈[−0.51,0.30]). (e) uh·∇ hZ(∈[−0.32,0.53]). (f) w∂Z/∂z (∈[−0.042,0.14]). Time t= 8 Tip. 34 CHAPTER 2 dP/dtdZ/dt w 0 10 20 30 40 50 −50 −40 −30 −20 −10 0 z (a) (T ) ip (m) time 0 10 20 30 40 50 −50 −40 −30 −20 −10 0 (b) (T ) ip time 0 10 20 30 40 50 −50 −40 −30 −20 −10 0 (c) (T ) ip time Figure 2.14: (a) σ{dP/dt}(z, t) (max = 0.11, ∆ = 0.01). (b) σ{dZ/dt}(z, t) (max = 5×10−6, ∆ = 5.4×10−5). Units are mmol N m−3d−1. (c) σ{w}(z, t) (max = 1.2×10−3, ∆ = 10−4). 2.5 Concluding Remarks We have first shown that 1D steady and continuously differentiable (in a numerical sense) solutions to the NPZ equations are possible. These solutions are potentially useful as initial steady ecosystem conditions to investigate the role of horizontal and vertical advection in 3D coupled physical-ecosystem numerical models. An example of ecological development due to vertical velocity enhancement during a baroclinic instability process has been presented. This example shows that once phytoplankton and zooplankton anomalies develop locally forced by balanced vertical velocity they are horizontally advected away from the upwelling or downwelling regions so that spatial distributions of vertical velocity and ecological fields become eventually uncorrelated (for experimental evidence of this process see e.g. Ruiz et al., 2001). Thus the biological distributions are more related to PV gradients than to PV itself. This fact, and the submesoscale vertical origin of NPZ anomalies, is consistent with (L´evy et al.,2001), who used a primitive equations model with vertical diffusion. However, the experimental work of (L´evy et al.,2005) questions the contribution of submesoscale total advection in the phytoplankton variability over large time scales. The physical-ecological model used here has several limitations. On the one hand, these NPZ solutions require very good vertical resolutions, with a grid size of few centimeters, to be properly discretisized. This imposes a severe handicap to the available random access memory of current computers running 3D coupled physical-ecosystem models. Though from a strict numerical perspective this fact is a serious modeling limitation, from a wider 2.5. CONCLUDING REMARKS 35 perspective other handicaps, for instance, errors in the mathematical parametrization of the different NPZ processes, are likely to be of larger relevance (see Anderson,2005;Flynn, 2005;Mitra et al.,2007). It is nevertheless important to know the degree at which ecosystem modeling solutions faithfully reproduce the underlaying ecosystem dynamics and that, even with poor vertical resolution, it is possible to obtain good qualitative results from these models. On the other hand, our results show that Pand Zapproximately behave as passive tracers while organisms are in fact active tracers. This is so because the simple initial steady NPZ profiles let little interaction between Pand Z. In future work we will address these interactions using both more realistic initial NPZ profiles and a more complex biological model. To conclude, these results are a first approximation towards a better understanding of biological processes forced by vertical velocity at the submesoscale. Many questions still remain to be answered in this context. How much does the submesoscale vertical velocity contribute to primary and secondary productivity in comparison to the mesoscale? Does the vertical advection induce different biological patterns in eutrophic and oligotrophic regimes? Is the submesoscale important in the seasonal biological variance? Chapter 3 Phytoplankton enhancement by oceanic dipoles Chapter submitted as M.Claret, A.Viúdez, and Y. H. Spitz, 2012: Phytoplankton enhancement by oceanic dipoles. Al cim d’un promontori que domina les ones de la mar, quan l’astre rei cap a ponent declina me’n pujo a meditar. Amb la claror d’aqueixa ll`antia encesa contemplo mon no-res; contemplo el mar i el cel, i llur grandesa m’aixafa com un pes. Vora la mar, Jacint Verdaguer 39 ABSTRACT Phytoplankton distributions associated to oceanic dipoles are investigated using a highresolution, non-hydrostatic, three-dimensional numerical model coupled to a NPZ (NutrientPhytoplankton-Zooplankton) oligotrophic model. Two scenarios are considered, a submesoscale surface and a mesoscale subsurface dipoles, in order to characterize both the effect of vortices in translation and the distant action of potential vorticity (PV) on phytoplankton dynamics. We observe that the dipole separatrix acts as an impermeable barrier dividing two different ecological niches. On the one hand, plankton trapped within vortices reaches a near steady state, which depends on the community entrained at the formation stage of vortices. On the other hand, phytoplankton advection by the dipole introduces ecological heterogeneity outside the separatrix. A balance between vertical advection (VA) and biological forcing generates a subsurface trail of phytoplankton at the wake of the translating dipole. The length of this trail is nearly constant and may be used to estimate the phytoplankton mortality rate for a known dipole speed. When negative gradients of PV exist in the upper ocean, horizontal advection accounts indirectly for a larger phytoplankton increase, though much localized, than that caused by VA. If a subsurface dipole is considered, an ecosystem patch is stirred such that a filament runs along its axis. The filament extent increases with depth and as a result phytoplankton self-shading decreases at the filament front in benefit of light irradiance, and thus of phytoplankton production. Finally, the horizontal extension of the filament is approximated to an analytical expression derived for a dipole of known geometry and intensity. 46 CHAPTER 3 since the dynamical ABmodel is adiabatic so that changes in water density due to solar radiation are ignored. Phytoplankton is removed by zooplankton grazing R(x, t), modeled with an Ivlev response (Parsons et al.,1967), and a linear term (MP) that represents death, exudation, excretion or other processes. A portion (1 −Γ0) of the Pgrazed is converted to Zproduction (GZ), which is balanced by a density dependent death rate (MZ). A quadratic closure term for Zhas a double purpose. On the one hand, it eliminates the oscillatory behavior of the NPZ model with u= 0 that appears when a linear closure term is instead considered (Edwards and Yool,2000). On the other hand, it introduces the effect of cannibalism (Pitchford and Brindley,1998;Ohman et al.,2002), and predation by higher trophic levels (Ohman and Hirche,2001). Finally, detrital matter (Γ0R+MP+MZ) is remineralized and nutrients are again available for Puptake. Therefore, the total nitrogen NT≡P+Z+Nis materially conserved, dNT dt= 0 .(3.11) The above equation lets us to define the original NPZ equations (3.5)–(3.7) depending on only two variables, Pand Z, which satisfy (3.5), (3.6), and (3.10). The ecological model is coupled to the physical one such that the latter provides the 3D velocity field to advect the variables of the former. Additionally, diffusion processes are neglected so that PV is materially conserved. 3.2.3 Numerical parameters We work in the quasigeostrophic (QG) space, where vertical dimension is stretched by c. Thus the numerical domain is isotropic with vertical extent LZ= 2π(which defines the unit of space) and horizontal extents LX=LY=cLZ. Triply-periodicity is also imposed but only in the dependent variable ϕ(x, t). For example, though D(x, t) is triply-periodic the total density field ρ(x, t) is not. The number of grid points is (nX, nY, nZ) = (128,128,128), and the number of isopycnals nL= 128. Equations (3.4), (3.5), and (3.6) are integrated forward in time using an explicit leapfrog scheme combined with a Robert-Asselin time filter to avoid the computational mode. Time-step is δt = 7×10−4d in all the cases considered. Finally, a biharmonic hyperdiffusion operator µ∇4 qfor Ah,P, and Zis added to their respective diagnostic equations in order 3.2. PHYSICAL-ECOLOGICAL COUPLED MODEL AND IMPLEMENTATION 47 to dump the amplitude of grid-size scale noise due to spatial discretization on a fixed grid. Above, ∇q≡c∇h+k∂zis the gradient operator in the QG space, and the hyperviscosity coefficient µis defined by specifying the e-folding time of the largest wave number in spectral space per inertial period, which corresponds to about 0.7 d for a mean latitude of 45o. In the AB-model ef= 100, while in the NPZ model ef= 10. 3.2.4 Initial conditions The AB-model is initialized using the so called PV initialization approach (Vi´udez and Dritschel,2003), which is unique to the PV conserving algorithm used in this numerical model. It consists in a gradual increase of PV in every fluid particle until a prescribed value is reached. This initialization technique largely avoids the generation of inertia–gravity waves due to the initial imbalance between density and velocity fields, which otherwise could contaminate the balanced vertical velocity. In the dipole cases here considered, an initialization time period of about 3.5 d is sufficient to avoid the appearance of the imbalance. The NPZ model is initialized with steady state solutions in order to isolate the processes that control plankton dynamics when the physical system is perturbed by an oceanic dipole. Stationary stable ecological solutions ˆ N(z),ˆ P(z), and ˆ Z(z) are obtained numerically by time integration of the NPZ model (3.5)–(3.7) in the state of rest (u= 0) with the following initial {NT, P, Z}profiles NT(z) =    T1cos h³z z1−1´2π 5i, z ∈[z1,0], 1, z ∈[zmin, z1), (3.12) (P, Z)(z) =    (P1, Z1) sin h³z z2−1´π 2i, z ∈[z2,0], 0, z ∈[zmin, z2), (3.13) where T1= 1 mmol N m−3,P1= 0.22 mmol N m−3, and Z1= 0.1 mmol N m−3. In contrast, the values of zmin,z1, and z2differ whether a surface or a subsurface dipole is considered. Independently of the {zmin, z1, z2}choice, time integration of (??)-(??) converges to a steady state (Fig. 3.1). The resulting stationary profiles {ˆ P, ˆ Z, ˆ N}, where ˆ Nis recovered by inversion of (3.10), are used to initialize the ecological model in the AB-NPZ coupled cases. Both models are coupled once the PV field is fully initialized at t=tisince, in the 48 CHAPTER 3 PV initialization approach, PV is not materially conserved during the initialization period (t < ti). Note that at initial time, t0= 0, isopycnals are flat and therefore d(x, t0) = z since D(x, t0) = 0. However, at the end of the initialization period doming of isopycnals implies D(x, ti)6= 0. As a result, two different ecological initial conditions can be further distinguished depending on whether NTis considered homogenous on isopycnals (constant d) or on horizontal layers (constant z). In this work, we have considered the former initialization when zP∼ =−61 m (section 3.3.1) and the latter when zP∼ =−90 m (section 3.3.2). 3.3 Numerical simulations The effects of both vortex translation and the distant PV action on ecological dynamics are investigated considering two oligotrophic scenarios typical of the open ocean. We first focus on how a surface dipole perturbs a stationary ecosystem at the submesoscale (section 3.3.1), where high-resolution remote sensing observations have recently unveiled many vortical structures (Munk et al.,2000). In this case phytoplankton subsurface maximum ˆ Pmax ≡ˆ P(zP) is initially placed underneath the dipole, so that zP≃ −61 m (Fig. 3.1a), by choosing (z1, z2, kw) = (−120 m,−150 m,0.03 m−1). Secondly, the distant influence of a subsurface dipole on a plankton 3D patch is investigated at the mesoscale (section 3.3.2), the scale at which subsurface long-lived vortices, such as meddies, are widespread in the ocean (Richardson et al.,2000). Since the shallowest dipole edge is at z∼ =−200 m, we have deepened ˆ Pmax to zP≃ −90 m (Fig. 3.1b) by setting (z1, z2, kw) = (−290 m,−187 m,0.023 m−1). 3.3.1 Submesoscale surface dipole The submesoscale is here introduced as c= 50, and only half of the vertical domain LZis considered with zmin =−150 m. Thus the horizontal domain is 15 km. The vortex dipole consists of a baroclinic cyclone (+) and anticyclone (−) with PV anomalies ± max = 0.5 (∼ =0.35 d−1). The maximum length of the horizontal semi-axes of the ellipsoids of constant PV in both vortices are aX∼ =1.9 km and aY∼ =2.9 km. The length of the vertical semiaxes a± Zare different in the initial configuration with flat isopycnals, being a+ Z∼ =62 m and a− Z∼ =52.5 m. However, during the initialization time the isopycnals stretch (shrink) in the anticyclone (cyclone), so that at ti∼ =3.5 d the vortices have similar vertical extent and the dipole describes a straight trajectory along the x-axis (Dubosq and Vi´udez,2007). At 3.3. NUMERICAL SIMULATIONS 49 0.0 0.2 0.4 0.6 0.8 1.0 −140 −120 −100 −80 −60 −40 −20 0 mmol N dm−3 (m)d N ^ P ^ Z ^NT ^ (a) 0.0 0.2 0.4 0.6 0.8 1.0 −300 −250 −200 −150 −100 −50 0 mmol N dm−3 (m)z N ^ P ^ Z ^ NT ^ (b) Figure 3.1: Vertical profiles of total nitrogen NT(z), and its associated stationary and stable profiles of {P, Z, N}obtained with the parameters given in Table 3.1 and without physical forcing. In the surface dipole case the profiles corresponding to (a) are initialized homogeneous on isopycnals (d), while in the subsurface dipole scenario those of (b) are homogeneous on horizontal levels (z). the end of the initialization time t=ti, the horizontal speed contours correspond to those of concentric deformed double tori intersecting an horizontal plane with a maximum speed |u|max = 4.76 cm s−1along the dipole axis (Figs. 3.2a,b). We note that submesoscale dipoles, as long as they remain isolated and do not interact with other submesoscale vortices, have a quadrupolar pattern of w(Fig. 3.2c) similar to mesoscale dipoles (Pall`as-Sanz and Vi´udez, 2007). The wmaximum absolute value |w|max = 71.23 cm d−1is placed at z=−35 m (iZ= 50, Fig. 3.2d). As the dipole moves forward, the isopycnals are displaced upwards (downwards) at the front (rear) of the cyclone. The opposite changes occur in the anticyclone. The dipole flow remains always statically and inertially stable since the Rossby number R≡ωh/N, where the squared total Brunt-V¨ais¨al¨a frequency N2(x, t)≡ −gα0∂ρ/∂z = N2[1−∂D/∂z(x, t)], minimum is Rmin =−0.43, and the Froude number F≡ζ/f maximum Fmax = 0.19. The NPZ model is coupled to the above described dynamical conditions at t=ti. In this case, total nitrogen is set homogeneous on isopycnal levels by projecting the initial profiles 50 CHAPTER 3 −7.5 −3.25 0 3.25 7.5 −7.5 −3.25 0 3.25 (km)y 0.0 0.8 1.6 2.4 3.2 4.0 4.8 | h| (x10−2 m s−1)u (a) −7.5 −3.25 0 3.25 7.5 −7.5 −3.25 0 3.25 (c) −7.5 −3.25 0 3.25 7.5 −100 −50 0 (km)x (b) −7.5 −3.25 0 3.25 7.5 −100 −50 0 (km)x (d) Figure 3.2: (a) Horizontal distribution of |uh|at z= 0 (iZ= 65). Only every other vector is plotted. A straight solid line is drawn where the vertical distribution of vis shown in (b) at y∼ =−2.6 km (iY= 43, v∈[−4.76,2.76] cm s−1, contour interval δv ≃0.4 cm s−1). Horizontal (c) and vertical (d) distributions of w(contour line range is |w|<71.23 cm d−1 with δw ≃1 cm d−1) at z∼ =−33 m (iZ= 50) and at y∼ =−0.82 km (iY= 58), respectively. Distributions correspond to the initial time ti∼ =3.5d. PV contours =±0.2 (∼ =0.28 d−1) are included for reference. Hereinafter, solid and dashed contours are used for positive values and negative values, respectively. We observe that the dipole translates southwards as a solid body since uhhas a surface maximum along the dipole axis. In contrast, extreme wbounds are reached at the subsurface. (Fig. 3.1a) on isopycnals as NT(x, t) = ˆ NT(d(x, t)) .(3.14) Note that this function does not depend explicitly on tbecause ˆ NTis materially conserved. In accordance to this NTinitialization, {P, Z, N}variables are also set constant on isopycnals at t=tisatisfying P(x, ti) = ˆ P(d(x, ti)) ,(3.15) with similar relations for Z(x, ti) and N(x, ti). Thus the initial distributions of these variables depend on their vertical gradient and the isopycnal vertical displacement D. For instance, Pshows a quadrupolar pattern at t=ti(Fig.3.3a) because Pihas a subsurface maximum at zP∼ =−61 m, and thus a dipolar pattern is observed within each vortex. This also explains 3.3. NUMERICAL SIMULATIONS 51 −7.5 −3.25 0 3.25 7.5 −150 −100 −50 0 (km) (m) x z 0.94 0.95 0.96 0.97 0.98 0.99 1.00 1.01 1.02 1.03 1.04 1.05 1.06 −7.5 −3.25 0 3.25 7.5 −150 −100 −50 0 (a) P’ D t=tip −7.5 −3.25 0 3.25 7.5 −150 −100 −50 0 (km) (m)z −7.5 −3.25 0 3.25 7.5 −150 −100 −50 0 (b) Z’ t=tip −7.5 −3.25 0 3.25 7.5 −150 −100 −50 0 (km) (m)z −7.5 −3.25 0 3.25 7.5 −150 −100 −50 0 (c) t=24.7 d −7.5 −3.25 0 3.25 7.5 −150 −100 −50 0 (km) (m) y z −7.5 −3.25 0 3.25 7.5 −150 −100 −50 0 (d) t=88.4 d Figure 3.3: (a) Vertical section crossing the vortex centers, at y∼ =−2.6km, (iY= 43), of P′(x, z) in grey scale and D(|D|<4.3 m, contour interval δD≃0.6 m) after PV initialization at t=ti. (b, c, d) Vertical sections parallel to the dipole axis of P′(y, z) and Z′(y, z) in contour lines (|Z′|<1.06, δZ′= 1%) at different times and (b) x∼ =2.1 km (iX= 83), (c) x∼ =2.3 km (iX= 85), and (d) x∼ =3.5 km (iX= 95). Phytoplankton subsurface maximum zP(straight line in (a) ) and contours of =±0.2 are indicated for reference. Note that P′increases about 2% at the cyclone wake. 52 CHAPTER 3 the quadrupolar initial distribution of Z, though of different magnitude (not shown). In contrast, Nhas a dipolar pattern (not shown) since Niincreases gradually with depth. The initial fields of P, and Zare advected with the 3D velocity at times t > ti, while N is obtained using (3.5), (3.6), and (3.10). In order to quantify the dipole impact on the ecosystem, ecological properties are adimensionalized such that they are normalized by a value outside the dipole where isopycnals are flat, that is, by the stationary solution at the isopycnal depth in the reference density configuration (d). So that, any ecological variable χvis characterized through the following enhancement factor, χ′ v(x, t)≡χv(x, t) ˆχv(d(x, t)) ,(3.16) while any ecological rate χris refered to the steady state Pproduction (GP) as χ′ r(x, t)≡χr(x, t) ˆ DP(d(x, t)) .(3.17) The P′maximum is located at the center of the cyclone, where isopycnal displacement is the largest, and moves with the translating dipole (Figs. 3.3b–d). It represents about a 6% increase referred to its steady state ˆ P(d). Zooplankton enhancement Z′is of the same order as P′and its distribution is positively correlated with that of P′. Based on the biological forcing balance (Eq. 3.5), we found that zooplankton grazing is much smaller than the other terms. In light of that, we shall focus on the phytoplankton response only. Two different phytoplankton dynamics are observed depending on a critical depth zc≃ −80.8 m, which corresponds to the maximum depth of the dipole. Above zc, the physical forcing dominates over the ecological forcing and Pbehaves mainly as a passive tracer. In contrast, below zc vertical advection triggers an ecological response of similar magnitude to physical terms and aP′trail develops at the dipole wake. Above zc, the dipole separatrix, which is defined by the minimum isosurface, acts as a dynamical barrier. On the one hand, the horizontal advection of Pis greater than its vertical advection (Fig.3.4a). On the other hand, vertical uplift of Pto lighter levels triggers aPincrease, however it is much smaller than the horizontal advection (not shown). Thus Pis translated at the dipole phase speed. Furthermore, Pcorrelates with D(Fig.3.4b), and ultimately with through (??). This fact suggests that stationary stable ecological solutions are reached within vortices, which are close to their initial conditions at tiwhen 3.3. NUMERICAL SIMULATIONS 53 −7.5 −3.25 0 3.25 7.5 −7.5 −3.25 0 3.25 7.5 (km)x (km)y −4.38 −3.12 −1.88 −0.62 1.25 2.50 3.75 −7.5 −3.25 0 3.25 7.5 −7.5 −3.25 0 3.25 7.5 (a) t=88.4 d (w P z)’ (%) u h hP( )’ −7.5 −3.25 0 3.25 7.5 −7.5 −3.25 0 3.25 7.5 (km)x 0.94 0.95 0.97 0.98 1.00 1.02 1.03 1.05 1.06 −7.5 −3.25 0 3.25 7.5 −7.5 −3.25 0 3.25 7.5 (b) P’ D Figure 3.4: Horizontal section at z≃ −34.6 m (iZ= 50) of (a) (w ∂P/∂z)′in grey scale and (uh· ∇hP)′in contour lines (|∆|<0.054, δ= 6.25 ×10−3), and (b) P′in grey scale and Din contour lines (D∈(−4.15,3.75) m, δD= 0.5 m). PV contours ±=±0.2 at z= 0 (thickest line) are included for reference. After δt = 88 d, P′correlates with Dwithin vortices as in the initial configuration at time ti. 20 40 60 80 0.0 0.5 1.0 1.5 2.0 2.5 -0.2 0.0 0.2 0.4 0.6 < P’>-1 (%) (x10-1 d-1)< ϖ > <( h h P)’> <(w P z)’> (%) u<( P/ t)’> d d days Figure 3.5: Time evolution of hP′i(solid line), hi(dash dotted line), h(dP/dt)′i(dotted line), h(uh·∇hP)′i(dash triple dotted line) and h(w ∂P/∂z)′i(dashed line), averaged over the cyclone volume. After an initial adjustment, hP′ireaches a steady state close to the initial configuration. 54 CHAPTER 3 NTis initialized homogeneous on isopycnals and so ecological variables. In order to unveil this steady state, any property χis integrated on the cyclone volume as hχi(t) = 1 nZΩ χ(x, t) dV,(3.18) where the spatial domain Ω comprises those grid points nwith  > 0. As expected, after an initial adjustment of about δt ∼ =20 d, hP′iremains nearly constant oscillating around a time averaged value of 1.9% (Fig. 3.5). The oscillation in hP′ihas a period T≃5.6 d different from that of advective terms. Actually, hP′iis in phase with hi, which oscillates due to vortex Rossby waves (VRWs) (Rodr´ıguez-Marroyo and Vi´udez,2009). VRWs cause a periodic azimuthal oscillation of the dipole geometry, relative to its time averaged configuration, with a period of 8 inertial periods, so that close to Tfor a midlatitude, and thus also alter Pdistributions. A distinct plankton dynamics occurs below the dipole at z < zc. In this case, P′increases at the cyclone wake as a combination of physical and biological processes, which have similar magnitude. Firstly, (w∂P/∂z)′is negatively correlated with P′−1 at the dipole rear (t= 7 d, Fig. 3.6a). This suggests that P′increases (decreases) due to subduction of shallower depths richer (poorer) in P′at the anticyclone (cyclone) wake. As dipole translates, the fluid parcels that experience these P′changes decelerate relative to vortices and reach a state of rest at later times. However, P′does not inmediately converge to its stationary initial condition (P′= 1) since the ecosystem responds with some inertia to a perturbation. As a result, changes in P′retain some memory of their causing mechanism and thus P′correlates in this case with the time integrated vertical advection R(w ∂P/∂z)′dt(Fig. 3.6b, t= 10.6 d). Later on, after δt ∼ =15.5 d of the increase in P′, biological forcing (dP/dt)′reaches the order of magnitude of (w∂P/∂z)′(Fig. 3.6c). Since ecological and advective terms are opposite, we conclude that P′is induced at the dipole wake due to a balance both and thus it becomes correlated with the time integrated local rate R(∂P/∂t)′dt(Fig. 3.6d, t= 24.7 d). The long term evolution of this P′trail is quantified by volume integration of properties using (3.11). In this case, the spatial domain Ω is enclosed within isosurface P′= 1.002 at the dipole wake. We observe that hP′ireaches a nearly stationary stable state after an initial adjustment of δt ∼ =35 d, however h(dP/dt)′i ≡ h(∂P/∂t)′i+h(u·∇P)′i+h(w∂P/∂z)′i 6= 0 (Fig. 3.7). In order to unveil what we are missing in the abovementioned balance, we take a Lagrangian approach. To this end, particles that experience a smaller P′increase than that 3.3. NUMERICAL SIMULATIONS 55 P’ 0.93 0.94 0.95 0.96 0.97 0.98 0.99 1.00 1.01 1.02 1.03 1.04 1.05 1.06 1.07 −7.5 −3.25 0 3.25 7.5 −7.5 −3.25 0 3.25 7.5 (km)y −7.5 −3.25 0 3.25 7.5 −7.5 −3.25 0 3.25 7.5 (a) t=7 d −7.5 −3.25 0 3.25 7.5 −7.5 −3.25 0 3.25 7.5 −7.5 −3.25 0 3.25 7.5 −7.5 −3.25 0 3.25 7.5 (b) t=10.6 d −7.5 −3.25 0 3.25 7.5 −7.5 −3.25 0 3.25 7.5 (km)x (km)y −7.5 −3.25 0 3.25 7.5 −7.5 −3.25 0 3.25 7.5 (c) t=15.5 d −7.5 −3.25 0 3.25 7.5 −7.5 −3.25 0 3.25 7.5 (km)x −7.5 −3.25 0 3.25 7.5 −7.5 −3.25 0 3.25 7.5 (d) t=24.7 d Figure 3.6: Horizontal sections of P′in grey-scale and contour lines of (a) (w ∂P/∂z)′ with maximum absolute value |∆|max ∼ =0.016 and contour interval δ∼ =3×10−3, (b) R(w ∂P/∂z)′dtwith |∆|max ∼ =0.24 and δ= 0.025, (c) (dP/dt)′with |∆|max = 0.019 and δ∼ =3×10−3, and (d) R(∂P/∂t)′dtwith |∆|max = 0.32 and δ∼ =0.04. Section layers are z≃ −71.5 m (iZ= 34) in (a) and z≃ −85.4 m (iZ= 28) in (b–d). Finally, PV contours = 0.2 are included for reference. Note that the P′increase results from a balance between vertical advection and biological forcing. 20 40 60 80 1.000 1.002 1.004 1.006 1.008 1.010 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 < P’> <( h h P)’> <(w P z)’> (%) u<( P/ t)’> d d days Figure 3.7: Time evolution of hP′i(solid line), h(dP/dt)′i (dotted line), h(w ∂P/∂z)′i (dashed line), and h(uh·∇hP)′i (dash dotted line) averaged within isosurface P′= 1.002 below the deepest PV edge of the dipole. Note that hP′i reaches a nearly steady state, though ecological terms do not balance out. 62 CHAPTER 3 center, respectively, is δt =Zr ri dr′ vg e(r′, z)(3.23) =γ0 3[R3(r, z)−R3 i(ri, z)] + γ0z2[R(r, z)−Ri(ri, z)] +γ0z3ln ¯¯¯¯ [R(r, z)−z][Ri(ri, z) + z] [R(r, z) + z][Ri(ri, z)−z]¯¯¯¯, where constant γ0= 3/(2N0R3 0). Since all the terms above are strictly increasing continuous functions of r,r(z, ti+δt) is easily obtained by numerical inversion. This approximation predicts a maximum separatrix depth zAclose to the observed zS(Fig.3.12b). Above zA the predicted NTfront translation is about half the observed one, a discrepancy that may be explained by the absence of non-linear terms in the QG PV dynamics. Below zA, the front line gradually converges to the separatrix curvature with a maximum displacement at z=−500 m close to twice the vortices radius. Though the differences observed, equation (3.15) may be used to obtain a first estimate of the length of the patch filamentation at any depth for a given dipole geometry and intensity. Finally, the ecological dynamics inside this NTpatch is investigated. To this end, ecological properties are quantified through the enhancement factors (??) and (??), which are referred to a location outside the dipole where isopycnals are flat. In this case, the ecosystem moves with a spatial and temporal dependent speed since the photic depth is shallower than that of the dipole separatrix, that is, z=−240 m < zS. As in the previous section, we focus only on Pdynamics since grazing is much smaller than the Pgrowth. We observe that two different mechanisms cause a Pincrease of different order of magnitude. On the one hand, a P′maximum increase of about 7% develops at the wake of the anticyclone, since this vortex type has isopycnals uplifted in the photic zone. As discussed in the previous section, this increase is generated by vertical advection (Fig. 3.13a) and diminishes afterwards by biological forcing (Fig. 3.13b). Note that the observed Pincrease is independent of whether NTis initialized homogeneous on horizontal layers or isopycnals, and thus of the ecological initial conditions. On the other hand, positive P′increases up to a 42% ˆ Pat the filament front. Since the patch translation speed has negative vertical shear, deep Players translate further southwards than shallow ones. As a result, Pself-shading decreases at the filament front in benefit of Pnet growth (dP/dt)′(Fig. 3.13c). 3.3. NUMERICAL SIMULATIONS 63 −20 −15 −10 −5 0 5 10 15 20 25 −500 −400 −300 −200 −100 0 (km) (m) zS y z (a) −20 −15 −10 −5 0 5 10 15 20 25 −500 −400 −300 −200 −100 0 (km) (m) zA rY z (b) Figure 3.12: Position of the NTfrontal edge along dipole axis, x= 0, centered at y1and relative to the dipole. Numerical results (a) are approximated analitically (b) using the exterior geostrophic velocity ˜vg e(r, t) of the QG approximation, equation (3.12), where rY is the distance between the front edge and the dipole center. In both cases, the vertically straight NTfront tilts gradually with time around the minimum separatrix depth zS. Each line corresponds to consecutive times ranging from t∈(ti,15.9) d with δt ∼ =0.7 d. In order to quantify the contribution of both mechanisms to the Pincrease, ecological properties are integrated on the patch volume above the photic depth using (3.11). At initial times, hP′irapidly increases up to 7% due to vertical advection (Fig. 3.13a), though (w∂P/∂z)′nearly balances out since whas a symmetric distribution within the patch (Fig. 3.14). At times t > 12 d, the patch has been decelerated relative to the dipole and thus Pvertical advection decreases while patch filamentation increases. The elongation of this filament enlarges with depth causing Pself-shading to diminish at its front in benefit 64 CHAPTER 3 −50 −25 0 25 50 −50 −25 0 25 50 (km)y 0.80 0.88 0.95 1.05 1.12 1.20 −50 −25 0 25 50 −50 −25 0 25 50 (a) P’ t=2.8 d −50 −25 0 25 50 −50 −25 0 25 50 (km)y −50 −25 0 25 50 −50 −25 0 25 50 (b) t=12.7 d −50 −25 0 25 50 −50 −25 0 25 50 (km)y (km)x 0.89 1.00 1.11 1.21 1.32 1.42 −50 −25 0 25 50 −50 −25 0 25 50 (c) P’ t=28.3 d Figure 3.13: Horizontal distributions at z=−125 m (iZ= 112) and at the indicated times of P′(grey scale) and (a) (w ∂P/∂z)′(contour lines range is (−0.28,0.48) with δ= 0.05), and (b)-(c) (dP/dt)′(contour line range is ∈(−0.14,0.4) with δ= 0.02 in b and δ= 0.1 in c). PV contours = 0.2 are included for reference. Two mechanisms are responsible for P′ increase, and biological forcing at the filament front. 3.4. CONCLUDING REMARKS 65 10 20 30 40 1.00 1.05 1.10 1.15 1.20 1.25 1.30 1.35 -2 -1 0 1 2 3 66 67 68 69 70 71 72 73 < P’> < (GP)’> < ( P/ t)’>d d -1 <I> (W m-2) <(w P z)’> x 10-2 days Figure 3.14: Time evolution of hP′i(solid line), h(dP/dt)′i(dotted line), hG′ Pi(long dashed line), h(w ∂P/∂z)′i(short dashed line), and hI′i(dash-dotted line) averaged on the patch above the photic depth z=−240 m (iZ= 97). The observed hP′iincrease is caused by two mechanisms. At initial times, vertical advection accounts for an enhancement of about 6%. Later on, a decrease in Pself-shading in benefit of solar radiation increases GP, which doubles hP′iat t= 45 d. of solar radiation I. As a consequence, GPconverges to a value around 1.23, which forces the ecosystem to be unsteady (hdP/dti>0) and ultimately doubles hP′iwithin the following t= 45 d. 3.4 Concluding remarks Three-dimensional plankton distributions were related to potential vorticity (PV) using a numerical physical-biological coupled model that explicitly conserved PV on isopycnals. In particular we investigated how vortex translation and the distant PV action perturbed an oligotrophic system in steady state. The former fact was addressed by considering a mature dipole embedded in the ecosystem, while the latter by initializing a subsurface dipole and a surface ecosystem patch. We stated that dipole separatrix, which was defined by PV edges, acted as an impermeable barrier, creating two different ecological niches. Firstly, phytoplankton and zooplankton trapped within vortices converged to an equilibrium state, which was sustained by regener- 66 CHAPTER 3 ated production. This may explain the subsistence of communities within long-lived vortices moving across the open ocean (Mackas et al.,2005;Whitney et al.,2005). The steady state observed depended on both the isopycnal vertical displacement, and thus on PV, and the ecological initial configuration. This was in agreement with previous numerical and experimental works which concluded that vortex geometry changes (L´evy,2003) and biogeochemical properties of source waters (Thompson et al.,2007), respectively, alter phytoplankton distributions. One caveat worth noting is that we assumed a fully-recycling ecosystem and slow mixing processes and sedimentation are likely to occur as vortices translate. The way these facts alter the ecological steadiness within vortices is left for future research. Secondly, plankton heterogeneity at the dipole outer was generated by vertical and horizontal advection, which were of the same order of magnitude as biological forcing. Plankton enhancement through vertical advection occurred at the wake of the dipole. As the dipole translated, fluid particles were displaced from the vortices front to their rear around PV edges. Therefore, they were upwelled to better lit depths when moving anticlockwise and clockwise in surface and subsurface dipoles, respectively. Phytoplankton response to this forcing had a time lag of about δt = 1.4 d causing a spatial uncoupling between the increase in phytoplankton and its causing mechanism. As a result, fluid particles with enhanced phytoplankton biomass accumulated at the vortices wake and since they moved slower than the dipole phase speed, a trail of phytoplankton was observed behind it. The extent of this trail was proportional to the phytoplankton material rate of change dP/dt. So that, for a given length of the phytoplankton trail and dipole speed, dP/dtcould be inferred. The formation of a phytoplankton trail behind vortices was first conceptualized by Olaizola et al. (1993). In this work, we took a step further by conceiving two upwelling mechanisms ocurring simultaneously in a single vortex. A permanent upwelling of trapped waters within the vortex separatrix and a transient upwelling along the vortical periphery (Flierl and McGillicuddy, 2002). Both upwellings may explain the chlorophyll trail developed below a shallow chlorophyll peak (Ning et al.,2004, Fig. 9d) at the wake of a pair of vortices of opposite sea level anomaly (Chang et al.,2010, Fig. 10). We also evidenced that horizontal advection was indirectly responsible for a phytoplankton increase when negative vertical gradients of PV existed. Since horizontal speed |u|h was maximum along the dipole axis and increased with depth, a filament of phytoplankton runned along the direction of the dipole trajectory, the extent of which was maximum at 3.4. CONCLUDING REMARKS 67 depth. The horizontal extension of this filament was qualitatively in agreement with the analytical relation between the front displacement and vortex PV derived using a simple theoretical quasigeostrophic dipole model composed of spheric vortices of fixed radius and depth. Therefore, the filament elongation was approximated for a PV structure of known geometry and intensity. As a result of the negative |u|hshear, phytoplankton increased at the filament front due to a decrease in self-shading in benefit of light irradiance. Though the light change was small, it accounted for significant enhancements of phytoplankton as was reported within an anticyclone over which waters with higher transparency flooded (Baird et al.,2011). Thus, mesoscale subsurface dipoles are responsible of off-shore transport of surface chlorophyll filaments Serra et al. (2010) and also of phytoplankton enhancement. The mechanism described may explain the subsistence of subsurface communities as filaments penetrate into the open ocean, which otherwise would decay through mortality and sinking. Finally, sharp phytoplankton fronts such as the one described here are often associated to well-defined density fronts (Jones et al.,1991), which induce upwelling along the filament edges (Moisan and Hofmann,1996). The simultaneous contribution of a change in the light regime and localized upwelling to the community subsistence in coastal filaments is left for further investigation. Chapter 4 Plankton resonant response to light and nutrients in mesoscale vortices Chapter submitted as M. Claret, A. Viúdez, and Y. H. Spitz, 2012: Plankton resonant response to light and nutrients in mesoscale vortices. Nom´es viu qui pregunta. P`ortic, Miquel Mart´ı i Pol 71 ABSTRACT The response of a fully-recycling ecosystem to light and nutrients variability within mesoscale cyclones is investigated. A simple physical-ecological coupled numerical model is used to study the horizontal and steady flow of ellipsoidal vortices. Ecological variables have only vertical spatial dependence through phytoplankton growth. This enables us to introduce mesoscale forcing in terms of vertical displacement of isopycnals in a NPZ (Nutrients-Phytoplankton-Zooplankton) model. Small-scale motions is parametrized as a Fickian-type diffusion along isopycnals. We observe that stationary stable ecological solutions are possible within mesoscale vortices with non-zero Pand Zprofiles. The dependence of these solutions on small-scale motion, vortex intensity, and trophic regime is explored. Firstly, small-scale motions increase the spatially integrated Pbiomass at a characteristic diffusion coefficient K∗, at which resonance between phytoplankton and diffusive timescales occurs. Two mechanisms are involved in this Pincrease, namely nutrients uplift to better lit levels, and Zgrazing pressure decrease. Secondly, vortex intensity exerts a double effect on the ecosystem. On one side, it enhances the irradiance affecting isopycnals, which induces a positive linear response in P. On the other side, it increases the diffusive fluxes. As a result, K∗is unaffected by vortex intensity, that is, Penhancements caused by isopycnal doming and small-scale motions are approximately additive. Finally, the trophic regime determines the magnitude of the ecosystem response. The increase in Pbiomass caused by small-scale motion is significant only in mesotrophic regimes, and it remains independent of the trophic condition when the vortex intensity varies. 78 CHAPTER 4 0.0 0.5 1.0 1.5 2.0 2.5 3.0 r −3.0 −2.5 −2.0 −1.5 −1.0 −0.5 0.0 z Figure 4.1: Vertical distribution of D(r, z) (dashed line) corresponding to the polynomial function (??) is compared with that obtained using the CASL algorithm (solid line) for a cyclone with Dmax = 0.39. PV contours with maximum = 1.8 and δ = 0.3 (thickest line) are included for reference. (b) Vertical distribution of D(s, d) corresponding to Dmax = 0.39 (upper case, λ0= 0.027, solid line), and Dmax = 0.16 (λ0= 0.011, dashed line). Contour interval δD= 0.025. places Dmax at z=−π/3≃zD. On the horizontal, Ddecreases exponentially with the radius, similarly to the PV distribution. However, a radial linear dependence of type π−ris only considered in the polynomial expression of D. This simplification allows a closed-form solution for the isopycnal radius, which facilitates the physical-ecological coupling (explained in section 4.2.3). The resulting polynomial expression for D(r, z) is D(r, z) = −λ0z(z+π)2(π−r),(4.11) where λ0is a dimensionless scaling factor. This parameter adjusts the magnitude of our analytical solution for Dto the numerical one ˆ Dthrough the relation λ0= 27/(4π4)ˆ Dmax. 4.2.2 The NPZ model The ecological model consists of three dependent variables, namely, phytoplankton (P), zooplankton (Z), and dissolved nutrients (N). Each variable represents a functional group and they are related through physiological processes. Functional groups are composed either by one dominant species or several species with similar physiological characteristics. This simplification of the food web minimizes the error associated to unconstrained free parameters, while maintains the ecological dynamics of our interest. The addition of complexity to the ecological model is found convenient when reproducing specific data sets (Friedrichs et al., 2006), rather than when used in a heuristic way as in this work. 4.2. PHYSICAL-ECOLOGICAL COUPLED MODEL 79 We consider first an oligotrophic environment typical of the open ocean (Lawson et al., 1996). In this case, the light-nutrients segregation implies an acute scarcity of resources, and species compete by optimizing their nutrient assimilation. This ecological system is fully-recycling and satisfies the set of equations dP dt=N K0+NL P |{z } GP −R0(1 −e−Λ0P)Z |{z } R −Ξ0P |{z} MP ,(4.12) dZ dt= (1 −Γ0)R |{z } GZ −Θ0Z2 |{z} MZ ,(4.13) dN dt=−GP+ Γ0R+MP+MZ,(4.14) where d( )/dtis the material time derivative. Phytoplankton production (GP) is nitrogen and light limited. The nutrient uptake follows a Michaelis-Menten kinetics (Dugdale,1967) governed by K0, an indicator of the Paffinity for N. The photosynthetic rate Lhas also a saturating response (Smith,1936;Jassby and Platt,1976) according to L(r, z, t)≡V0Ψ0I(r, z, t) pV2 0+ Ψ2 0I2(r, z, t),(4.15) where the available radiation I(r, z, t)≡I0exp ½Awz−ApZ0 z P(r, z′, t)dz′¾(4.16) is the surface solar radiation attenuated by both sea water and Pself-shading. Phytoplankton is removed by zooplankton grazing R, which has an Ivlev response (Parsons et al.,1967), and by a linear term (MP) including death and exudation. Only a fraction (1 −Γ0) of the Pgrazed is converted into Zproduction GZ, while the remainder is transferred to the N compartment. In contrast to P,Zhas a quadratic closure term for two purposes. Firstly, it gives stability to the model (Edwards and Yool,2000), and secondly, it introduces implicitly the effect of intratrophic cannibalism (Ohman et al.,2002;Pitchford and Brindley,1998) and predation by higher trophic levels (Ohman and Hirche,2001). Finally, detrital matter (Γ0R+MP+MZ) is remineralized making nutrients again available for Puptake. 80 CHAPTER 4 4.2.3 Physical ecological coupling Ecological steady-state solutions are found here numerically by considering material conservation of total nitrogen NT. The small-scale and sub-grid ecological dynamics are parametrized through Fickian diffusion. We focus next on the spatial heterogeneity introduced at the mesoscale, which depends on three main factors. Firstly, Pgrowth is affected by the vertical light regime, which is determined by the isopycnal depth. Secondly, the diffusive flux depends on the isopycnals slope. And thirdly, NTis homogeneous on isopycnals but we assume a negative diapycnal gradient of NT. Consequently, the perturbation in the ecosystem caused by the mesoscale forcing can be characterized through the spatial distribution of isopycnals, which is given by D(r, z). We focus next on expressing the NPZ equations in terms of D(r, z). In order to simplify the derivatives along isopycnals we switch from cylindrical (r, z) to isopycnal (s, d) coordinates. An expression for the isopycnal radius s(r, z) is obtained by combination of (4.7) and the line element definition ds2≡dr2+ dz2, that is, s(r, z) = Zr 0s1 + µ∂D ∂r′(r′, z)¶2 dr′.(4.17) Considering (??), the above integral has the following exact solution s(r, z) =   p1 + β2r, γ = 0 , −1 γ³τ+β 2ln ¯¯¯τ−β τ+β¯¯¯´+C(z), γ 6= 0 ,(4.18) where β(z)≡λ0z(z+π)2,(4.19) γ(z)≡3λ0(z+π)(z+π/3) ,(4.20) τ(r, z)≡p(1 + γ(z)(π−r))2+β(z)2,(4.21) and C(z) is obtained from the boundary condition s(r= 0, z) = 0. The vertical displacement of isopycnals D(s, d) as a function of sand d, that is D(s(r, z), d(r, z)) ≡D(r, z), is found using a numerical iterative procedure. 4.3. INITIAL CONDITIONS 81 Finally, parametrizing the small-scale isopycnal mixing as a Fickian diffusion, the concluding NPZ equations to solve are GP−R−MP+κ∂2P ∂s2= 0 ,(4.22) GZ−MZ+κ∂2Z ∂s2= 0 ,(4.23) −GP+ Γ0R+MP+MZ+κ∂2N ∂s2= 0 ,(4.24) where κis the isopycnal mixing coefficient, and the spatial dependence of the available radiation is computed using z(s, d) = d+D(s, d). Above, and hereinafter, unless otherwise specified, all the dependent variables are assumed to be function of the independent variables (s, d), e. g. , P=P(s, d), etc. Note that the material rate of change of the the total nitrogen NT≡P+Z+Nbalances its isopycnal diffusion dNT dt=κ∂2NT ∂s2.(4.25) Equations (4.16)–(4.18) are used in the next section to obtain stationary ecosystem solutions with the fluid at rest. 4.3 Initial conditions 4.3.1 Numerical parameters 4.3.1.1 Spatio-temporal domain The unit of space in the QG domain is defined by the vertical extent LZ=π. Thus the radial extent LR=cLZ. The number of grid points is (nR, nZ) = (64,128), and the number of isopycnals nL= 128. Variables can be adimensionalized by specifying two dimensional parameters. Firstly, the maximum depth Zmin ≡ −150 m determines the length conversion factor according to L≡ |Zmin|/LZ. Secondly, the maximum photosynthetic rate Lmax, corresponding to that at the sea surface Lmax =α0I0V0/pV2 0+ (α0I0)2∼ =3 d−1with constant sea surface irradiance I0, gives the time conversion factor as T≡L−1 max. Consequently, ecological variables are expressed in terms of mmol N. Dimensional physical values can therefore be 82 CHAPTER 4 recovered by multiplying the adimensional values by the spatio-temporal conversion factors Land Televated to the appropiate powers to match physical dimensions. For exemple, the dimensional equivalent of a Pproduction GP= 700 mmol N is 700/(L3×T) = 0.02 mmol N m−3d−1. 4.3.1.2 Parameter λ0 The intensity of cyclones is modeled through the magnitude of the vertical displacement of isopycnals D, and thus through the dimensionless scaling factor λ0. We have considered spherical cyclones from = 0.75 to = 2, which is at the limit of their stability. This range corresponds to Dmax ∈[0.16,0.39], using the CASL algorithm aforementioned (section 4.2.1), and thus to λ0∈[0.11,0.27]. An arbitrary δλ0= 1.7×10−3is chosen. Note that λ0does not alter the pattern of D(r, z) but that of D(s, d). As Dmax increases, the corresponding isopycnal level deepens (Fig. 4.3). 4.3.1.3 Isopycnal diffusion coefficient κ Suggestive eddy diffusivities are extracted from the North Atlantic Tracer Release Experiment (Ledwell et al.,1998;Polzin and Ferrari,2004). Scale-dependent measurements of lateral dispersion were performed at 300 m and estimated κfor scales ranging 1–30 km at about 2 m2s−1. Subinertial vortical structures had a major contribution to the mentioned 0.0 0.5 1.0 1.5 2.0 2.5 3.0 s −3.0 −2.5 −2.0 −1.5 −1.0 −0.5 0.0 d Figure 4.2: Vertical distribution of D(s, d) corresponding to Dmax = 0.39 (upper case, λ0= 0.027, solid line), and Dmax = 0.16 (λ0= 0.011, dashed line). Contour interval δD= 0.025. 4.3. INITIAL CONDITIONS 83 isopycnal coefficient with κ∈(0.5,2.5) m2s−1, followed by mesoscale-submesoscale interactions (κ≃0.1 m2s−1), ultimately developing into inertia–gravity waves (κ≃0.01 m2s−1) (Polzin and Ferrari,2004). The estimates above are based on mesureaments of an inert tracer. However, the dispersion rate of κof a reactive tracer, such as plankton, depends on its reactive timescale (Pasquero,2005). Furthermore, spectral slopes between inert or active tracers and chlorophyll diverge at a critical lengthscale 2–7 km, where the flow and ecological timescales have the same order of magnitude (see review Martin,2003). Many theories may account for this divergence such as the scale of the tracer supply (Bracco et al.,2009), and thin layers (Franks,2005). As a result of the uncertainity about bounding an isopycnal difussion coefficient for ecological variables, we constrain κwithin the widest plausible range κ∈[0.01,10] m2s−1. According to the above mentioned scaling, the isopycnal diffusion coefficient is adimensionalized such that K≡κT (Lc)2.(4.26) Note the squared dependence on cdue to the fact that the QG space is vertically streched. It must be kept in mind also that the parameter cspecifies QG mesoscale vortices as well as submesoscale or small scale vortices. 4.3.1.4 Trophic regime The nutrients availability of the system is modified through the maximum of total nitrogen T1. We have considered a range of T1from 1L3, which is representative of an oligotrophic condition, to 8L3, with increments δT1= 1L3mmol N (see Fig. 4.4 for associated vertical profiles). 4.3.2 Stationary ecosystem with the fluid at rest A first step in solving the time evolution of the ecosystem is to define the appropriate initial conditions for the NPZ dynamics. Stable steady-state solutions of the coupled model are obtained by time integration of (4.16)–(4.18), which converge to an equilibrium state when using an initial solution close enough to that state. To this end, we first seek for ecological steady state profiles with the fluid at rest, that is, outside the vortex. In this case, isopycnals 84 CHAPTER 4 0.0 0.1 0.2 0.3 0.4 0.5 0.6 −3.0 −2.5 −2.0 −1.5 −1.0 −0.5 0.0 d P ^ x 3(mmol N)L 01234 −3.0 −2.5 −2.0 −1.5 −1.0 −0.5 0.0 d Z ^ x 3(mmol N)L 02468 −3.0 −2.5 −2.0 −1.5 −1.0 −0.5 0.0 d N ^ x 3(mmol N)L 02468 −3.0 −2.5 −2.0 −1.5 −1.0 −0.5 0.0 d NT x 3(mmol N)L Figure 4.3: Steady-state profiles ˆ P,ˆ Z, and ˆ Nfor a given profile of NTwith the fluid at rest. Ecological variables are adimensionalized and expressed in terms of mmol N. Hereinafter, progressive line thickness indicates an increase in NTmaximum δT1= 1 L3mmol N (solid lines) and δT1= 0.2L3mmol N (dashed lines). 4.4. RESULTS 85 are flat (D= 0), so that light availability is horizontally homogeneous and the ecosystem has only vertical variability. The vertical dependence of {NPZ}is obtained by time integration of (4.16)–(4.18), with K= 0 and using the oligotrophic non-dimensional parametrization (Table 4.1) of Lawson et al. (1996) and the following initial profiles NT(d) =    T1cos h³1−d d1´2π 5i, d ∈[d1,0], 1, d ∈[−π, d1), (4.27) (Pi, Zi)(d) = S1µ1−d dmin ¶,(4.28) where T1∈(1,8) L3mmol N, S1= 0.1L3mmol N, d1∼ =−0.91 (isopycnal index, iL= 91), and dmin =Zmin. Time integration is performed using an explicit leap-frog scheme, together with a Robert-Asselin time filter to avoid the computational mode. The initial profiles smoothly converge to steady-state solutions (Fig. 4.4), which are hereinafter referred to as {NT,ˆ P, ˆ Z}. Phytoplankton and zooplankton have similar vertical distributions (Fig. 4.4), though the latter is slightly larger than the former, with a subsurface maximum at d∼ =−0.86 (iL= 92). The obtained {NT,ˆ P, ˆ Z}distributions have a double purpose. On the one hand, they are used as a first guess in the coupled model with given distribution of D(s, d) by projecting them on isopycnals. On the other hand, they define the boundary condition at smax, that is, outside the vortex where isopycnals are flat. 4.4 Results In this section we investigate the phytoplankton response to light and nutrients. To this end, different cases are explored considering the isopycnal vertical displacement D, the nondimensional diffusion coefficient K, and the maximum of total nitrogen T1(see section 4.3.1 for parameters range). The distribution of Dmodifies the phytoplankton light regime of a given isopycnal, while diffusion tends to homogenize ecological variables along isopycnals. The factor T1defines the trophic regime of the system. The ecosystem dynamics without diffusion is first analyzed considering vortices of variable isopycnal vertical displacement embedded in various trophic regimes (section 4.4.1). Next, diffusion is considered in these 86 CHAPTER 4 symbol description value units Awlight attenuation by sea water 0.03 m−1 Aplight attenuation by phytoplankton 9.5×10−3m2mmol N−1 I0surface available radiation 234.11 W m−2 Ψ0initial slope of the P-Icurve 0.025 m2W−1d−1 V0phytoplankton maximum uptake rate 3.6 d−1 K0half-saturation for phytoplankton uptake 0.5 mmol N m−3 Ξ0phytoplankton specific mortality rate 0.1 d−1 R0zooplankton maximum grazing rate 0.8 d−1 Λ0Ivlev constant 0.4 mmol N−1m3 Γ0fraction of zooplankton grazing egested 0.25 Θ0zooplankton excretion/mortality rate 0.05 d−1 Table 4.1: List of ecological constants corresponding to an oligotrophic environment as extracted from Lawson et al. (1996). The Zspecific excretion/mortality rate Φ0is computed normalizing Θ0by the spatially averaged Z. All constants are made adimensional in the physical-ecological coupled model. different trophic conditions (section 4.4.2) and vortex types (section 4.4.3). One general result is that ecological stationary profiles numerically stable are reached in all considered cases. 4.4.1 Isopycnal vertical displacement and trophic regime Isopycnal doming deepens the isolines of N,P,and Ztowards the vortex center (Fig. 4.5) because the irradiance received by an isopycnal increases when it shallows. Thus Pincreases in detritment of N, and so does the grazer community. In addition, when T1increases from 1L3mmol N to 8L3mmol N, the maximum of Pis doubled (Figs. 4.5a,b), while that of Zis enhanced by one order of magnitude (Figs. 4.5c,d). In fact, this is not an original contribution since phytoplankton growth is limited by light and nutrients. Our aim is to quantify this limitation when the light regime is modulated by cyclones. To this end, we define the anomaly of any variable χ(s, d, K,Dmax, T1) with respect to a reference state where isopycnals are flat (dotted lines in Fig. 4.5) as χ′(s, d, K0,0.39, T1)≡χ(s, d, K0,0.39, T1)−χ(s, d, K0,0, T1),(4.29) where K0= 0. The largest isopycnal vertical displacement, and thus irradiance change, occurs at the vortex center at isopycnal dI=−1 (Fig. 4.3). Accordingly, in the T1= 8L3 mmol N case, the maximum of P′is located close to dI, but it turns out to be much deeper 4.4. RESULTS 87 (d∼ =−1.9) when T1= 1 L3mmol N (Figs. 4.6a,b). This paradox is addressed by subsequently splitting the anomalies of the physiological processes. Anomalies of phytoplankton growth GPand grazing Rare positively coupled as are P′and Z′in both trophic conditions. However, they are one order of magnitude larger in the T1= 8 L3mmol N case than the lowest T1case (Figs. 4.6c,d). The production of Pdepends upon the light but also upon nutrients. While the former attenuates with depth, the latter increases. When T1= 1 L3mmol N, the distribution of P′resembles that of the nutrient uptake. In contrast, when T1= 8 L3mmol N the pattern of P′correlates better with the photosynthetic rate L (Figs. 4.6e,f). Thus we have characterized two different ecological dynamics depending on the light regime and the trophic condition. We have focused so far on one vortex type and on two trophic regimes. In order to NT= 1 L3mmol N NT= 8 L3mmol N 0.0 0.5 1.0 1.5 2.0 2.5 3.0 s −3.0 −2.5 −2.0 −1.5 −1.0 −0.5 0.0 (a) d 0.0 0.5 1.0 1.5 2.0 2.5 3.0 s −3.0 −2.5 −2.0 −1.5 −1.0 −0.5 0.0 0.00 11.41 22.81 34.22 45.63 57.04 68.44 79.85 (b) d x 103 0.0 0.5 1.0 1.5 2.0 2.5 3.0 s −3.0 −2.5 −2.0 −1.5 −1.0 −0.5 0.0 2.75 17.91 33.07 48.22 63.38 78.54 93.69 108.85 (c) d x 103 0.0 0.5 1.0 1.5 2.0 2.5 3.0 s −3.0 −2.5 −2.0 −1.5 −1.0 −0.5 0.0 16.94 138.92 260.90 382.88 504.85 626.83 748.81 870.79 (d) d x 103 Figure 4.4: Vertical distributions with K= 0 and Dmax = 0.39 of (a,b) P(s, d) in shaded contours and ˆ P(d) (dotted line, max ∼ =70.5×103mmol N, δ≃105mmol N); and (c,d) N(s, d) and Z(s, d) (solid line, Zmax ∼ =52.79 (592.57) ×103mmol N, δZ ≃104(105) mmol N in c (d)). Left panels correspond to NTmaximum T1= 1 L3mmol N, while right panels to T1= 8 L3mmol N. 94 CHAPTER 4 4.4.3 Isopycnal vertical displacement and diffusion We analyze next whether the characteristic resonant diffusion coefficient depends on the isopycnal vertical displacement. To this end, we fix the trophic condition to an arbitrary value T1= 1.4L3mmol N within the range 1 < T1<2L3mmol N, where K∗shows the highest sensitivity (Fig. ??a). The maximum Presponse to diffusion occurs when Kis O(0.1). Thus we explore a range of Kcomprised within (0,0.4) considering a subset of vortices with Dmax ∈(0.16,0.39). Each case is characterized by the enhancement factor with respect to the state without diffusion, that is, ´ P(K,Dmax,1.4L3mmol N). We observe that Dmax affects the magnitude of ´ Pmaximum but not the value of K(Fig. 4.9a). Firstly, as stated in the first section, Pincreases linearly with Dmax and accordingly a linear increase is also expected in ´ P. Secondly, K∗shifts from 0.08 to 0.1 at Dmax = 0.16 and Dmax = 0.4, respectively. The resonant behaviour is reached when Pand diffusive timescales are comparable. So that, the factors determining the location of ´ Pmaximum in the (K,Dmax) space should be related to the Pdiffusive rate and the increase in the Pgrowth rate caused by diffusion. In this regard, ´ Pmaximum is reached when η≡K∂2P/∂s2|Ω (G′′ P−M′′ P)|Ω∼ =1.08 ,(4.36) which corresponds to the characteristic non-dimensional diffusion coefficient ˜ Kat which ´ Pis maximum for this trophic condition. Isopycnal doming increases diffusion, since it steepens isopycnals slopes, and also phytoplankton growth, because it enhances the irradiance received by the uplifted phytoplankton. As a result, K∗is nearly unaffected by Dmax. However, phytoplankton growth and diffusion do not account totally for ´ P(Fig. 4.9b). We observed that diffusion uncouples Pand Z, and thus decreases the grazing pressure on P. Consequently, when Zgrazing is included, the resulting balance correlates with ´ P(not shown). 4.5 Concluding remarks We constructed a simple numerical model to investigate the dynamics of a fully-recycling ecosystem within mesoscale vortices. Particularly, we focused on surface cyclones and subsurface anticyclones since their isopycnals are domed in the upper layer. The novelty of 4.5. CONCLUDING REMARKS 95 0.0 0.1 0.2 0.3 0.4 0.16 0.19 0.22 0.25 0.27 0.30 0.33 0.36 0.39 99.7 100.0 100.2 100.4 100.6 100.9 101.1 101.3 101.5 101.7 102.0 102.2 P (x 10−2) ′ (a) K Dmax 0.0 0.1 0.2 0.3 0.4 0.16 0.19 0.22 0.25 0.27 0.30 0.33 0.36 0.39 97.3 97.6 97.9 98.2 98.5 98.9 99.2 99.5 99.8 100.1 100.5 100.8 GP’’+ 2P/ s2 (x 10−2) _______________ K′ ∂ ∂ (b) K Dmax Figure 4.9: Distributions of (a) ´ P(K,Dmax) and η′=η−1.08 (contour line, η′∈(−0.4,0.6), δη′∼ =0.035), (b) ´ G′′ P+K∂2P/∂s2(K,Dmax) and ´ R′′ (K,Dmax) (contour line, ∈(−4,0.2) × 10−2,δ∼ =1.5×10−3). The trophic condition is T1= 1.4L3mmol N. this model was that the mesoscale forcing was introduced in terms of vertical displacement of isopycnals in the NPZ (Nutrients-Phytoplankton-Zooplankton). Using this physicalecological coupled model, the impact of the trophic regime, vortex intensity, and small-scale motions, parametrized as a Fickian-type isopycnal diffusion, on the ecosystem was explored. In all considered cases, stationary ecological solutions numerically stable, with coexistence of phytoplankton (P) and zooplankton (Z), were obtained. The trophic regime, determined by the maximum of total nitrogen, characterized two different ecological dynamics. In the oligotrophic regime case, Pand Zanomalies (with respect to a reference state where isopycnals are flat as outside the vortex) had similar magnitude. 96 CHAPTER 4 These anomalies were localized at the vortex bottom edge, where the nutrients (N) flux was the largest. When a mesotrophic regime was instead initialized, Zanomalies were one order of magnitude larger than those of P. Both anomalies had their maximum localized at the isopycnal with the largest vertical displacement, and thus the one experiencing the largest change in irradiance. Small-scale motions, parametrized as a Fickian-type diffusion along isopycnals, increased the spatially integrated Pbiomass up to 1% or 15% depending on the trophic condition. This increase occurred at a characteristic diffusion coefficient (K∗) due to a resonance mechanism between diffusive and phytoplankton timescales. That is, when the ratio between diffusive and phytoplankton growth rates was close to one. Two main mechanisms were involved in this Pincrease. Firstly, diffusion uplifted Nto better lit levels through a Ninflux that counterbalances the Poutflux. Secondly, Zuncoupled to its prey increase because the growth rate of Pwas larger than that of Z. As a result, the Zoutflux was larger than that of Pand the grazing pressure diminished in benefit of the latter. The largest contribution of diffusion in maximizing Pwas observed in mesotrophic regimes. However, large K∗values were required and some caution should be taken when interpreting them. Considering K∗= 0.3 and the upper bound of the small-scale motion parametrization (κ= 10 mmol N m3), results in a Prandtl ratio c∼ =20, which is comprised within the submesoscale. Since our basic equations were inappropiate to model submesoscale viscid dynamics, we constrain instead cto the lowest mesoscale bound (c= 35). In this case, κ∼ =30 m2s−1, which approaches to the parametrization of mesoscale motions (Martin et al.,2001;Ledwell et al.,2008). This suggests that resonance between plankton growth and advective rates occurs within non-axisymmetric vortices, which has been already reported by Pasquero (2005). Finally, the vortex intensity, modeled through the maximum isopycnal vertical displacement, exerted an antagonistic effect on the ecosystem. On the one side, it induced a positive linear Presponse since Pand Nwere uplifted to better lit depths. Isopycnal doming accounted for an increase about 8% or 15% of the Pbiomass corresponding to a scenario with flat isopycnals depending on the trophic condition. On other side, it also increased diffusive fluxes since steepened the isopycnal slopes. As a result, K∗was nearly independent of the isopycnal vertical displacement. This suggests that the Penhancement factor may be approximated by additive contributions of isopycnal doming and diffusion at K∗. Chapter 5 Vertical velocity in the interaction between inertia-gravity waves and submesoscale baroclinic structures This chapter has been published as: Claret, M., and A. Viúdez, 2010: Vertical velocity in the interaction between inertia-gravity waves and submesoscale baroclinic structures. J. Geophys. Res., 115, C12060, doi:10.1029/2009JC005921. A la vida o al cor quelcom li prenen les ones que se’n van; si no tinc res, les ones que ara v´enen dieu-me qu`e voldran? Vora la mar, Jacint Verdaguer 99 ABSTRACT The interaction between submesoscale baroclinic vortical structures and large amplitude inertia–gravity waves (IGWs), with emphasis on the vertical velocity, is numerically investigated using a high-resolution three-dimensional non-hydrostatic model. A rich variety of vortex-wave interactions are possible depending on the potential vorticity (PV) content and length scale of the submesoscale monopoles or dipoles, and on the amplitude and wavenumber of the IGWs. On the one hand the large amplitude IGWs cause horizontal and vertical advection of the vortices, which conserve their stability though their geometry is largely modified by the wave motion. On the other hand the horizontal vortical motion Doppler shifts the local frequency of the IGWs. The vortical angular velocity and vortex density stratification lead to a wave dispersion relation involving the effective Coriolis frequency (Coriolis frequency plus the vortical angular velocity) and the total Brunt-V¨ais¨al¨a frequency. This inhomogeneous change in the local wave frequency causes the IGWs depart from their initial plane geometry. In the particular case of inertial waves, the non-linear vortex-wave interaction generates spiral IGWs, having vertical velocities one order of magnitude larger than the submesoscale vortical flow in the absence of waves. 5.1. INTRODUCTION 101 5.1 Introduction Recent numerical works have reproduced the fully three-dimensional nature of submesoscale flows (e. g., Capet et al.,2008), where vertical velocity can reach values one order of magnitude greater than those at the mesoscale (Mahadevan,2006). Submesoscale structures have been reported both in the upper ocean (Rudnick and Luyten,1996;Shay et al.,2003;Capet et al.,2008) and deep ocean (McWilliams,1985;Testor and Gascard,2003;Steffen and D’Asaro,2004;Kasajima et al.,2006). On the other hand, inertia–gravity waves (IGWs) are also ubiquitous in the ocean (Garrett and Munk,1979;Miropol’sky,2001;Pedlosky, 2003), and consequently interaction between submesoscale flows and IGWs is a frequent phenomenon. Here we address this interaction, focusing on the vertical velocity, in the special case where vortical and wave flows have similar amplitudes. The submesoscale refers to flows with horizontal scales Lof order 1, and Rossby Rand Froude Fnumbers also of order 1. It plays an important role in the ocean because it facilitates the energy transfer from the mesoscale to smaller scales (Molemaker et al.,2005), and the vertical flux of momentum, buoyancy, potential vorticity (PV), and biogeochemical properties (L´evy et al.,2001;Thomas et al.,2008). In the deep ocean, long-lived submesoscale vortices are also responsible for both deep convection (Gascard et al.,2002), and horizontal transport, as they are advected away from their origin by mean currents (Testor and Gascard,2003). In the particular case of near-inertial oscillations, anisotropy of the wave field caused by the mesoscale geostrophic motions has been extensively reported (Mooers,1975a,b;Perkins, 1976;Weller,1982;Kunze,1984;van Meurs,1998;Niwa and Hibiya,1999). Several vortexwave interactions have been proposed to explain this wave heterogeneity, like wave trapping of IGWs inside vortices (Kunze,1985), wave capture (B¨uhler and McIntyre,2005), dispersion of near-inertial energy by geostrophic eddies (Young and Jelloul,1997;van Meurs,1998), inertial pumping (Rubenstein and Roberts,1986), or resonance mechanisms (Niwa and Hibiya,1999;Danioux and Klein,2008b). Here we address both the vortex-wave mechanisms that explain the wave frequency shift by submesoscale vortices and the PV structures that remain coherent after being advected by large amplitude IGWs. Our results extend the works mentioned above by using a non-hydrostatic numerical model, which considers the fully nonlinear three-dimensional momentum equations and resolves the vertical velocity 102 CHAPTER 5 with high accuracy, to simulate baroclinic assymmetric PV flows of length scales similar to those of the pre-existent wave field. The first vortex-wave interaction we introduce implies vortical motion affecting plane IGWs. This occurs both through the wave frequency Doppler-shift (Kh·ubh), where Khis the horizontal wavenumber and ubh the horizontal vortical (balanced) velocity, and through the vortical angular velocity and vortex density stratification anomaly, which lead to effective Coriolis and Brunt-V¨ais¨al¨a frequencies in the dispersion relationship of Kunze (1985). The second interaction is the advection of the vortical flow by large amplitude IGWs. In this case, the vortex geometry is largely deformed giving rise to new circulation generated in the process towards geostrophic adjustment. Finally, nonlinear vortex-waves interactions trigger a spiral IGW when a pure inertial wave is present in a vortex flow. In this case horizontal gradients of the vertical vorticity ζgenerate gradients in the effective Coriolis frequency, through the ζ/2 shift (Mooers,1975a;Kunze,1985;Rubenstein and Roberts, 1986) originating divergence of the wave field from which vertical velocity develops. In this work we use a triply periodic non-hydrostatic numerical model under the Boussinesq and f-plane approximations (section 5.2) to examine the interaction between submesoscale baroclinic vortex structures with Rossby number R.1 and pre-existent IGW fields. We particularly focus on the generation of spiral patterns of vertical velocity. The flow has constant background Prandtl ratio N/f = 10, where Nand fare constant background BruntV¨ais¨al¨a and Coriolis frequencies, respectively. We consider two types of vortical structures and waves, namely the monopolar vortex (cyclonic and anticyclonic) and the vortex dipole, and two types of waves, namely pure inertial and gravity plane waves (section 5.3). Next we investigate the flow of a monopole embedded in initially plane inertial and gravity wave fields (sections 5.4.1 and 5.4.2). The vortex, although no longer homogeneous nor steady, remains always stable despite the substantial advection by the large IGWs. Interactions between a submesoscale dipole, the simplest vortical structure having linear momentum, and large amplitude IGWs of different wavenumbers are addressed in section 5.5. The baroclinic dipole remains coherent despite the presence of large amplitude wave fields. The balanced and unbalanced components of the flow are extracted from the total flow and are separately analyzed. Finally, conclusions are given in section 5.6. 5.2. NUMERICAL MODEL AND PARAMETERS 103 5.2 Numerical Model And Parameters 5.2.1 AB-model The non-hydrostatic numerical model (hereinafter referred to as the AB-model) simulates the isochoric (volume-preserving) flow of a stratified rotating fluid under the Boussinesq and f-plane approximations (Dritschel and Vi´udez,2003). Here the flow is initialized with 1) a localized vortical flow specified by the potential vorticity PV using the PV initialization approach (Vi´udez and Dritschel,2003), and 2) a plane IGW background field (described in section 5.3.1). The theoretical basis of the numerical model is explained in detail in the references above, succinctly here in appendix Aand only a brief definition of the physical quantities is given next. The Froude number F≡ωh/Nand the Rossby number R≡ζ/f, where ωhand ζ are the horizontal and vertical components of the relative vorticity ω≡ωh+ζk, and N is the total Brunt-V¨ais¨al¨a frequency. The vertical displacement of isopycnals is defined as D(x, t)≡z−d(x, t), where d≡(ρ−ρ0)/zis the depth that an isopycnal located at xat time thas in the reference density configuration defined by ρ0+zz. Above ρ(x, t) is the mass density, and ρ0>0 and z<0 are constant values that do not need to be specified in the Boussinesq approximation. The squared total Brunt-V¨ais¨al¨a frequency is therefore N2(x, t) = N2µ1−∂D ∂z (x, t)¶.(5.1) Static instability occurs when the stratification number Dz≡∂D/∂z > 1, and inertial instability when R<−1. The AB-model integrates the dimensionless ageostrophic horizontal vorticity Ah= (A,B)≡˜ ωh−c2∇ hD, dAh dt =−fk×Ah+ (1 −c2)∇ hw+˜ ω·∇uh+c2∇ hu·∇D,(5.2) where N2≡ −gz/ρ0, the Prandtl ratio c≡N/f, the relative vorticity ω≡∇×u, the velocity u=uh+wk,∇is the gradient operator, subscript hdenotes the horizontal component, and ˜χ≡χ/f, for any quantity χ. The material derivative dχ/dt ≡∂χ/∂t +u· ∇χ. The third prognostic equation is the explicit conservation of PV anomaly through 110 CHAPTER 5 V¨ais¨al¨a frequency. In the cyclone case (C1) an effective wave frequency ωl∼ =0.13 is predicted from (??). This is confirmed by numerical results which show that the wave frequency peak evolves from inertial ωl= 0.1 to near-inertial frequency ωl= 0.12 ±0.01 (Fig. 5.3). So that the local frequency shift ζ/2 (Mooers,1975a;Kunze,1985;Rubenstein and Roberts,1986) is also caused by baroclinic PV structures that remain no longer axisymmetric after the vortexwave interaction. In order to analyze how the balanced flow affects the wave motion and viceversa, we have extracted the balanced flow from the total flow. The balanced vector potential ϕb= (ϕb, ψb, φb) is here diagnosed using the Optimal PV Balance (OPVB) approach (Vi´udez and Dritschel,2004a), and the balanced quantities are derived therefrom. From a given PV field anomaly (x, y, z), the OPVB approach diagnoses a flow having only those IGWs that have been spontaneously generated during the process of acquiring its own PV (that is, during a time interval set equal to the initialization time ti= 5 Tip). The OPVB flow does not contain most of the IGWs, which remain, almost entirely, in the unbalanced vector potential ϕi≡ϕ−ϕb. The unbalanced velocity and vertical displacement of isopycnals, are obtained directly from ϕithrough the usual relations ui=−f∇×ϕiand Di=−ǫ2∇·ϕi. An alternative way to obtain the interaction between the inertial waves and the vortical flow in this case is using the near-inertial oscillation (NIO) equation of Young and Jelloul (1997), which is valid for small Rossby numbers. We note that their geostrophic streamfunction Ψ Figure 5.3: Clockwise rotatory spectrogram of u+iv. The spectrogram comprises 194 spectra from t= 0 to t= 14.25 Tip using a window of 5 Tip and a time lag of δt = 0.5Tbp. The distribution shows the Fourier transform magnitude of the components pˆu2(ωF) + ˆv2(ωF) for the Fourier frequencies ωF< 0. The initial inertial peak evolved to near-inertial, ωl= 0.12 ±0.01. 5.4. VORTEX-WAVE INTERACTION 111 is similar to our vertical potential φ. However, we use here the OPVB because is valid for largely ageostrophic flows. The wave frequency shift mentioned before is noticeable because the vertical wave phase velocity at the vortex center (σZ= (f+ Ω)/m) is larger than outside the vortex (σZ= f/m). Therefore, phase lines of uiat x∼ =y∼ =0 accelerate inside the vortex (Figs. 5.4a,b). As a result, the originally straight phase lines of uih are broken by the cyclone, and the vertical distribution of the speed anomaly of the unbalanced horizontal velocity is distorted (Fig. 5.4c), reaching negative values at the vortex center. As we have seen, the initial inertial wave field is strongly modified by the vortex, but at the same time the vortex is also deformed by the wave velocity, which causes the PV contours to depart from the spherical geometry, modifying the vertical distribution of D(Fig. 5.5a). Since Dis related to the geostrophic velocity shear ug hz ≡∂ug h/∂z by the thermal-wind relation ug hz =−N2 fk×∇hD,(5.11) the |ug h|contours (Fig. 5.5b) depart from the circular geometry typical of a vortex in the absence of a wave field. This is also confirmed when extracting the balanced flow from the total flow (Fig. 5.5c). An important result of the vortex-wave interaction is the generation of win the form of spiral waves (Fig. 5.6). Since the motion of the isolated spherical vortex on the one hand, and the motion of the isolated inertial waves on the other hand, are purely horizontal, the development of win the vortex-wave system is a clear result of non-linear vortex-wave interaction. The maximum wamplitude reaches |w|max = 5 ×10−2, that is 5% of the horizontal inertial wave speed, from t= 5 Tip to t= 6 Tip (Figs. 5.6a,b). This spiral wpattern seems to be related to wave motion rather than to balanced motion since the QG vertical velocity wq obtained by solving the QG omega equation (Hoskins et al.,1978) c2∇2 hwq+∂2wq ∂z2= 2∇ h·Qg h,(5.12) where Qg h≡c2∇ hug h·∇ hDis the geostrophic Q-vector and ug his the geostrophic velocity. Though having a spiral pattern as well, is about one order of magnitude smaller than the total w(not shown). The unbalanced origin of the total wis confirmed by splitting it into wband wiusing the OPVB approach. The vertical distribution of w(Fig. 5.6b) follows the 112 CHAPTER 5 (a) (b) (c) Figure 5.4: Vertical distributions in the x-z plane iY= 65 (y= 0) at t= 5 Tip of (a) ui (ui∈[−1.3,1.28], δui= 0.15), (b) vi(vi∈ [−1.17,1.19], δvi= 0.15), and (c) the speed anomaly of the unbalanced horizontal velocity U′ i=|ui| − 1 (U′ i∈[−0.82,0.31], δU′ i= 0.09). Domain extend is x∈[−π, π]c,z∈ [−3/4π, 0]. The PV contour = 0.2 (thick line) is included for reference. Hereinafter, solid and dashed lines indicate positive and negative values, respectively. Figure 5.5: Vertical distributions in the xzplane at iY= 65 (y= 0) and t= 5 Tip of the (a) isopycnal displacement D(D∈ [−0.06,0.17], δD= 0.017), (b) horizontal geostrophic speed anomaly Ug′=|Ug h| − 1 (Ug′∈[−0.15,0.82], δUg′= 0.07). and (c) balanced v(vb∈[−1.55,1.18], δvb= 0.15). Domain extent is x∈[−π, π]c,z∈[−π, 0]. The PV contour = 0.2 is included. (a) (b) (c) 5.4. VORTEX-WAVE INTERACTION 113 pattern of wiand both have the same order of magnitude, which is two times larger than that of wb(Fig. 5.6c). The motion of the spiral IGWs is that of a right-handed helix (the height increasing with increasing phase, Fig. 5.7), rotating anticyclonically so that the phases propagate upwards. The wave packet propagates downward and horizontally leaving the vortical region in a few inertial periods (not shown). This spiral IGW has a local frequency ranging from fto feat initial times and extending to 2fand 3ffrequencies afterwards (Fig. 5.8). Near-inertial wis generated by divergence of the uifield, which becomes horizontally inhomogeneous because ζshifts the frequency of pure inertial waves. When separating unbalanced from balanced flows we observe that ζiand ζbhave the same order of magnitude after t > ti. On the one hand, ζiis in phase with w(not shown), as predicted from (5.6). However pure inertial waves have ζi= 0 at t= 0. On the other hand, we observe that wcorrelates with |∇ζb|maxima after t > tiand in deeper layers, where horizontal advection is minimum (Fig. 5.9a). As a result, widevelops at vortex edges (note that F>1, Fmax = 1.37, occurs once the spiral wave has been already generated). Since the vortex geometry is largely horizontally advected by an initially pure inertial wave, |∇ζb|isosurfaces become spiralized with depth (Fig. 5.9b) generating an helical IGW. Thus, while the frequency of the total wis directly related to ζ, its 3D structure is explained by |∇ζb|, in accordance with the stated correlation between wand the eddy relative vorticity (Danioux and Klein,2008a). Finally, superinertial wobserved at later times is due to resonance mechanisms, in agreement with the results of Niwa and Hibiya (1999) and Danioux and Klein (2008b), that occur when PV structures and IGWs have similar length scales. Analogous results were obtained with an axisymmetric anticyclone (case C2), having min =−0.5 and semi-axes ah/c =aZ= 1, initialized with an inertial wave field of |uih|= 0.1 and m= 6 (wavelength λZ∼ =1). In this case the vortex has Ω <0 and the local frequency ωl∼ =0.79 < f. Consequently, σZat the vortex center is smaller than that far away from it, which is the opposite effect to that described in the cyclonic case, and wave phase lines accumulate at eddy edges (Fig. 5.10a). The total walso shows a right handed helical structure (not shown), consistent with the anticyclonic rotation with time of the inertial wave velocity, which propagates horizontally and downwards at initial times but it is trapped at the eddy bottom later on (Figs. 5.10b–d). This spiral IGW has a subinertial frequency ωl≃0.08±0.01 (Fig. 5.12) and therefore close to the predicted fe. 114 CHAPTER 5 (a) (b) (c) Figure 5.6: Distributions of total vertical velocity w(w∈[−4.81,4.23]×10−2,δw = 5×10−3) at t= 5 Tip (a) in the x-yplane at iZ= 45 (z=−0.98) and (b) in the x-zplane at iY= 65 (y= 0). (c) Vertical distribution at the same time of wb(wb∈[−5.52,3.6] ×10−4, δwb= 5×10−5). Domain extend is x, y ∈[−π, π]c,z∈[−3/4π, 0]. PV contour = 0.2 (thick line) is included. Straight lines mark the horizontal (a) and vertical (b,c) sections plotted. x/c z Figure 5.7: Isosurfaces of total vertical velocity w(w=±0.02) at t= 5 Tip. The view is from the south. 5.4. VORTEX-WAVE INTERACTION 115 Figure 5.8: Domain averaged spectrogram w(ω, tk)≡ 1 nPn j=1 ˆw(xj, ω;tk) from t= 0 to t= 15 Tip, where ˆw(xj, ω;tk) is the Fourier transform of the time series w(xj, t) with t∈[tk− ∆t/2, tk+ ∆t/2]. The spatial average comprises n= 83time series equally distributed in the 3D domain. The spectrogram window is ∆t= 5Tip and the time lag δt= 0.5Tbp. The vertical dashed lines mark the frequencies f,f+ζ/2∼ =0.12, 2f, and 3f. 0.0 0.9 1.7 2.6 3.5 4.3 5.2 (x100) (a) (b) z x/c Figure 5.9: (a) As in figure 5.6 but at iZ= 43 (z=−1.03) (w∈[−4.18,3.48] ×10−2, δw = 5 ×10−3). Shaded contours are |∇hζb|(max = 5.1×10−2, ∆ = 8.7×10−3). Domain extend is x, y ∈[−3/4π, 3/4π]. (b) Isosurface |∇hζb|= 0.04 at the same time t= 5 Tip. The view is from the south. 116 CHAPTER 5 Though we do not consider in detail the long term vortex-wave interaction we note that, starting at t= 8.25 Tip and during the next 16 Tip, the inertial waves may cause the vortex become unstable in the sense that the vortex losses PV by PV filamentation (not shown). This PV filamentation increases the horizontal PV gradients remaining in the vortex and as a result, later on at t= 13.43 Tip, the flow becomes inertially unstable (R<−1). This long term instability is left for future research. (a) (b) (c) (d) Figure 5.10: Vertical distributions at iY= 65 (y= 0) of (a) vi(x, z) (vi∈[−0.24,0.2], δvi= 0.025) at t= 5 Tip and of total w(x, z) (w∈[−8.5,8.7] ×10−3,δw = 8.3×10−4) at (b) t= 5 Tip, (c) t= 5 Tip, and (d) t= 15 Tip. PV contour =−0.2 is included (solid thick line). The domain extend is x∈[−π, π]c,z∈[−2/3π, 0]. 5.4. VORTEX-WAVE INTERACTION 117 Figure 5.11: Domain averaged spectrogram w(ω, tk) as in Fig. 5.8. Spatial average comprises n= 92horizontal points equally distributed over 17 vertical levels. In this case vertical dashed lines mark frequencies fand f− ζ/2∼ =0.08. 5.4.2 Vortex And Gravity Waves Interaction In this section a spherical (in the QG space) anticyclone, having min =−0.75 and semiaxes ah/c =aZ= 1.5 in the initial configuration, is initialized (Fig. 5.13a) embedded in a gravity wave field with k= 8/c (wavelength λX/c = 2π/(ck)∼ =0.78) and Di= 10−2(case C3, Fig. 5.13b). Thus, the wave spatial scale is smaller, though of the same order, than the vortical flow scale. Note that the vortical D(due to the otherwise balanced vortex) is zero at z= 0. The maximum amplitude of Dcaused by the vortical motion is |Db|max = 0.20 at the end of the initialization period (Fig. 5.13c), which is about 20 times larger than Di. However, the amplitude of the vertical wave motion is wi= 0.06, which is 10 times larger than the typical mesoscale QG vertical velocity wq(that is, wqis about 10−3times the horizontal vortical speed). The most noticeable result is the deformation of the initially straight phase lines of wof the gravity waves (Fig. 5.14). This occurs because the oscillating fluid particles are horizontally advected by the vortex giving a new local (absolute) wave frequency ωlwhich is the Doppler shifted particle (intrinsic) wave frequency ωpby the vortex motion, according to ωl=ωp+Kh·ubh .(5.13) Since |Dz|max = 0.2, we have N∼ =N, and thus ωpis approximately homogeneous. Consequently, ωlis affected mainly by the Doppler shift Kh·ubh. The anticyclonic vortical motion ubh =ubi+vbjhas ub>0 (ub<0) at y > 0 (y < 0), implying a positive (negative) 118 CHAPTER 5 (a) (b) Figure 5.12: Distributions at t= 4 Tbp of the vertical displacement Dat (a) iZ= 40 (z=−1.2,D∈[−18,4.2] ×10−2, δD= 1 ×10−2), and (b) iZ= 65 (z= 0,D∈ [−1,1] ×10−2, δD= 2.5×10−3). The PV contour =−0.2 at z= 0 is included. The domain extent is δx =δy = 2πc in (a), and δx =δy = 3.62cin (b). frequency Doppler shift. Hence ωlincreases (decreases) in the northern (southern) region of the vortex, so that the initially straight phase lines acquire an anticyclonic pattern. The deformation of wave phase lines caused by the vortex is not confined to the vortical region but is transferred through all the water column (Fig. 5.14). As a first approximation we assume that the vortical vertical motion can be neglected, so that ∇ h·ubh = 0, since uih =0for gravity waves, and the non-divergence condition yields ∂w/∂z = 0. Therefore the horizontal phase velocity σh≡(−∂w/∂t)/|∇hw|is constant along the water column, ∂σh/∂z = 0. Similar results were obtained for a cyclone with max = 0.75 and semi-axes ah/c =aZ= 1.5 in a gravity wave field identical to the case above (case C4, not shown). In this case the initial configuration has |D|max = 0.18, that is, 18 times larger than Di. Contrary to the anticyclonic case the initially straight phase lines acquire a cyclonic pattern because now ωl decreases (increases) in the northern (southern) side of the vortex due to the Doppler shift frequency (??). 5.5 Dipole-Wave Interaction We address here the interaction between a vortex dipole, which, unlike the monopolar vortex, possesses a net linear momentum, and large amplitude IGWs. With that purpose we first 5.5. DIPOLE-WAVE INTERACTION 119 (a) (b) Figure 5.13: Distributions of w(w∈[−6.8,8.5] ×10−2, δw = 0.13) at t= 26.8Tbp (a) in the x-yplane at iZ= 65 (z= 0) and (b) in the x-zplane at iY= 65 (y= 0). PV contours =−0.2 at (a) z= 0 and (b) y= 0 are shown (dashed thick line). Domain extent is x, y ∈[−π, π]c, and z∈[−π, 0]. describe the flow characteristics of the dipole (section 5.5.1) initialized when pure inertial wave (section 5.5.2) or gravity wave (section 5.5.3) fields are included. 5.5.1 The dipole A submesoscale baroclinic dipole is here initialized, in the absence of waves, as two ellipsoidal PV distributions with max = 0.75 and min =−0.75, horizontal semi-axes a± X= 0.6cand a± Y= 0.4c, and vertical semi-axes a+ Z= 0.4 and a− Z= 0.27 for the cyclone (+) and the anticyclone (−), respectively (case C5, Fig. 5.15a). The initial asymmetry in the prescribed a± Zis due to the fact that these vortices are defined in the initial (reference) configuration which has flat isopycnals. During the initialization time the isopycnals stretch (shrink) in the anticyclone (cyclone), so that at the end of the initialization period (ti= 5 Tip) the 126 CHAPTER 5 (a) (b) Figure 5.20: Horizontal distributions of wat iZ= 65 (z= 0) and at t= 2.08 Tip with different wavenumbers (a) (k, l) = (8/c, 0) (w∈[−8.2,7.3] ×10−2,δw = 18 ×10−2), and (b) (k, l) = (0,8/c) (w∈[−7.5,7.5]×10−2,δw = 18×10−2). Domain extent is δx =δy = 4.32c. 5.6 Concluding Remarks In this work we have numerically investigated the interaction between idealized baroclinic vortical structures and pre-existent plane inertia–gravity waves with similar horizontal velocity or isopycnal vertical displacement amplitudes at the submesoscale. There is a large number of different possible interactions depending on the initial parameters of the vortical structures and IGWs, and we have not attempted to exhaust the very large parameter space. Two main mechanisms are usually involved in this vortex-wave interaction. The first mechanism is the advection of PV by the waves, which makes the vortical structure unsteady and forces it to be permanently in a state of geostrophic adjustment, at the same time that it modifies the upper and lower limits of the IGW frequency wave band. The second mechanism is the advection of waves by the vortices, which changes the local wave frequency through the Doppler-shift frequency relation. These mechanisms operate on submesoscale vortical structures with Rossby numbers close to, but smaller than 1, which remain always stable despite the large amplitude waves. A remarkable result is the enhancement of the total vertical velocity by an order of magnitude when inertial waves are present in vortical flows. This is a clear example of a 5.6. CONCLUDING REMARKS 127 non-linear vortex-wave interaction, which results in the generation of right-handed helical waves.Therefore, the wave frequency ranges at initial times from the Coriolis frquency fto an effective frequency feand aftewards reaches also suprainertial frequencies due to resonance mechanisms. Finally, we have considered only interactions between two kinds of submesoscale vortical structures (monopolar and dipole vortices) and two kinds of plane waves (inertial and gravity waves) and many other interactions remain still unexplored. Some examples are the interaction between localized wave packets of IGWs and submesoscale vortical structures, and the long term vortex instability of these vortical flows in presence of an inertia–gravity wave field. We also leave for further research the catalytic behaviour of vortical structures triggering IGWs. Chapter 6 Discussion This thesis aims to characterize plankton patterns associated to unstable jets and long-lived vortices at mesoscales and submesoscales. Specifically, the effects of horizontal advection, vertical advection, and ecological isopycnal mixing are numerically investigated. Additionally, the flow resulting from a particular kind of vortex-wave interaction is also analyzed due to its likely role at introducing plankton heterogeneity. We summarize here the main results raised in the context of the thesis purpose. 6.1 Ecological initialization First of all we addressed the non-trivial problem of ecological initialization using a threevariable NPZ (Nutrients-Phytoplankton-Zooplankton) ecological model (chapter 2). We sought vertical profiles in stable equilibrium with the fluid at rest because the purpose of this work was to quantify the ecological response to physical disturbances. To this end, analytical steady solutions with a constant profile of total nitrogen NTwere found. We observed that these profiles are non continuously differentiable, which implies an error source when computing vertical gradients in the advection term of the Eulerian physical-ecological coupled equations. To overcome this error, we found numerical solutions through time integration of a continuously differentiable vertical profile, which converged to a stable stationary state. In contrast to their analytical counterparts, these numerical solutions are continuously differentiable but require vertical resolutions of a few cm, which are computationally unfeasible. However, the numerical error introduced when using larger resolutions is smaller than those associated to uncertainties in ecological parametrization. Thus, numerically stable steadystate vertical profiles are suitable to initialize one dimensional NPZ models, but how are they implemented into three-dimensions? Is NThomogeneous on horizontal or isopycnal levels? Is 129 130 CHAPTER 6 NTpatchy or continuously distributed in the domain? We observed that these initial conditions led to different plankton distributions (chapter 3). Thus care must be taken in choosing proper initial conditions pertinent to the objectives of the specific physical-ecological coupled reasearch. 6.2 Horizontal and vertical advection The role of horizontal and vertical advection on plankton dynamics was investigated through different three-dimensional initializations of NT. When NTwas assumed homogeneous on horizontal levels we were able to unveil the plankton heterogeneity caused by vertical advection. Instead, when NTwas homogeneous on isopycnals of fully-developed vortices, an initial seed of horizontal heterogeneity was created. Thus we were able to compare the effects of vertical and horizontal advection on this initial heterogeneity. We considered first the former initialization to quantify the ecological impact of submesoscale vertical velocities associated to a baroclinic unstable jet (chapter 2). We observed that phytoplankton anomalies develop by vertical advection. However, these anomalies are uncorrelated from vertical velocity because phytoplankton responds to the upwelling slower than the timescale of horizontal advection. Based on that, we next considered a submesoscale surface vortex dipole where NTwas initialized constant on isopycnals (chapter 3). This scenario could be conceived as a consequence of an eddy pumping event and led us to investigate its long term evolution. We observed two different plankton dynamics spatially divided by the vortex separatrix, and thus by outer PV isosurfaces. Within vortices plankton distribution is dominated by horizontal processes. The generation of plankton anomalies due to vertical advection was insignificant compared to preexisting plankton anomalies, which were trapped inside vortices and thus translated at dipole phase speed. In contrast, outside vortices plankton heterogeneity was introduced by vertical advection, though plankton biomass was dispersed from upwelling regions, analogously to the jet case. As a result, a trail of phytoplankton developed at the cyclone wake. When frontal fluid particles moved anticlockwise to the vortex rear they were uplifted. Since phytoplankton response to the perturbation had some time lag, biomass increased at the cyclone wake. Once there, it decayed at constant mortality rate. This mortality rate was approximated to the parametric value considering a constant trail extent and dipole speed. Firstly, the dominance of horizontal over vertical advection on plankton dynamics when horizontal gradients of plankton exist was in agreement with L´evy 6.3. ECOLOGICAL ISOPYCNAL MIXING 131 (2003). Secondly, spatial uncorrelation between vertical velocities and the plankton increase caused by them had been already stated at mesoscales and submesoscales (L´evy et al.,2001; Lima et al.,2002). In this work, we took a step further by representing an scenario where both processes occurred, and suggested that potential vorticity contours may act as a spatial divider between them. The above mentioned results gave us insight on the dynamics of a plankton patch perturbed by a subsurface mesoscale dipole (chapter 3). We observed that the distant action of PV deformed this patch in such a way that a filament runned along the dipole axis, where horizontal speed was maximum. Since this speed increased exponentially with depth, so did the filament elongation. An analytical approximation to this filament extent was derived using a quasi-geostrophic model with vortices of given radius and constant PV. As a result of the negative vertical shear, a given phytoplankton layer was advanced from the layer above, phytoplankton self-shading decreased at its front, and plankton anomalies developed. Horizontal advection has been often considered incapable of introducing heterogeneities unless they already exist (see review Martin,2003), as in the previous dipole case. We developed further this concept, showing an indirect way through which vertical shear of horizontal speed creates phytoplankton heterogeneity. 6.3 Ecological isopycnal mixing We have observed that ecosystems get trapped at the interior enclosed by the vortex separatrix. This initial seed of heterogeneity is often introduced by isopycnal doming when vortex formation, such as occurs in eddy pumping. Since its long-term evolution remains largely unknown we investigated whether it reaches a steady-state considering three factors: isopycnal vertical displacement D, trophic condition, and isopycnal diffusion (chapter 4). In order to explore a wide range of the above mentioned factors, an isopycnic physicalecological coupled model, less complex than the AB-NPZ model, was first constructed. We considered a spherical mesoscale vortex in the QG space, where the flow is in gradient balance, and thus, is steady and horizontal. As a result, ecological dynamics is only vertically dependent, which let us to introduce the mesoscale vortex forcing into the NPZ model in terms of D. To this end, an analytical expression for Dwas found by approximating a polynomial function to a known Ddistribution. Finally, this physical-ecological coupled model was initialized using stationary vertical profiles numerically stable with the fluid at 132 CHAPTER 6 rest (section 6.1). We observed that these profiles always converged to a steady-state with coexistence of phytoplankton and zooplankton. When isopycnal mixing was not considered, phytoplankton biomass increased linearly with D, due to an enhancement of light irradiance, and logarithmically with total nitrogen, suggesting a saturating response with nutrients. Though these results are restricted to QG spherical vortices, phytoplankton biomass reached also a nearly steady-state in the submesoscale dipole case (Chapter 3). Thus these results may give us some insight in more complex cases that consider three-dimensional fluid motion. Another caveat worth noting is that we assumed a fully-recycling system, while in fact nitrogen is lost across vortex boundaries due to phytoplankton sedimentation. Sinking rates of phytoplankton range from 0.5 m d−1to 10 m d−1in laboratory experiments (Smayda,1970). This explains why vortices become nutrient depleted after a few months. When isopycnal mixing was considered, phytoplankton biomass was maximized at a characteristic isopycnal diffusion coefficient K∗due to a resonance mechanism between diffusive and plankton timescales. Two main effects were involved. Isopycnal doming induced phytoplankton growth, and thus nutrient consumption. This created isopycnal gradients of phytoplankton and zooplankton opposite to those of nutrients, which resulted in an outward diffusive flux of the formers and an inward of the latter. As a result, phytoplankton biomass increased directly through an upward flux of nutrients, and indirectly through a grazing decrease, which was a consequence of a diffusion of zooplankton faster than that of phytoplankton. Phytoplankton resonant response to nutrients and light has already been stated when including horizontal (Pasquero,2005;McKiver and Neufeld,2011) and vertical (Huisman et al.,1999;Ghosal and Mandre,2003) diffusions, individually. A step further was taken here by exploring a range of nutrient conditions and vortex types considering both diffusions. We observed that the trophic regime determined the magnitude of the phytoplankton response, being meaningful in mesotrophic conditions, in which the largest K∗was obtained. Conversely, K∗was nearly unaffected by isopycnal doming. In this case, the phytoplankton growth caused by isopycnal uplift to better lit levels was balanced by an enhancement of phytoplankton outward flux. We have state therefore that isopycnal mixing upwells nutrients to better lit levels, but can it balance the long-term nutrient depletion within vortices? To answer this question, we 6.4. THE ROLE OF PV IN VORTEX-WAVE INTERACTIONS 133 can compare the timescales of the vertical component of the isopycnal diffusion coefficient with a constant plankton sinking rate. The observed nondimensional resonant isopycnal diffusion coefficient ranged from 0.06 to 0.3, for oligotrophic to mesotrophic conditions, respectively. Considering the isopycnal level with the greatest slope, this interval corresponds to a vertical diffusion coefficient Kz∈[0.35,1.8] ×10−2. Estimates of phytoplankton sinking rates are about wP=−0.65 m d−1(Spitz et al.,2003). This value has been calibrated for a coastal ecosystem and is chosen as an upper bound of sinking rates in oligotrophic environments, where phytoplankton sedimentates slower than in coastal systems. Appropiately adimensionalizing wPinto WPwe obtain a |Kz/WP|-ratio comprised within [0.78,4]. Thus isopycnal mixing could account for long-term plankton subsistence in mesotrophic conditions. The validation of this gross approximation with the simple physical-ecological model constructed is left for future research. 6.4 The role of PV in vortex-wave interactions In general terms, we have stated that plankton distributions are related with PV when horizontal processes are dominant, which occurs if plankton heterogeneities already exist. Within vortices plankton and PV distributions are in phase since both propagate at vortex translation speed. Outside vortices, the distant action of PV deforms plankton patches in a non-linear way. In contrast, plankton is correlated to horizontal gradients of PV when vertical processes are relevant. However, this correlation is only observed at initial times since ecological anomalies are fastly transported by horizontal advection far from the upwelling location. Additionally, isopycnal mixing accounted for an increase in phytoplankton biomass close to the vortex center. Specifically, this occurs where isopycnal vertical displacement reached its maximum, and thus where vertical gradients of PV are the largest. In order to get an insight on how vortex-wave interaction may alter plankton dynamics we investigated the relation between PV and the flow resulting from a particular type of interaction (chapter 5). We considered spherical vortex monopoles, in the QG space, and vortex dipoles embedded in an initial pure inertial and gravity wavefield. When inertial waves were involved, a near-inertial right-handed helical wave was developed through a non-linear vortex-wave interaction. Firstly, vortices shifts the inertial frequency to an effective frequency fe=f+ζ/2 (Mooers,1975a;Kunze,1985;Rubenstein and Roberts,1986), where ζis the vertical component of the relative vorticity. Since the initial configuration of ζwas spatially dependent, 134 CHAPTER 6 horizontal gradients of fewere generated. Secondly, waves forced the vortical structure to be in permanent geostrophic adjustment through PV horizontal advection. As a result, PV contours remained no longer axysimmetric and unbalanced vertical velocities develop, which correlated with horizontal gradients of PV. When fast gravity waves were involved, the wave advection by vortices caused a Doppler shift of the local wavefrequency ωl. Thus the largest ωlshifts were observed where horizontal speeds were maximum, and hence where horizontal gradients of PV were large. Inertia–gravity waves (IGWs) have timescales from minutes to hours, often smaller than the phytoplankton growth timescale, which ranges from half a day to a couple of days. However, a combination of horizontal and vertical advection may allow plankton coupling to IGWs (Franks,1995b). In fact, some in-situ observations manifest the role of IGWs creating heterogeneities (Franks,1995a;Granata et al.,1995). The impact of the observed near-inertial spiral wave on plankton dynamics is left for future research. Overall, the present work contributes to characterize the three-dimensional plankton structure associated to mesoscale and submesoscale vortices through PV. It is thought as a process investigation rather than an attempt to simulate any particular ecosystem. Our results are the consequence of some initial numerical assumptions, which have been designed to widen the comprehension plankton dynamics in the stratified and oligotrophic open ocean. Chapter 7 Conclusions 1. Ecological steady-state vertical profiles numerically stable were found suitable to initialize physical-ecological coupled models in the Eulerian description. 2. Initialization of the above mentioned profiles, homogeneous on horizontal or isopycnal levels in the whole domain or in spatial patches, led to different plankton distributions. Thus, three-dimensional initial conditions should be carefully chosen based on the proposed research objectives. 3. Vortex separatrix divided two different plankton dynamics. Inside the separatrix, plankton distribution was dominated by horizontal advection and was in phase with potential vorticity (PV) since both translated at vortex phase speed. In contrast, outside the separatrix, plankton vertical and horizontal advections were of the same order of magnitude than the ecological forcing. Vertical advection generated plankton anomalies, which were immediately transported far away from the upwelling location by horizontal advection. As a result, plankton was initially correlated with vertical velocity, and hence with horizontal gradients of PV. 4. A phytoplankton trail was developed at the wake of a translating surface cyclone due to horizontal advection, vertical advection, and phytoplankton intrinsic timescale. The horizontal length of this trail was nearly stationary and depended linearly on the vortex propagation speed and phytoplankton mortality rate. 5. A baroclinic subsurface vortex dipole deformed surface ecosystem patches such that a filament ran along its axis. The elongation of this filament increased with depth because the vertical shear of the horizontal speed was negative. As a consequence, phytoplankton self-shading decreased at the filament front and phytoplankton biomass increased. 135