Alexandre Pereira

Dos Santos

As difficult as it may be to put every person that helped me during my path from undergraduate education until the end of this PhD I will try. On a personal level I would like to thank my girlfriend Aroa Monterde Altaba, for being always present during this years of PhD in France, for being constantly supportive and loving, and for being the best gift I received for relocating to France for this project. I thank my family, my parents Paulo Afonso Tergolina and Rosa Maria Beltran Tergolina, for always giving me the best education they could afford, for supporting my choices and for the love. I thank my sisters Camila, Silvana and Letícia for the love and support. My friends, both from Brazil and from France, who were also supportive and made my life happier. Finally, I would like to thank my colleagues and professors from UML, many who had a difficult time during these years of pandemic, but maintained their perseverance and good humour, specially Dario Canossi who became my friend and was very kind to me when I arrived in France.

On a professional level I would like to thank my co-advisors: Stefano Berti, for trusting me to work on the idea of this project, for the valuable lessons on turbulence and for the commitment and hard work during these years. Enrico Calzavarini, for providing valuable ideas, for the support with numerical problems and teaching, and as well for the commitment to this project. I thank my advisor Gilmar Mompean for his professional advice, teaching tips, for receiving me when I arrived in France

v "We have no right to assume that any physical laws exist, or if they have existed up to now, that they will continue to exist in a similar manner in the future." Plankton organisms are of paramount importance for aquatic ecology and climate due to their participation in biogeochemical cycles. They supply the planet with more than half of the total oxygen production, taking part in the carbon-dioxide exchanges with the atmosphere, and they are at the base of the marine food web. Turbulence has been recognized as a crucial factor for the survival or extinction of phytoplankton species, however, dealing with its multiscale nature in models of coupled fluid and biological dynamics is a complex task. Turbulent transport in aquatic media can manifest in different ways, involving a broad range of processes, going from large-scale eddy transport to smaller scale wave breaking and bubble formation phenomena, which are all thought to be relevant for plankton dynamics In this thesis I revisit a theoretically appealing model for phytoplankton vertical dynamics, and numerically investigate how different fluid motions affect the survival conditions and the spatial distribution of the biological population. For this purpose, I first present a literature review of the main existing theories previously developed to explain phytoplankton survival conditions and blooms. The numerical work, subsequently reported, relies on realistic parameter values and on a kinematic flow model, allowing to account for different spatial and temporal scales of turbulent motions. The dynamics of the population density are described by an advectionreaction-diffusion model with a vertically heterogeneous growth term proportional to sunlight availability. The numerical methodology is thoroughly documented, before the results are illustrated.

Therefore I explored the role of fluid transport by progressively increasing the complexity of the flow in terms of spatial and temporal scales. I found that, due to advection by large and persistent structures, phytoplankton accumulates in downwelling regions and its growth is reduced. An explanation of the observed phenomenology is provided in terms of a plankton filament model. Moreover, by contrasting the results in the different flow cases that I examined, I could show that the large-scale coherent structures have an overwhelming importance. Indeed, smallerscale motions were found to only quite weakly affect the dynamics, without altering the general mechanism identified.

The results reported in the present thesis bring indications that are relevant for parameterizations in numerical models of phytoplankton life cycles in realistic oceanic flow conditions. The thesis ends with an illustration of preliminary results from an extension of the previous model taking into account also horizontally heterogeneous light conditions, and with an exposition of different possible directions for future research in the field. Le phytoplancton et le plancton jouent un rôle primordial dans le maintien de la vie sur Terre à travers les cycles biogéochimiques. Ils fournissent à la planète plus de la moitié de la production totale d'oxygène, ils régulent le climat et sont à la base de la chaîne alimentaire marine. La turbulence a été reconnue comme un facteur d'importance cruciale pour la survie ou l'extinction des espèces de phytoplancton. Cependant, traiter de sa nature multi-échelle dans les modèles de dynamique couplée des fluides et de la biologie est une tâche complexe. Le transport turbulent dans les milieux aquatiques peut se manifester sous de nombreuses formes, des courants d'advection à grande échelle aux déferlements de vagues et formations de bulles à plus petite échelle, tous reconnus pour jouer un rôle dans la survie du phytoplancton.

Dans cette thèse, je revisite un modèle pour la dynamique verticale du phytoplancton et j'étudie numériquement comment différents mouvements de fluide affectent les conditions de survie et la distribution spatiale de la population biologique. A cet effet, je présente dans un premier temps une revue de la littérature des principales théories existantes précédemment développées pour expliquer les conditions de survie et les blooms du phytoplancton. Le travail numérique, rapporté par la suite, repose sur des valeurs de paramètres réalistes et sur un modèle d'écoulement "kinematic", permettant de prendre en compte différentes échelles spatiales et temporelles des mouvements turbulents. La dynamique de la densité de population est décrite par un modèle d'advection-réaction-diffusion avec un terme de croissance verticalement hétérogène proportionnel à la disponibilité de la lumière solaire. La méthodologie numérique est minutieusement documentée avant que les résultats soient illustrés. J'ai donc exploré le rôle du transport des fluides en augmentant progressivement la complexité de l'écoulement en termes d'échelles spatiales et temporelles. J'ai découvert qu'en raison de l'advection par les grandes structures persistantes, le phytoplancton s'accumule dans les régions de flux descendant et sa croissance est réduite. Une explication de la phénoménologie observée est fournie en termes d'un modèle de filament de plancton. De plus, en contrastant les résultats dans les différents cas d'écoulement que j'ai examinés, j'ai pu montrer que les structures cohérentes à grande échelle ont une importance primordiale. En effet, les mouvements à plus petite échelle n'affectent que très faiblement la dynamique, sans altérer le mécanisme général identifié.

x Les résultats rapportés dans la présente thèse apportent des indications pertinentes pour les paramétrisations dans les modèles numériques des cycles de vie du phytoplancton dans des conditions réalistes d'écoulement océanique. La thèse se termine par une illustration des résultats préliminaires d'une extension du modèle précédent prenant également en compte les conditions d'éclairage horizontalement hétérogènes, et par un exposé des différentes directions possibles pour la recherche future dans le domaine. Lattice-Boltzmann method [START_REF] Calzavarini | Eulerian-Lagrangian fluid dynamics platform: The ch4-project[END_REF] with a forcing of the form of Eq. 2.23 for a grid size of 600X300 cells. Here the color code represent u intensity and the black arrows the the direction of the velocity field. . . . . . 15 2.9 Okubo measures properties of ocean diffusion through the use of dye experiments. He thus verified that the 4/3 law is locally valid, given that the rate of turbulent energy dissipation varies with the time and length scales in which one is interested. Figure extracted from [START_REF] Okubo | Oceanic diffusion diagrams[END_REF] the left two dinoflagellates (photosynthetic primary producers) and a tintinnid ciliate (phytoplankton predator) on the right, both are active swimmers, Thau Lagoon of Sète, France -Uwe Kils. (II) A copepod from Antarctica, planktonic crustaceans that can be 1mm to 2mm long, Uwe Kils. (III) Antarctic krill, feeds primarily on phytoplankton and and it is probably the largest biomass of a single species on the planet, Professor Dr. Habil, Uwe Kils. (IV) Amphipod a scavenger species, Uwe kils. (V) A cyanobacterium Prochlorococcus, a photosynthetic organism and major contributor to atmospheric oxygen, Luke Thompson from Chisholm Lab and Nikki Watson from Whitehead, MIT. (VI) The sea sparkle dinoflagellate, an active swimmer species that exhibits bioluminescence when disturbed producing the mareel (Milky seas) effect, Maria Antónia Sampayo, Instituto de Oceanografia, Faculdade Ciências da Universidade de Lisboa -http://planktonnet.awi.de. . . . . . . . . . . . . . . 3.2 (I) A few hundred thousands square kilometers of phytoplankton patches in the Southeastern Pacific Ocean, filaments appear to follow the flow pattern of mesoscale ocean vortices. (II) Algae patches drawn into thin ribbons by a turbulent eddy field in the northern area of the Arabic Sea. (III) An image with seven separate spectral bands to highlight differences in plankton communities on the coast of Patagonia.

(IV) The southeast part of Rio de la Plata shows large floating populations of phytoplankton favored by the proximity to the metropolis of Buenos Aires and Montevideo, huge suppliers of nutrients such as nitrate and phosphate. https://oceancolor.gsfc.nasa.gov/. . . . . . . . . . . 3.3 Representation of the elements responsible for phytoplankton survival dynamics in aquatic environments. The euphotic zone is the sunlit layer of the ocean where photosynthesis is viable, the mixed layer is a region with active movement that mixes waters in an almost homogeneous way. In an ideal scenario for reproduction the mixed layer allows phytoplankton to have access to both the euphotic zone and nutrients. . ficient that allow a phytoplankton bloom and the combinations that prevent a phytoplankton bloom. Extracted from [START_REF] Huisman | Critical depth and critical turbulence: two different mechanisms for the development of phytoplankton blooms[END_REF]. . . . . . . . . . . . 3.8 Results analog to Fig. 2 but with the inclusion of sinking speed in the simulations. Above v = 0.04m/h and below v = 0.4m/h. Above we see the concepts of maximal turbulence and minimal turbulence represented by the horizontal lines. The critical depth and the compensation depth are equally represented by vertical lines. Extracted from [START_REF] Huisman | How do sinking phytoplankton species manage to persist?[END_REF]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Illustration of the concept behind the "Dilution-Recoupling Hypothesis". Dilution impact of the winter increase in mixed layer depth is higher on grazers than on phytoplankton (in loss rates), resulting in the bloom starting at winter. Extracted from [START_REF] Behrenfeld | Resurrecting the ecological underpinnings of ocean plankton blooms[END_REF] . . . . . . . . . . . . .

Introduction

In phytoplankton ecology, biological and fluid dynamics are tightly linked, ruling how planktonic populations interact with their environment and other organisms. Interdisciplinary work aiming to establish relations between plankton biology and the underlying physical factors has refined our understanding of plankton, and phytoplankton, ecology [START_REF] J C Prairie | Biophysical interactions in the plankton: A cross-scale review[END_REF]. Theoretical progress in the field goes in parallel with advancements in methodology and instrumentation, with the later allowing for higher resolution measurements and sampling of smaller scales than previously possible. The studies in the domain are fed by both scientific and societal interests, as it can be understood, e.g., considering the relevance of plankton for both water resource quality [START_REF] Jungo | Artificial mixing to reduce growth of the blue-green alga Microcystis in Lake Nieuwe Meer, Amsterdam: an evaluation of 7 years of experience[END_REF], and climate through its role in the carbon pump [START_REF] Severin | Openocean convection process: A driver of the winter nutrient supply and the spring phytoplankton distribution in the Northwestern Mediterranean Sea[END_REF]. Within this broad framework, this thesis focuses on the influence of physical effects, such as turbulence and light absorption in aquatic environments, on phytoplankton growth dynamics. More specifically, studying phytoplankton affected by light limitation implies the need of information as a function of depth, and which is not possible to obtain at a global scale from satellite observations. Similarly, data from ship measurement campaigns have limited coverage. Furthermore, the link between the surface and the deep ocean is still poorly known, which critically limits the current understanding of carbon transport to and from the interior [START_REF] De Lavergne | Getting to the bottom of the ocean[END_REF]. While a considerable quantity of data of horizontal properties of oceans is becoming available thanks to satellites images (Figure 1.1) and drifters with good spatial and temporal resolution, at least to some extent, this is still not the case for measurements along the vertical. The latter are sparse in space, as already mentioned, but also often discontinuous in time. Collecting ocean measurements (Fig. 1.2), for example, mandates adequate weather conditions, as if the latter are severe one deals with limited possibilities to perform research campaigns from instrumented vessels. Theoretical and numerical studies can therefore reveal particularly useful. As a matter of fact, in some situations they may represent the only viable possibility for investigation. A major question in plankton dynamics concerns the development and evolution of blooms, meaning the phases characterized by a rapid increase of the population. A notable example is the North Atlantic spring bloom. The problem is still a subject of debate. The seminal work by Riley [START_REF] Riley | The relationship of vertical turbulence and spring diatom flowerings[END_REF] was perhaps one of the first theoretical studies addressing the mathematical modeling of blooms and proposing a mechanistic explanation for their occurrence. Among the many different subsequent developments, I would like to cite, here, those by Huisman and collaborators (see, e.g., [START_REF] Ebert | Critical conditions for phytoplankton blooms[END_REF][START_REF] Huisman | Critical depth and critical turbulence: two different mechanisms for the development of phytoplankton blooms[END_REF][START_REF] Huisman | Changes in turbulent mixing shift competition for light between phytoplankton species[END_REF]), which provided a more complete picture in terms of the relevant physical parameters, particularly addressing the role of turbulence intensity, and which have been instructive as a term of comparison for the work performed in this thesis.

It is in this context of public and scientific interest that I developed this work on phytoplankton population dynamics in turbulent aquatic environments, hoping to contribute some of the unanswered questions present in the field. I present this thesis in the subsequent manner:

• In Chapter 2 I introduce turbulence and the main fluid transport processes relevant to phytoplankton populations, from a general perspective. In particular, I will discuss the basic properties of diffusive motion, turbulent transport and the concept of eddy diffusivity. The focus is put on statistical features. Both Eulerian and Lagrangian approaches are shortly presented, for future reference in subsequent chapters. The chapter ends with a concise overview of mixing processes, driven by either wind or buoyancy effects, in real aquatic environments.

• Chapter 3 is devoted to biological dynamics. Its aim is to present the modeling framework, in terms of advection-reaction-diffusion equations, that led to the main predictions about phytoplankton survival conditions, and which constitutes the basis for the numerical work I performed. A historical review of the most influential theories is also presented, from the pioneering work of Riley [START_REF] Riley | The relationship of vertical turbulence and spring diatom flowerings[END_REF], to more recent studies and currently open questions.

• Chapter 4 illustrates the numerical methodology utilized through this thesis.

It starts with an overview on advection-reaction-diffusion equation's solution techniques, that is then presented in a more specific manner to what was chosen for the thesis problem. There I will specifically detail how the advective and diffusive terms are dealt with in the advection-reaction-diffusion equation describing vertical phytoplankton dynamics in a stirred, two-dimensional, vertical fluid layer. The results of validation tests are also reported at the end of the chapter.

• In Chapter 5 I report on the main part of the research work I conducted. The primary goal of the latter was to add a more explicit representation of fluid motions to the description provided by the previous theoretical models mentioned above, in relation to the conditions allowing for a bloom. Here, I focus on the interplay between biological dynamics and fluid transport by a flow characterized by the presence of both small-scale turbulence and large-scale coherent structures. In spite of the idealized nature of the flow model employed, both these features are encountered in turbulent flows in realistic aquatic environments. The complexity of the flow is added incrementally, in terms of temporal and spatial scales, in the model I adopted. Such an approach allowed me to quantify, by means of numerical simulations, the relative importance of flow features on different scales, and to provide indications about the relevant role, hindering phytoplankton growth, of large-scale persistent motions.

• Lastly, Chapter 6 presents the work's general conclusions and perspectives. The latter include both a preliminary study of an extension of the previous model, taking into account also horizontally heterogeneous light conditions, and a discussion of other possible extensions of the work reported in this thesis.

Chapter 2

Turbulence and Turbulent Transport

Collectively turbulence and turbulent transport present one of the greatest challenges in physics, as a problem it is often considered as the last open one in classical mechanics. The difficulty is based on its non-linear multiscale characteristics, as one can glimpse by examining the sketches of Leonardo da Vinci (Figure 2.1) describing turbulent phenomena in a cascade, resulting in vortices of multiple sizes dispersed in a seemingly disorganized manner. Despite the introduction of the Navier-Stokes equation ruling the evolution of fluid flows almost two century ago, its understanding is still far from complete. A general analytic solutions of the Navier-Stokes for turbulent flows is still not achieved to date, the only available predictive theory for fully developed turbulence has been proposed by A. N. Kolmogorov in 1941 (which we will latter discuss in this chapter). Even though not all fluid transport seen in oceans, rivers and lakes is turbulent, mathematically it is the turbulent transport that present the unpredictability, that is the high sensitivity of the solution to very small perturbations that are always present in real physical systems or numerical simulations. Hence the fascination and difficulty around this concept. In this chapter we will present the basic concepts of fluid transport and turbulence which were necessary for the development of this thesis. We start by laying the fundamentals of diffusive transport in the form of a review of stochastic processes and the equation of diffusion. Then we discuss turbulence in a compact manner considering some fundamental concepts such as the Navier-Stokes equations, Kolmogorov Chapter 2. Turbulence and Turbulent Transport theory and the energy cascade picture. Then we appraise the topic of fluid transport by considering both the Eulerian and the Lagrangian approaches as we characterize tracer dispersion. We then come back to the topic of diffusion when we discuss turbulent diffusivity and its usefulness in the modelling framework. Finally, we end the chapter by presenting some basic real fluid flow phenomena in oceans and lakes.

Diffusion

In day to day life we associate diffusion to the tendency of a initially closely concentrated group of something (particles, animals, people...) to spread out in time, slowly occupying a larger area around the initial point. Still, this process is not necessarily controlled by the same rules as the formal theory of diffusion. Formally diffusion represents the transport due to the irregular thermal motion at molecular scales, in spite of that, other types of irregular motion may also be modelled by diffusion [START_REF] Neufeld | Chemical and biological processes in fluid flows: a dynamical systems approach[END_REF]. In this sense, the theory can be used to represent a phenomenon by which an initial group as a whole spreads according to the irregular motion of each individual [START_REF] Okubo | Diffusion and ecological problems: modern perspectives[END_REF]. In this picture, diffusion arises when we look at the total group or the statistical majority and we observe regularity. This way diffusion is easily associated with the macroscopic behaviour of stochastic variables.

Stochastic Processes

A stochastic or random variable refers to a variable whose values are specified only in terms of probabilities. The classical example is the result of a coin flip, it is not possible to predict weather we will get heads or tails if we are using a "fair" coin, but the probability of getting either is 50%. The same can logic can be stretched to study a variety of processes in fluids, such as fluid molecule's trajectories under thermal motion. As random variables each individual jump of a molecule is unpredictable, but their bulk behaviour is not, and as we will see mathematically, the probability that a initially close group of molecules will remain that way is much smaller than the probability it will disperse. The theory that dictates the unpredictable trajectories of particles or molecules in a diffusion process is referred to as Brownian motion or random walk.

The Random Walk theory is based on the irregular motion of the individual pollen particles, studied by botanist, Robert Brown in 1828 (hence Brownian motion). Although the theory rely on randomness there are other forms of random walk that do not behave in the same way as the pollen particles. Trajectories in a random walk can be correlated in the sense that each step has a higher probability of maintaining its direction, this is a process that has been used to analyze the movement of insects in nature [START_REF] Kareiva | Analyzing insect movement as a correlated random walk[END_REF]. Trajectories can also be biased, meaning that the probability of traveling in a certain direction is higher, this can be a factor of space variability, like chemical gradients, or an individual decision for each step, in figure 2.2 a poor drunk man tries to arrive home despite his precarious state. In here we focus on the isotropic random walk as the fluid particles we will be modelling later do not follow any biased or correlated direction.

If we consider a line (Figure 2.3) where a particle is able to perform a 1D isotropic random walk, every singular step over the line is a coin flip. In mathematical terms at every time step τ the particle has a probability p(τ) = 1/2 of going to the left or FIGURE 2.2: A drunk man returning home is an effective metaphor of a biased random walk, although his preferred direction is the one pointing to his house his inebriated state is likely to throw him off the path in an unpredictable manner.

https://prakhartechviz.blogspot.com/2019/09/random-walk-termweighting-for-text.html.

to the right (x ± 1). After one time step, the particle can either be at a distance x = 1 or x = -1, with probability 1/2 each. After the next time step (2τ), the particle will either be at a distance x = 2 or x = -2 (with probability 1/4 each) or will have returned to the origin (with probability 1/2). Therefore the probability that a particle arrives at x = n or x = -n after a time t = nτ gets lower as n increases, being that the number of possible trajectories is 2 n . FIGURE 2.3: A 1D lattice with 2n steps of size 1. In an isotropic random walk the probability p(τ) that a particle will go to the left or the right, in a time step τ, is 50%.

The probability that the particle arrives at a point m after n steps is described by a binomial distribution:

p(m, n) = 1 2 n n! {(n + m)/2}!{(n -m)/2}! . (2.1)
If the above variables are considered continuous, i.e., nτ = t and m = x, the probability density function random walk in the continuous limit is:

p(x, t) = 1 2(πDt) 1/2 exp(-x 2 /4Dt). (2.2)
Where D is called the coefficient of diffusion and it is related to the size of the random walk length steps δx and time steps τ. D is obtained as δx and τ approach zero in a way that their ratio becomes constant: (2.3) Equation 2.2 describes a Gaussian or Normal distribution which implies the random walk has finite variance and central mean (in the center of the distribution). This means that, both in our lattice-particle problem, and in a multiple coin flip situation, we would have the same result in a probabilistic frame since the random walk is a good description for a wide variety of random variable problems. The problem here was presented in a Lagrangian frame, but the interest in the random walk can be further explored when we go to the Eulerian frame and study the transport of matter by diffusion.

The Equation of Diffusion

Diffusion represents the transport due to the irregular thermal motion at molecular scales [START_REF] Neufeld | Chemical and biological processes in fluid flows: a dynamical systems approach[END_REF]. In spite of that, there are other processes that can be modelled by diffusion, the lattice-particle problem we were discussing for example can be modelled by 1D diffusion in the continuous limit. Later in this chapter we will take advantage of this assertion to explain how sometimes turbulent transport can be modeled by a diffusive process, but before we have to understand the simplest problem in diffusive transport. For this reason we start with Fick's equation of diffusion. According do Fick's law the transport amount of matter across a unit of length x in a unit of time t (the flux J x ) is proportional to the gradient of the concentration of matter [START_REF] Okubo | Diffusion and ecological problems: modern perspectives[END_REF]:

∂θ(x, t) ∂t = - ∂J x ∂x = ∂ ∂x D ∂θ(x, t) ∂x . (2.4)
Where θ(x, t) is the concentration of matter and D the diffusivity. Considering an initial condition θ(x, t = 0) = δ(xx 0 ) (where all matter is concentrated in a point x 0 at the initial time) and boundary conditions of the form θ(x → ±∞, t) → 0, the solution for D = constant has the form:

θ(t, x) = exp -(x -x 0 ) 2 /4Dt √ 4πDt . (2.5)
Not surprisingly it is very similar to equation 2.2 as it was discussed that diffusion transport behaves as a random walk. This shows that the result of the diffusion process starting from a localized concentration or density patch is a growing isotropic concentration cloud [START_REF] Neufeld | Chemical and biological processes in fluid flows: a dynamical systems approach[END_REF] whose width, or second moment is

σ 2 ≡ x 2 θ(x, t)dx θ(x, t)dx , ( 2.6)
and expands as,

σ 2 ≈ 2Dt. (2.7)
The same result as equation 2.3 as expected from a random walk in the continuous limit. So as we can see that diffusion transport is not that far from a drunk walk after all. Using a computer code with a semi-Lagrangian scheme (which we will further explore in the next chapter), we can simulate for example the diffusion of a initially concentrated patch in one dimension in an one dimensional grid. In Figure 2.4 we show that by using a sufficiently high number of grid points i.e., by decreasing δx, When thermal motion is the only mechanism contributing to diffusion, D is proportional to the temperature and is ruled by Einstein's relation D = K b T/(6πrη), where T is the temperature, K b is Boltzmann's constant, r is the diffused particle radius and η the dynamic viscosity of the fluid. In addition, relation 2.7 can be generalized for d dimensions as D ≈ σ 2 /2td. As stated, in the following sections we will take advantage of the knowledge presented here to refer to a simple but effective way to model the turbulent transport of plankton cells in oceans and lakes. But first we shall explain the basis of the concept of turbulence.

Turbulence

Da Vinci's drawing description of turbulent waters as a highly unorganized multiscale phenomenon is not far from what we can verify by opening a tap to the maximum in a full bathtub or blowing a cloud of vapor from a hot cup of tea. Disregarding his artistic genius, and in the simplest of forms: When a fluid flows in an orderly fashion, the flow is called laminar when the fluid flows in an irregular fashion with mixing, the flow is called turbulent (figure 2.5).

Chapter 2. Turbulence and Turbulent Transport

In fact this is not far from what observed on the first systematic work on turbulence carried out by Osborne Reynolds in 1883 [START_REF] Kundu | Fluid mechanics[END_REF]. When injecting a thin stream of dye into a flow of water (through a tube) Reynolds observed that at low flow rates the dye would follow a well-defined straight path in parallel lines (laminae) with no macroscopic motion across layers, a laminar flow. As the flow rate was increased the dye stream broke up into irregular motion and spread throughout the section of the tube, indicating the presence of macroscopic mixing motions perpendicular to the direction of the flow, a turbulent flow. Then he demonstrated that the transition from laminar to turbulent flow always occurred when the non dimensional rate Ud/ν ∼ 3000, where U was the flow velocity, d the tube diameter and ν the fluid's kinematic viscosity. The ratio would be later baptized by Sommerfeld as the Reynolds number:

Re = UL ν (2.8)
with L as the flow's characteristic size (in the case of Reynolds the tube diameter). This allows us to present the following characterization: Turbulence is high Reynolds number flow, dominated by non-linearity, with both spatial and temporal disorder [START_REF] Vallis | Atmospheric and oceanic fluid dynamics[END_REF]. The mathematical characterisation of turbulence is ultimately shaped by the Navier-Stokes equation.

Navier Stokes Equations

In here we summarize the importance of the Navier-Stokes set of equations to the theory of turbulence. These, if solved, can exactly describe the motion of any incompressible Newtonian fluid:

∂ ∂t u + u • ∇u = - 1 ρ ∇p + ν∇ 2 u + f (2.9) ∇ • u = 0, (2.10)
where u is the fluid speed, p its pressure and f represents the external body or volume forces acting on the fluid. Equation 2.10 is the continuity equation representing the mass conservation of the incompressible fluids and equation 2.9 is the fluid's momentum equation. In there we have represented the inertial or non-linear term (u • ∇u), responsible for the transfer of kinetic energy, the pressure gradients (∇p/ρ) and the viscous dissipation of energy (ν∇ 2 u). One can take advantage of the Reynolds number (Eq. 2.8) to make equation 2.9 non-dimensional:

∂ ∂t u + u • ∇u = -∇p + 1 Re ∇ 2 u + f. (2.11)
The Navier-Stokes (NS) equations require boundary and initial conditions to be solved, as any partial differential equation. An idealized but practical choice is to adopt periodic boundary conditions to simulate the NS equations in an infinite space (such as u(x, y, z) = u(x + aL, y + bL, z + cL), where L is the simulation box size and (a, b, c) integers, this comes with the advantage of maximum symmetry. Despite of the choice of boundaries, it is in symmetry that relies the onset of turbulence from the NS equations point of view, more specifically, in the breaking of symmetries. NS equations are characterized by multiple symmetries; space translations, time translations, scale invariance, between others. As the Reynolds number of a fluid increases symmetries are broken, laminar flows break into vortices (Fig. 2.5), fluctuations in velocity increase as we look into smaller scales (becoming scale dependent), until we have fully developed turbulence. That said, for large enough Re symmetries are statistically restored at sufficiently small scales far from boundaries. So a good definition is: turbulence is characterized by the presence of significant scale separation between energy injection and dissipation, with injection at large scales and dissipation at small scales for three dimensional (3D) turbulence.

Energy Cascade and K41 Theory

The statement about energy injection and dissipation above was coined by Lewis Richardson in his first book in 1922 [START_REF] Richardson | Weather prediction by numerical process[END_REF], there he proposed that turbulent kinetic energy is transferred from large to small eddies, until it is consumed by viscous dissipation (figure 2.6). Even though figure 2.6 provides us a good illustration of the process, in reality eddies can vary a lot in shape and the smaller ones are actually embedded in the larger ones. In spite of that, the idea of spectral energy cascade is at the heart of our present understanding of 3D turbulent flows. We thus have to think about turbulence as an intrinsically multi-scale system where energy is contained in the large scale eddies that, via non-linear processes, generate smaller eddies that keep reducing in size until they are dissipated by viscous processes.

The work of Richardson was almost ignored until 20 years after when the idea of a spectral cascade took quantitative shape in the hands of Kolmogorov and Obukhov in Russia. Kolmogorov, generally regarded as the greatest probabilist of the 20th century took advantage of the unpredictable behaviour of turbulence to express its general attributes in a statistical frame. The Russian mathematician hypothesized that the statistics of small scales are isotropic and depend only on two parameters, viscosity ν and the energy flux [START_REF] Nikolaevich | The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers[END_REF]. Following [START_REF] Vallis | Atmospheric and oceanic fluid dynamics[END_REF] we can take advantage of Kolmogorov's work to better understand the energy cascade. If we decompose a flow velocity field into Fourier components in a finite domain:

u(x, y, z, t) = ∑ k x ,k y ,k z ũ(k x , k y , k z , t)e i(k x x+k y y+k z z) .
(2.12)

Where ũ is the Fourier transformation of the velocity field u, k(k x , k y , k z ) = 2πn x,y,z /L x,y,z , with n x,y,z as integers and L x,y,z the domain size. The kinetic energy in the fluid (per unit mass) is given by:

Ê = 1 2V |u| 2 dV = 1 2 ∑ | ũ| 2 dk,
(2.13) using Parseval's theorem, where Ê is the total energy and E is the energy per unit of volume V. This leads to:

Ê ≡ E(k)dk, (2.14)
where E(k) is the energy spectral density, or energy spectrum. Here we represent homogeneous isotropic turbulence (statistically invariant under translations and rotations of the reference frame). The statistical characteristics of the flow are described by the spectrum E(k):

E(k) = 1 2 < | ũ(k)| 2 > |k|=k . (2.15)
Here <> |k|=k is averaged over spherical shells of radius |k| = k.

The relation 2.15 is valid in a range where forcing and dissipation are negligible, in this range one can then expect that E(k) = f ( , k). Also here energy is passed from large scales to small scales unaltered, so the energy injection rate i (at large scales) is equal to the energy flux (from large scales to small scales) and the energy dissipation rate ν (at small scales). Dimensional reasoning can then be made here, if E(k) has dimensions (lenght 3 /time 2 ) and has dimensions (length 2 /time 3 ) then

E(k) = C k 2/3 k -5/3 , (2.16)
where C k is a dimensionless constant. The range where this relation is valid would be later baptized as the inertial range. The result was exposed by Kolmogorov in 1941 and would present such importance in the following era in the turbulence field that even the abbreviation K41 became of common use. The inertial range came to complete the energy spectrum cascade picture in between the energy injection and the viscous dissipation (figure 2.7 ). We can also derive the scale that separates the inertial range from the dissipation range, the Kolmogorv scale η, by considering the eddy turnover time τ k . This is the time needed for a parcel with velocity u k to move a distance 1/k, where u k is the velocity associated with the scale k = 1/L. On dimensional consideration u k = (kE(k)) 1/2 so that:

τ k = (k 3 E(k)) -1/2
(2.17) and so, Energy is injected at rates i and transferred at rates from the injection range to the dissipation range (where it is dissipated at rates ν ). According to Kolmogorov i = = ν , this leads to the 'Kolmogorov's -5/3 spectrum'.

τ k ∼ -1/3 k -2/3 . (2.18)
Representing the inertial time scale. At some small length scale viscosity becomes important and the inertial range scaling we presented will fail to predict the proper energy spectrum. The reason for that is because in the inertial range viscosity is ignored since the timescales on which it acts are too long for consideration, it does becomes important when τ k ∼ 1/(k 2 ν) and in combination with Eq. 2.18 we get:

k η ∼ ν 3 1/4 , (2.19
) .20) where Eq. 2.19 represents the dissipation wavenumber. Therefore for k << k η inertial effects dominate and for k >> k η viscous dissipation dominates. The number of degrees of freedom (number of scales) is given by:

τ η ∼ ν 1/2 . ( 2
L 0 /η = (Re) 3 4 , (2.21)
Re → ∞ illustrates the multiscale nature of the phenomenon.

These results were extensively checked and reproduced by laboratory and numerical experiments [START_REF] Frisch | Turbulence: the legacy of A N Kolmogorov[END_REF]. They thus constitute a kind of 'boundary condition' on theories of turbulence: such theories, to be acceptable, must either satisfy them, or explicitly violate the assumptions made in deriving them. Moreover, this so found scaling characteristics of turbulence are also of great value for the understanding of Chapter 2. Turbulence and Turbulent Transport dispersion processes as we will see in the next session. This allows us to make assumptions that permit us to rely on Lagrangian dynamics, which are simpler to deal with in contrast to the Navier-Stokes equation.

Cellular Flow

A cellular flow is described by the presence of cellular motion usually in the form of vortices of similar size, a visualization is available in Fig. 2.8. The cellular flow mathematical form was originally conceived to describe the trajectory of particles of different settling velocities in thermal convection cells [START_REF] Stommel | Trajectories of small bodies sinking slowly through convection cells[END_REF], and it has the form (in 2 dimensions):

Ψ(x, y) = Ψ 0 sin(kx)sin(ky), (2.22)
where k = 2π/L, with L as the vortices sizes. As a consequence of incompressibility the velocity field is:

u x = -∂ y Ψ ; u y = ∂ x Ψ (2.23)
Later in this chapter we will explain how a flow field of this form allows for a simplified description of buoyancy and wind-driven currents [START_REF] Thorpe | The turbulent ocean[END_REF] such as convective currents and Langmuir circulations, often encountered in the upper layers of oceans and lakes. But first it is interesting to notice that if we force the NS equations with a forcing term f (x, y) ∼ sin(kx)sin(ky) we obtain a velocity field corresponding to the streamfunction of Eq. 2.22. The procedure can be demonstrated by expressing the NS equations through vorticity (ω):

ω = ∂ x u y -∂ y u x (2.24)
So taking the curl (z component) of Eq. 2.9 and considering the 2D incompressible flow, ω(∇ • u) = 0, leaving us with:

∂ ∂t ω + (u • ∇)ω = ν∇ 2 ω + f ω (2.25)
where

f ω = (∇ × f) z = ∂ x f y -∂ y f
x is the vorticity forcing. Finally, if the velocity field has the form of Eq. 2.23, u = (u x , u y ) and we see that (u x ∂ x ω + u y ∂ y ω) = 0, we also see that the vorticity is not time dependent ∂ t ω ∼ 0 and we are left with:

f ω = -ν∇ 2 ω = 2νk 2 ω. (2.26)
The result shows that the vorticity forcing has the same functional form as the vorticity field. The same type of proportionality relation holds for the forcing f and the velocity field u from the original Navier-Stokes equation. We obtain the numerical solution of the NS equations forced by a cellular forcing via Lattice-Boltzmann method [START_REF] Calzavarini | Eulerian-Lagrangian fluid dynamics platform: The ch4-project[END_REF] and plot the velocity field in Fig. 2.8 .

The velocity field in this case was stable and time independent. It is known from multiscale analysis that the 2D cellular flow is linearly unstable with respect to the large scale perturbation, meaning much larger than 1/k [START_REF] Sivashinsky | Negative viscosity effect in large-scale flows[END_REF]. This solution is in agreement with what we obtained later on in Chapter 5 with a kinematic method for the solution of our model. As we will further discuss, our approach did not rely on the solution of the NS equations but on a Lagrangian procedure for the calculations of phytoplankton concentrations in flow fields. This was due to the high numerical cost of adding small scales in direct numerical simulations using a realistic domain size, relevant for phytoplankton in oceans and lakes.

Transport and Dispersion in Fluids Flows

Here we discuss transport in fluids in regards to the choice of approach (Eulerian/Lagrangian), and we disclose the characteristics of the different transport regimes for passive tracers inside turbulence . Passive tracers, in opposition to active ones, have no influence in the properties of the fluid flow velocity field.

Eulerian and Lagrangian Approaches

The study of scalar fields is a topic of great interest in physics and engineering. From the study of aerosol dispersion in the atmosphere [START_REF] S P Arya | Air pollution meteorology and dispersion[END_REF] to the vertical stratification of temperature in lakes and oceans [START_REF] Vallis | Atmospheric and oceanic fluid dynamics[END_REF], the number of problems that fit the classification is enormous. These passive fields are described by defined values at every point in a region. The Eulerian approach concerns with the the evolution of a scalar field θ(x, t) at each point in space and at each instant in time t (Eq. 2.27), this typically contain a molecular diffusion term D and a velocity field u(x, t) [START_REF] Bohr | Dynamical systems approach to turbulence[END_REF] which may be given or dynamically determined by the Navier-Stokes equations. On the other hand, the Lagrangian approach deals with the trajectory of particles released into the fluid by following individual particles positions x i in time. Analogous to the Eurelian frame these particles are subject to advection, but this time by a velocity field u(x i , t) measured at the particle position. The Langevin equation 2.28 representing this approach also accounts for diffusion by using a random vector ζ(t) (δ-correlated in time, with unity variance and zero mean) modelling the distribution of equation 2.5.

∂θ(x, t) ∂t = D∇ 2 θ(x, t) -∇.(u(x, t)θ(x, t)), (2.27)
Chapter 2. Turbulence and Turbulent Transport

dx i dt = u(x i , t) + √ 2Dζ(t).
(2.28)

Note that although the equations presented are useful as an illustration of the approaches used in fluid mechanics, the Eulerian counterpart for the Langevin equation corresponds to a Fokker-Plank equation for the probability density function of finding particle i at position x i and time t. One can study scalar fields using Eulerian and Lagrangian approaches, choosing between one or both according to precision and convenience. In summary, either they deal with a fluid flow field in the space and time domains (Eulerian) or they deal with the trajectory of each fluid particle in time (Lagrangian). Furthermore, it is common in both experimental and numerical approaches to study fluid flow properties trough tracers. Tracers are point-like particles with density and velocity equal to the fluid. They may also be divided in: Non-reactive tracers which are not subject to any chemical or biological reaction and are usually employed in the study of fluid motions, as they behave in the exact same manner as the medium [START_REF] Sandulescu | Kinematic studies of transport across an island wake, with application to the Canary islands[END_REF], or reactive tracers, which are tracers that undergo some type of chemical or biological reaction during their dynamics in a flow-field. Examples go from the study of chemical reactions to the phytoplankton dynamics we are presenting in this thesis. The reaction characteristic here also means that, in contrast to the passive tracers, the interest in the reactive ones is in the interplay between transport and reaction dynamics, presenting a further degree of complexity.

Relative Dispersion in Turbulent Flows

Our interest for the Lagrangian approach comes from the fact that relative dispersion in turbulent flows can be conveniently used to explore transport and mixing properties of turbulent flows. This approach has been widely utilized by authors in the fields of mixing, transport, oceanography and meteorology [START_REF] Lacorata | Lagrangian predictability characteristics of an Ocean Model[END_REF][START_REF] Prants | Chaotic Lagrangian transport and mixing in the ocean[END_REF][START_REF] Boffetta | Nonasymptotic properties of transport and mixing[END_REF][START_REF] Lacorata | 3D chaotic model for subgrid turbulent dispersion in large eddy simulations[END_REF], particularly through the study of tracers dispersion. In the following the main theoretical expectations for relative dispersion will be discussed with reference to a turbulent flow in homogeneous isotropic conditions. The separation distance at time t between a pair of particle is given by R(t) = |x 2 (t)x 1 (t)|, where x 1,2 represents the position of particles 1 and 2. Chaos is a common manifestation of nonlinear dynamics associated with a strong sensitivity to infinitesimal errors on the initial conditions [START_REF] Lacorata | 3D chaotic model for subgrid turbulent dispersion in large eddy simulations[END_REF]. Therefore, chaotic trajectories are to be expected in a turbulent flow assuming an initially infinitesimally small difference in the positions of two particles, relative dispersion should be exponentially felt in time:

R 2 (t) ∼ R 2 (0) e 2λt , (2.29)
where λ is known as the maximum Lyapunov exponent (MLE). However, turbulence is characterized by a multiplicity of scales, and the expected behavior for R 2 (t) depends on how the initial separation compares to the turbulent length scale.

Recalling that the length scales η (the Kolmogorov scale below where energy dissipation is dominating and the flow becomes smooth) and L 0 (where energy is injected into the system) are respectively the smallest and the largest ones, three regimes are expected. The first regime is the exponential one described by Eq. 2.29, which is valid when typical separation is small, R 2 (t)

1/2 η. 2.5. Turbulent Diffusion If R 2 (t) 1/2
L 0 , instead, we expect a diffusive regime, particle velocities are practically uncorrelated and their trajectories perform diffusive motion:

R 2 (t) 2D E t, (2.30)
where D E is the effective diffusivity, dominated by the turbulent contribution, which is much larger than the coefficient of molecular diffusivity D. In other words, only when the distance between two particles is sufficiently larger than the correlation length of the velocity field, the relative velocity can be approximated by a stochastic process. Between the diffusive and the chaotic regime, in the range of scales corresponding to the turbulent inertial range, relative dispersion is expected to have a power law dependence on time in the presence of a direct energy cascade from large to small eddies [START_REF] Bohr | Dynamical systems approach to turbulence[END_REF], hence:

R 2 (t) = C R t 3 , (2.31)
where is the energy dissipation rate and C R is the non-dimensional Richardson constant. This equation is known as Richardson law and it describe the superdiffusive behaviour of turbulent dispersion. The regimes presented here are not by coincidence analogous to the ones presented in the past section as dispersion is another way of expressing the Richardson cascade [START_REF] Frisch | Turbulence: the legacy of A N Kolmogorov[END_REF].

Turbulent Diffusion

As we saw in the previous section, when we examine the tracer dispersion phenomenon at sufficiently large scales the general behavior is equivalent to isotropic diffusive motion. Indeed this relation has been widely used in geophysical and biological models under the name of turbulent diffusion or eddy diffusion to approximate the complex dispersion generated by turbulence eddies over large scales of space and time. In spite of that, the key factor here is that molecular diffusion is a process of random motion that affects particles at very small scales, typical mean time scales for molecular motion are 10 -10 s for air and 10 -12 s for water, and typical length scales are 10 -7 m for air and 10 -9 m for water [START_REF] Okubo | Diffusion and ecological problems: modern perspectives[END_REF]. Since these values are considerably smaller than our typical time and length scales of observation, molecular diffusion can be safely considered to correspond to a random motion. However, the typical time and length scales of turbulence in oceans and lakes are of the same order of magnitude as typical time and length scales of observation. For this reason we have no guarantee that a simple analogy with the random walk process is applicable to turbulent flows.

In addition to that, in geophysical or biological fluid models one wants to represent the macroscopic effect of turbulence as a diffusion process through the use of a diffusivity D E = D T + D (that can be space dependent as well), in the same manner as in the diffusion equation. Nonetheless, in terms of macroscopic dispersion turbulent diffusion D T is much more efficient than the molecular diffusion as it is associated with a flow of much larger scales. Due to the wide range of eddy scales involved, D E should vary with the considered scale . Estimating the required eddy diffusivity K a may not be obvious [START_REF] Gargett | Vertical eddy diffusivity in the ocean interior[END_REF][START_REF] Haskell | Estimates of vertical turbulent mixing used to determine a vertical gradient in net and gross oxygen production in the oligotrophic South Pacific Gyre[END_REF], particularly for what concerns its scale dependence [START_REF] Richardson | Weather prediction by numerical process[END_REF]. Perhaps the most commonly accepted measure of effective diffusivity comes in terms of tracer concentration variance σ 2 over time K a ≡ (1/2)dσ 2 /dt.

Chapter 2. Turbulence and Turbulent Transport

This was first measured by Richardson in 1926 and repeated by Okubo [START_REF] Okubo | Oceanic diffusion diagrams[END_REF] (Fig. 2.9) and it provides evidence that K a scales with the scale of interest as 4/3 . In other terms, K a , as estimated from Eulerian measurements, scales as the relative diffusivity that can be estimated from Lagrangian measurements, (1/2)d R 2 /dt. This result means that as becomes larger more and more eddies participate in diffusion and the effective diffusivity increases. The 4/3 law does not always apply but it is a good general approximation. FIGURE 2.9: Okubo measures properties of ocean diffusion through the use of dye experiments. He thus verified that the 4/3 law is locally valid, given that the rate of turbulent energy dissipation varies with the time and length scales in which one is interested. Figure extracted from [START_REF] Okubo | Oceanic diffusion diagrams[END_REF].

Another way to look into an expression for the diffusivity is to examine eddy fluxes in a tracer concentration field c [START_REF] Williams | Ocean dynamics and the carbon cycle: Principles and mechanisms[END_REF]. The eddy flux is the average of the product of fluctuations in both the concentration field c and the flow velocity field u, and it relates to the eddy/turbulent diffusivity tensor as:

D (i,j) T ≡ -F (i) ∇ (j) < c > = - < |u (i) c | > ∇ (j) < c > , ( 2.32)
where F (i) represents the i th component of the eddy flux and the primes fluctuations with respect to the mean in their respective fields. This approach can be used for estimations of D T but it presents clear disadvantages. First one has to have access to sufficiently resolved measurements of the tracer concentration and the velocity field, which may not be easy to obtain, and second one has to count on the concentration field being sufficiently stratified in space as if not the measurement of D diverges. Alternatively historically scientists estimate diffusivity via dye dispersion experiments as K a = (1/2)dσ 2 /dt [START_REF] Okubo | Oceanic diffusion diagrams[END_REF] where an initial portion of dye can be measured multiple times as it spreads, which provide good estimations for small scale experiments but it is obviously unfeasible for large length scales. Still the most applied technique in measuring diffusivity are experiments with drifters where one can estimate diffusivity either by measuring their dispersion or eddy kinetic energy for example. The advantage is that one can follow drifters for many kilometers, allowing large domains to be explored, on the other hand these are usually floaters and therefore unusable for vertical diffusivity measurements. In fig. 2.10 we see an example of an apparent horizontal diffusivity map made from eddy kinetic energy measurements through drifters. Moreover, the study of vertical diffusivity in the oceans is considerably less explored than lateral diffusivity given the increased difficulty to perform such experiments. These usually require scientific boat expeditions for many hours in multiple periods of time for proper water probing with depth. In addition to the arguments presented in the introduction where we expressed the convenience of ocean surface monitoring through satellite images, numerous experiments with floating drifters deployed by research institutes are available for scientists, whereas constant time series depth probing experiments are significantly more difficult to find. Nevertheless, we can mention some comparisons, for example, horizontal diffusivity is several orders of magnitude larger than the vertical diffusivity. In many cases the influence of the fluid density distribution is more important in determining the vertical diffusivity than the scale of the phenomenon. In the next section we will discuss some examples of physical vertical mixing processes in oceans and lakes, there we can better clarify the problem of anisotropic turbulence in vertical mixing processes.

Vertical Mixing Processes in Oceans and Lakes

It is beyond the scope of this work to describe in detail the intricate complexity of vertical mixing processes in aquatic environments. Instead, we will present a short summary of the theory behind some of these mechanisms with the goal of better relating the numerical simulations and observations we will later expose in this thesis to what is known for these processes from the theoretical and experimental literature and observations. The number of processes that govern vertical movement in the ocean is very numerous. Fig. 2.11 provides a representation of some small and large scale ocean flow processes. The physical phenomena that originates them is also vast, to cite some; thermal processes, topological features, mechanical forcing (wind for example), salinity gradients... In this thesis we are concerned with the ones that occur in length scales < 1km as this scale coincides with what can be expected for the mixed layer depth. Even though large geophysical currents can affect ocean dynamics locally [START_REF] Williams | Ocean dynamics and the carbon cycle: Principles and mechanisms[END_REF], we will abstain from talking about them as it is not the objective of this work. That said, we can divide the physical process that originates the fluid flow in the aforementioned range into two categories: wind driven currents, and buoyancy driven currents (both effects are often concomitant, so the category is defined by the dominant one).

Buoyancy driven currents

The density of fluids has fundamental effects on their dynamics and also on the genesis and nature of turbulence. Forces can originate from the variations of density in fluids lying within gravitational fields. Differences in density lead to pressure differences that, if unopposed, drive movement. These forces are named buoyancy forces and are linked to temperature or salinity gradients, or driven by heat fluxes at the boundaries of fluids, they are the cause of a wide range of turbulence related phenomena [START_REF] Thorpe | The turbulent ocean[END_REF]. The buoyancy force on a water parcel is determined by its buoyancy anomaly:

b = -g (ρ -ρ amb ) ρ 0 , (2.33)
where g is the acceleration of gravity, ρ the parcel density, ρ amb the surroundings density and ρ 0 is a constant reference density equal to 1000kg/m 3 [START_REF] Marshall | Open-ocean convection: Observations, theory, and models[END_REF]. In fluid mechanics one usually associates the buoyancy term to the Boussinesq approximation to the Navier-Stokes equation, which allows for the calculation of the dynamics of a fluid velocity field by exploring smallness of density variations in many liquids.

In aquatic environments complications arise because the density of water can depend on temperature T, salinity S and, although very weakly, on pressure p, i.e. ρ = ρ(T, S, p). However, a simplified equation of state can be adopted:

ρ = ρ 0 [1 -α T (T -T 0 ) + β S (S -S 0 )], (2.34)
with α T and β S being respectively the thermal expansion and haline contraction coefficients, and T 0 and S 0 are reference temperature and salinities. In here pressure can be neglected due to the almost incompressible nature of water. The equation is not valid for low temperatures as these coefficients will not be constant anymore and will vary with T and p [START_REF] Marshall | Open-ocean convection: Observations, theory, and models[END_REF] (for example the sign of the thermal expansion coefficients reverts for temperatures below 4 o C in unsalted water). Density gradients will typically be the origin of convection currents, although wind forcing also has a important role in this.

An example of deep sea convection can be seen in Fig. 2.12, here one can see that when an initially stratified fluid is subject to buoyancy loss from the surface its surface layer will become denser sinking and driving convective currents in the form of circular motions. These are similar to eddies, sinking denser water parcels from the surface and restoring its buoyancy by bringing lighter waters from the bottom. This persistent structures are one of the motivations of the numerical simulations that will be exposed in this thesis. This phenomenon can be observed in oceans [START_REF] Backhaus | Phyto-convection': the role of oceanic convection in primary production[END_REF] and lakes [START_REF] Fer | Winter cascading of cold water in Lake Geneva[END_REF] and its origin is mostly dominated by meteorological forcing via temperature gradients generated by seasonal cooling. As atmosphere temperature drops the surface of water cools and loses buoyancy in relation to deeper layers, driving convective currents in the process. It has been even hypothesized that winter convection is crucial for the survival of phytoplankton in mid and high latitudes by allowing resting spores and vegetative cells from deep ocean waters to be uplifted [START_REF] Backhaus | Phyto-convection': the role of oceanic convection in primary production[END_REF]. This comes to add to the variety of ways in which advection by currents and dispersion by turbulence can affect biological dynamics in the sea.

Wind driven currents

Rarely you will find a aquatic environment in which wind does not plays a role in mixing dynamics (Fig. 2.13 as example), that said, the quid pro quo in between wind forcing and fluid dispersion is not at all obvious. Starting from what the eye can see, when touching the first layer of water turbulent winds tend to produce breaking surface waves (Fig. 2.11). Even tough these play a role in surface dispersion, At low wind speeds the sea/lake surface may behave, over much of its area, like a viscous-conductive boundary layer above a rigid surface, and molecular processes may dominate the exchanges of some properties across the air-water interface in sufficiently low wind speeds (typically those less than about 3m/s) [START_REF] Thorpe | The turbulent ocean[END_REF]. Laboratory experiments confirm that as much as 40% of the energy lost by breakers may go into the production of bubbles [START_REF] Lamarre | Void-fraction measurements and sound-speed fields in bubble plumes generated by breaking waves[END_REF], these are mostly associated with higher turbulence zones than non bubble zones. But as much as breaking waves are associated with turbulence, and help in the maintenance of the mixed layer in combination with buoyancy processes, their penetration in depth is limited by the viscous nature of the aquatic environment. Between many processes that participate in the production and regulation of the mixed layer we can cite also shear, convection and Langmuir circulation.

While flying in an aircraft over the sea, or while one observes a lake from surrounding mountains one can probably observe parallel streaks of foamy water, 30-50 meters apart, extending hundreds of meters across the water (Fig. 2.14). These are a direct result of Langmuir circulation, a feature first observed by Langmuir [START_REF] Langmuir | Surface motion of water induced by wind[END_REF]. The phenomenon originates from the drag between the surface water and steady wind, giving rise to a converging cellular motion in a plane perpendicular to the wind direction (Fig. 2.11). Beneath them, the downward vertical current speed increases with depth, reaching a maximum of about 1-20 cm/s at a depth of 0.2-0.5 times the mixed layer depth, before decreasing towards the base of the mixed layer [START_REF] Thorpe | The turbulent ocean[END_REF]. The beginning of cellular circulation in a previously quiescent fluid was observed to occur at Re = 530 ± 20 and, after the first appearance of surface waves, at Re = 370 ± 10, where Re = U 0 (t/ν) 1/2 , U 0 is the surface water speed and t is the time since the onset of wind [START_REF] Veron | Experiments on the stability and transition of wind-driven water surfaces[END_REF]. Langmuir circulation and its associated turbulence, as well as the mentioned buoyancy and wind driven currents, play a pivotal role in both climate and biological processes, regulating the mixed layer temperature In this chapter we presented a summary for the theories of diffusion, turbulence, turbulent, diffusive and advective transport, as well as some real flow phenomena in aquatic environments. We proceed this way because later on, when we present our advection-reaction-diffusion model for vertical phytoplankton dynamics, we hope that the reader can understand the origin of the necessary assumptions needed to build the fluid transport part of the model. As well as the importance of the different time and length scales that one tries to replicate in order to bring as much realism as possible to our idealized system.

Chapter 3 Population Dynamics, Plankton and Phytoplankton

Population dynamics problems date from even before the the sketches of Da Vinci seen in Fig. 2.1. Origins can be tracked to a book by Fibonacci in 1202 entitled Liber abaci (Book of Calculation) in which he describes a recurrence problem where a man wants to know how many rabbits he will have in a year after buying a pair [START_REF] Bacaër | A short history of mathematical population dynamics[END_REF]. Population dynamics is at the intersection of various fields: mathematics, physics, social sciences (demography), biology (population genetics and ecology) and medicine (epidemiology).

In this chapter we well briefly review basic biological concepts of plankton and phytoplankton which are necessary for understanding, both the importance of studying these organisms, and how to properly model their life cycle dynamics in a mathematical/numerical framework. First we introduce some basic understanding regarding plankton and phytoplankton biology, then we briefly review the basic knowledge behind population modelling problems and situate light-limited growth phytoplankton models in this vast field. We finish the chapter by exposing the mathematical framework and history behind phytoplankton bloom studies, concluding the chapter with an overview of some of the knowledge acquired in the past decades about the subject trough modelling.

Plankton and Phytoplankton

Plankton

The word Plankton is a classification based on the Greek word planktos meaning errant (drifter). This refers to the fact that planktonic life is considered to be composed of organisms that are not able to actively control their location but are rather advected by the surrounding medium (unable to overtake the currents) [START_REF] Schmitt | Turbulence et écologie marine[END_REF]. The classification can be debated as in fact many forms of planktonic life like the zooplankton or the phytoplanktonic dinoflagellates are motile species, able to swim (mostly vertically) hundreds of meters daily. Regardless of the motility of some species, their movement is nonetheless dominated by current flows and turbulence. What is remarkable about plankton is their diversity (Figure 3.1), the classification embraces organisms from marine viruses smaller than 0.2 µm to the amphipod of tenths of centimeters. In this size range we also find bacteria, algae, crustaceans and jellyfish. Each possessing different habits and roles in the maintenance of aquatic biodiversity. We usually separate the larger species that predate on the primary producers as grazers or zooplankton (zoo standing for "animal") these will Chapter 3. Population Dynamics, Plankton and Phytoplankton ultimately sustain fish and mammals at the top of the food chain. The remaining planktonic species either receive their energy from chemical reactions through organic/inorganic compounds and/or photosynthesis, these being referred to as phytoplankton (the prefix phyto meaning "plant").

Phytoplankton

Phytoplankton is a generic name for photosynthesizing microscopic organisms that inhabit the upper sunlit layer of almost all oceans and bodies of freshwater. These tiny plants, some as small as a micron, use chlorophyll to capture energy from sunlight to create the organic molecules that make up their bodies, a reaction that turns them into the fuel that is imperative for sustaining of the aquatic food web. In other words they are known to control biogeochemical cycles through primary production and are believed to contribute for about half of the global primary production [START_REF] Tréguer | Silica control of carbon dioxide[END_REF]. They are also refereed to as algae (all algae is phytoplankton but not the inverse). Phytoplankton are mostly too small to be seen by the naked eye, this means that in a sunlit surface of water every millilitre may contain hundreds of thousands of phytoplankton [START_REF] Williams | Ocean dynamics and the carbon cycle: Principles and mechanisms[END_REF], and therefore collectively their chlorophyll colors the water in a way that it is easily spotted from cruise ships to satellite images. Figure 3.2 shows some satellite images from NASA's OceanColor website, in there we see a set of different physical conditions to which phytoplankton may be subjected, these range from advection by mesoscale ocean vortices and mixing by turbulent eddy fields, to increase in population due to eutrophication (enriched nutrient zones).

The number of phytoplankton species that exist in the ocean are very numerous (hundreds of thousands), between those the most common are the diatoms, then the cyanobacteria and dinoflagellates. Diatoms refer to a group of phytoplankton that cannot swim and thus rely on vertical turbulence to oscillate between sunlit and nutrient rich zones. The group nonetheless is composed of hundreds of thousands of species of different shapes. Cyanobacteria are bacteria found in almost every aquatic environment on earth, from moistened rocks in the desert to rocks in the Antarctica. They sometimes are described as dangerous because their algae blooms can release cyanotoxins that may poison and even kill animals. Both groups rely solely on inorganic nutrients and photosynthesis to survive. Dinoflagellates on the other hand may combine photosynthesis with grazing, though not all follow this rule. These organisms differ from the past groups also in the fact that they posses two flagella, an appendage that allows them to swim according to survival strategies (flee from predators for example). It is far from the scope of this work to cover the many features of phytoplankton species like sizes, shapes and motion techniques, characteristics that show importance in very small length and time scale scale dynamics [START_REF] Borgnino | Scale-dependent colocalization in a population of gyrotactic swimmers[END_REF][START_REF] Borgnino | Gyrotactic swimmers in turbulence: shape effects and role of the large-scale flow[END_REF], therefore we focus on describing the general properties of these organisms survival dynamics.

Phytoplankton Life Cycle and Turbulence

Phytoplankton are at the bottom of the food chain: they create fresh organic matter from dissolved nutrients, carbon dioxide, and energy from sunlight. Unfortunately for them nutrients and sunlight are in opposite places. While the sunlit layer of the water (euphotic zone) can extend up to 200m, depending on how clear the aquatic medium is, usually it only provides enough energy for survival on a fraction of this (first 20 meters for clear water for example). Nutrients on the other hand (with the exception of coastal zones) are localized in deep ocean (Figure 3.3), and so it comes that what allows these organisms to thrive in the way we see in satellite images is turbulence, and more specifically the fluid mixing associated with it. Phytoplankton survival in aquatic environments depends mostly in being trapped in a layer of active movement where the almost homogeneous mixing in the water column will allow them to access both the euphotic zone and nutrients, either by bringing the nutrients to them or the inverse, this is known as the mixed layer [START_REF] Williams | Ocean dynamics and the carbon cycle: Principles and mechanisms[END_REF]. The definition of mixed layer is independent on the existence of organisms or sediments, the mixed layer is present in all oceans and most bodies of water and it possesses a seasonal evolution that is associated to its genesis mechanisms due to thermal and mechanical processes. In the north Atlantic, for example, the mixed layer known to shoal in spring, an event historically associated with phytoplankton spring blooms, and then regain depth in winter, when phytoplankton populations get smaller. Turbulence can nonetheless have positive or negative influence over these survival dynamics as we will see. Apart from what it was mentioned, the influence of turbulence on plankton can be summarized according to [START_REF] J C Prairie | Biophysical interactions in the plankton: A cross-scale review[END_REF]:

• Turbulence increases turbulent diffusion of nutrients, hence contributing to population growth.

• According to the intensity of turbulence one species or another is favored, defining the composition of plankton communities.

• Turbulence influences in sedimentation defining the types of particles in suspension.

• It transports zooplankton having a positive or negative effect on the predator/prey and reproductive dynamics according to the rate of encounters it generates.

In addition to these, many, if not most plankton species are denser than water [START_REF] Huisman | How do sinking phytoplankton species manage to persist?[END_REF], and so they tend to sink. This means that these species depend on turbulence to rise to the sunlit zones (balance their weight), resulting that in the absence of turbulence they would not have formed a successful part of the phytoplankton community in lakes and oceans for millions of years. The role of turbulence will be further explored in the following sections.

Chapter 3. Population Dynamics, Plankton and Phytoplankton

Biologically the two factors that play a role in phytoplankton survival are photosynthesis and respiration. Photosynthesis is the conversion of electromagnetic energy into chemical energy through the formation of organic molecules [START_REF] Williams | Ocean dynamics and the carbon cycle: Principles and mechanisms[END_REF]:

6CO 2 + 6H 2 O photons ----→ C 6 H 12 O 6 + 6O 2 (3.1)
The chemical reaction turns carbon dioxide and water into glucose and molecular oxygen. The creation of organic matter from inorganic compounds is called primary production. The reaction is exactly what makes phytoplankton so valuable for life in planet earth, they not only contribute to about half of the oxygen that allows us to breath but also assimilate the carbon dioxide in a process that helps regulate the climate, and, as a bonus, they generate huge quantities of organic matter that compose the base of the ocean's food chain. The reverse of this process is the respiration:

C 6 H 12 O 6 + 6O 2 → 6CO 2 + 6H 2 O + energy (3.2)
By consuming part of the glucose and oxygen generated by photosynthesis, plankton (and not only phytoplankton) can produce the energy that fuels biochemical reactions and mechanical work. This process is one of which that can make plankton blooms dangerous as the extensive consumption of oxygen by a plankton patch can generate "dead" anoxic (depleted of oxygen) regions in oceans and lakes [START_REF] Pitcher | Suffocating Phytoplankton, Suffocating Waters-Red Tides and Anoxia[END_REF]. In terms of phytoplankton blooms photosynthesis is considered a production process, since it generates organic matter, whereas respiration is considered a loss process, since it consumes. In the ensemble of loss processes we can also add the natural degradation of organic matter and predating as processes that reduce and can extinguish populations. Later in this chapter, when we describe the modelling of phytoplankton population dynamics, we will further explore the processes that affect the life cycles of plankton and phytoplankton. In the last section of this chapter we will make a summary of some of the models and theories that were coined to represent phytoplankton dynamics in oceans and lakes. These models approach the representation of plankton's life cycles from a number of different parameters; from nutrients and radiation availability to turbulence and water column depth. We start from the most minimalist of models and from that we build the complexity of the current and most refined models.

Modelling of Phytoplankton Populations

The value of mathematical treatment in ecosystem dynamics is the same as in classical mechanics or electromagnetism for example, quantitative analysis and predictions are indispensable. The difference arises in the fact that population models cannot be compared to Newton's laws or Maxwell's equations, the natural environment has not yet found a firm mathematical basis akin to these. In the words of Okubo and Levin [START_REF] Okubo | Diffusion and ecological problems: modern perspectives[END_REF]:

"The mere fact that a mathematical model agrees well with a small amount of data does not suffice, as the agreement could be coincidental. Moreover, models should not be confused with fundamental equations or laws. Only those hypotheses that have withstood large amounts of critical scrutiny can be elevated to the status reserved for laws."

That said, it is expected that a complex field such as ecology will take a very long time to reach such a goal, making mathematical models an important step in that direction. In this sense one can crudely separate models into stochastic models or deterministic models, the first tend to better capture the probabilistic and discrete nature of biological processes, and the second are useful in the creation of a general understanding, particularly at the collective level, i.e. as far as populations are concerned (here we will only focus on deterministic models). A straightforward way to evaluate this is to construct a mathematical model that describes a generic population increase/decrease. According to deterministic theory, the rate of population change, dN/dt, is proportional to the number of individuals at time t, N(t), and to a coefficient of net growth (birth minus death) r:

dN dt = rN. (3.3) If N(t = 0) = N 0 , than N = N 0 e rt , (3.4)
meaning that if r is constant we will get the same result for each experiment in a deterministic fashion. Regardless of how good of a solution Eq. 3.4 may turn out to be for the macroscopic behavior of said population, the process of population growth is usually stochastic, and at a given instant the growth will most likely not depend solely on rN. Instead we can express our solution in terms of the expected value, N:

N = ∞ ∑ N =0 N p(N , t, a, b, ...) = N 0 e rt , (3.5)
where, p is a growth probability density function that may depend on numerous factors.

If p represents the growth of phytoplankton for example, a may be related to nutrient availability and b to the position of phytoplankton in the sunlit water column. In this case the deterministic method expresses the average state (or mean field approximation) of the actual stochastic process, this statement is true for a large number of individuals with small enough fluctuations in dynamics, but it is not true in general. Equation 3.4 represents the Malthusian exponential population growth (as long as birth > death), an example of a mechanism that can stop the Malthusian explosion is the finiteness of resources, in the biological context it is common to write it as a logistic model [START_REF] Tsoularis | Analysis of logistic growth models[END_REF]:

dN dt = rN 1 - N K ≡ g(N)N. (3.6)
Here r, the difference between birth and death rates is the maximum growth rate, since the effective population-dependent growth rate g(N) = r(1 -N/K) decreases with N. K is the carrying capacity, identified with the maximum population that can be sustained by the available resources, that is the population value at which growth stops. Fig. 3.4 illustrates the form of the population growth described by the integration of equation 3.6. Also one could cite population growth dynamics known as Allee effects [START_REF] Neufeld | Chemical and biological processes in fluid flows: a dynamical systems approach[END_REF], which reduce individual fitness when population density is small, this effect is introduced to deal with the codependence of individuals in certain situations such as winter, when animals need to pack close together to avoid death by cold. One could also enter the next step of modelling by adding predator-prey dynamics through the coupling of two equations, one representing the prey and another the predators, as introduced by Volterra in 1926 [START_REF] Volterra | Fluctuations in the abundance of a species considered mathematically[END_REF]. But in this work we chose to follow a single population approach for convenience since our main interest is in the interplay between phytoplankton dynamics and advection diffusion features.

The Advection-Reaction-Diffusion Framework

Before entering the modelling methodology specific to phytoplankton it is interesting to present the mathematical framework that allows to couple the aquatic transport seen in Chapter 2 with the biological conditions seen in this chapter. To represent the survival dynamics of these organisms it is common to lay the different parts standing for the transport and biology of phytoplankton organisms into an advection-reaction-diffusion equation. The standard advection-reaction-diffusion model deals with the time evolution of a chemical or biological species in a flowing medium such as water or air. The mathematical equations describing the evolution of the model are partial differential equations (PDEs) that can be derived from the mass action law [START_REF] Murray | Mathematical biology: I. An introduction[END_REF] (for the biological part). One can separate the processes that compose the general equation as to better explain their function. A classical representation of a advection equation is given by the continuity equation:

∂ ∂t θ = -∇ • (vθ) (3.7)
θ can represent either a fluid density or an inert chemical/biological substance concentration, v represents the velocity field affecting both the fluid and the inert substance. θ can be easily generalized as any scalar field, if it was temperature for example, Eq. 3.7 would be part of the dynamics involved, e.g, in thermal convection.

The equation of diffusion was already presented in section 2.1.2. We can add here that if D represents turbulent diffusivity, as seen in Chapter 2, it can be chosen to have different constant values depending on different regions for example, something that is used to represent different layers of mixing in the ocean [START_REF] Huisman | Population dynamics of sinking phytoplankton in light-limited environments: simulation techniques and critical parameters[END_REF]. We will talk about grid discretisation later on, but it is worth remarking that it can also be used to represent sub-grid turbulence [START_REF] Lindemann | Dynamics of phytoplankton blooms in turbulent vortex cells[END_REF]. The final piece of the model is the reaction equation, in here for simplicity we express it with a function that is responsible for expressing local changes due to biological or chemical reactions:

∂ ∂t θ = g (θ) (3.8)
where as said, g can represent any type of biological or chemical reaction kinetics, from light-limited phytoplankton growth dynamics [START_REF] Platt | Mathematical formulation of the relationship between photosynthesis and light for phytoplankton[END_REF], to combustion reactions [START_REF] Abel | Front propagation in laminar flows[END_REF]. Meaning that g controls the value of θ accordingly to its definition and in conjunction to the advection-diffusion processes. In our case, the complete advectionreaction-diffusion equation is used to describe a passive-reactive scalar, phytoplankton density, in two dimensions:

∂θ(x, z, t) ∂t = ∇ • (D∇θ(x, z, t)) -∇.(v(x, z, t)θ(x, z, t)) + g(θ)θ(x, z, t). (3.9)
Hence θ(x, t) represents the density of phytoplankton cells at point x = (x, z) and time t. The novelty comes in the form of the reaction term g(θ)θ(x, z, t), this part of the equation represents a biological reaction (parallel to the growth term of, e.g., Eq. 3.6). Eq. 3.9 no longer obeys the mass conservation principle since now cells are produced or lost due to the form of the reaction term, this means that the equation relies solely on its initial condition and on its boundary conditions to be solved for each time t and coordinates [x, z]. Consequently an advection-reaction-diffusion equation of this type can be fairly more complicated to solve than the advection-diffusion equations previously discussed, even for an incompressible velocity flow field. In reality, they are analytically solvable for phytoplankton dynamics only when a sensible amount of approximations or assumptions are made [START_REF] Ebert | Critical conditions for phytoplankton blooms[END_REF][START_REF] Huisman | How do sinking phytoplankton species manage to persist?[END_REF]. Nonetheless, they have been widely utilized in plankton mathematical modelling studies for more than 50 years [START_REF] Shigesada | Analysis of the self-shading effect on algal vertical distribution in natural waters[END_REF][START_REF] Okubo | Oceanic diffusion diagrams[END_REF][START_REF] Hernández | Sustained plankton blooms under open chaotic flows[END_REF][START_REF] Taylor | Shutdown of turbulent convection as a new criterion for the onset of spring phytoplankton blooms[END_REF], but instead of being approached analytically they are mostly dealt with numerically. This can be done in many ways, the methods utilised in this thesis will be further discussed in Chapter 4. What remains for us in this chapter is to understand the form of the reaction growth term g, this can depend on the spatial position [START_REF] Lindemann | Dynamics of phytoplankton blooms in turbulent vortex cells[END_REF], on one or more population densities [START_REF] Sandulescu | Biological activity in the wake of an island close to a coastal upwelling[END_REF] and subsequently on time.

In the next subsection we will show how g can be used to portray the biological dependence phytoplankton has on light.

Light limited growth

The history of phytoplankton light limited growth models is quite rich, starting from the seminal work of Sverdrup [START_REF] Sverdrup | On conditions for the vernal blooming of phytoplankton[END_REF], it gained much impetus through the introduction of specific studies on the importance of light limitation for phytoplankton dynamics [START_REF] Platt | Mathematical formulation of the relationship between photosynthesis and light for phytoplankton[END_REF][START_REF] Shigesada | Analysis of the self-shading effect on algal vertical distribution in natural waters[END_REF], reaching the status of a well established methodology in the modern days [START_REF] Huisman | Critical depth and critical turbulence: two different mechanisms for the development of phytoplankton blooms[END_REF][START_REF] Ebert | Critical conditions for phytoplankton blooms[END_REF][START_REF] Huisman | Population dynamics of sinking phytoplankton in light-limited environments: simulation techniques and critical parameters[END_REF][START_REF] Taylor | Shutdown of turbulent convection as a new criterion for the onset of spring phytoplankton blooms[END_REF]. According to Huisman [START_REF] Huisman | Light-limited growth and competition for light in well-mixed aquatic environments: an elementary model[END_REF] a variety of competition models lead to the following prediction: if one resource is the main limiting factor for the survival of multiple species, the species that is able to survive at the lowest resource density equilibrium value (R * ) competitively excludes all other species. In here he Chapter 3. Population Dynamics, Plankton and Phytoplankton also quotes Tilman [START_REF] Tilman | Mechanisms of plant competition for nutrients: the elements of a predictive theory of competition[END_REF], "light competition is conceptually more complex than nutrient competition. We do not yet have either rigorous theoretical predictions or experimental results that indicate a single number, analogous to R * , that can predict the outcome of competition for light." This naturally arises from the fact that in relation to nutrients, light as a resource present further modelling complications, it can never be assumed to be distributed homogeneously vertically due to absorption/shading effects, and the response of micro organisms to such resource are still not fully understood. With that said, one can start the model construction by considering the Beer-Lambert law for the attenuation of light due to the water medium with depth:

dI dz = -κI, (3.10)
where κ is the extinction coefficient, I is the light intensity, and z the depth. The extinction coefficient can consist of two parts:

κ = K bg + kθ, (3.11)
K bg is the background extinction coefficient associated to the water turbidity, without phytoplankton, and k is the specific light attenuation coefficient of the phytoplankton pigments, representing the effect known as self-shading. The introduction of the self-shading effect is not mandatory for vertical phytoplankton numerical studies [START_REF] Lindemann | Dynamics of phytoplankton blooms in turbulent vortex cells[END_REF], specially when phytoplankton population density is low. In spite of that, it is known from numerical models [START_REF] Shigesada | Analysis of the self-shading effect on algal vertical distribution in natural waters[END_REF][START_REF] Platt | Mathematical formulation of the relationship between photosynthesis and light for phytoplankton[END_REF] that the self-shading effect of phytoplankton plays an important role, especially in the dynamics of spring algal blooms in coastal waters. Without self-shading, the plankton concentration would increase rapidly, making the nutrient limitation apparently too strong. The final form of the light intensity modulation throughout the water column is then given by:

I(z, t) = I in exp - z 0 [K bg + kθ(z , t)]dz , (3.12)
with I in as the incident light intensity, and θ(z, t) the phytoplankton concentration.

From the form of Eq. 4.22 one can see why the the self-shading effect is not always taken into account by numerical models. Its inclusion in the complete population dynamics with advection and diffusion leads to an integro-differential equation, making the solution non-trivial for numerical methods and computational resources. The final step is to introduce a growth term analogous to r in Eq. 3.3. This has to consider the light limitation and death dynamics so it can be later used as a biological reaction term in a complete population dynamics simulation. Therefore we can consider a growth rate:

g(I(z, t)) = aI(z, t) -l. (3.13)
with a reproduction rate aI proportional to the local light intensity and constant a, and a loss term l representing population losses due to respiration, sedimentation, excretion, between other processes. Eq. 3.13 is a good approximation at low light intensities, however, due to the physiological limits of the maximal reproduction rates of organisms, the reproduction rate generally saturates for high light intensities [START_REF] Ebert | Critical conditions for phytoplankton blooms[END_REF]. Such a behaviour can be modelled by Michaelis-Menten kinetics, originally introduced to represent the rate of enzymatic reactions [START_REF] English | Ever-fluctuating single enzyme molecules: Michaelis-Menten equation revisited[END_REF], in the phytoplankton modelling field it has the form [START_REF] Monod | Technique, theory and applications of continuous culture[END_REF]:

g(I(z, t)) = aI 1 + cI -l, (3.14)
where c is a reproduction saturation constant. For cI 1, the function reduces to the approximation of equation 3.13. The mathematical expression from Eq. 3.14 with parameters taken from [START_REF] Huisman | How do sinking phytoplankton species manage to persist?[END_REF] can be seen in Fig. 3.5. The presented curve is independent of phytoplankton population and therefore does not take into account the self-shading effect (explained in the next section). Even so we can still evaluate that the depth at which production is equal to loss is quite shallow (17.8m), this means that for considerably clear water, if the mixing layer has a size of 100 meters (quite usual for deep lakes and mid latitude oceans [START_REF] Williams | Ocean dynamics and the carbon cycle: Principles and mechanisms[END_REF] in winter), only about 20% of the water column on which algae are circulating will lead to production large enough to surpass death. This gets even more complicated with the introduction of the selfshading effect and sinking dynamics, showing how complex the modelization of the life of phytoplankton can be.

Phytoplankton Bloom Dynamics

Scientific observations relating the spatial distribution of algae patches with ocean currents date from the beginning of the last century. In 1927 for instance, in a cruise between the United States and England, Langmuir reported that large quantities of seaweed were arranged parallel to the wind direction [START_REF] Langmuir | Surface motion of water induced by wind[END_REF] in a phenomena that would later be baptized with his name, the Langmuir circulation. From then to the current days a lot has been done in the efforts of explaining the behavior of plankton blooms. Here we describe the critical depth hypothesis, a theory based on a article [START_REF] Sverdrup | On conditions for the vernal blooming of phytoplankton[END_REF] that achieved much impetus in oceanography and aquatic ecology [START_REF] Lucas | Does the Sverdrup critical depth model explain bloom dynamics in estuaries?[END_REF][START_REF] Huisman | Critical depth and critical turbulence: two different mechanisms for the development of phytoplankton blooms[END_REF]. The concept presented there was later thoroughly debated and so refinements Chapter 3. Population Dynamics, Plankton and Phytoplankton and alternatives to the theory appeared, the critical turbulence hypothesis is one of these spawns that we will also discuss in this section.

Critical Depth Hypothesis

In 1935 Gran and Braarud [START_REF] H H Gran | A quantitative study of the phytoplankton in the Bay of Fundy and the Gulf of Maine (including observations on hydrography, chemistry and turbidity)[END_REF] pointed out that phytoplankton production in a water column could not exceed losses if there existed a deep mixed top layer. In this one can think of a compensation depth, defined as the depth at which the energy intensity is such that the production by photosynthesis balances losses by respiration. Therefore the reasoning was that within a well mixed layer the plankton organisms are about evenly distributed, but a net production takes place only above the compensation depth (Fig. 3.6), whereas below the compensation depth there is a net loss of organic matter since the light depletion makes production by photosynthesis smaller than the loss by respiration. The total population cannot increase if this loss exceeds the net production. This implies that there must exist a critical depth such that blooming can occur only if the depth of the mixed layer (the water column slice where a plankton bloom may occur) is less than the critical value. Gran and Braarud concluded that the critical depth was 5 to 10 times the compensation depth. From this idea Sverdrup writes his classical article in 1953 [START_REF] Sverdrup | On conditions for the vernal blooming of phytoplankton[END_REF] creating a model that will appear in many oceanography text books and will be considered a milestone for many other models in the field. To justify the conclusions Sverdrup arrives in his model he makes the following assumptions:

1. Within the top layer the turbulence is strong enough to distribute the plankton organisms evenly through the layer.

2. Production is not limited by lack of plant-nutrition salts.

3. κ, the extinction coefficient of radiation energy is constant in the water column.

4. The production of organic matter by photosynthesis is proportional to the energy of the radiation at the level under consideration.

5. I c , the energy at the compensation depth, is known.

Many of these assumptions were also made before by Riley in 1942 [START_REF] Riley | The relationship of vertical turbulence and spring diatom flowerings[END_REF] in a model we will discuss later. In Sverdrup's case the result was a model for the total production and destruction of phytoplankton in a water column based on Beer-Lambert's law. Starting from the expression for the total production:

P = m T 0 0 -D z I in e κz dtdz = m k (1 -e -κD z ) T 0 I in .dt (3.15)
Where m is the production rate constant, I in is the effective energy at the surface, κ light depletion constant, T the interval of interest, D z the mixed layer depth and z the vertical coordinate. Following the expression for the total losses by respiration:

R = n T 0 0 -D z dtdz = nTD z , (3.16)
where n is the death rate constant. The condition for an increase in phytoplankton population is P > R, this happens when the area bcD c in Fig. 3.6 is greater than the area dabD c . Introducing the average incident light:

Īin = 1 T T 0 I in dt. (3.17)
Also I c = n/m, at the compensation depth, the production and loss rates are equal, m = n, the critical depth results by equating Eq. 3.15 and Eq. 3.16. From these Sverdrup gave the following expression for the critical depth:

D cr 1 -e -kD cr = Īe κI c . (3.18)
Sverdrup remarked that the speed in which the plankton grows is proportional to how much the thickness of the top layer deviates from the critical value. The more it deviates, the faster it grows. Also he stated that the presence of weaker turbulence may affect the critical value as a higher concentration of plankton near the surface may compensate the destruction under the compensation depth (as we will see later). Data from the months of March, April and May in Norway in 1949 further supported his reasoning. In months where the mixed layer was deeper than the critical depth no increase in population was detected, in months where the mixed layer was shallower than the critical depth a bloom occurred with intensity proportional to the deviation from the critical value.

Critical Turbulence Hypothesis

Here we discuss the work used as basis for our phytoplankton model [START_REF] Huisman | How do sinking phytoplankton species manage to persist?[END_REF]. Even before Sverdrup presented his study scientists in the field already discussed the importance of phytoplankton vertical stability for production. Riley [START_REF] Riley | The relationship of vertical turbulence and spring diatom flowerings[END_REF] was the first to propose a relation for the loss rate of diatoms considering the water column depth:

r = r p D r p D r c -r c -r s (3.19)
Chapter 3. Population Dynamics, Plankton and Phytoplankton

With r as the phytoplankton loss rate, D r p the thickness of the euphotic zone, D r c the thickness of the vertical turbulence layer, r c the loss by respiration and r s the loss by grazing and sinking.

Riley coined this based on data collected on four cruises made to Georges Bank in the United States during the spring of 1941. By analyzing the number of plankton cells against the vertical stability (density difference between different depths) he realized that during early spring when the bloom was peaking there was no relation between the parameters, but in April when the peak was gone there was a positive relation between them. He finally concluded that if the euphotic zone and the rates of photosynthesis and respiration are constant, the rate of increase in the plankton population is a linear function of the reciprocal of the thickness of the zone of vertical turbulence (mixed layer). Later in 1949 Riley et al. [START_REF] Riley | Quantitative ecology of the plankton of the western North Atlantic[END_REF] delivered a one-dimensional (1D) model that represented the vertical processes of diffusion and settling combined with the simple growth-death processes of phytoplankton represented by a two-layer model:

∂θ ∂t = p(I in )θ + D ∂ 2 θ ∂z 2 -v ∂θ ∂z (3.20)
for phytoplankton in the euphotic zone and,

∂θ ∂t = -lθ + D ∂ 2 θ ∂z 2 -v ∂θ ∂z (3.21)
for phytoplankton below the euphotic zone.

Here θ represents the plankton concentration, p(I in ) the growth rate as a function of the incident radiation, l the death rate, D the turbulent diffusivity constant, and v the sinking velocity. Through this Riley was able to derive a simple expression for the minimal turbulence [START_REF] Riley | Quantitative ecology of the plankton of the western North Atlantic[END_REF][START_REF] Ryabov | Population growth and persistence in a heterogeneous environment: the role of diffusion and advection[END_REF][START_REF] A B Ryabov | Vertical distribution and composition of phytoplankton under the influence of an upper mixed layer[END_REF], an important concept that represents the turbulence necessary to maintain a sinking species of phytoplankton enough time in the euphotic zone as to have their growth superior to death:

D min ≈ v 2 4(p(I in ) -l) . ( 3.22)
The modelling evolution continued with Shigesada and Okubo [START_REF] Shigesada | Analysis of the self-shading effect on algal vertical distribution in natural waters[END_REF] who in 1981 built a model which incorporated light-dependent growth rates but neglected light absorption by the water column. This would be refined a posteriori by [START_REF] Huisman | Critical depth and critical turbulence: two different mechanisms for the development of phytoplankton blooms[END_REF] through the development of a model [START_REF] Huisman | Critical depth and critical turbulence: two different mechanisms for the development of phytoplankton blooms[END_REF] that provided a theoretical underpinning for observations of Townsend et al. [START_REF] Townsend | Spring phytoplankton blooms in the absence of vertical water column stratification[END_REF] and Eilertsen [START_REF] H C Ellertsen | Spring blooms and stratification[END_REF], which challenged the critical depth hypothesis. The model was a growth-diffusion model in which the population dynamics of phytoplankton were governed by light-limited growth, local phytoplankton losses, and local transport of the phytoplankton by turbulent diffusion. Following Sverdrup [START_REF] Sverdrup | On conditions for the vernal blooming of phytoplankton[END_REF], Huisman et al. kept the model as simple as possible in order to have a representation that was close to Sverdrup, but with the consideration of turbulence. Thus, the model neglected many additional complexities like nutrient limitation, photo-inhibition (decreased efficiency in photosynthesis), buoyancy, and sinking: Where θ(z, t) is the population density (cells per unit volume) at depth z and time t, p(I(z, t)) is the specific production rate as an increasing function of light intensity I(z, t), l is the specific loss rate, p max is the maximal specific production rate and D is the turbulent diffusion coefficient (also known as "vertical eddy diffusivity"). The light intensity is modeled by the incident light intensity I in , K bg is the total background turbidity (responsible for modelling the water depletion of light), H is the half-saturation constant of light-limited growth and k is the specific light attenuation, responsible for accounting for the self shading of the plankton. The model also had no flux boundary conditions:

∂θ(z, t) ∂t = [p(I(z, t)) -l]θ(z, t) + D ∂ 2 θ(z,
D ∂θ(z, t) ∂z = 0, at z = 0 and z = z max . (3.25)
Where z max is the domain depth. From simulations utilizing Eq. 3.7 Huisman builds a diagram (Fig. 3.7) for the existence or non existence of plankton blooms depending only on the depth of the simulated water column and the diffusivity coefficient.

Here the bloom condition is defined as the state where the phytoplankton population grows from an initial state of small uniform density and then reaches a stationary state. The no bloom condition is defined as the condition in which the phytoplankton population disappears from an initial state equal to the aforementioned one, after enough simulation time has elapsed.

Chapter 3. Population Dynamics, Plankton and Phytoplankton

The line perpendicular to the Turbulent diffusion axis that delimits the bloom and no bloom regions represents the critical turbulence (or maximal turbulence), a concept that complements the idea of minimal turbulence developed by Riley et. al [START_REF] Riley | Quantitative ecology of the plankton of the western North Atlantic[END_REF]. This is the key concept of his work, the critical turbulence delimits the maximal turbulence a water environment of defined depth can have to support phytoplankton blooms.

The turbulence importance in phytoplankton life-cycles, an idea that has been discussed through theories [START_REF] Ebert | Critical conditions for phytoplankton blooms[END_REF][START_REF] J R Koseff | Coupled effects of vertical mixing and benthic grazing on phytoplankton populations in shallow, turbid estuaries[END_REF], observations [START_REF] J E Cloern | Tidal stirring and phytoplankton bloom dynamics in an estuary[END_REF][START_REF] Lauria | Contrasting phytoplankton distributions controlled by tidal turbulence in an estuary[END_REF] and experimentation [START_REF] Visser | Artificial mixing prevents nuisance blooms of the cyanobacterium Microcystis in Lake Nieuwe Meer, the Netherlands[END_REF][START_REF] Reynolds | Some effects of artificial mixing on the dynamics of phytoplankton populations in large limnetic enclosures[END_REF] gained much momentum through the advancement of phytoplankton dynamics modelling. The next logical step was for Huisman to include a sinking/buoyancy velocity for the plankton in the model and that he did in 2002 [START_REF] Huisman | How do sinking phytoplankton species manage to persist?[END_REF]. By adding a velocity therm to Eq. 3.23 we have the following partial differential equation (PDE): In Fig. 3.8 wee see that the inclusion of sinking speed in the model creates two thresholds for both turbulence and water column depth. For turbulence we have the minimal turbulence, the minimal amount of turbulence required to compensate sinking in order to remain in the euphotic zone enough time to outgrow death, and the maximum turbulence, the maximum amount of turbulence a water environment can provide above which the mixing is too strong for the bloom to exist. The critical depth is the same as described by Sverdrup's theory and the compensation depth is the same mentioned before (the depth below which death by respiration surpasses growth by photosynthesis). The second graph in Fig. 3.8 shows what happens when the sinking speed is considerably increased. In this scenario the maximal and minimal turbulence limits merge and the system is only weakly affected by the amount of turbulence, being dominated by the depth of the water column. These results are important because they suggest clear boundaries for important physical quantities that affect the phytoplankton life cycle and will be a solid starting point for the study we carry on in this thesis.

∂θ(z, t) ∂t = [p(I(z, t)) -l]θ(z, t) + D ∂ 2 θ(z, t) ∂z 2 -v ∂θ(z,

Current Bloom Scientific Debate

We saw that most past theories relied on a handful of approximations but nowadays modelling in this field has achieved a much more complex state. As we will later discuss in the next chapter it has incorporated further dynamics that are very difficult to include explicitly in these theories, from physical quantities as water temperature and buoyancy [START_REF] Dube | Modelling the effects of variable salinity on the temporal distribution of plankton in shallow coastal lagoons[END_REF][START_REF] Wroblewski | An Ocean Basin Scale Model of plankton dynamics in the North Atlantic: 1. Solutions For the climatological oceanographic conditions in May[END_REF], to biological factors as predators and nutrients [START_REF] Franks | NPZ models of plankton dynamics: their construction, coupling to physics, and application[END_REF]. When we include into these the considerable amount of field data, experiments and satellite observations that were and are still being performed in the last decades, the room for discussion has increased. On the topic of plankton bloom dynamics the main debate still is around the spring bloom observed every year in the north Atlantic [START_REF] Siegel | The North Atlantic spring phytoplankton bloom and Sverdrup's critical depth hypothesis[END_REF]. A large part of what inspired the phytoplankton bloom discussion comes from observations, mostly in the north Atlantic, where a strong plankton bloom is seen every spring. And so a debate has formed around the phenomenon that triggers this event. Historically, it is believed that the phenomenon is a consequence of the shoaling of the mixed layer during spring, in agreement with the critical depth hypothesis. From the recent theories that have followed the critical depth theory in hopes to better explain the spring bloom phenomenon we can highlight the hypothesis raised by an article by Behrenfeld [START_REF] Behrenfeld | Abandoning Sverdrup's critical depth hypothesis on phytoplankton blooms[END_REF]. In this it is suggested that the critical depth hypothesis has to be abandoned based on its incompatibility with a nine-year satellite record of phytoplankton biomass in the subarctic Atlantic that shows that bloom initiation occurs in the winter when mixed layer depths are maximum, not in spring. Further it is argued that the coupling between phytoplankton growth and losses increases during spring stratification rather than decreases, and maxima in net population Chapter 3. Population Dynamics, Plankton and Phytoplankton growth rates are as likely to occur in midwinter as in spring. As an alternative a "Dilution-Recoupling Hypothesis" is reported where rather than being triggered by mixed layer shoaling, the spring bloom is triggered by the dilution of phytoplankton grazers (predators) when the inverse occurs, the mixed layer increases Fig. 3.9.

FIGURE 3.9:

Illustration of the concept behind the "Dilution-Recoupling Hypothesis". Dilution impact of the winter increase in mixed layer depth is higher on grazers than on phytoplankton (in loss rates), resulting in the bloom starting at winter. Extracted from [START_REF] Behrenfeld | Resurrecting the ecological underpinnings of ocean plankton blooms[END_REF] .

In another important work Taylor and Ferrari [START_REF] Taylor | Shutdown of turbulent convection as a new criterion for the onset of spring phytoplankton blooms[END_REF] postulate that the onset of a bloom can also be triggered by a reduction in air-sea fluxes at the end of winter. The hypothesis is supported by numerical simulations and satellite remote sensing data. They suggest that when net cooling subsides at the end of winter, turbulent mixing becomes weak, increasing the time phytoplankton spends in the euphotic zone. The hypothesis is compatible with the critical depth hypothesis but more complete as it does not rely simply on the mixing layer depth shoaling but on reduction of atmospheric forcing as a trigger. As we can see, the field is far from reaching a consensus on what phenomenon causes phytoplankton blooms or extinction. This constitutes an important motivation for the work presented in this thesis. We hope that we can provide some interesting insight on the basic features related to the interplay between biological and fluid dynamics, at the hearth of such a complex problem, in the next chapters.

Chapter 4

Numerical Methods

In this chapter I describe the numerical methods adopted in this thesis and report technical points related to the space and time convergence of the algorithms applied to the investigated problem (see Sec. 4.4). I finish with the presentation of the model considered in the numerical work and provide comparisons with literature results for validation purposes. In the first section I review advection-reaction-diffusion equations from a numerical point of view and describe the discretization and time marching schemes. I then describe the specific semi-Lagrangian method chosen for advection and finally present the method adopted for diffusion solution scheme.

Advection-Reaction-Diffusion Equations

Here we discuss the advection-reaction-diffusion equation introduced in Chapter 3 from a numerical point of view. The complete equation (see also Eq. 3.9) has the form:

∂θ ∂t = D∇ 2 θ -∇ • (vθ) + g(θ)θ. (4.1)
The equation in the above form is quite daunting as one should expect from any equation representing the time evolution of a scalar field under these different processes. In addition to that, one can rarely find an analytical solution for equations of this type without many assumptions. Indeed it is known that already the passive scalar problem is highly nontrivial in the presence of advection, even for simple flows [START_REF] Falkovich | Particles and fields in fluid turbulence[END_REF]. Here, the presence of a biological reaction adds to the complexity through its typically nonlinear character, giving rise to very rich dynamics [START_REF] Neufeld | Chemical and biological processes in fluid flows: a dynamical systems approach[END_REF]. Nonetheless one can count on numerical methods to deal with this problem.

One typically starts by discretizing space into cells of size h > 0 (a small number), this way our once continuous density distribution is represented into cell intervals. The procedure is described in Fig. 4.1, starting from a continuous 2D scalar density field θ(x, z, t), one has to make the equivalence between the continuous coordinates (x, z) and the discrete grid positions given by the indices i and j respectively. The cell size (in our case h for both directions) and the resolution in terms of the number of grid points, N x and N z , will be chosen based on the value of the largest and smallest scales of interest. Fluid mechanics numerical methods that rely on grid discretizations such as DNS (Direct Numerical Simulations), LES (Large Eddy Simulations), semi-Lagrangian, Lattice-Boltzmann, Finite-Differences, between others, have strengths and weaknesses related to both their accuracy and their computational costs. DNS, for example, is usually regarded as a method that provides very realistic solutions since it is a direct solution of the Navier-Stokes equations, but Chapter 4. Numerical Methods is also computationally costly, often depending on Fourier transforms and no turbulent viscosity approximation, differently from LES. In this thesis the choice for a semi-Lagrangian method to account for the effect of the velocity flow field was partially justified by the necessity to simulate large physical domains extending over several meters in order to reproduce realistic setups, something that would have been impossible in simulation time with DNS. Finally, in addition to the spatial discretization, one has to discretize the time evolution in time steps ∆t. The advection-reaction-diffusion equation is integrated in time at each grid point using standard discretisation algorithms (Euler, Velocity-Verlet, Runge-Kutta). The method chosen for this thesis was the 4th order Runge-Kutta, and it was used both for the integration of the reaction term and advection terms, diffusion was integrated via Crank-Nicholson finite differences method (see Sec. 4.3). In the next sections we discuss in depth the methodology for the step-bystep solution of the advection-reaction-diffusion equation of this project.

Semi-Lagrangian Scheme

During my thesis I worked with a semi-Lagrangian scheme through the use of a discrete-time map approach for the solution of the advection part of the mentioned model equation. This involves both the flow field and the sinking velocity of phytoplankton. In the first year of the thesis I relied on a so called forward scheme, that will be explained here, for the solution of both the diffusion and the advection parts, in which it was easier to implement the no-flux boundary conditions at the top and bottom sides of the domain. During the second year of PhD, due to the limited accuracy of the resulting advection-diffusion solution (related to extrapolation on the grid points and number of noise realizations needed for diffusion), I changed the algorithm for a backward advection scheme with diffusion by finite differences. The backward advection method is also explained in this section but the diffusion solution is presented in the next section. With the new methodology we were able to ensure field smoothness, which is clearly an essential prerequisite. The semi-Lagrangian scheme is very practical to work with for its simplicity and its relative flexibility to adopt complex boundary conditions. Also, as advection is integrated in a continuous space, an advantage of this approach is that it does not require sophisticated implementations and limiting constraints to deal with transport by fluid flows.

Discrete-Time Map

Discrete-time maps are widely used in the numerical solution of advection-reactiondiffusion equations [START_REF] Abel | Front speed enhancement in cellular flows[END_REF][START_REF] Berti | Mixing and reaction efficiency in closed domains[END_REF], here we can follow a similar procedure to the one described by Abel et al. [START_REF] Abel | Front propagation in laminar flows[END_REF] to explain the methodology. The numerical solution of equation 4.1 for scalar parcels θ(x, z, t) can be separated, we can first calculate the trajectory followed by the parcel due to diffusive and advective transport in our simulation box, and then integrate the reaction part according to the reaction term g(θ)θ, fulfilling the complete scalar time evolution. Hence we establish F t as the formal evolution operator for a 2D trajectory described by the position vector x(t) = (x(t), z(t)), as x(t) = F t x(0) and G t as the time evolution operator according to θ(t) = G t θ(0). The total dynamics are represented by:

θ(x, z, t) = G t θ[F t (x(0), z(0)), 0]. (4.2)
The next step is to introduce the time discretization to the operators in the form of a Lagrangian map x(t + ∆t) = F ∆t (x(t)). This follows a discrete time march of intervals ∆t. Then Eq. 4.2 can be generalized in this discrete time form in the two ways that were relevant for my work. The forward advection with diffusion in a Lagrangian way:

θ(x(t + ∆t), t + ∆t) = G ∆t [θ(F ∆t (x(t)) + √ 2D∆tζ(t), t)] ζ (4.3)
where ζ(t) is a random vector of zero mean and unity variance (as seen in Chapter 2), and the process has to pass trough an ensemble average <> ζ to represent the diffusion process. In practice this means multiple implementations of ζ have to be accounted for and then averaged. The other method I applied is represented by the backwards advection scheme without diffusion (that is implemented separately via finite differences):

θ(x(t + ∆t), t + ∆t) = G ∆t [θ(F -1 ∆t (x(t + ∆t), t)]. (4.4)
The main differences between the schemes represented by equations 4.3 and 4.4 are illustrated in Fig. 4.2: In the Forward case we integrate a scalar parcel trajectory (via the operator F ∆t ) starting from a grid point and ending at a different position (typically not a grid point). Consequently the value of the parcel has to be spread (or extrapolated) surrounding grid points, the procedure is repeated until all scalar parcels, each starting from a grid point, have been evolved following their trajectories. We then have the advected-diffused population field resulting from the sum of all scalar parcels that were transported and then proportionally divided according to grid discretization. The final step is to implement reaction (via the operator G ∆t ) over the entire population field in accordance to the reaction term g(θ)θ. While in this case diffusion trajectories are implemented by the ensemble average of multiple calculations of the random kick √ 2D∆tζ(t) using first order discretisation, both the advection trajectory and the reactive evolution of the scalar parcels are calculated via the 4th order Runge Kutta method. On the other hand, the Backwards scheme relies on the calculation of trajectories backwards in time. Starting from the future grid point we calculate the past trajectory that arrived at that point (via the operator F -1 ∆t ). The origin of such trajectory is off grid too, therefore interpolation using the nearest grid points is needed to compute the scalar density at the Lagrangian origin. The so-determined θ values is then reacted (via the operator G ∆t ) and the procedure is repeated for the entire grid. After we have advected the θ field and let it react, we calculate diffusion via finite differences (explained in the next section). Also here, both backwards advection and reaction are implemented via the 4th order Runge Kutta method.

Interpolation/Extrapolation

The interpolation used for the backwards advection method was a bilinear one, meaning employing the nearest neighbours on the grid. The Lagrangian origin F -1 ∆t (x(t + ∆t)) is used to obtain the scalar density θ(F -1 ∆t (x(t + ∆t), t) that will be advanced in time according to the given reaction dynamic. It is not difficult to see that the coefficients needed for the interpolation are:

c 1 = (x old -x grid )(z old -z grid ) h 2 , ( 4.5)
c 2 = x old -x grid h -c 1 , (4.6
)

c 3 = z old -z grid h -c 1 , (4.7
)

c 4 = 1 -(c 1 + c 2 + c 3 ), (4.8)
where (x old , z old ) is the Lagrangian origin, (x grid , z grid ) are the coordinates of one of the lower left vertex of the square formed by the 4 points on the grid nearest to (x old , z old ) (as depicted by Fig. 4.2), and h is the cell size. The interpolated, prereaction, density field is:

θ pre = c 4 θ old (x grid , z grid ) + c 2 θ old (x grid + h, z grid ) +c 3 θ old (x grid , z grid + h) + c 1 θ old (x grid + h, z grid + h), (4.9)
where θ old is the gridded density field at the previous time step. Boundary conditions can be dealt with by reducing the order of interpolation, or they can be solved separately via finite differences as well.

The extrapolation process needed for the forward advection scheme follows the same logic of coefficient calculation as equations 4.5, 4.6, 4.7 and 4.8, but instead of relying on a scalar parcel interpolated from 4 grid contributions the scalar parcel is divided between the grid points according to: θ pre (x grid , z grid ) = c 4 θ old , (4.10)

θ pre (x grid + h x , z grid ) = c 2 θ old , (4.11)
θ pre (x grid , z grid + h z ) = c 3 θ old , (4.12)

θ pre (x grid + h x , z grid + h z ) = c 1 θ old , (4.13)
where the now gridded θ pre can receive contributions from multiple particle based values θ old , transported starting from different grid points. This means that the entire field has to be advected before the reactive dynamics can be evolved. something that is not necessary in the backwards case. Boundary conditions are imposed in the same manner as in the backwards case.

When implementing the backwards advection-diffusion technique I was faced with a problem. Since transport had to be reverted at the no-flux boundaries, reverse trajectories that were reflected resulted in the θ old parcel being taken from a slightly higher density region at the top boundary due to the phytoplankton sinking. One possible solution for this problem would be to rely on very small time steps for the simulation dynamics, but due to the constraints related to the implementation of the stochastic term for diffusion, the grid would have to be very fine as well, making the solution impracticable. The first solution for this problem found in this work was the discussed forward advection-diffusion scheme, this solved the issue of the Lagrangian origins of grid points close to the top boundary. However a better solution was implementing diffusion via finite differences in the backwards method, as the stencil scheme chosen for diffusion with the backwards methodology for advection also resulted in smoother density fields. The latter was then chosen as the reference implementation adopted for the simulation code.

An advantage of the semi-Lagrangian implementation of diffusion is that it does not impose an upper bound on the time step ∆t for stability reasons. To simulate diffusion on scale h is then enough to respect the condition h ≤ √ 2D∆t. Nonetheless, a strong limitation of the method with forwards or backwards diffusion (not through the Eulerian scheme) is the need to perform ensemble averages, which require tuning for the parameters, namely the number of noise realizations. This can make computations costly (the needed number of realizations in our case was ∼ 50).

Eulerian Scheme for Diffusion

As we saw in the previous sections the advection-reaction-diffusion equation can be solved numerically by a splitting approach, in Sec. 4.2 we have shown the solution of the advection and diffusion parts with a semi-Lagrangian scheme with subsequent reaction integration. Here we show an alternative that is using the backwards-Lagrangian scheme only for the advection (without changing the reaction integration method), but diffusion dynamics are calculated through a pseudo-Crank-Nicholson method with spatial finite differences. Starting from the Crank-Nicholson method, let us introduce the function f representing the right side of the diffusion equation as:

∂θ ∂t = D∇ 2 θ = f (θ, D) (4.14)
Then, on our 2D grid indexes (i, j), numerical integration is calculated similarly to the Euler method as:

θ t+∆t i,j -θ t i,j ∆t = 1 2 f t+∆t i,j + f t i,j . (4.15)
The Crank-Nicholson method is a semi-implicit scheme based on the trapezoidal rule and has second order convergence in time. One could go further and use the RK4 method for better precision, but the coupling of the method with a finitedifference scheme for the spatial discretization increase the number of grid elements to be computed as O(2d 2 ) where d is the number of dimensions. This means that one has to weight the necessity of better precision, with two or more dimensions calculations can become quite heavy since the RK4 method requires even more operations. In our case we found that the a pseudo-Crank-Nicholson method with a second order central finite differences scheme presented satisfactory convergence. So for the spatial discretization in 2D:

f i,j = D θ i+1,j + θ i,j+1 + θ i-1,j + θ i,j-1 -4θ i,
θ t+∆t i,j -θ t i,j ∆t = D 2h 4 θ t i+1,j + θ t i,j+1 + θ t i-1,j + θ t i,j-1 -4θ t i,j + θ t+∆t i+1,j + θ t+∆t i,j+1 + θ t+∆t i-1,j + θ t+∆t i,j-1 -4θ t+∆t i,j ].
(4.17)

Here a problem arises as the terms originated from the f t+∆t i,j function are not known making this an implicit method, these result in the need to solve a system of linear equations that are highly time consuming. As an alternative I realized that I could obtain good accuracy by assuming that θ t+∆t i,j is given by the advected-reacted density field at time t + ∆t (as diffusion is implemented last in each time step), hence one could dub this numerical scheme a pseudo-Crank-Nicholson method. The Eulerian scheme shown here has to respect some constraints [START_REF] Sandulescu | Biological activity in the wake of an island close to a coastal upwelling[END_REF]; interpolation induces a numerical diffusion of the order D n ∝ h 2 /∆t, it is therefore desirable to have D > D n . In addition, the stability condition for the Eulerian diffusion step requires that D∆t d /h 2 < 1, where ∆t d denotes the time step for the diffusion calculations. Both conditions imply that ∆t d < ∆t, so for example, if ∆t d = ∆t/10 then after every advection and reaction steps in the algorithm (that solves the advectionreaction-diffusion equation) one should run 10 steps of diffusion (our choice in this work). Regarding the previous boundary problem, here no-flux boundary conditions ([vθ -D∂ z θ] z=boundaries = 0) can be better dealt with by being discretized directly with a 3-points forward (top boundary) or backwards (bottom boundary) finite-differences scheme (here the terms backwards and forward refers to space, different from the schemes of the previous section). So after the rest of the grid was updated we calculate the scalar values at the boundary nodes, this way:

θ t+∆t i,j=0 = D 4θ t+∆t i,j+1 -θ t+∆t i,j+2 2vh + 3D (top) (4.18)
θ t+∆t i,j=N z = D θ t+∆t i,j-2 -4θ t+∆t i,j-1 2vh -3D (bottom). (4.19)
The difference between methods is that the forward-Lagrangian scheme is more easily applied to non-divergent flows; otherwise, accumulation of Lagrangian fluid parcels into convergence regions and rarefaction of the particle concentration in divergence regions can lead to a spatially variable resolution in the scalar concentration fields [START_REF] Pasquero | Impact of the spatiotemporal variability of the nutrient flux on primary productivity in the ocean[END_REF]. Also the extrapolation adopted for the time evolution in the forward scheme contributes for the generation of less smooth scalar fields. I extensively tested both the forward-Lagrangian and the Eulerian scheme for the solution of our problem, confirming the above statement. Fig. 4.3 shows an extreme case for the model described in Sec. 4.4, where dynamics are calculated with the forward semi-Lagrangian scheme (a) and the backwards scheme with the Eulerian method for diffusion (b) (both with the same grid resolution). We see that forward method results in a much more granular density field after interpolation, whereas in the latter the Eulerian diffusion scheme produces a smooth field. In terms of general results such as average total population density or growth measurements in time, both methods behave very similarly, (a) being fairly quicker (around 2 times faster) in computational time. In spite of that, the necessity for precise spatial analysis makes the field smoothness indispensable. We finish the numerical methods description with a flowchart that summarizes the two numerical method implementations described so far (Fig. 4.4).

Phytoplankton Model Dynamics and Numerical Validation

As expressed in chapter 3 and chapter 2 the history of research on vertical phytoplankton dynamics models is quite recent (less than a century) but rather rich. What is common to many previous studies is that they mainly focused on biological factors, only quite crudely representing fluid motions. The question of plankton blooms was particularly addressed. The concept of critical depth was introduced and investigations explored the role of the mixed-layer depth (MLD) [START_REF] Sverdrup | On conditions for the vernal blooming of phytoplankton[END_REF][START_REF] Huisman | Critical depth and critical turbulence: two different mechanisms for the development of phytoplankton blooms[END_REF], turbulence intensity [START_REF] Riley | The relationship of vertical turbulence and spring diatom flowerings[END_REF][START_REF] Huisman | How do sinking phytoplankton species manage to persist?[END_REF], light-limited growth [START_REF] Platt | Mathematical formulation of the relationship between photosynthesis and light for phytoplankton[END_REF][START_REF] Shigesada | Analysis of the self-shading effect on algal vertical distribution in natural waters[END_REF], to cite a few. In addition, the biological modelling was debated and it still is to some extent when dealing with phytoplankton vertical dynamics. The problem is settled on how to express complex biological relations into a compact mathematical framework. Between those relations are the ones referring to phytoplankton survival directly; nutrient limitations, photosynthesis and biomass creation/depletion and predation by zooplankton. In addition, it is worth stressing once more that mathematical models of biological and physical dynamics are typically difficult to solve analytically because of their nonlinear character and coupling. In the following I revise the 1D model [START_REF] Huisman | How do sinking phytoplankton species manage to persist?[END_REF] that was the basis for our study and perform a validation study with the numerical methodology presented so far. In there I also discuss the findings of this 1D study that are relevant for the comprehension of the 2D work we perform in the next chapter, which can be seen as an extension of former. I end this chapter with the 2D extension of the 1D model and further validation results for the 2D numerical approach.

Light-Limited Sinking Phytoplankton Model

In this extent the work of Huisman and Weissing [START_REF] Huisman | Light-limited growth and competition for light in well-mixed aquatic environments: an elementary model[END_REF] was of paramount importance, they assumed that light-limited phytoplankton dynamics could be examined Phytoplankton sinking velocity in a mono-culture environment with homogeneous nutrient availability. So by neglecting much of the complexity of real-world systems they [START_REF] Huisman | Critical depth and critical turbulence: two different mechanisms for the development of phytoplankton blooms[END_REF] later set of parameters based on measurements from "average" freshwater phytoplankton species in the culture collection of the Laboratory for Microbiology, at the University of Amsterdam. With this methodology they were able to derive a relatively simple model. This was done by taking advantage of the well known advection-reaction-diffusion models.

Our approaches are the ones described in chapter 3, at first we relied solely on the forward advection technique (Sec. 4.2) to reproduce Huisman 1D results [START_REF] Huisman | How do sinking phytoplankton species manage to persist?[END_REF] and to extend them to 2D in the absence of a flow. So recalling the model 1D evolution equation:

∂θ(z, t) ∂t = [p(I(z, θ, t)) -l]θ(z, t) + D ∂ 2 θ(z, t) ∂z 2 -v ∂θ(z, t) ∂z , ( 4.20)
where;

p(I(z, θ, t)) = p max I(z, θ, t) H + I(z, θ, t) , ( 4.21)
I(z, θ, t) = I in e -z 0 kθ(σ,t)dσ-K bg z , (4.22)
and where v is the phytoplankton sinking speed, θ(z, t) as the population density (cells per unit volume) at depth z and time t, p(I(z, t)) is the specific production rate as an increasing function of light intensity I(z, θ, t), l is the specific loss rate, p max is the maximal specific production rate and D is the turbulent diffusion coefficient. The light intensity is modeled by the Lambert-Beer law with the incident light intensity denoted I in , K bg is the total background turbidity (responsible for light absorption by the water medium) [START_REF] Platt | Mathematical formulation of the relationship between photosynthesis and light for phytoplankton[END_REF], H is the half-saturation constant of light-limited growth and k is the specific light attenuation, accounting for the self shading of plankton [START_REF] Shen | Improved efficiency of the biological pump as a trigger for the Late Ordovician glaciation[END_REF]. Additionally the system is subject to no flux boundary conditions:

vθ(z, t) -D ∂θ(z, t) ∂z = 0, at z = 0 and z = z max , (4.23)
where z max refers to the lower boundary of the water column that represents the mixed layer depth. Furthermore we take advantage of the biological parameters mentioned at the beginning of the section, reported in table 4.1, these will be used throughout all studies in this and in the next chapters.

For illustration purposes a surface representing the specific net growth rate (p(I)l) as a function of θ and z) for the parameters of table 4.1, can be appraised in Fig. 4.5 where an arbitrary population distribution is used. There, when the population density is small, we see that the reaction rate has a form imposed by Beer-Lambert law, but when the population density increases also the effect of the self-shading becomes apparent. For higher values of θ, the depth at which growth is positive shrinks, a factor that may play a role when the total population (biomass) (P tot (t))

P tot (t) = L z 0 θ(z, t)dz (4.24)
is considered.

Model Exploration in 1D

We move on to investigate the time evolution of equation 4.20 with 4.21 through computer simulations. We start from a homogeneously distributed population corresponding to θ ≈ 10 6 cells/m 3 , as we expect the asymptotic dynamics not to depend on the choice of the initial population density (see also [START_REF] Huisman | How do sinking phytoplankton species manage to persist?[END_REF]), and this choice allowed us to accelerate the approach to the later time behavior of the system. In Figures 4.6 and 4.7 we see that after a period of exponential growth (for bloom cases) P tot saturates and maintains a constant value, a consequence of the self-shading presence (without it population would grow indefinitely). We also performed tests to verify the convergence of the method for the 1D simulations by measuring P tot for different time steps dt, in Fig. 4.6 we can see an example of such tests. For a depth (z max ) of 20m we found that dz = 0.133m was an appropriate spatial discretization, given that the length over which light decreases is 1/K bg = 5m, we were able resolve well the light-limited dynamics. We then measured the largest time step we could use the ones of the original work [START_REF] Huisman | How do sinking phytoplankton species manage to persist?[END_REF]. These are here presented in the exact same form as in [START_REF] Huisman | How do sinking phytoplankton species manage to persist?[END_REF] to appreciate the quantitative agreement with those previous results. They are useful to observe the effect of mixing (turbulent diffusivity) over the population through depth: as D decreases the population distribution becomes less homogeneous and more concentrated in the upper part of the water column, the zone with the highest growth rate. We can also observe surface population density maxima for the smaller D cases. They become apparent when growth rates exceed both sinking rates and mixing rates, and are a manifestation of the interplay between mechanisms (light limitation and sinking) acting in opposite directions along the vertical.

2D Phytoplankton Model Extension and Comparison

It is worth remarking again that our approach led to less smooth vertical profiles than the original one [START_REF] Huisman | How do sinking phytoplankton species manage to persist?[END_REF], this is clear from Fig. 4.9 where the curves are somehow "noisy". We solved this problem when validating the 2D code using the approach previously described, namely by resorting to backwards advection and to finite-differences discretization of the diffusion term. Equation 4.20 in 2D reads:

∂θ(x, z, t) ∂t = [p(I) -l] θ(x, z, t) -v • ∇θ(x, z, t) + D∇ 2 θ(x, z, t), (4.25)
with v = v ẑ, where ẑ is the unitary vector pointing downwards in the vertical direction. We solve equation 4.25 through numerical simulations starting from a homogeneous θ field, in the same manner as equation 4.20 and, since dynamics are independent of the horizontal direction x, we expect the same results. In addition to the previous vertical boundary conditions, we chose periodic boundary conditions for the horizontal direction.

Repeating the same process specified in the section for the 1D dynamics we again computed vertical profiles and compared them to the previous results. For this, we choose the highest and the lowest value of D in Fig. 4.9 and we present the comparison with the 1D results in Fig. 4.10. Results show good agreement, which validates also the 2D code. In addition, it is clear that the chosen method for the 2D numerical simulations produces much smoother profiles, which revealed essential for the analysis of the spatial structure of the population density field in the presence of a flow, as it will be shown in the next chapter. In this chapter I described the numerical methodology utilised during this thesis to solve the advection-reaction-diffusion equation representing the phytoplankton dynamics vertical model. First I introduced the way of implementing advection through a Lagrangian scheme. Then I explained two different ways of implementing diffusion, the first through stochastic kicks and the second through spatial finite differences. The second method proved more reliable for the production of smooth scalar fields, necessary for the scalar spatial analysis that will be presented in the following study. I finished the chapter with a review of the 1D model [START_REF] Huisman | How do sinking phytoplankton species manage to persist?[END_REF] that served as basis for the work in this thesis, there I also validated the solution provided by numerical methods I presented previously through results comparison, both for 1D and 2D. Now that I established that the chosen numerical methods are appropriate for the presented problem, I will present an extension to the study last shown in this chapter [START_REF] Huisman | How do sinking phytoplankton species manage to persist?[END_REF] in the next chapter.

.

Chapter 5

Effects of large-scale advection and small-scale turbulent diffusion on vertical phytoplankton dynamics

Introduction

In this chapter we present a study of the different transport mechanisms that can affect phytoplankton growth and spatial distribution. The chapter is organized as follows: In the three forthcoming sections we will discuss the same biological model presented in Chapter 4 but with the inclusion of a kinematic flow field of different degrees of complexity [START_REF] Lacorata | Chaotic Lagrangian models for turbulent relative dispersion[END_REF][START_REF] Boffetta | Nonasymptotic properties of transport and mixing[END_REF]. This has the objective of investigating both the impact of different scale transport mechanisms on the modelling of turbulence, as well as the influence of large scale advective coherent structures in the previous purely turbulent diffusive problem.

Therefore the three upcoming sections will refer to a stationary (Sec. 5.2) flow case, a time-dependent (Sec. 5.3) flow case and a multiscale (Sec. 5.4) flow case. These cases are based on a flow configuration with two large scale vortices, that will be better explained later on, and are respectively characterized by a case where only the stationary coherent structures (vortices) are present, a case where these structure are able to oscillate horizontally with time, and a case where we study the mentioned oscillating structures with the addition of smaller scale structures (in space), which are able to oscillate horizontally but also vertically. These worked as an idealized model of flow circulation structures that are seen both in the ocean and lakes such as the Langmuir circulation and winter convection [START_REF] Thorpe | The turbulent ocean[END_REF][START_REF] Vallis | Atmospheric and oceanic fluid dynamics[END_REF]. With this study we aim at contributing to a better understanding of the modelling of phytoplankton vertical dynamics and the importance of the representation of different transport mechanisms. This improves previous studies where the quantities of interest were estimated by models that accounted only for turbulence by turbulent diffusivity [START_REF] Huisman | Population dynamics of sinking phytoplankton in light-limited environments: simulation techniques and critical parameters[END_REF][START_REF] Ebert | Critical conditions for phytoplankton blooms[END_REF].

Steady Large-Scale Flow

Introduction

Studies discussing the influence of turbulence and horizontal advection over plankton cycles, patchiness and survival through modelling have shown that advection has considerable impact over the population evolution and distribution [START_REF] Bracco | Horizontal advection, diffusion, and plankton spectra at the sea surface[END_REF][START_REF] Koszalka | Plankton cycles disguised by turbulent advection[END_REF]. Chapter 5. Effects of large-scale advection and small-scale turbulent diffusion on vertical phytoplankton dynamics Moreover both simulations through Lagrangian tracers transport [START_REF] Backhaus | Phyto-convection': the role of oceanic convection in primary production[END_REF] and largeeddy simulations (LES) [START_REF] Taylor | Shutdown of turbulent convection as a new criterion for the onset of spring phytoplankton blooms[END_REF] have obtained success in replicating experimental and remote sensing data with the aid of proper meteorological modelling. A common feature of the upper ocean is that many buoyancy/wind driven currents such as Langmuir circulation and convection can be crudely approximated by cellular flow motions [START_REF] Thorpe | The turbulent ocean[END_REF], providing a valuable simplification on the study of the general impact of these currents on phytoplankton populations. With this in mind, in this section we implement a kinematic flow-field model over the base 2D model described in the past section.

Phytoplankton dynamics in a Flow Field

Beginning from equation 4.25 we include the flow-field mathematical form inside the velocity term v, meaning now that advection is realized by a 2D incompressible flow u = (u x , u z ) and thus v = u + v sink ẑ. More precisely, in this section we consider advection by a prescribed cellular flow, which on the following sections is intended to model the presence of eddy fluid motions on different scales, but in this sections is studied solely for its large scale properties. The velocity field is then obtained as u = -∂ z Ψ st L , ∂ x Ψ st L from a streamfunction that is written as:

Ψ st L (x, z) = - U k sin(kx) sin(kz), (5.1)
where the superscript and subscript in Ψ st L (x, z) represent respectively the stationary/steady (no time dependency) nature of the flow, and the large scale representation. This is done to better separate the flow from this study from the ones from the remaining studies since they will be directly compared. Also k = 2π/L x = π/L z , where L x and L z refer respectively to the horizontal and vertical dimensions of the system and U will be referred to as the flow's velocity. The plotted streamfunction can be assessed in Fig. 5.1. Furthermore, in this study and in the remaining ones in this chapter L x = 2L z as to maintain the two large-scale vortices configuration studied here. Additionally we consider a depth for which the turbulence window exists for the no-flow system (L z ≥ 60 m), as documented in [START_REF] Huisman | How do sinking phytoplankton species manage to persist?[END_REF], and we fix the turbulent diffusivity to a value that is intermediate between the minimum and maximum critical ones for blooming. Due to the computational cost of simulations in larger spatial domains, we choose a depth value close to the minimum possible one, namely L z = 60 m.

A relevant feature of the original model is the existence of a turbulence window allowing for phytoplankton bloom, for large enough system depths. As cited in chapter 3, determining an analytical expressions for the critical conditions for population survival (i.e. blooming) or extinction is not an easy task, even in such a simple model. This difficulty is due to the heterogeneity of the environment and is common to different population dynamics' models (see, e.g., [START_REF] Speirs | Population persistence in rivers and estuaries[END_REF][START_REF] Ryabov | Population growth and persistence in a heterogeneous environment: the role of diffusion and advection[END_REF][START_REF] Vergni | Invasions in heterogeneous habitats in the presence of advection[END_REF] for other 1D systems). Adopting some simplifying assumptions, it is possible to obtain an approximate estimate of the minimum turbulent diffusivity (the lower bound of the turbulence window) required to compensate the sinking of phytoplankton, and hence to let the population survive [START_REF] Riley | Quantitative ecology of the plankton of the western North Atlantic[END_REF][START_REF] Shigesada | Analysis of the self-shading effect on algal vertical distribution in natural waters[END_REF][START_REF] Speirs | Population persistence in rivers and estuaries[END_REF][START_REF] Huisman | How do sinking phytoplankton species manage to persist?[END_REF]. Nevertheless, for the maximum turbulent diffusivity (the upper bound of the window), beyond which the population cannot outgrow the turbulent mixing rate to sustain the bloom in the upper part of the water column, no simple analytical expression is known [START_REF] Ebert | Critical conditions for phytoplankton blooms[END_REF][START_REF] Huisman | How do sinking phytoplankton species manage to persist?[END_REF]. Hence we focus on the upper bound of the turbulence window (taking D ≤ 20 cm 2 s -1 ), for which numerical simulations reveal more useful.

Effects of Steady Large-Scale Advection

Starting from a homogeneous θ field initial condition we begin our investigation by varying the flow intensity U in a broad range to examine possible changes of behavior due to advection by the coherent flow. We started by measuring the 2D total biomass (per unit of length) P tot ,

P tot (t) = L x 0 L z 0 θ(x, z, t)dzdx = L x L z θ (t). (5.2)
where θ (t) is the average biomass density. In Fig. 5.2 we first notice that the presence of the flow field monotonically hinders population growth proportionally to its intensity. We also notice that simulations did not achieve a stationary converged values (as in the previous no-flow simulations) for P tot at reasonable times, with the exception of the weakest flow intensity U = 0.493m/h. Moving on we decided to base our measurements on the per-capita growth rate (see, e.g., [START_REF] Murray | Mathematical biology: I. An introduction[END_REF]):

r p (t) = 1 θ ∂ θ ∂t , (5.3)
This was possible for the reason that after an initial transient, the latter quantity attained a statistically constant value r p , corresponding to exponential growth (r p > 0) or decay (r p < 0) in the early regime before the onset of nonlinear dynamical effects (due to self-shading). Note that we also use r p to represent time averages of r p (t) over the entire simulation (this was useful in the time dependent flow cases from next sections), in simulations of duration T such time average can be expressed as: In Fig. 5.3 we expose the spatial influence of the flow-field for six different cases, in the first column we have a) no-flow, b) U = 1.24m/h and c) U = 4.93m/h. The first shows the dynamics of the base model where population accumulates in the top layer of the system where growth rates surpass sinking and mixing rates. In b) we see that the spatial effect of the flow is to accumulate population as an elongated central filament in the convergence area of the rolls. In c) we see that as U increases the density field appears more well mixed, filling the outskirt zone of the flow's vortices, this was a no-bloom case (we will later see that the bloom/no bloom transition is controlled by the advection-reaction time scale ratio). So for the majority of bloom cases the snapshots of the density field are similar to the b) case, although diffusivity plays a role on this as well. In d), e) and f) we expose the effect of increasing D on cases with the same U as b); from top to bottom on the second column the result of increasing D is shown as more homogenized fields that resemble more the no-flow case as D increases. One can say that diffusivity partially counterbalances the flow's effect in this sense.

r p = 1 T ln θ (T) θ (0) . ( 5
Before showcasing the r p behavior for the different study cases, we investigate the influence of the initial condition on the dynamics we saw in Fig. 5.2. In the original base model, results are shown to converge to the same quantities independently of the initial condition as they depend solely on the biological parameters, D and the MLD [START_REF] Huisman | How do sinking phytoplankton species manage to persist?[END_REF]. In Fig. 5.4 results were obtained using a localized initial condition corresponding to a small patch of population density located in the central-upper part of the domain (close to z = 0, x = L x /2), here we verified that the overall phenomenology stays unchanged when considering a uniformly spread initial population (Fig. 5.2). The coherent flow reduces the growth of θ (t) and eventually causes an extinction when its intensity is large enough. The growing or decaying temporal behavior is already established after one large eddy turnover time (here estimated as 2πL z /U, approximating streamlines with perfectly circular orbits of radius L z , see vertical lines in Fig. 5.4. Indeed, after this time interval, the initial population density results to be either increased or decreased, depending on the magnitude of U. At later times, the total population appears to grow exponentially at a constant rate, as it is also noticeable in Fig. 5.2. Results then show that the critical flow intensity U responsible for turning a bloom case into a no-bloom case is measurable with good

accuracy.

In order to characterize the bloom to no-bloom transition induced by advection, we measure the per-capita growth rate r p (t) (see Eq. (5.3)). This quantity, normalized by the intrinsic total (birth minus death) growth rate at the surface r b = I in /(H + I in )p maxl, versus time normalized by L z /U is shown in Fig. 5.5a. Here, a uniform initial population density was chosen. As it can be seen, at large enough times, for all U, r p (t) approaches an asymptotic constant value r p , confirming the exponential character of growth or decay of θ (t). Furthermore, the large-time limit r p decreases from positive values (for low U) to negative ones (at larger U), therefore allowing a robust estimate of the critical flow intensity at the transition. We remark that we could not detect a transition to a no-bloom regime for D < 5 cm 2 s -1 , even with very large values of U, as it can be seen in Fig. 5.5b. There we show that for D = 1 cm 2 s -1 the dynamics take longer to achieve a stationary growth regime, a proof that higher diffusivity simulations homogenize quicker, achieving such state much faster. Nevertheless it is clear that increasing the flow speed U in this case only inhibit growth up to a certain value, after which r p /r b tends to converge to similar values. With this in mind we proceed to run multiple simulations, again, starting from a homogeneously spread θ initial condition and L z = 60 m, with 5 cm 2 s -1 ≤ D ≤ 20 cm 2 s -1 and multiple U values near the transition zone as in Fig. 5.5. As first indicated in [START_REF] Abraham | The generation of plankton patchiness by turbulent stirring[END_REF], where horizontal patchiness was numerically studied adopting an NPZ (for nutrient-phytoplankton-zooplankton) model in a turbulent flow, we expect that also in the present case the dynamics are primarily controlled by the interplay between advection and reaction mechanisms. To quantify the relative weight of the latter, we consider the ratio of the biological time scale r -1 b to the flow time scale L z /U, i.e.:

γ = U r b L z . ( 5.5)
Figure 5.6 reports the (asymptotic) per-capita growth rate as a function of γ. From this plot, one can clearly see that the survival/extinction transition caused by the flow occurs for γ = O(1), in correspondence with r p turning from positive to negative. Essentially, a bloom can take place (r p > 0) when the biological growth is faster than the advective transport (γ < 1) to the less favorable deeper part of the domain. The proximity of the data obtained with different values of D highlights the generality of this mechanism and confirms the weak effect of the turbulent diffusivity in this picture. In Fig. 5.7 we display the values for the transition velocity U transition from bloom to no-bloom for multiple D, taken from the same simulations as in Fig. 5.6. There it is easier to verify the previous statement as it is clear that the flow velocity threshold value, that triggers the transition, changes very little as we increase D, even though it increases monotonically. So for the range of D values studied we found that both the bloom/no-bloom transition advection intensity U transition (Fig. 5.7) and the spatial distribution of the population density (Fig. 5.3) were only weakly affected by varying the diffusivity. This leads us to search for a better explanation on how to compare the slight changes reported on the density spatial distribution of different D and γ cases, and their relation to the transition bloom/no-bloom.

The Filament Model

Further insight on the mentioned topic comes from inspection of the spatial structure of the population density field θ(x, z, t)/ θ at a given time (Fig. 5.3). While in the absence of flow the population is uniformly distributed along the horizontal and decreases with depth, nonzero advection causes an increase of θ in the downwelling region (at x = L z ). This feature gets accentuated by increasing U, with the population accumulating in thinner and thinner filaments outside vortices, and particularly in the one located at x = L z . Such a behavior points to the relevance of strain-dominated flow regions for the spatial organization of the population and the formation of fine structures. In our flow, as it can also be easily seen in Fig. 5.3 (b to f panels), the latter regions are close to the hyperbolic points corresponding to the vertices of the squares of side L z containing the rolls. Among such points, clearly, a prominent role is played by the point (x, z) = (L x /2, 0), where the flow locally compresses the scalar field θ along the x-axis (and stretches it in the z-direction), in the region of highest growth rate (i.e. at the surface). Relying on the above picture, a useful interpretation of the dynamics observed in our simulations is offered by an appropriate adaptation of the plankton filament model [START_REF] Martin | On filament width in oceanic plankton distributions[END_REF], originally introduced to describe the formation of fine structures in 2D flows. To apply this reasoning, we neglect the sinking speed, which is considerably smaller than the advecting velocity close to the transition to extinction, as well as self-shading, as close to an extinction the population density is low everywhere and because our main point of interest is at the surface. Under these hypotheses, Eq. (4.25) (with the flow) becomes:

∂ t θ + u • ∇θ = I in H + I in p max e -κ bg z -l θ + D∇ 2 θ. ( 5.6)
Since, as argued above, we are interested in the dynamics at the surface, close to the hyperbolic point at x = L z , the net growth rate will be r b . Moreover, we can write the population density as θ(x, z, t) = θ back exp r b t + θ (x, z, t), where θ back is the background population density and θ represents the perturbation determined by the flow. Using this decomposition, it is not difficult to see that the equation for θ dynamics is the same as Eq. (5.6). Following [START_REF] Martin | On filament width in oceanic plankton distributions[END_REF], we can then consider only the 1D dynamics for the population fluctuation in the compressing (or cross-filament) direction, because along the filament θ should vary less due to the stretching operated by the flow. In this region, the flow can be locally approximated as u = (-λ(x -L x /2), λz), where λ = kU ≈ U/L z is the strain rate. Therefore, from Eq. (5.6), for the cross-filament dynamics one has:

∂ t θ -λ (x -L x /2) ∂ x θ = r b θ + D∂ 2
x θ .

(5.7)

The solution of the above equation is (see also [START_REF] Martin | On filament width in oceanic plankton distributions[END_REF]):

θ (x, t) = θ 0 e -(x-Lx /2) 2 λ 2D e (r b -λ)t , (5.8)
with θ 0 a constant. From this expression we can see that, in the x-direction, the population density field keeps the same (Gaussian) shape at different times. The filament width σ = √ D/λ does not depend on time and is only determined by the physical parameters associated with fluid transport. As the flow intensity increases, so does the strain rate, which explains the thinning of filaments and the more important localization of the population for higher values of U. Growth or decay over time, instead, depends on whether r b is larger or smaller than λ, respectively. This simple model thus provides theoretical support to the survival/extinction criterion based on the ratio between the biological and flow time scales, γ.

The above model accounts for the dynamics at the surface and, strictly speaking, it is only there that its predictions should apply. If the population cannot survive at the surface, however, it should not deeper below either due to the reduced growth rate, which makes the conclusion appear more general. Considering that, differently from the 1D filament case, in our fully 2D model both the strain rate and the growth rate vary with depth, and that sinking and self-shading might also play a minor role, the comparison between our previous estimate of the control parameter, γ = U/(L z r b ), and that from Eq. (5.8), λ/r b , seems to us reasonable also from a quantitative point of view. Regarding the dependence on the vertical coordinate, we further note that the biological growth rate monotonously decays with z, and that the strain rate, in absolute value, decreases until half the total depth, before growing again in the lower half of the domain, but now acting in the opposite way (stretching instead of compressing the scalar in the x-direction). The combination of these effects, impacting both the width and the intensity of the filament, can then explain, in a qualitative way, the tendency, visible in panels b), d), e) and f) of Fig. 5.3, of this localized downwelling structure to fade around z = L z /2. Such effect can be better appraised in Fig. 5.8, there we show multiple population density horizontal profiles θ(x, z * , t * ) normalized by their mean value θ(x, z * , t * ) , for multiple D and U = 1.24m/h (these are the same cases exposed in Fig. 5.3 b, d, e and f). It is clear that, as we approach the flow's maximum vertical speed region (z = L z /2), the profiles present higher population density concentration at x = L x /2, and then the filament configuration tends to fade as profiles homogenize for z > L z /2. Here we also notice the role of diffusivity in counterbalancing the filament concentration by promoting more homogeneous horizontal profiles in the higher diffusivity cases. To test the validity of the above argument for our system, we examined the horizontal profiles of population density at z = 0 from simulations with different values of D and U, once r p (t) converged to its asymptotic value r p . We found that such profiles are to good extent time independent and that their shape is reasonably described by a Gaussian function. In Fig. 5.9a we show an example of the time invariance of these profiles after the stationary r p regime is achieved, there for D = 20cm 2 s -1 and U = 1.01mh -1 we see that the normalized population density horizontal profiles remain unaltered for multiple non-dimensional times. After verifying the stationary nature of these profiles, we decided to proceed by performing a fit through a Gaussian function in a subregion centered around x = L z of each profile. We defined a window of x = [50; 70]m for the surface (z = 0) population density horizontal profiles data and then we estimated the standard deviation of the Gaussian curves fitted to the mentioned data by a non-linear least squares method (see Fig. 5.9b for an example). These provided a measure of the filaments widths.

The results are compared to the theoretical prediction in Fig. 5.10, which indicates a strong correlation between the numerical and theoretical estimations of σ. As one can observe in the figure, we actually detect a tendency of the numerically estimated σ to grow slightly faster than the theoretical one. However, such a small difference seems quite reasonable, taking into account the assumptions made for the theoretical prediction with respect to the details of our numerical setup. We further mention that multiple fitting intervals for the horizontal profiles were tried apart from the final choice of x = (50-70)m ( i.e. x = (0-120)m, (20-100)m and (40-80)m), all present in Fig. 5.9b. It is worth remarking that at the bottom, due to the similar structure of the deep and surface flow, the spatial organization of the reactive scalar parallels that found at the surface but with a more homogeneous distribution as seen in Fig. 5.8. An analogue reasoning in the straining regions close to x = 0 or x = 2L z would always give extinction locally, however, as the prefactor in the exponent of the exponential involving time would be -lλ < 0, assuming that the growth rate is negligibly small there. Hence, the relatively high values of population density at the bottom appear to be due to fluid transport (including sinking) only and the zero-flux boundary conditions. Finally, according to Eq. (5.8), in the limit of very small diffusivity, the filament width approaches zero while its density amplitude grows exponentially. Consequently, it becomes more and more difficult to observe an advection-driven extinction. From a practical point of view this is essentially impossible in numerical simulations as it would require an infinite spatial resolution in order to resolve the crossfilament structure. These are likely the reasons why we could not detect the transition to no-bloom at sufficiently small values of D.

Summary

In this section we expanded the 2D phytoplankton dynamics model presented in the past section by including explicit flow field advection in the form of two large scale vortical coherent structures. The addition of such structures into the previously purely diffusive/sinking transport was shown to inhibit growth proportionally to the flow speed of the flow rolls. The effect presented by these structures was to monotonically decreased per capita growth rates to the point where phytoplankton bloom cases (in the no-flow original dynamics) could be turned into no bloom cases for intense enough flows (for multiple D in the range 5cm 2 s -1 ≤ D ≤ 20cm 2 s -1 ). The bloom/no bloom transition results were quantified and presented in the form of the reactive to advective time scales ratio γ, exhibiting that the transition occurred at γ = O(1) independently of the choice of D, as the diffusivity showed weak influence over the transition dynamics, mostly dominated by γ.

The influence of the flow structures was also assessed spatially, first through simulation snapshots and then through horizontal population density profiles. With these it was possible to realize that the presence of the flow field triggered the formation of a persistent vertically elongated central plankton filament, localized on the separatrix between the flow rolls. By means of Gaussian fits we were able to measure the filament widths and compare them to a theory presented by Martin [START_REF] Martin | On filament width in oceanic plankton distributions[END_REF]. Considering that the original theory accounted for 1D dynamics without sinking, our numerical results presented a good match. This unveiled a more robust explanation for why the bloom/no bloom transition is controlled by γ, presenting also an argument for why the bloom/no bloom transition could not be measured for smaller diffusivities (i.e. D < 5cm 2 s -1 ). This now set foot for the next stage of the work which was to further add complexity to the flow field, in hopes to increase the realism of the simulations.

Unsteady Large-Scale Flow

Here we extend our large scale advection study over the phytoplankton lightlimited growth dynamics by including explicit time dependency over the largescale advective transport. This is done by allowing for lateral oscillations of the flow pattern considered in the previous section, which is enough to produce chaotic Lagrangian trajectories of fluid particles [START_REF] Vulpiani | Chaos: from simple models to complex systems[END_REF][START_REF] Neufeld | Chemical and biological processes in fluid flows: a dynamical systems approach[END_REF]. Numerical studies of plankton predator-prey dynamics under chaotic flows [START_REF] Hernández | Sustained plankton blooms under open chaotic flows[END_REF][START_REF] Neufeld | Excitable media in open and closed chaotic flows[END_REF] shown that when time scales associated to fluid stirring become slower than the ones associated to biological growth a phytoplankton bloom is sustained. Even though conditions are fairly different from our study (multiple biological species, excitability) it is interesting that we have shown dependency on the equivalent ratio (for our study) for light-limited dynamics in the previous section (steady large-scale flow). Consequently, in this section we hope to investigate if the phenomena presented in the last section matches what is observed when chaotic trajectories are added to the advective transport or if substantial changes are observed.

Adding Flow Time Dependency

Here we follow the same procedure of last section for the inclusion of the largescale incompressible flow-field into the base model 2D equations. The difference here is that the prescribed cellular flow of equation 5.1 now operates under periodic lateral oscillations, vertical oscillations are not accounted for here to respect the noflow boundary conditions of equation 4.23. Hence the velocity field is still described by a streamfunction as u = (-∂ z Ψ L , ∂ x Ψ L ), but now it represents a periodic flow field:

Ψ L (x, z, t) = - U k sin {k[x -s sin(ω 1 t)]} sin(kz), (5.9)
where one can clearly see that when t = 0 we have the same streamfunction of equation 5.1, so the same visualization as in Fig. 5.1, whereas when t > 0 the flow rolls start oscillating with amplitude and pulsation respectively set to s = L z /5 and ω = πU/L z , corresponding to a fraction of the roll size and a period comparable to the advective time scale L z /U, a choice that has been shown to be optimal to enhance chaotic diffusion [START_REF] Solomon | Chaotic particle transport in time-dependent Rayleigh-Bénard convection[END_REF][START_REF] Vulpiani | Chaos: from simple models to complex systems[END_REF][START_REF] Lacorata | Chaotic Lagrangian models for turbulent relative dispersion[END_REF][START_REF] Lacorata | FSLE analysis and validation of Lagrangian simulations based on satellite-derived GlobCurrent velocity data[END_REF]. Note that we took away the superscript of Ψ st L as to differentiate the streamfunction of this section with the one from the last section. A visualisation of the flow under the periodical oscillations can be seen in Fig. 5.11. Chapter 5. Effects of large-scale advection and small-scale turbulent diffusion on vertical phytoplankton dynamics

Steady vs Time-Dependent Flow Comparison

In order to compare the results of the stationary flow with the ones from this section we repeat the analysis made for the first. In Fig. 5.11 we show two snapshots of the population field at different times in the asymptotic per capita growth rate regime. These visualizations suggest that the dynamics are fairly similar to the stationary flow case although horizontal symmetry is now broken due to the lateral oscillations of the flow. We can notice that the central surface plankton patch remains, but now it is subject to strain in a region that oscillates horizontally as the flow varies with time. In Sec. 5.2 we showed that P tot did not achieve a stationary value at feasible simulation times like in the purely diffusive transport cases, nevertheless, the per capita growth rate r p (t) did. Consequently we repeat the same measurements done in Fig. 5.5a, for the same U cases and D = 20cm 2 s -1 and expose them in Fig. 5.12. There we see that for all cases r p (t)/r b reaches a statistically stationary value characterised by small oscillations around a constant value. These oscillations have periodicity proportional to L z /U and tend to decrease in amplitude as U increases, they are clearly a result of the lateral oscillations present in the unsteady flow. The same type of oscillations around a stationary value were reported on numerical simulations for phytoplankton predator-prey dynamics on open chaotic flows in the form of biomass measurements [START_REF] Hernández | Sustained plankton blooms under open chaotic flows[END_REF]. That said, we can see that the overall effect of the flow over the per capita growth rates is still quite similar to the stationary flow case. Additionally we do not achieve extinction cases for D < 5cm 2 s -1 .

The critical advection intensity U transition determining the bloom/no bloom transition is found to be slightly higher in the present time-dependent case for all D. The increase with respect to the previous, steady, case depends on the value of the small-scale diffusivity (about 6% for D = 20 cm 2 s -1 and 18% for D = 10 cm 2 s -1 ), but the dependency of U transition on D remains weak. Considering that the explicit time dependency of Ψ in Eq. (5.9) now gives rise to chaotic diffusion of Lagrangian fluid particles, and hence to an effective diffusivity larger than D, such an increase of U transition seems to us reasonable, from a qualitative point of view. A more quantitative assessment of the comparison between the unsteady and steady flow cases is illustrated in Fig. 5.13. Here we show rp /r b of the steady flow case as a function of rp /r b in the unsteady case, for several values of U and D. As exposed previously, for this comparison we use rp (the time average of r p (t) over the entire simulations) rather than r p (the asymptotic stationary value achieved by r p (t) in the stationary flow simulations) since we have shown in Fig. 5.12 that the growth rate r p (t) now oscillates around a stationary value. Finally, Fig. 5.13 illustrates that over the range of values of D and U explored, the two quantities are almost perfectly correlated, corroborating the idea that the lateral oscillations do not produce any major modifications. The analysis indicates that time dependency of such large scale flow does not have a major impact on the resulting biological dynamics.

Summary

We extended our previous study with advection by two large scale vortex structures by including time dependency in their dynamics. This was made by allowing the flow rolls to oscillate in the horizontal direction. Results showed that the oscillation was reflected in the per capita growth rate dynamics as, in comparison with the past study case, the unsteady flow simulations presented r p (t) curves that still oscillated around a stationary value after achieving an asymptotic regime (as seen in Fig. 5.12). Both flow study cases displayed were characterized by their large scale vortex structures dominating transport dynamics as their weak growth dependency on D indicates. In addition, their dynamics were quite similar, both in terms of quantifiers (see Fig. 5.13) and spatial distribution (see Fig. 5.11), showcasing that the addition of time dependency to the flow in the present study case did not substantially alter the large-scale influence over the biological dynamics in comparison to Sec. 5.2. Hence, we decided to explore the inclusion of additional flow structures into the transport dynamics of the simulations with the purpose of understanding if the picture reflected by these results would remain unaltered.

Multiscale Flow

Introduction

As specified through this chapter, the assumption of past phytoplankton growth models [START_REF] Shigesada | Analysis of the self-shading effect on algal vertical distribution in natural waters[END_REF][START_REF] Huisman | Critical depth and critical turbulence: two different mechanisms for the development of phytoplankton blooms[END_REF][START_REF] Huisman | How do sinking phytoplankton species manage to persist?[END_REF][START_REF] Huisman | Light-limited growth and competition for light in well-mixed aquatic environments: an elementary model[END_REF][START_REF] Ebert | Critical conditions for phytoplankton blooms[END_REF] was that the multi-scale characteristics of turbulence could be approximated by a diffusion process in the form of a turbulent diffusivity term D. In sections 5.2 and 5.3 of this chapter we investigated the influence of an extra element to this modelling, this was the inclusion of the mentioned large scale structures into the advective transport. In this section we implement the last step of this chapter's goal which was to fully understand the influence of turbulent multiscale transport into the assumptions of our dummy turbulence model.

The kinematic model behind our multiscale flow field has already been utilized in phytoplankton [START_REF] Pasquero | Coherent vortices, Lagrangian particles and the marine ecosystem[END_REF] and Lagrangian tracer dynamics [START_REF] Lacorata | FSLE analysis and validation of Lagrangian simulations based on satellite-derived GlobCurrent velocity data[END_REF], mostly to study the transport of these by submesoscale vortices on the ocean surface. In the latter it was shown that the method is quite accurate in comparison to real drifter trajectories measurements. Nonetheless, to our knowledge, the application of this technique to light limited phytoplankton growth models is new. Therefore we present our model with the final form of the stream function Ψ, representative of our flow field, and make a series of comparative studies with the objective of disclosing the importance of the inclusion of multiple smaller scales vortical structures into the turbulence dynamics modelling so far presented.

Model

Here again we depart from equation 4.25 with the inclusion of the flow velocity term to the advection part as u = (-∂ z Ψ, ∂ x Ψ). Since in this section the flow field is time dependent as in Sec. 5.3 we take advantage of the streamfunction defined by equation 5.9 as the first part of our stream function:

Ψ(x, z, t) = Ψ L (x, z, t) + Ψ s (x, z, t), (5.10)
where

Ψ s (x, z, t) = - n k ∑ i=2 U i k i sin {k i [x -s i sin(ω i t)]} × sin{k i [z -β(z)s i sin(ω i t)]}. (5.11)
Here Ψ L represents the large scale persistent structure that is allowed to oscillate in the horizontal, as in Sec. 5.3, and Ψ s represents n k smaller scale vortices that oscillate both in the horizontal and in the vertical. In Eq. 5.11, n k is the number of modes selected, L i = L 1 η 1-i is the typical length scale of mode i, with η > 1 a scale separation factor and L 1 the largest flow scale, k i = 2π/L i the corresponding wavenumber, and U i the typical flow intensity at scale L i . Even though equation 5.9 presents no indexes, it is equivalent to the first mode, e.g. i = 1 (hence why the sum in equation 5.11 starts at i = 2). In addition, we recall that the first mode is not allowed to oscillate in the vertical (whereas the subsequent are), oscillation have amplitudes s i = L i /10 and pulsations ω i = πU i /L i (a choice that is analogous to that of Sec. 5.3, considering that L 1 = L x ). It is important to clarify that we will use U to refer to the flow velocity in spite that it describes only the first mode amplitude (i.e. U = U 1 ). We choose to account for the explicit time dependency of the flow field in the form of oscillations with amplitudes s i and pulsations ω i , parallel to the notations of Sec. 5.3. To respect the no-flux boundary conditions in Eq. 4.23 we use the function:

β(z) = 1 2 tanh z -z 1 ξ -tanh z -z 2 ξ , (5.12)
to dump vertical oscillations near the vertical boundaries and therefore guarantee u z = 0 when z = 0, L z . Parameters values z 1 = 2 m, z 2 = L zz 1 and ξ = 1 m were shown to be adequate for this purpose, the use of such artifact was already proposed for chemical reaction simulations in closed domains [START_REF] López | Efficiency of a stirred chemical reaction in a closed vessel[END_REF]. We choose a number of modes that allows spanning the scale range going from the domain Chapter 5. Effects of large-scale advection and small-scale turbulent diffusion on vertical phytoplankton dynamics size L 1 = L x = 120 m to the smallest length scale L n k = 1 m, corresponding to ≈ 1/(5κ bg ), where κ -1 bg is related to the growth dynamics, as it is the typical length over which light is absorbed by the medium. Such small length scale also roughly corresponds to the scale that can be estimated from Richardson scaling of diffusivity with length, ∼ (2/3) 3/4 -1/4 D( ) 3/4 [START_REF] Lacorata | Chaotic Lagrangian models for turbulent relative dispersion[END_REF][START_REF] Boffetta | Relative dispersion in fully developed turbulence: the Richardson's law and intermittency corrections[END_REF], using the values of diffusivity explored in the previous sections, 5 cm 2 s -1 < D < 20 cm 2 s -1 , and values of the kinetic energy dissipation rate ≈ (10 -8 -10 -6 ) m 2 s -3 that appear reasonable for oceanic turbulence [START_REF] Kiørboe | Planktivorous feeding in calm and turbulent environments, with emphasis on copepods[END_REF][START_REF] Barton | The impact of fine-scale turbulence on phytoplankton community structure[END_REF][START_REF] Lindemann | Dynamics of phytoplankton blooms in turbulent vortex cells[END_REF]. We then set the scale separation factor to η = 2 and the number of modes to n k = 7. Finally, we assume a Kolmogorov scaling of velocity,

U i = U 1 (L i /L 1 ) 1/3 .
A visualization of the of the multiscale flow field streamlines is available at Fig. 5.14, the isocontours of the streamfunction at the same time (black lines in the figure), allow to appreciate the presence of eddies of different sizes and the more disordered spatial structure of the velocity field. There we also display a simulation where, starting from a localized phytoplankton patch in the top central area of the system, we let the system evolve until t * = 845h (in the asymptotic r p regime) for U = 1.01m/h and D = 0.1cm 2 s -1 . Both the initial condition and the small diffusivity setup are chosen for exposition, to allow the reactive scalar to separate into the multiscale streamlines, exhibiting the smaller scale structures of the flow. Notice that even in such a setup the evolved system present, to some extent, the localized central patch configuration reminiscent of the past sections (Figs. 5.3 and 5.11), although we can clearly see the scalar being better spread by the smaller scale modes of the flow. We expose this snapshot in here and not as a result because, maintaining the coherence of past sections, we only obtain quantifiers for 5cm 2 s -1 ≤ D ≤ 20cm 2 s -1 .

Results

We start our results section with another snapshot of the population density field in Fig. 5.15, for a higher diffusivity than what was presented in the model section, D = 5cm 2 s -1 and U = 1.021m/h (streamlines isocontours in black). We see here that in spite of the population density field distribution appearing more irregular than in To further investigate the impact of small-scale fluid motions, we also analyze vertical profiles θ x (z) of the phytoplankton distribution (similarly to what was done in Sec. 4.4 and in [START_REF] Taylor | Shutdown of turbulent convection as a new criterion for the onset of spring phytoplankton blooms[END_REF][START_REF] Huisman | Population dynamics of sinking phytoplankton in light-limited environments: simulation techniques and critical parameters[END_REF][START_REF] A B Ryabov | Vertical distribution and composition of phytoplankton under the influence of an upper mixed layer[END_REF]), obtained by averaging θ(x, z, t) over the horizontal coordinate x at fixed instants of time. Such profiles, normalized by the corresponding global spatial averages θ , are shown in Fig. 5.18 for all the flow cases studied (1-mode steady-flow, 1-mode oscillatory flow, multiscale time-dependent flow) at common given times. Independently of the considered flow or value of D, their shape is always characterized by a maximum at small, but finite, depth and a decrease deeper below the surface, plus a second inflection point close to the bottom boundary. These features are typical for sinking phytoplankton species [START_REF] Huisman | How do sinking phytoplankton species manage to persist?[END_REF] (as discussed in Sec. 4.4), whereas non-sinking ones would display a maximum at the surface [START_REF] Huisman | Critical depth and critical turbulence: two different mechanisms for the development of phytoplankton blooms[END_REF].

The similarity of the profiles obtained in different configurations (Fig. 5.18) points to the dominance of advection by the large-scale coherent flow, as including its time dependence or smaller scales does not alter the general picture substantially. Note, however, that for sufficiently low D, the addition of small scales favors, to small but measurable extent, a localization of the population close to the surface, akin to the vertically non-homogeneous distribution typical of the no-flow case (see Fig. 5.3a) and increased possibility of survival.

The importance of the large-scale flow can be better appreciated by inspecting Fig. 5.19. Here, again for a common fixed time (tU 1 L -1 z = 10) in the asymptotic growth-rate regime, we show the normalized vertical profiles θ x / θ , for the smallest and largest value of diffusivity used (D = 5 cm 2 s -1 and 20 cm 2 s -1 in panels (a) and (b), respectively), for different flow types. Specifically, we examine the following different combinations: Ψ = 0 (no flow), Ψ = Ψ st L (large-scale steady flow), Ψ = Ψ L (large-scale time-dependent flow), Ψ = Ψ L + Ψ s (multiscale time-dependent flow), Ψ = Ψ s (time-dependent flow without the large-scale contribution provided by Ψ L , according to equation 5.11). The last case was explicitly added to test the relevance of the large-scale advection. It is apparent that whenever Ψ L is present the population gets homogenized in the vertical direction, with respect to the no-flow case. The addition of time dependency and small scales to the flow turns out to play only a minor role, as the corresponding profiles are essentially indistinguishable from the one obtained with only Ψ L . When the latter contribution is removed and the flow only possesses smaller scales, instead, the population distribution retrieves the vertically non-homogeneous character typical of the Ψ = 0 case. In such a case, in fact, the vertical profile θ x approaches the one obtained without flow, as it is particularly evident in Fig. 5.19b (where D = 20 cm 2 s -1 ). Finally, it seems to us that Fig. 5.19 summarizes in an effective way the main outcome of this work, meaning the outstanding relevance of advection by the large-scale coherent flow, as the dominant mechanism controlling phytoplankton dynamics in the present setting.

Summary

In the final section of this chapter we presented the last extension to our series of flow field studies over the modelling of phytoplankton light-limited dynamics in turbulent flows. This was carried out by adding smaller scale structures (vortices) to the time dependent large-scale flow field previously presented in Sec. 5.3, this way achieving a better representation of turbulence (as described by Richardson cascade, discussed in chapter 2). In a qualitative way, simulation snapshots showed that the increment of small scales to the flow field was reflected in the spatial spreading of the population density θ as irregular distributions (not seen in the past studies presented in this chapter) were achieved. Nonetheless, quantitative measurements displayed in the form of the normalized average per-capita growth rates rp /r b exposed that dynamics were quite similar to the past steady and unsteady flow field cases. The overall phenomenology remained untouched as the transition bloom/no bloom was still dominated by the large-scale flow field velocity U ( rp decreasing as U increases) with a weak dependence on D. And again we did not report bloom/no bloom transitions for cases with D < 5cm 2 s -1 . At the end of the section we exposed multiple population density depth profiles that furthermore support the argument that dynamics are predominantly controlled by the flow field first mode, the largest scale structure. This was ultimately clarified by the exposition of the θ x depth profiles of cases where the large-scale vortices were not present (but the remainder were) juxtaposed against the no-flow profiles.

Conclusion

In this chapter we numerically investigated the dynamics of sinking phytoplankton in a stirred 2D fluid layer where the vertically decreasing light availability is the only limiting factor for biological growth. For this purpose we extended a previous theoretical 1D model [START_REF] Huisman | Critical depth and critical turbulence: two different mechanisms for the development of phytoplankton blooms[END_REF][START_REF] Ebert | Critical conditions for phytoplankton blooms[END_REF][START_REF] Huisman | How do sinking phytoplankton species manage to persist?[END_REF], where turbulent motions were only described in terms of an effective diffusivity, by taking into account in an explicit way the transport operated by a structured fluid flow. The choice to neglect possible heterogeneities in the nutrient distribution was motivated by our goal to focus on the role of transport mechanisms. While clearly this poses some limitations in relation to real natural environments, where nutrients can also affect biological growth, such a configuration still appears reasonable for, nutrient-rich, eutrophic habitats, namely shallow warm lakes or high-latitude oceans.

A major outcome of the simplified theoretical model [START_REF] Huisman | How do sinking phytoplankton species manage to persist?[END_REF] mentioned above was was to provide evidence of two transitions between extinction and survival, blooming, of the population, depending on the turbulent intensity (for deep enough fluid layers). Our aim, here, was to explore the impact of a more realistic representation of the advecting velocity field on the survival-to-extinction transition, for which no analytical prediction is available, occurring at large turbulent intensity when biological production cannot compensate turbulent mixing to sustain the bloom in the well-lit region close to the surface. Using realistic parameter values for the biological dynamics [START_REF] Huisman | How do sinking phytoplankton species manage to persist?[END_REF], we then considered a domain with a fixed depth representative of the Chapter 5. Effects of large-scale advection and small-scale turbulent diffusion on vertical phytoplankton dynamics mixed layer, in the presence of flows of progressively increasing complexity, relying on a kinematic-flow approach [START_REF] Lacorata | Chaotic Lagrangian models for turbulent relative dispersion[END_REF]. We first examined a velocity field possessing a single large-scale stationary mode, in the form of two recirculating cells spanning the horizontal extent of the system [START_REF] Stommel | Trajectories of small bodies sinking slowly through convection cells[END_REF][START_REF] Solomon | Chaotic particle transport in time-dependent Rayleigh-Bénard convection[END_REF]. Such a spatial structure was intended to mimic the large-scale features observed in realistic flows, as those arising from buoyancy driven convection [START_REF] Solomon | Chaotic particle transport in time-dependent Rayleigh-Bénard convection[END_REF][START_REF] Vallis | Atmospheric and oceanic fluid dynamics[END_REF] or wind-driven Langmuir circulation [START_REF] Stommel | Trajectories of small bodies sinking slowly through convection cells[END_REF][START_REF] Okubo | Diffusion and ecological problems: modern perspectives[END_REF][START_REF] Denman | Time and space scales of vertical mixing and advection of phytoplankton in the upper ocean[END_REF]. We then added time dependence in the form of lateral oscillations of such a flow pattern, and finally included spatially and temporally varying smaller scales.

Our results indicate that advection plays a relevant role on the biological dynamics. Indeed, persistent large-scale motions reduce the per-capita growth rate and can eventually lead to the suppression of the bloom, when the flow is intense enough. This effect is found to be controlled by the ratio between the characteristic biological and flow timescales, similarly to what occurs for plankton horizontal dynamics stirred by mesoscale ocean eddies [START_REF] Mckiver | Plankton bloom controlled by horizontal stirring[END_REF]. From a general perspective, a similar harmful role of the advecting flow was also put in evidence in previous large-eddy simulations of turbulent thermal convection [START_REF] Taylor | Shutdown of turbulent convection as a new criterion for the onset of spring phytoplankton blooms[END_REF], and in a study considering a steady cellular flow and a matrix-based approach to compute the biological growth rate [START_REF] Lindemann | Dynamics of phytoplankton blooms in turbulent vortex cells[END_REF]. However, those studies neglected the phytoplankton self-shading [START_REF] Lindemann | Dynamics of phytoplankton blooms in turbulent vortex cells[END_REF] and also sinking [START_REF] Taylor | Shutdown of turbulent convection as a new criterion for the onset of spring phytoplankton blooms[END_REF]. Moreover, in both of them it is less straightforward than in our work to disentangle contributions from large and small flow scales, either because the latter are essentially absent [START_REF] Lindemann | Dynamics of phytoplankton blooms in turbulent vortex cells[END_REF], or because they dynamically interact with the large-scale ones [START_REF] Taylor | Shutdown of turbulent convection as a new criterion for the onset of spring phytoplankton blooms[END_REF].

The main finding of the present study is that the large-scale flow dominates the dynamics, which is only weakly affected by (temporally and/or spatially) more complex fluid motions. This is revealed by both the strong correlation found for the critical flow intensities (for the transition), and the similar vertical population profiles, in the different flow cases. Even in the presence of a multiscale flow, the velocity field at the largest scale has a strong signature on the dynamics, as it drives the localization of the population in a patch at the center (with respect to the horizontal coordinate) of the domain. This is a region of phytoplankton downwelling, corresponding to the location of the straining point associated with the largest-scale flow mode ( the separatrix between the two largest rolls). A similar accumulation was observed in [START_REF] Lindemann | Dynamics of phytoplankton blooms in turbulent vortex cells[END_REF], and in suspensions of neutrally buoyant particles (compatible with several phytoplankton species that are only slightly heavier than sea water) in vortex-like turbulent flows such as Langmuir circulation ones [START_REF] Barstow | The ecology of Langmuir circulation: a review[END_REF][START_REF] Pershing | Evidence for vertical circulation cells in the well-mixed area of Georges Bank and their biological implications[END_REF].

In the (large-scale) one-mode, steady, flow case, we have been able to rationalize the picture by adapting the plankton-filament model originally introduced in [START_REF] Martin | On filament width in oceanic plankton distributions[END_REF]. This allowed us to provide a quantitative justification for the control exerted by the biological-to-flow timescale ratio on the transition to extinction. As shown by our analysis of the multiscale-flow case, the presence of smaller-scale fluid motions tends to partially disrupt the regular spatial distribution of the population due to the flow at large scale, and the associated central downwelling filament. This was further confirmed by the comparison of phytoplankton-density vertical profiles in two multiscale flows, one of which does not possess the largest-scale mode. Indeed, in the strain region between the largest eddies, the combined action of the flow and of small-scale diffusivity vertically homogenizes the population, thus hindering survival. When only smaller eddies are present, however, the planktonic population localizes closer to the surface, and spreads more over the horizontal, giving rise to a situation resembling that of the no-flow case, which is less prone to extinction.

We hope that the analysis reported here can contribute to the understanding of the basic mechanisms at the heart of the interplay between fluid transport and phytoplankton growth dynamics. The favorable comparison of some of our results with those obtained in the framework of more realistic fluid models [START_REF] Taylor | Shutdown of turbulent convection as a new criterion for the onset of spring phytoplankton blooms[END_REF] seems to us interesting in light of the reduced computational cost of our approach. Several extensions can be envisaged, in a rather natural way. On one side, it would be interesting to consider a 3D extension of our kinematic flow [START_REF] A S Lanotte | Effects of vertical shear in modelling horizontal oceanic dispersion[END_REF][START_REF] Lacorata | 3D chaotic model for subgrid turbulent dispersion in large eddy simulations[END_REF], to explore possible links between the vertical organization of phytoplankton and its horizontal patchiness. On the other, we believe that accounting for vertical variations of the turbulent intensity, could provide a more realistic representation of real aquatic environments under stirring. Finally, introducing nutrient dynamics could allow investigating new and more realistic biological phenomena as, e.g., the response of the system to a sudden increase of the turbulent intensity [START_REF] Liccardo | Response of the deep chlorophyll maximum to fluctuations in vertical mixing intensity[END_REF], as the one generated by the arrival of a storm.

Chapter 6

Conclusions and Perspectives

Conclusions

This thesis presented an overview of some physical processes that affect phytoplankton life cycles. More specifically, I focused on vertical dynamics by addressing light limited growth conditions, advection and turbulent transport, factors that were shown to be relevant for the mentioned dynamics [START_REF] Huisman | Critical depth and critical turbulence: two different mechanisms for the development of phytoplankton blooms[END_REF][START_REF] Backhaus | Phyto-convection': the role of oceanic convection in primary production[END_REF][START_REF] Shigesada | Analysis of the self-shading effect on algal vertical distribution in natural waters[END_REF]. Chapter 2 presented an introduction to the concepts of turbulent transport [START_REF] Frisch | Turbulence: the legacy of A N Kolmogorov[END_REF][START_REF] Neufeld | Chemical and biological processes in fluid flows: a dynamical systems approach[END_REF], large scale advection [START_REF] Vallis | Atmospheric and oceanic fluid dynamics[END_REF] and some important effects relevant for oceans and lakes [START_REF] Thorpe | The turbulent ocean[END_REF]. Subsequently, in chapter 3 I introduced the basic features of phytoplankton dynamics [START_REF] Williams | Ocean dynamics and the carbon cycle: Principles and mechanisms[END_REF], its relation to turbulence [START_REF] Huisman | Population dynamics of sinking phytoplankton in light-limited environments: simulation techniques and critical parameters[END_REF][START_REF] Schmitt | Turbulence et écologie marine[END_REF], modelling frameworks [START_REF] Sandulescu | Biological activity in the wake of an island close to a coastal upwelling[END_REF][START_REF] Ebert | Critical conditions for phytoplankton blooms[END_REF][START_REF] Hernández | Sustained plankton blooms under open chaotic flows[END_REF] and I discussed different views on algal blooms [START_REF] Riley | Quantitative ecology of the plankton of the western North Atlantic[END_REF][START_REF] Huisman | Critical depth and critical turbulence: two different mechanisms for the development of phytoplankton blooms[END_REF]. The description of the numerical methodology in chapter 4 had the objective not only of clarifying the author's approach, but also of documenting how one can proceed to simulate advection-reaction-diffusion dynamics in a more general framework.That chapter also addresses the validation of the numerical algorithm adopted, by comparison with literature results [START_REF] Huisman | Population dynamics of sinking phytoplankton in light-limited environments: simulation techniques and critical parameters[END_REF]. Finally in chapter 5 I presented the thesis' main results relying on the content presented in the previous chapters.

The present work aimed to contribute to the understanding of vertical phytoplankton dynamics by including a further mechanism in previous models [START_REF] Huisman | Population dynamics of sinking phytoplankton in light-limited environments: simulation techniques and critical parameters[END_REF][START_REF] Shigesada | Analysis of the self-shading effect on algal vertical distribution in natural waters[END_REF][START_REF] Huisman | Critical depth and critical turbulence: two different mechanisms for the development of phytoplankton blooms[END_REF]. A discussion on the relevance of diagrams such as the ones displayed in Fig. 3.8 in chapter 3, arises. An attempt to extend the latter with the inclusion of advection was already made [START_REF] Lindemann | Dynamics of phytoplankton blooms in turbulent vortex cells[END_REF], although in a very simple flow setting (essentially the steady flow case discussed in chapter 5 of this thesis). The representation of those phase diagrams, obtained by Huisman and coauthors through multiple numerical simulations, was unique in the sense that it brought a general picture illustrating the influence of mixing over phytoplankton dynamics in an extremely idealized, but highly relevant system (see more in chapter 3). The interest to include a flow field, here was not aimed to extend the mentioned diagram, but rather in understanding how advection by coherent structures can be important to phytoplankton survival dynamics. In this sense the filament model [START_REF] Martin | On filament width in oceanic plankton distributions[END_REF] prediction revealed particularly useful to gain insight on the role of advection by large-scale flow features, as specified in the conclusions of chapter 5. In particular, it allowed an understanding of the phytoplankton survival/extinction transition in the presence of both mixing by turbulent diffusion and stirring by coherent structures.

Moreover, one can note that thermal/mechanical exchanges between the atmosphere and oceans/lakes can force convection in the upper layer of the latter, with limited areas of downward flow (downwelling) and of upward currents (upwelling).

At sufficiently small scales, we underline as well that persistent wind conditions can drive Langmuir circulation, as seen in chapter 2. Both processes, on different scales (convection < 1km, Langmuir circulation < 200m [START_REF] Vallis | Atmospheric and oceanic fluid dynamics[END_REF]), are examples of large-scale advection by vortex-like structures that can be accounted for in an idealized manner by the cellular flow configuration adopted in this thesis. Hence I expect that the work here displayed can be fruitful not only for my future research, in possible extensions of this project (explained in the next section), but also to other scientists interested in the fundamental understanding of phytoplankton vertical dynamics. It is worth stressing that the problem is still open. Even though there is some understanding of the importance of certain phenomena like convection [START_REF] Macias | Deep winter convection and phytoplankton dynamics in the NW Mediterranean Sea under present climate and future (horizon 2030) scenarios[END_REF][START_REF] Backhaus | Convection and primary production in winter[END_REF], the interplay between the multiple factors affecting phytoplankton dynamics (nutrients availability, predation, light availability, to cite a few) and the transport by complex fluid flows is still not completely understood. A major difficulty is here due to the broad range of relevant scales. This means that one has to make choices when studying these systems, it is impossible to take into account all the relevant scales and parameters of the problem. Therefore mathematical biological models will always contain gaps of knowledge, but they can be used to explore basic mechanisms. Furthermore, considering the difficulties related to obtaining field data in aquatic regions such models may also help to fill gaps in knowledge, by providing general expectations and new visions on the problem.

Perspectives

Here I introduce possible perspectives for the work presented in this thesis. They contain both studies that I would like to pursue in the future of my research and a work that is underway.

Regarding what I believe would be interesting for the continuation of this thesis, I start with the extension of the model studied in chapter 5 to a 3D system. First, one would like to understand how the phytoplankton patchiness on the horizontal relates to the one on the vertical given the growth dynamics and flow conditions specified on chapter 5. The idea is very appealing considering the availability of data about phytoplankton blooms at the ocean's surface, for comparison purposes. In practical terms the kinematic cellular flow model can be extended to three dimensions, without major conceptual difficulties, as it was done in previous works [START_REF] Lacorata | FSLE analysis and validation of Lagrangian simulations based on satellite-derived GlobCurrent velocity data[END_REF][START_REF] Lacorata | Lagrangian predictability characteristics of an Ocean Model[END_REF]. The main difficulties in our case are the interpolation and finite-difference scheme, as these can become computationally costly in 3D, possibly requiring a shift in the numerical methodology. One can cite as well the inclusion of nutrient dynamics, in order to relate blooms with areas of upwelling, which may be located inside large eddies or outside of them. Previous studies [START_REF] Sandulescu | Biological activity in the wake of an island close to a coastal upwelling[END_REF][START_REF] Pasquero | Impact of the spatiotemporal variability of the nutrient flux on primary productivity in the ocean[END_REF] have shown that under certain conditions the interplay between upwelling structures and biological growth leads to plankton blooms inside vortices. It would be interesting to verify if the same applies to light-limited growth dynamics with a more explicit representation of vertical transport.

Another continuation that was foreseen during this project was a comparison between the results of the present study, relying on a kinematic flow model, and those from a dynamical turbulent flow model. For the latter case, one possibility would be to couple our biological dynamics to a velocity flow field generated via 6.2. Perspectives 87 DNS, which provides an accurate description of the multiscale character of turbulent flows. However, this is unfeasible for realistic domain sizes (meaning for mixedlayer depths of O(10 -100)m). As an alternative, to cope with this limitation, one could resort to LES. This was already done by Taylor et al. [START_REF] Taylor | Shutdown of turbulent convection as a new criterion for the onset of spring phytoplankton blooms[END_REF], who showed that meteorological forcing plays a role in the onset of phytoplankton blooms as it drives convective fluid motions. Another compelling prospect is to move to a turbulent system that is more refined than our model but still idealized to some extent, namely to a 2D model of a convective mixed layer [START_REF] Bhamidipati | Turbulent mixing of a passive scalar in the ocean mixed layer[END_REF]. In the latter a buoyancy-driven flow is considered by solving the Boussinesq equations. While also in this case it is not possible to adopt realistic values of the mixed-layer depth, an interesting feature is the possibility to use boundary conditions such that the system reaches a dynamical equilibrium which could allow investigating statistical features in a reliable way. After some time, the flow attains a statistically steady state characterized by a vertically inhomogeneous turbulent intensity (higher in the mixed layer and lower below it). A preliminary, promising, attempt to reproduce such a velocity field was made by using the same Lattice-Boltzmann code [START_REF] Calzavarini | Eulerian-Lagrangian fluid dynamics platform: The ch4-project[END_REF] used for the solution of the Navier-Stokes equations in Fig. Finally, as explained in chapter 1, it was part of the planning of this project to participate in a campaign of data collection and analysis in the English Channel. Given the unfavorable conditions faced during the pandemic of 2020/2021 we could not proceed with this plan. Even though some measurements of salinity, PAR (radiation), temperature and chlorophyll were obtained, measurements of the vertical profiles of velocity and of the kinetic energy dissipation rate, could not be obtained. Consequently there were insufficient data to proceed with this study. Nevertheless, Chapter 6. Conclusions and Perspectives it is still in the perspectives of this project to compare the numerical results presented to real experimental data.

Phytoplankton Dynamics Under Ice

Multiple studies in the Arctic ocean reported under-ice phytoplankton blooms that were initiated by the onset of ice melt [START_REF] Assmy | Leads in Arctic pack ice enable early phytoplankton blooms below snow-covered sea ice[END_REF][START_REF] Boles | Under-ice phytoplankton blooms in the Central Arctic Ocean: Insights from the first biogeochemical IAOOS platform drift in 2017[END_REF]. This phenomenon increases light transmittance through the ice and induces a shoaling of the mixed layer. An example of phytoplankton communities near ice being shaped by turbulent eddies can be seen in Fig. 6.2, there one can also visualize openings in the ice that may benefit these organisms. It is not known whether under-ice blooms are controlled by the increase of light transmittance, leads (openings in the ice), upwelling from springtime convective mixing, biological factors, or other factors (most likely a combination of more than one). Nonetheless, it was observed that waters beneath sea ice with leads, had significantly lower phytoplankton biomass, despite high nutrient availability, than fully consolidated sea ice regions (without leads) [START_REF] Lowry | Under-ice phytoplankton blooms inhibited by spring convective mixing in refreezing leads[END_REF]. Phytoplankton blooms beneath snow-covered ice might become more common and widespread in the future Arctic Ocean, as lead formation becomes more frequent due to increasingly thinner and dynamic ice coverage [START_REF] Kwok | Arctic sea ice circulation and drift speed: Decadal trends and ocean currents[END_REF]. This could alter primary production and consequently the local marine food webs. With this in mind we decided to reformulate the model presented in chapters 4 and 5 to account for light obstruction by obstacles. The latter may be thought as idealized ice coverage, and can be altered both in size, location and in the intensity of light transmitted. An illustration of our first numerical setup can be seen in Fig 6 .3. We proceed to obtain results by taking advantage of our advection-reaction-diffusion model exposed in chapters 4 and 5. As a preliminary study we maintain the same For this study we planned to vary the obstacle position, size and the fraction of light intensity transmitted by the obstacle. Furthermore, we will increase the flow field complexity, maintaining the cellular flow configuration with distinct upwelling and downwelling regions, in order to see the effect of having upwelling under ice or under a lead/thinner ice.

biological parameters shown in table 4.1 and the same stationary cellular flow configuration through the use of the streamfunction of Eq. 5.1. Consequently the only change in the system compared to the previous chapter is the presence of an obstacle of size L x /2, centered at x = L x /2, where L x is the system's horizontal size. The obstacle does not transmit light and hence the entire region underneath it is biologically ruled by a null specific production rate p(I) = 0 since I in = 0. The production underneath the lead is still given by Eq. 4.25, so dynamics are controlled by: ∂θ(x, z, t) ∂t = -lθ(x, z, t)v • ∇θ(x, z, t) + D∇ 2 θ(x, z, t), (6.1)

under the obstacle, and ∂θ(x, z, t) ∂t = [p(I)l] θ(x, z, t)v • ∇θ(x, z, t) + D∇ 2 θ(x, z, t), (6.2)

under the unobstructed space.

Furthermore, the obstacle is located at the surface z = 0 and does not have any physical influence over the scalar field or the flow field other than the obstruction of light. The velocity term v still accounts for phytoplankton sinking and the advective flow field. Boundary conditions are maintained as before, no flux on the vertical and periodic on the horizontal. Given that, even in the absence of the flow field, bloom conditions are not controlled anymore by the diagram of Fig. 3.8, we decided to consider shallower systems. Our interest is still in the bloom/no bloom transition, but since production is now globally lower as a result of the blockage of light, the transition can occur also in shallower systems. Thus one can save computational resources while still simulating a depth close to the critical depth of the system. So we begin simulations in a rectangular domain of sides L x = 60m and L z = 30m. We proceed to make qualitative analyses (some quantitative measurements will be also displayed) on this preliminary results; these are shown in the regime of stationary per-capita growth rate r p (t) = r p , according to Eq. 5.3 (as in chapter 5). The first results are presented in Fig. 6.4 in the form of density field snapshots. There, the influence of the obstacle (represented by the black rectangle) is tested for a no-flow case and two different flow speeds U. In addition, two different values of D were considered (corresponding to the different columns in Fig. 6.4). In the left column we show the D = 1cm 2 /s set of results. There, the presence of the obstacle is made clear by the accumulation of phytoplankton in the unobstructed zones in the absence of advection. Also, for this diffusivity the plankton patches are localized closer to the surface, in contrast to the D = 10cm 2 /s no-flow case (right column). In the higher D no-flow case, the phytoplankton spreads more evenly through the water column, while horizontally it is still mostly restricted to the lit zones. In the D = 1cm 2 /s cases, the flow appears to have a strong influence in the mixing, allowing the phytoplankton to populate more evenly the entire simulation domain as U increases. In the D = 10cm 2 /s cases the increase in flow speed has a much weaker influence. We notice that curiously the D = 1cm 2 /s, U = 0 case has a very similar configuration to the D = 10cm 2 /s, U = 0.7m/h case. Both are bloom cases, although in the first population growth is much higher (see Fig. 6.5). Another remarkable feature seen in Fig. 6.4 is that the phytoplankton patches in the smaller diffusivity cases break into more irregular shapes, while in the higher D cases they maintain more regular configurations. Also in this case of light intensity modulation at the surface, the phytoplankton survival hindering effect of the flow field remains as U increases, in parallel with the results of chapter 5. Here the flow proves particularly harmful for survival in the D = 1cm 2 /s cases; when it is present, net growth rates are always negative, for the values tested here (Fig. 6.5). This is in contrast with our findings in chapter 5, where for smaller D the transition bloom/no bloom was never possible. Even though in our previous study the system was deeper (L z = 60m), it is only with the presence of the obstacle that we see the transition. This seems to indicate that the flow field killing action is more effective in the present setting. This finding is in agreement with a study [START_REF] Lowry | Under-ice phytoplankton blooms inhibited by spring convective mixing in refreezing leads[END_REF] that shows that waters beneath loosely consolidated sea ice (with leads) had weak stratification and were frequently mixed below the critical depth, resulting in these locations having lower phytoplankton biomass than regions covered by fully consolidated ice.

We conclude by showing more density field snapshots obtained in simulations with advection by the multiscale flow field. We remind that in this case the flow corresponds to the streamfunction of equation 5.10. In Fig. 6.6 the signature of the inclusion of smaller vortices is identifiable in the irregular shapes taken by the phytoplankton patches. This is specially true for the D = 1cm 2 /s case where the advection action is more dominant, but in both cases we now have horizontal asymmetry. Nonetheless, the phytoplankton patches are still localized in the unobstructed regions, as in the stationary flow case. What is more interesting in the multiscale flow cases is that the D = 1cm 2 /s case now gives rise to a bloom event (see Fig. 6.7). To understand this, we can recall one of the relevant features identified in chapter 5, namely that the multiscale flow produces an increase in diffusivity, that counteracts the hindering effect of advection in phytoplankton growth. This result illustrates the beneficial effect of the small turbulent scales on the survival of phytoplankton in the presence of obstacles to light transmission. The multiscale D = 10cm 2 /s case however, presents almost no difference with respect to the corresponding stationary flow case, again a result that agrees with the study of the previous chapter.

The preliminary results here exposed show a promising continuation for the work presented in this thesis. The immediate next steps will deal with the study of different obstacle conditions, such as: different localization and size, as well as partial transmittance of light. We believe that this work could prove complementary to some of the studies cited in this section [START_REF] Lowry | Under-ice phytoplankton blooms inhibited by spring convective mixing in refreezing leads[END_REF][START_REF] Boles | Under-ice phytoplankton blooms in the Central Arctic Ocean: Insights from the first biogeochemical IAOOS platform drift in 2017[END_REF], providing further evidence of the importance of fluid motion on phytoplankton survival under ice.

  et Technologies Doctor of Philosophy in Fluid Mechanics Phytoplankton Vertical Dynamics in the Presence of Advection by Vinicius BELTRAM TERGOLINA

1 . 1 2 2. 1

 1121 Lagoa dos Patos, Rio Grande do Sul, Brazil. Suspended sediment, phytoplankton, submerged aquatic vegetation, and colored dissolved organic matter all contribute to the many colors observed, https://oceancolor.gsfc.nasa.gov/. 1 1.2 Measuring equipment being deployed in the English Channel from the CNRS research vessel Sepia (based in Boulogne-sur-Mer) during a campaign to which the author participated in 2020. The objective of the campaign was to sample temperature, salinity, chlorophyll, PAR (radiation), ocean kinetic energy dissipation rate and velocity in multiple locations. Unfortunately due to the restricted number of sampling days under the pandemic of 2020/2021 the data obtained was insufficient for analysis. . . . . . . . . . . . . . . . . . . . . . . . . . . . Studies of turbulent water shows Leonardo's understanding of how vortices move [9]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 A drunk man returning home is an effective metaphor of a biased random walk, although his preferred direction is the one pointing to his house his inebriated state is likely to throw him off the path in an unpredictable manner. https://prakhartechviz.blogspot.com/2019/09/randomwalk-term-weighting-for-text.html. . . . . . . . . . . . . . . . . . . . . . . 7 2.3 A 1D lattice with 2n steps of size 1. In an isotropic random walk the probability p(τ) that a particle will go to the left or the right, in a time step τ, is 50%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.4 Example of discrete vs. continuous solution of Fick's law for one dimension, here D = 0.01m 2 h -1 . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.5 A schematic representation of laminar flow (a) and a turbulent flow (b). 9 2.6 A schematic representation of the turbulence cascade proposed by Richardson. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.7 Schema of energy spectrum in three-dimensional turbulence, according to Kolmogorov's theory. The inertial range is situated between the large scale L 0 and the Kolmogorov length scale η. Energy is injected at rates i and transferred at rates from the injection range to the dissipation range (where it is dissipated at rates ν ). According to Kolmogorov i = = ν , this leads to the 'Kolmogorov's -5/3 spectrum'. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.8 Solution of the Navier-Stokes equation (in non-dimensional units) via

3 . 4

 34 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvii Integration of equation 3.6 in non dimensional units, K = 100000, r = 3 • 10 -4 and dt = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Growth term vs. depth from Eq. 3.14, without the self-shading effect. The red line indicates the compensation depth (when net growth becomes negative, i.e. when l > aI/(1 + cI)) at z = 17.8m for K bg = 0.2m -1 , a coefficient of absorption equivalent to clear lakes or coastal areas. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Visualization of the concept of critical depth. Here dr is the loss rate, dp is the production rate and D c the compensation depth. If the area bcD c is larger than the area abD c , production surpasses losses. Extracted from [58]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Combinations of water-column depth and turbulent diffusion coef-

4. 1

 1 Illustration of the discretization process of a continuous 2D density field θ(x, z, t) into a grid of size N x × N z , indexes i, j and cell size h. . . 4.2 Illustration of the backwards and forward schemes for diffusion and advection, respectively corresponding to moving fictitious fluid particle in the past and the future. Here a single noise realization (needed for diffusion) is represented. . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Example of the comparison between population density field snapshots calculated for advection-reaction-diffusion dynamics corresponding to the model of Sec. 4.4, by the forward semi-Lagrangian method (a) and the backwards semi-Lagrangian method with diffusion by the Eulerian scheme (b). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Flowchart portraying the order of implementation of the numerical methodology explained in the past sections. The final method chosen for this thesis is on the right side. . . . . . . . . . . . . . . . . . . . . . . 4.5 The specific growth rate (p(I)l) for a 20 meter deep water column and multiple values of θ. The surface is obtained with the simulation parameters of table 4.1 and equation 4.20. . . . . . . . . . . . . . . . . . 4.6 Total population versus time for dz = 0.133m, z max = 50m, D = 10cm 2 /s and multiple time steps dt. . . . . . . . . . . . . . . . . . . . . 4.7 P tot vs. time (log scale on the vertical axis) for z max = 20m and multiple D values. The inset shows a zoom of the top area. . . . . . . . . . . 4.8 Reproduction of the bloom diagram in Fig. 3.8 (top), from [45]. This was done by running multiple simulations for different D and z max values, and measuring the value of P tot after a sufficiently large time. .

xviii 4 . 9

 49 Population density vertical profiles for different values of D and z max = 20m. Results are presented in the same units as in the original [45]: a) D = 10cm 2 s -1 , b) D = 1cm 2 s -1 , c) D = 0.5cm 2 s -1 and d) D = 0.1cm 2 s -1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Comparison of population density vertical profiles for the 1D approach and the 2D approach. a) D = 0.1cm 2 s -1 , b) D = 10cm 2 s -1 ; z max = 20m in both cases. One can observe that curves in the 1D case are more noisy due to the stochastic way of calculating diffusion in the methodology utilized for the 1D simulations. . . . . . . . . . . . . . . 5.1 Streamlines of the steady large-scale flow. The color bar |u(x, z)|/U represents the flow's velocity absolute value normalized by its maximum value. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 P tot dynamics for D = 10cm 2 s -1 , L z = 60 m and multiple values of U (log-scale on the vertical axis). . . . . . . . . . . . . . . . . . . . . . . . 5.3 Instantaneous normalized population density fields θ(x, z, t * )/ θ at a fixed instant of time t * = 640h, where θ stands for the spatial average. The white lines are the isolines θ/ θ = 1, black arrows represent the stream function of Eq. 5.1 and t * is in the regime of stationary per-capita growth rate (r p (t) = r p ). In a), b) and c) U = (0, 1.24, 4.93) m h -1 (respectively) and D = 5 cm 2 s -1 . In d), e) and f) D = (10, 15, 20) cm 2 s -1 (respectively) and U = 1.24 m h -1 . . . . . . . 5.4 Average biomass density θ (t), on a logarithmic scale, versus time, for D = 20 cm 2 /s, L z = 60 m and different values of the advection intensity U in the steady-flow case. Vertical lines indicate t = 2πL z /U, the time of one flow roll revolution. . . . . . . . . . . . . . . . . . . . . 5.5 Population per-capita growth rate r p (t), normalized by the intrinsic net growth rate r b , as a function of time, normalized by the advective time scale L z /U of the steady-flow case, for various values of U and L z = 60 m. a) D = 20 cm 2 s -1 and b) D = 1 cm 2 s -1 . . . . . . . . . . . 5.6 Asymptotic per-capita growth rate r p , normalized by the intrinsic net growth rate r b , versus the ratio of reactive to advective time scales γ, in the steady-flow case, for different values of D. . . . . . . . . . . . . 5.7 Measurements of the transition flow intensity U transition (from bloom to no-bloom) for different D obtained from the simulations of Fig. 5.6. 5.8 Instantaneous normalized population density horizontal profiles θ(x, z * , t * )/ θ(x, z * , t * ) at a fixed instant of time t * = 640h. In (a), (b) (c) and (d) D = (5, 10, 15, 20) cm 2 s -1 (respectively) and U = 1.24 m h -1 . . . . . . . . . 5.9 (a) Horizontal profiles of surface population density at different times (in non-dimensional units), normalized by their average values, θ(x, z = 0, t)/ θ(x, z = 0, t), with U = 1.001 m/h and D = 20 cm 2 s -2 . (b) Gaussian fitting examples on the data of (a) (here we use θ = θθ back e r b t ) for different interval slices of the profile. . . . . . . . . . . 5.10 Filament width, estimated from a fit (in the interval 50 m < x < 70 m) with a Gaussian function of horizontal profiles of population density from simulations with different values of D and U, versus its theoretical prediction σ = D/(kU), with k = π/L z . The dashed black line corresponds to σ numerical = 0.085 + 1.276σ. . . . . . . . . . . . . . . . . xix 5.11 Instantaneous population density field, normalized by its spatial average, in the unsteady-flow case, θ(x, z, t)/ θ for U = 1.4m/h, D = 5 cm 2 s -1 , at two different times in panels (a) and (b). The white line is the isoline θ/ θ = 1 and the black lines are the streamlines of the flow field at the considered instants of time. . . . . . . . . . . . . . . . . 5.12 Population per-capita growth rate r p (t), normalized by the intrinsic net growth rate r b , as a function of time, normalized by the advective time scale L z /U of the steady-flow case, for various values of U, D = 20 cm 2 s -1 and L z = 60 m. . . . . . . . . . . . . . . . . . . . . . . . . . 5.13 Time averaged normalized per-capita growth rate rp /r b of the steadyflow case vs. the corresponding quantity from the unsteady-flow case. In here we see several values of U and D are considered. Fitting the data corresponding to a given value of D [(5, 10, 15, 20) cm 2 s -1 ] with a linear function, we obtain slopes that are always quite close to 1 (slopes between 0.84 and 0.93). The black dashed line has unitary slope. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.14 Instantaneous population density field, normalized by its spatial average, in the multiscale flow case, θ(x, z, t)/ θ for U = 1.001m/h, D = 0.1 cm 2 s -1 , at t * = 845h. The initial condition was set as a uniform plankton patch in the central top area of the simulation box. The black lines are the streamlines of the flow field at the considered instant of time according to Eq. 5.10. . . . . . . . . . . . . . . . . . . . . 5.15 Instantaneous population density field, normalized by its spatial average, θ(x, z, t)/ θ , in the multiscale-flow case, for U = 1.021 m/h and D = 5 cm 2 s -1 , in the asymptotic regime r p (t) = const. The white line is the isoline θ/ θ = 1 and the black lines are the streamlines of the flow field at the considered instant of time. . . . . . . . . . . . . . . 5.16 (a) Population per-capita growth rate r p (t), normalized by the intrinsic net growth rate r b , as a function of time, normalized by the advective time scale L z /U of the multiscale flow case, for various values of U and D=20 cm 2 s -1 . (b) Total biomass curve vs. time for U = 2.466m/h and D = 20cm 2 s -1 (in black), in red we see a curve corresponding to the expression P tot (t = 0)exp(r p t) where P tot (t = 0) corresponds to the first point of the black curve and rp to the time averaged per capita growth rate extracted from the black curve. . . . . 5.17 Time averaged normalized per-capita growth rate rp /r b of the multiscaleflow case vs. the corresponding quantity from the unsteady-flow case. In here several values of U and D are considered. Fitting the data corresponding to a given value of D [(5, 10, 15, 20) cm 2 s -1 ] with a linear function we obtain slopes between 0.86 and 1.21. The black dashed line has unitary slope. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.18 Vertical population density profiles θ x , normalized by the global spatial average θ , for the different stream functions considered, for D = 5 and 20 cm 2 s -1 . Panels a-d) correspond to different instants of time, in the asymptotic growth-rate regime r p (t) = const, as specified in the plot titles (in units of the advective time scale L z /U 1 ). . . . . . . xx 5.19 Normalized vertical population density profiles θ x / θ for Ψ = 0 (no flow), Ψ = Ψ st L (large-scale steady flow), Ψ = Ψ L (large-scale timedependent flow), Ψ = Ψ L + Ψ s (multiscale time-dependent flow), Ψ = Ψ s (small-scale time-dependent flow, without Ψ L ). Panels (a) and (b) respectively refer to D = 5 cm 2 s -1 and D = 20 cm 2 s -1 . Note the different value ranges on the horizontal axes in (a) and (b). All the profiles here reported are computed at a common fixed time tU 1 L -1 z = 10, for which r p (t) = const. . . . . . . . . . . . . . . . . . . . 6.1 Temperature (color) and velocity (arrows) fields at a common instant of time from preliminary simulations of mixed layer model [112] obtained using the Lattice-Boltzmann method [20]. Here one can see the presence of convective plumes in the turbulent flow field. . . . . . . . . 6.2 Satellite image from Suomi-NPP/VIIRS collected on June 14, 2015 (https://oceancolor.gsfc.nasa.gov/gallery/454/). Phytoplankton communities and sea ice limn are shaped by the turbulent flow field between Greenland and Iceland. Openings in the ice (leads) near the coast, and their influence over vertical phytoplankton dynamics, are a motivation for this study. . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Illustration of the phytoplankton under-ice model setup. For this study we planned to vary the obstacle position, size and the fraction of light intensity transmitted by the obstacle. Furthermore, we will increase the flow field complexity, maintaining the cellular flow configuration with distinct upwelling and downwelling regions, in order to see the effect of having upwelling under ice or under a lead/thinner ice. . . . 6.4 Instantaneous normalized population density fields θ(x, z, t * )/ θ at a fixed instant of time t * = 600h, where θ stands for the spatial average. The obstacle is represented by the black rectangle and the dashed lines delimit the light obstruction zone. The white lines are the isolines θ/ θ = 1. The solid black lines represent flow streamlines from Eq. 5.1, with arrows indicating the circulation direction. Also, t * is in the regime of stationary per-capita growth rate (r p (t) = r p ). In a), b) and c) U = (0, 0.71, 1.19) m h -1 (respectively) and D = 1 cm 2 s -1 . In d), e) and f) U = (0, 0.71, 1.19) m h -1 (respectively) and D = 10 cm 2 s -1 . a), d) and e) are bloom cases while b), c) and f) are no bloom cases. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.5 Population per-capita growth rate r p (t), normalized by the intrinsic net growth rate r b , as a function of time in the steady-flow case, for various values of U and D, with a light obstruction of size L x /2 and zero transmitted light, centered at L x /2, localized at z = 0. The inset shows the same graph with the addition of the no-flow cases. . . . . . 6.6 Instantaneous normalized population density fields θ(x, z, t * )/ θ at a fixed instant of time t * = 600h, where θ stands for the spatial average. The obstacle is represented by the black rectangle and the dashed lines delimit the light obstruction zone. The white lines are the isolines θ/ θ = 1. The solid black lines represent flow streamlines from Eq. 5.10, with arrows indicating the circulation direction. Also, t * is in the regime of stationary per-capita growth rate (r p (t) = r p ). In a) and b) U = 0.71 m h -1 and D = (10, 1) cm

FIGURE 1 . 1 :

 11 FIGURE 1.1: Lagoa dos Patos, Rio Grande do Sul, Brazil. Suspended sediment, phytoplankton, submerged aquatic vegetation, and colored dissolved organic matter all contribute to the many colors observed, https://oceancolor.gsfc.nasa.gov/.

FIGURE 1 . 2 :

 12 FIGURE 1.2: Measuring equipment being deployed in the English Channel from the CNRS research vessel Sepia (based in Boulognesur-Mer) during a campaign to which the author participated in 2020. The objective of the campaign was to sample temperature, salinity, chlorophyll, PAR (radiation), ocean kinetic energy dissipation rate and velocity in multiple locations. Unfortunately due to the restricted number of sampling days under the pandemic of 2020/2021 the data obtained was insufficient for analysis.

FIGURE 2 . 1 :

 21 FIGURE 2.1: Studies of turbulent water shows Leonardo's understanding of how vortices move [9].

Chapter 2 .

 2 Turbulence and Turbulent Transport

FIGURE 2 . 4 :

 24 FIGURE 2.4: Example of discrete vs. continuous solution of Fick's law for one dimension, here D = 0.01m 2 h -1 .

FIGURE 2 . 5 :

 25 FIGURE 2.5: A schematic representation of laminar flow (a) and a turbulent flow (b).

FIGURE 2 . 6 :

 26 FIGURE 2.6: A schematic representation of the turbulence cascade proposed by Richardson.

FIGURE 2 . 7 :

 27 FIGURE 2.7: Schema of energy spectrum in three-dimensional turbulence, according to Kolmogorov's theory. The inertial range is situated between the large scale L 0 and the Kolmogorov length scale η.Energy is injected at rates i and transferred at rates from the injection range to the dissipation range (where it is dissipated at rates

FIGURE 2 . 8 :

 28 FIGURE 2.8: Solution of the Navier-Stokes equation (in nondimensional units) via Lattice-Boltzmann method[START_REF] Calzavarini | Eulerian-Lagrangian fluid dynamics platform: The ch4-project[END_REF] with a forcing of the form of Eq. 2.23 for a grid size of 600X300 cells. Here the color code represent u intensity and the black arrows the the direction of the velocity field.

FIGURE 2 .

 2 FIGURE 2.10: A map of the apparent horizontal diffusivity (K) in the Atlantic obtained through measurements of eddy kinetic energy by drifters. Figure extracted from [33].

FIGURE 2 . 11 :

 211 FIGURE 2.11: Pictorial representation of ocean flows at various scales, increased mixing is represented by stippling. Figure extracted from [19].

FIGURE 2 . 12 :

 212 FIGURE 2.12: Pictorial representation of open ocean deep convection, Buoyancy flux through sea surface is represented by curly arrows, stratification is represented by continuous lines. (a) preconditioning (buoyancy loss at the surface), (b) deep convection (currents are driven by density gradients), and (c) lateral exchange and spreading (the phenomenon slowly dissipates). Figure extracted from [34].

FIGURE 2 .

 2 FIGURE 2.13: Photo of the surface of the Garonne river taken by the author during a stay in Bordeaux, France. A constant flux of wind over the water surface produced a series of vortical motions with a life span averaging from a couple of seconds to a minute.

FIGURE 2 . 14 :

 214 FIGURE 2.14: Langmuir streaks on Quake Lake, United States. Image by Wayne Wurtsbaugh -www.aslo.org.

FIGURE 3 . 1 :

 31 FIGURE 3.1: Examples of planktonic microorganisms and macro organisms. (I) On the left two dinoflagellates (photosynthetic primary producers) and a tintinnid ciliate (phytoplankton predator) on the right, both are active swimmers, Thau Lagoon of Sète, France -Uwe Kils. (II) A copepod from Antarctica, planktonic crustaceans that can be 1mm to 2mm long, Uwe Kils. (III) Antarctic krill, feeds primarily on phytoplankton and and it is probably the largest biomass of a single species on the planet, Professor Dr. Habil, Uwe Kils. (IV) Amphipod a scavenger species, Uwe kils. (V) A cyanobacterium Prochlorococcus, a photosynthetic organism and major contributor to atmospheric oxygen, Luke Thompson from Chisholm Lab and Nikki Watson from Whitehead, MIT. (VI) The sea sparkle dinoflagellate, an active swimmer species that exhibits bioluminescence when disturbed producing the mareel (Milky seas) effect, Maria Antónia Sampayo, Instituto de Oceanografia, Faculdade Ciências da Universidade de Lisboa -http://planktonnet.awi.de.

FIGURE 3 . 2 :

 32 FIGURE 3.2: (I) A few hundred thousands square kilometers of phytoplankton patches in the Southeastern Pacific Ocean, filaments appear to follow the flow pattern of mesoscale ocean vortices. (II) Algae patches drawn into thin ribbons by a turbulent eddy field in the northern area of the Arabic Sea. (III) An image with seven separate spectral bands to highlight differences in plankton communities on the coast of Patagonia. (IV) The southeast part of Rio de la Plata shows large floating populations of phytoplankton favored by the proximity to the metropolis of Buenos Aires and Montevideo, huge suppliers of nutrients such as nitrate and phosphate. https://oceancolor.gsfc.nasa.gov/.

FIGURE 3 . 3 :

 33 FIGURE 3.3: Representation of the elements responsible for phytoplankton survival dynamics in aquatic environments. The euphotic zone is the sunlit layer of the ocean where photosynthesis is viable, the mixed layer is a region with active movement that mixes waters in an almost homogeneous way. In an ideal scenario for reproduction the mixed layer allows phytoplankton to have access to both the euphotic zone and nutrients.

FIGURE 3 . 4 :

 34 FIGURE 3.4: Integration of equation 3.6 in non dimensional units, K = 100000, r = 3 • 10 -4 and dt = 1.

FIGURE 3 . 5 :

 35 FIGURE 3.5: Growth term vs. depth from Eq. 3.14, without the selfshading effect. The red line indicates the compensation depth (when net growth becomes negative, i.e. when l > aI/(1 + cI)) at z = 17.8m for K bg = 0.2m -1 , a coefficient of absorption equivalent to clear lakes or coastal areas.

FIGURE 3 . 6 :

 36 FIGURE 3.6: Visualization of the concept of critical depth. Here dr is the loss rate, dp is the production rate and D c the compensation depth. If the area bcD c is larger than the area abD c , production surpasses losses. Extracted from [58].

FIGURE 3 . 7 :

 37 FIGURE 3.7: Combinations of water-column depth and turbulent diffusion coefficient that allow a phytoplankton bloom and the combinations that prevent a phytoplankton bloom. Extracted from [7].

FIGURE 3 . 8 :

 38 FIGURE 3.8: Results analog to Fig. 2 but with the inclusion of sinking speed in the simulations. Above v = 0.04m/h and below v = 0.4m/h. Above we see the concepts of maximal turbulence and minimal turbulence represented by the horizontal lines. The critical depth and the compensation depth are equally represented by vertical lines. Extracted from [45].

FIGURE 4 . 1 :

 41 FIGURE 4.1: Illustration of the discretization process of a continuous 2D density field θ(x, z, t) into a grid of size N x × N z , indexes i, j and cell size h.

FIGURE 4 . 2 :

 42 FIGURE 4.2: Illustration of the backwards and forward schemes for diffusion and advection, respectively corresponding to moving fictitious fluid particle in the past and the future. Here a single noise realization (needed for diffusion) is represented.

FIGURE 4 . 3 :

 43 FIGURE 4.3: Example of the comparison between population density field snapshots calculated for advection-reaction-diffusion dynamics corresponding to the model of Sec. 4.4, by the forward semi-Lagrangian method (a) and the backwards semi-Lagrangian method with diffusion by the Eulerian scheme (b).

FIGURE 4 . 4 :

 44 FIGURE 4.4: Flowchart portraying the order of implementation of the numerical methodology explained in the past sections. The final method chosen for this thesis is on the right side.

FIGURE 4 . 5 :

 45 FIGURE 4.5: The specific growth rate (p(I)l) for a 20 meter deep water column and multiple values of θ. The surface is obtained with the simulation parameters of table 4.1 and equation 4.20.

2 )FIGURE 4 . 7 :

 247 FIGURE 4.7: P tot vs. time (log scale on the vertical axis) for z max = 20m and multiple D values. The inset shows a zoom of the top area.

FIGURE 4 . 8 :

 48 FIGURE 4.8: Reproduction of the bloom diagram in Fig. 3.8 (top), from [45]. This was done by running multiple simulations for different D and z max values, and measuring the value of P tot after a sufficiently large time.

FIGURE 4 . 9 :

 49 FIGURE 4.9: Population density vertical profiles for different values of D and z max = 20m. Results are presented in the same units as in the original [45]: a) D = 10cm 2 s -1 , b) D = 1cm 2 s -1 , c) D = 0.5cm 2 s -1 and d) D = 0.1cm 2 s -1 .

FIGURE 4 . 10 :

 410 FIGURE 4.10: Comparison of population density vertical profiles for the 1D approach and the 2D approach. a) D = 0.1cm 2 s -1 , b) D = 10cm 2 s -1 ; z max = 20m in both cases. One can observe that curves in the 1D case are more noisy due to the stochastic way of calculating diffusion in the methodology utilized for the 1D simulations.

FIGURE 5 . 1 :

 51 FIGURE 5.1: Streamlines of the steady large-scale flow. The color bar |u(x, z)|/U represents the flow's velocity absolute value normalized by its maximum value.

. 4 )FIGURE 5 . 2 :

 452 FIGURE 5.2: P tot dynamics for D = 10cm 2 s -1 , L z = 60 m and multiple values of U (log-scale on the vertical axis).

FIGURE 5 . 3 :

 53 FIGURE 5.3: Instantaneous normalized population density fields θ(x, z, t * )/ θ at a fixed instant of time t * = 640h, where θ stands for the spatial average. The white lines are the isolines θ/ θ = 1, black arrows represent the stream function of Eq. 5.1 and t * is in the regime of stationary per-capita growth rate (r p (t) = r p ). In a), b) and c) U = (0, 1.24, 4.93) m h -1 (respectively) and D = 5 cm 2 s -1 . In d), e) and f) D = (10, 15, 20) cm 2 s -1 (respectively) and U = 1.24 m h -1 .

FIGURE 5 . 4 :

 54 FIGURE 5.4: Average biomass density θ (t), on a logarithmic scale, versus time, for D = 20 cm 2 /s, L z = 60 m and different values of the advection intensity U in the steady-flow case. Vertical lines indicate t = 2πL z /U, the time of one flow roll revolution.

FIGURE 5 . 5 :

 55 FIGURE 5.5: Population per-capita growth rate r p (t), normalized by the intrinsic net growth rate r b , as a function of time, normalized by the advective time scale L z /U of the steady-flow case, for various values of U and L z = 60 m. a) D = 20 cm 2 s -1 and b) D = 1 cm 2 s -1 .

FIGURE 5 . 6 :

 56 FIGURE 5.6: Asymptotic per-capita growth rate r p , normalized by the intrinsic net growth rate r b , versus the ratio of reactive to advective time scales γ, in the steady-flow case, for different values of D.

FIGURE 5 . 7 :

 57 FIGURE 5.7: Measurements of the transition flow intensity U transition (from bloom to no-bloom) for different D obtained from the simulations of Fig. 5.6.

FIGURE 5 . 8 :

 58 FIGURE 5.8: Instantaneous normalized population density horizontal profiles θ(x, z * , t * )/ θ(x, z * , t * ) at a fixed instant of time t * = 640h. In (a), (b) (c) and (d) D = (5, 10, 15, 20) cm 2 s -1 (respectively) and U = 1.24 m h -1 .

FIGURE 5 . 9 :

 59 FIGURE 5.9: (a) Horizontal profiles of surface population density at different times (in non-dimensional units), normalized by their average values, θ(x, z = 0, t)/ θ(x, z = 0, t), with U = 1.001 m/h and D = 20 cm 2 s -2 . (b) Gaussian fitting examples on the data of (a) (here we use θ = θθ back e r b t ) for different interval slices of the profile.

FIGURE 5 . 10 :

 510 FIGURE 5.10: Filament width, estimated from a fit (in the interval 50 m < x < 70 m) with a Gaussian function of horizontal profiles of population density from simulations with different values of D and U, versus its theoretical prediction σ = D/(kU), with k = π/L z . The dashed black line corresponds to σ numerical = 0.085 + 1.276σ.

FIGURE 5 . 11 :

 511 FIGURE 5.11: Instantaneous population density field, normalized by its spatial average, in the unsteady-flow case, θ(x, z, t)/ θ for U = 1.4m/h, D = 5 cm 2 s -1 , at two different times in panels (a) and (b). The white line is the isoline θ/ θ = 1 and the black lines are the streamlines of the flow field at the considered instants of time.

FIGURE 5 . 12 :

 512 FIGURE 5.12: Population per-capita growth rate r p (t), normalized by the intrinsic net growth rate r b , as a function of time, normalized by the advective time scale L z /U of the steady-flow case, for various values of U, D = 20 cm 2 s -1 and L z = 60 m.

FIGURE 5 . 13 :

 513 FIGURE 5.13: Time averaged normalized per-capita growth rate rp /r b of the steady-flow case vs. the corresponding quantity from the unsteady-flow case. In here we see several values of U and D are considered. Fitting the data corresponding to a given value of D [(5, 10, 15, 20) cm 2 s -1 ] with a linear function, we obtain slopes that are always quite close to 1 (slopes between 0.84 and 0.93). The black dashed line has unitary slope.

FIGURE 5 . 14 :

 514 FIGURE 5.14: Instantaneous population density field, normalized by its spatial average, in the multiscale flow case, θ(x, z, t)/ θ for U = 1.001m/h, D = 0.1 cm 2 s -1 , at t * = 845h. The initial condition was set as a uniform plankton patch in the central top area of the simulation box. The black lines are the streamlines of the flow field at the considered instant of time according to Eq. 5.10.

FIGURE 5 . 16 :

 516 FIGURE 5.16: (a) Population per-capita growth rate r p (t), normalized by the intrinsic net growth rate r b , as a function of time, normalized by the advective time scale L z /U of the multiscale flow case, for various values of U and D=20 cm 2 s -1 . (b) Total biomass curve vs. time for U = 2.466m/h and D = 20cm 2 s -1 (in black), in red we see a curve corresponding to the expression P tot (t = 0)exp(r p t) where P tot (t = 0) corresponds to the first point of the black curve and rp to the time averaged per capita growth rate extracted from the black curve.

FIGURE 5 . 17 :

 517 FIGURE 5.17: Time averaged normalized per-capita growth rate rp /r b of the multiscale-flow case vs. the corresponding quantity from the unsteady-flow case. In here several values of U and D are considered. Fitting the data corresponding to a given value of D [(5, 10, 15, 20) cm 2 s -1 ] with a linear function we obtain slopes between 0.86 and 1.21. The black dashed line has unitary slope.

FIGURE 5 . 18 :

 518 FIGURE 5.18: Vertical population density profiles θ x , normalized by the global spatial average θ , for the different stream functions considered, for D = 5 and 20 cm 2 s -1 . Panels a-d) correspond to different instants of time, in the asymptotic growth-rate regime r p (t) = const, as specified in the plot titles (in units of the advective time scale L z /U 1 ).

FIGURE 5 . 19 :

 519 FIGURE 5.19: Normalized vertical population density profiles θ x / θ for Ψ = 0 (no flow), Ψ = Ψ st L (large-scale steady flow), Ψ = Ψ L (large-scale time-dependent flow), Ψ = Ψ L + Ψ s (multiscale time-dependent flow), Ψ = Ψ s (small-scale time-dependent flow, without Ψ L ). Panels (a) and (b) respectively refer to D = 5 cm 2 s -1 and D = 20 cm 2 s -1 . Note the different value ranges on the horizontal axes in (a) and (b). All the profiles here reported are computed at a common fixed time tU 1 L -1 z = 10, for which r p (t) = const.

  2.8. Simultaneous snapshots of temperature (color coded) and velocity (arrows) fields are shown in Fig. 6.1.

FIGURE 6 . 1 :

 61 FIGURE 6.1: Temperature (color) and velocity (arrows) fields at a common instant of time from preliminary simulations of mixed layer model[START_REF] Bhamidipati | Turbulent mixing of a passive scalar in the ocean mixed layer[END_REF] obtained using the Lattice-Boltzmann method[START_REF] Calzavarini | Eulerian-Lagrangian fluid dynamics platform: The ch4-project[END_REF]. Here one can see the presence of convective plumes in the turbulent flow field.

FIGURE 6 . 2 :

 62 FIGURE 6.2: Satellite image from Suomi-NPP/VIIRS collected on June 14, 2015 (https://oceancolor.gsfc.nasa.gov/gallery/454/). Phytoplankton communities and sea ice limn are shaped by the turbulent flow field between Greenland and Iceland. Openings in the ice (leads) near the coast, and their influence over vertical phytoplankton dynamics, are a motivation for this study.

FIGURE 6 . 3 :

 63 FIGURE 6.3: Illustration of the phytoplankton under-ice model setup.For this study we planned to vary the obstacle position, size and the fraction of light intensity transmitted by the obstacle. Furthermore, we will increase the flow field complexity, maintaining the cellular flow configuration with distinct upwelling and downwelling regions, in order to see the effect of having upwelling under ice or under a lead/thinner ice.

FIGURE 6 . 4 :

 64 FIGURE 6.4: Instantaneous normalized population density fields θ(x, z, t * )/ θ at a fixed instant of time t * = 600h, where θ stands for the spatial average. The obstacle is represented by the black rectangle and the dashed lines delimit the light obstruction zone. The white lines are the isolines θ/ θ = 1. The solid black lines represent flow streamlines from Eq. 5.1, with arrows indicating the circulation direction. Also, t * is in the regime of stationary per-capita growth rate (r p (t) = r p ). In a), b) and c) U = (0, 0.71, 1.19) m h -1 (respectively) and D = 1 cm 2 s -1 . In d), e) and f) U = (0, 0.71, 1.19) m h -1 (respectively) and D = 10 cm 2 s -1 . a), d) and e) are bloom cases while b), c) and f) are no bloom cases.

FIGURE 6 . 5 :

 65 FIGURE 6.5: Population per-capita growth rate r p (t), normalized by the intrinsic net growth rate r b , as a function of time in the steadyflow case, for various values of U and D, with a light obstruction of size L x /2 and zero transmitted light, centered at L x /2, localized at z = 0. The inset shows the same graph with the addition of the noflow cases.

FIGURE 6 . 6 :

 66 FIGURE 6.6: Instantaneous normalized population density fields θ(x, z, t * )/ θ at a fixed instant of time t * = 600h, where θ stands for the spatial average. The obstacle is represented by the black rectangle and the dashed lines delimit the light obstruction zone. The white lines are the isolines θ/ θ = 1. The solid black lines represent flow streamlines from Eq. 5.10, with arrows indicating the circulation direction. Also, t * is in the regime of stationary per-capita growth rate (r p (t) = r p ). In a) and b) U = 0.71 m h -1 and D = (10, 1) cm 2 s -1(respectively). Both are bloom cases.

FIGURE 6 . 7 :

 67 FIGURE 6.7: Population per-capita growth rate r p (t), normalized by the intrinsic net growth rate r b , as a function of time in the multiscale flow case, with a light obstruction of size L x /2 and zero transmitted light, centered at L x /2, localized at z = 0.

  . . . 18 xvi 2.10 A map of the apparent horizontal diffusivity (K) in the Atlantic obtained through measurements of eddy kinetic energy by drifters. Figure extracted from [33]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.11 Pictorial representation of ocean flows at various scales, increased mixing is represented by stippling. Figure extracted from [19]. . . . . . 2.12 Pictorial representation of open ocean deep convection, Buoyancy flux through sea surface is represented by curly arrows, stratification is represented by continuous lines. (a) preconditioning (buoyancy loss at the surface), (b) deep convection (currents are driven by density gradients), and (c) lateral exchange and spreading (the phenomenon slowly dissipates). Figure extracted from [34]. . . . . . . . . . . . . . . . 2.13 Photo of the surface of the Garonne river taken by the author during a stay in Bordeaux, France. A constant flux of wind over the water surface produced a series of vortical motions with a life span averaging from a couple of seconds to a minute. . . . . . . . . . . . . . . . . . . . 2.14 Langmuir streaks on Quake Lake, United States. Image by Wayne Wurtsbaugh -www.aslo.org. . . . . . . . . . . . . . . . . . . . . . . . .

3.1 Examples of planktonic microorganisms and macro organisms. (I) On

  Population per-capita growth rate r p (t), normalized by the intrinsic net growth rate r b , as a function of time in the multiscale flow case, with a light obstruction of size L x /2 and zero transmitted light, centered at L x /2, localized at z = 0. . . . . . . . . . . . . . . . . . . . . . . .[START_REF] Abraham | The generation of plankton patchiness by turbulent stirring[END_REF]

2 s -1 (respectively). Both are bloom cases. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxi 6.7 xxiii List of Tables 4.1 Parameters of the biological dynamics. . . . . . . . . . . . . . . . . . . . 52 xxv Dedicated to my parents, Rosa Maria Beltran Tergolina and Paulo Afonso Tergolina.

TABLE 4 . 1 :

 41 Parameters of the biological dynamics.

	Parameter	Value	Meaning
	κ bg	0.2 m -1	Background turbidity
	κ	1.5 • 10 -11 m 2 cell -1	Specific light attenuation
			of phytoplankton
	p max	0.04 h -1	Maximal specific production rate
	l	0.01 h -1	Specific loss rate
	H	30 µmol photons m -2 s -1	Half-saturation constant of
			light-limited growth
	I in	350 µmol photons m -2 s -1	Incident light intensity
	v	0.04 m h -1

Acknowledgements

Acknowledgements xi

thanks to the Université de Lille, for the financial support for this project and for the professional opportunity, both in scientific research and teaching. I thank as well the people of Brazil and France

Code T = 15.5h Code T = 93h Theoretical Result T = 93.h Theoretical Result T = 15.5h

Chapter 4. Numerical Methods while still obtaining good convergence. Again, as described in this chapter, the condition required for the stochastic implementation of diffusion is that dz 2 /D < dt, so we started from dt = 0.01h and found that until dt = 0.5h we obtained good convergence for the D = 10cm 2 /s case (small enough to resolve well the reaction time given by 1/p max = 25h). In addition, the figure shows that D = 10cm 2 /s is a bloom case, in agreement with the diagram of Fig. 3.8 (top), extracted from [START_REF] Huisman | How do sinking phytoplankton species manage to persist?[END_REF]. We then performed multiple simulations with several values of D. The results are reported in Fig. 4.7 for a depth of 20 meters, close to the smallest one for which the minimal turbulence threshold exists (see Fig. 3.8, top). These results are also in agreement with the diagram of Fig. 3.8, as we see that when D = 10 -2 cm 2 /s we arrive at the first no bloom case since P tot decreases after enough time has elapsed. In the inset we can also see that D = 0.5cm 2 /s is the optimal case for growth for the chosen z max as the steady state value of P tot is the highest. To complement the code validation we decided to compare our results to the original ones by reproducing the full bloom diagram of Fig. 3.8 (top) with an appropriate number of simulations. In order to efficiently explore the parameter space and limit the use of computational resources, we performed 25 simulations spanning a a wide range of D and z max values as shown in figure 4.8. These simulations were carried out until P tot reached a constant steady-state value or a very small one. The first corresponds to a bloom case where population survives, and the latter to a no bloom case, where population vanishes. With these results we were able to obtain the same "phase" diagram as in the original work [START_REF] Huisman | How do sinking phytoplankton species manage to persist?[END_REF], which proves that our code behaved properly. The horizontal dashed line in the top part of A further comparison can be made through population density vertical profiles, which are straightforward to obtain in 1D simulations as they correspond to θ(z, t) at a given time. In Fig. 4.9 we show our vertical profiles in the same conditions as the precedent flow cases (Fig. 5.3 and 5.11), a result of the smaller vortices present in the flow field, the distribution of θ is more organized in here than in Fig. 5.14. In fact the signature of the largest-scale flow is much more apparent here, particularly in the surface θ patch at the center of the domain (x Lz = 60m), close to the surface. This hinted us again that dynamics can be quite similar as they appear to remain under preeminence of the large scale structures described by Ψ L . We return to the per capita growth rate analysis in Fig. 5.16a where we expose the same U, D cases as in the previous sections (Figs. 5.5 and 5.12). Here the dynamics are quite clearly more chaotic than in the mentioned cases, we can see that growth rate oscillation amplitudes increase considerably in relation to the unsteady flow case, although their periods are still proportional to to L z /U. Also one can question weather rp is still a robust quantifier for these dynamics as we cannot spot as easily as in the unsteady case whether r p (t) curves are oscillating around a constant value. In Fig. 5.16b we take the case with the highest oscillation amplitude (U = 2.466) from Fig. 5.16a and plot its total biomass dynamics vs. time (black), we then use the first value of the curve as a base point and plot this value multiplied by exp(tr p ) with rp being extracted from Fig. 5.16a. The result shows that rp is still a robust quantifier for describing the growth dynamics of the simulations.

Moreover we can also identify in Fig. 5.16a that rp dependence on U is still that of monotonic decreasing function, as in Secs. 5.2 and 5.3. Consequently we perform the same type of comparison of the last section involving rp for the different cases, so from the quantitative point of view, it is interesting to compare the asymptotic growth rate r p /r b measured in this multiscale setting and in the previous ones, e.g. in the 1-mode unsteady-flow case. As it can be seen in Fig. 5.17, the estimates from the two cases are still quite correlated, and diffusivity plays a rather weak role (similarly to the previous comparison, reported in Fig. 5.13). Still we notice that in Fig. 5.17, for small enough D, it is possible to observe that the multiscale estimate of the growth rate r p tends to be slightly larger than its counterpart in the absence of small eddies. We then argue that the latter flow features provide an effect that partially compensates the deadly action of the large-scale coherent flow, but that is only measurable for D ≤ 10cm 2 s -1 .