Modelling the exchange of microplastics between the inner surf zone and the open sea under different beach configurations

1 TREBALL DE FI DE GRAU TITLE: Modelling the exchange of microplastics between the inner surf zone and the open sea under different beach configurations DEGREE: Enginyeria de Sistemes Aeroespacials AUTHORS: David Balboa Mato DIRECTORS: Francesca Ribas Prats Albert Falqu´ es i Serra DATA: November 23, 2020 2 3 T´ıtol: Modelitzaci´ o de l’intercanvi de micropl` astics entre la zona de rompents interior i mar obert sota diferents configuracions de platja Autor: David Balboa Mato Directors: Francesca Ribas Prats Albert Falqu´ es i Serra Data: 23 de novembre de 2020 Resum En aquest treball s’ha modelitzat l’intercanvi de micropl` astics entre la zona de rompents interior i mar obert sota diferents configuracions de platja. El dany que causen els pl` astics a l’oce` a i la fauna que aquest cont´ e´ es un problema emergent i cada vegada m´ es estudiat. Per arribar a comprendre en la seva totalitat la din` amica dels micropl` astics als mars i oceans s’ha de con` eixer l’efecte de la morfodin` amica de les platges. La platja ´ es un lloc complex i dif´ ıcil de predir perqu` e ocorren multitud d’esdeveniments simult` aniament a diferents escales que cont´ ınuament la deformen. Utilitzant un model anomenat morfo55 desenvolupat pel grup Nonlinear Fluid Dynamics de la UPC s’han simulat diferents estats morfodin` amics en qu` e pot estar una platja. Amb aquestes dades, i utilitzant el programari Matlab, s’ha desenvolupat un algoritme capac¸ de modelitzar la din` amica de la concentraci´ o de micropl` astics sota les diferents configuracions. L’objectiu principal d’aquest projecte ´ es quantificar com influeix la formaci´ o d’un patr´ o de sorra r´ ıtmic al llarg de la costa anomenat barra cresc` entica, amb el corresponent sistema de corrents, en la din` amica d’aquests micropl` astics en la superf´ ıcie de l’aigua. Els resultats indiquen que els micropl´ astics es mouen influenciats principalment per Stokes drift. L’aparici´ o d’una barra crescentica que ocasiona un sistema de deep-currents pot accelerar el moviment cap a la costa. Incl` os pot ocasionar esquinc¸aments on el pl` astic trobi impediment per avanc¸ar, fins i tot zones on quedi estancat. 4 5 6 Title : Modelling the exchange of microplastics between the inner surf zone and the open sea under different beach configurations Author: David Balboa Mato Advisors: Francesca Ribas Prats Albert Falqu´ es i Serra Date: November 23, 2020 Overview In this project we have modelled the exchange of microplastics between the inner surf zone and the open sea under different beach configurations. The damage caused by plastics to the ocean and the fauna it contains is an emerging and increasingly studied problem. To fully understand the dynamics of microplastics in the seas and the oceans, the effect of the morphodynamics of the beaches must be studied. The beach is a complex and unpredictable place because many events occur simultaneously at different scales that continuously deform it. Using a model called morfo55 developed by the Group of Nonlinear Fluid Dynamics of the UPC, we have simulated different morphodynamic states that can occur on a beach. With this data, and using Matlab software, we have developed an algorithm capable of modelling the dynamics of microplastic concentration under the different configurations. The main objective of this project is to quantify how the formation of a rhythmic sand pattern along the coast called crescentic bar and the corresponding rip-current system influences the dynamics of these plastics on the surface of the water. The results indicate that the microplastics moves influenced mainly by Stokes drift. The appearance of a crescent bar causes a system of deep-currents, which accelerates the onshore movement. In turn, it can cause tears where the plastic finds an impediment to advance, even areas where it is stagnant. CONTENTS CHAPTER 1. Introduction ........................... 1 1.1. The plastic ocean. ............................... 1 1.2. Nearshore morphodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3. Outline of the thesis .............................. 5 CHAPTER 2. Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1. Governing equations and boundary conditions ............... 7 2.1.1. Governing equation for concentration . . . . . . . . . . . . . . . . . 7 2.1.2. Parametrization of γand~v....................... 8 2.1.3. Boundary and initial conditions . . . . . . . . . . . . . . . . . . . . 9 2.2. Discretization of the problem ......................... 10 2.2.1. Physical system and mesh . . . . . . . . . . . . . . . . . . . . . . 10 2.2.2. Discretization of the equation . . . . . . . . . . . . . . . . . . . . . 12 2.2.3. Discretization of the boundary conditions . . . . . . . . . . . . . . . 13 2.3. Study cases ................................... 15 2.3.1. Analytical cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3.2. Morfo55 cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 CHAPTER 3. Validation and sensitivity to numerical method . . . 19 3.1. Analytical case ................................. 19 3.1.1. CFL stability condition . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.1.2. Sensitivity of the analytical case to the numerical method . . . . . . 20 3.2. Morfo55 cases .................................. 23 CHAPTER 4. Results of plastic dynamics . . . . . . . . . . . . . . . . 25 4.1. Input data obtained from the morfo55 model ................. 25 4.2. Results for shore-normal wave incidence . . . . . . . . . . . . . . . . . . 28 4.3. Results for shore-oblique wave incidence . . . . . . . . . . . . . . . . . . 36 CHAPTER 5. Discussion and conclusions . . . . . . . . . . . . . . . . 41 8CONTENTS 5.1. Discussion of the results ............................ 41 5.2. Final conclusions ................................ 43 Bibliography .................................... 45 LIST OF FIGURES 1.1 Open sea plastic soup known as Garbage Island (Source: bbc.com). Center of the North Pacific Ocean. Between the coasts of California and Hawaii. . . . . . 1 1.2 Nearshore zone with its different parts (Source: Garnier, 2006) . . . . . . . . . 3 1.3 Time exposure images of a straight bar configuration (a) and (b) a crescentic bar configuration, in Duck, North Carolina, USA. The coast is at the top of the images. Courtesy of Prof.R.Holman, Oregon State University. Figure adapted from Garnier et al. 2003. ............................. 4 2.1 Physical system and frame of reference (Source: Garnier, 2006). . . . . . . . . 10 2.2 Staggered grid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.3 Sketch of the beach with he location of zones 1 and 2, where initial plastic concentration is placed. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.1 Comparison between the Euler method and the Adams - Bashforth method in the analytical case by varying the parameter cofor (a) co=1, and (b) co=0.52 after 1 second with a velocity of 0.25 m s. . . . . . . . . . . . . . . . . . . . . . 21 3.2 Result of he numerical 2D solution in the analytical case with Adams-Bashfort scheme using, (a) co=5(unstable solution), and (b) co=0.25 (stable solution) after 0.15 second with a velocity of 0.5 m son both axes. . . . . . . . . . . . . . 21 3.3 Comparison of he 1D analytical and numerical solutions with Adams - Bashforth scheme for (a) co=5, and (b) co=0.1after 1 second with a velocity of 0.25 m s.22 3.4 Comparison between (a) analytical solution and (b) numerical solution using the Adams-Bashforth scheme with co=0.1after 0.15 second with a velocity of 0.5 m son both axes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.5 Result of the 2D numerical solution corresponding to a morfo55 case with Dirichlet boundary condition for (a) co=0.7, and (b) co=0.6after 1 hour. . . . . . . 24 3.6 Result of the 2D numerical solution corresponding to a morfo55 case with mixed boundary condition for (a) co=0.7, and (b) co=0.6after 1 hour. . . . . . . . . 24 4.1 Data obtained with morfo55 for θ=0o. Panel (a) shows the alongshore-averaged sea bed level, wave height and alongshore current; panel (b) represents the bottom topographic perturbation along with the velocity field; panel (c) shows the diffusivity factor; and panel (d) displays the intensity of the velocity field where the Stokes drift has been included. . . . . . . . . . . . . . . . . . . . . . . . . 25 4.2 Data obtained with morfo55 for θ=5o. Panel (a) shows the alongshore-averaged sea bed level, wave height and alongshore current; panel (b) represents the bottom topographic perturbation along with the velocity field; panel (c) shows the diffusivity factor; and panel (d) displays the intensity of the velocity field where the Stokes drift has been included. . . . . . . . . . . . . . . . . . . . . . . . . 26 4.3 Data obtained with morfo55 for θ=10o. Panel (a) shows the alongshoreaveraged sea bed level, wave height and alongshore current; panel (b) represents the bottom topographic perturbation along with the velocity field; panel (c) shows the diffusivity factor; and panel (d) displays the intensity of the velocity field where the Stokes drift has been included. . . . . . . . . . . . . . . . . . . 27 CHAPTER 2. METHODOLOGY In this chapter, the methodology used for the analytical and numerical development of the project is explained. It is divided into three sections. The first section presents the equation that governs the phenomenon of advection and diffusion of plastic concentration and the corresponding boundary conditions. The second section describes the discretization applied to these equations. Finally, the third section introduces the different study cases, starting by the comparison between the numerical results and an analytical case that will allow us to validate the discretization, and followed by the different cases analyzed using results of the morfo55 model. 2.1. Governing equations and boundary conditions 2.1.1. Governing equation for concentration The volumetric concentration of a contaminant in a turbulent fluid in motion can be described, in the framework of the shallow-water approximation. by the depth-integrated advection-diffusion equation, ∂C ∂t+(~v·∇)C | {z } Advection =∇·(γ ∇C) | {z } Di f f usion , (2.1) where C is the concentration, which can be interpreted as grams of plastic per m2of water surface, ~v= (u,v)is the horizontal velocity field, γis the diffusion coefficient, ∇= (∂ ∂x,∂ ∂y)represents the horizontal nabla operator, x is the cross-shore coordinate, y is the alongshore coordinate and t is time. On the one hand, we have the advective term that describes the variation of the concentration location due to the presence of the flow. For example, if in a certain instant of time our concentration is located in a certain place, some time later it will have shifted to a different position due to the presence of the velocity field. On the other hand, the diffusive term describes the dynamics of the concentration due to the spatial differences in the amount of concentration. For instance, if a droplet of ink is placed on a clean water surface, it will tend to diffuse in time due to the random component of molecular motion. Notice that equation (2.1) is linear and, thereby, the scale of C is arbitrary. The only exception would be to consider a flux of plastic at the boundary but we will never use such boundary condition. Expanding the nabla operator in Cartesian coordinates we get the following equation: ∂C ∂t=−u∂C ∂x−v∂C ∂y+∂γ ∂x ∂C ∂x+∂γ ∂y ∂C ∂y+γ∂2C ∂x2+∂2C ∂y2, (2.2) 7 8 Modelling the exchange of microplastics between the inner surf zone and the open sea under different beach configurations 2.1.2. Parametrization of γand~v As can be seen in equation (2.2), it is necessary to parameterize some variables. The diffusion constant γin the nearshore is related to turbulence, which occurs because part of the water momentum is transferred into small turbulent eddies. A good parameterisation for this constant in the nearshore uses the fact that wave breaking is the most important source of turbulence (Ribas, 2004): γ=MD ρ1 3 D, (2.3) where M is a parameter of O(1) that characterizes the turbulence, Dis the rate of wave energy dissipation due to breaking water per unit area, ρis the water density and Dis the water depth. The dissipation of wave energy is due to the transfer of energy from the ordered wave motion to the turbulent eddies (Battjes et al., 1990). Given the statistical description of both broken and non-broken waves (Thornton and Guza, 1983), the rate of energy dissipation per unit area can be estimated as follows (Ribas, 2004): D=3√π 16 ρgB3fp H5 λ2 cD3 1− 1+H λcD2!−5 2 , (2.4) where gis the gravity, Bindicates a parameter describing the type of breaking, fpis the frequency peak of wave field, His the wave height and λcis a parameter giving the expected saturation value of H/D. The velocity field ~va complicated 3D field. As explained in section 1.3 the governing equations and variables used by the morfo55 model are timeand depth-averaged with the intention of simplifying the physical system. The velocity ~vaffecting the plastic contains is in principle the depth-averaged currents that emerge in the surf zone.However, an extra component must be added related with the onshore mass transport that is produced due to wave orbital motion in the upper part of the water column (Stokes drift). Due to this second component, when a particle is floating experiences a speed in the direction of wave propagation even in the absence of mean currents. The final velocity field affecting the plastic particles is supposed to be: vi(x1,x2,t) = Ui(x1,x2,t)+ ˜vi(x1,x2,t), i=1,2 , (2.5) There, Uirepresents the depth-averaged ’mean’ currents that may exists, which are given directly from the morfo55 model. and ˜viis the ’Stokes drift’ component. To calculate the second term (’Stokes drift’) we use the following equation: ˜vi=a2kD tanh(2kD)·H2 8Dsgk tanh(kD)· ¯ ki k, (2.6) CHAPTER 2. METHODOLOGY 9 where~ kis the wave vector. The constant a, taken equal to 0.5 in the present study, accounts for the unknowns behind the role of Stokes drift. If the microplastic is completely floating and there is no mean current, then a = 1. If the microplastic concentration is found along the whole water column and/or there are strong mean currents, then a = 0. It should be noted that this parameterization has been used with the morfo55 beach configurations. Before using these realistic parameterizations, initial tests are performed, with viand γtaken as constants. 2.1.3. Boundary and initial conditions The boundary conditions used are called Dirichlet, Neumann and Periodic. Some of these conditions can also be combined, obtaining what is known as mixed conditions. The corresponding equations allow us to model different situations that may be interesting for the study. •Dirichlet boundary conditions: Such conditions consist of assigning specific values to the variables in the domain boundary, and are also known as first-class conditions. They can be applied to ordinary differential equations (ODE) and partial differential equations (PDE). Being Ωthe solution domain: C(x,y,t) = φ(x,y,t)∀x,y∈∂Ω Where φis a known value. (2.7) •Neumann boundary conditions: Such conditions consist of defining a known value of the variable derivative at the boundary of the domain. They are also called second-class conditions and can be applied to ODEs and PDEs. ∂C(x,y,t) ∂x=φ(x,y,t)∀x,y∈∂Ω Where φis a known value. (2.8) •Mixed boundary conditions: They are a linear combination of the Dirichlet and Neumann conditions, typically used to impose a known value of a total flux. For a problem with flux equal to zero (insulating boundary), must fulfill that: 0=−γ(x,y,t)∇C(x,y,t)+C(x,y,t)¯v(x,y,t)∀x,y∈∂Ω,(2.9) •Periodic boundary conditions: These conditions consist in forcing that the values of the solution are equal at both sides of the domain, thereby imposing a certain alongshore periodicity. In the present work, this condition is imposed at the lateral boundaries, hence: C(x,0,t) = C(x,Ly,t)∀x∈∂Ω,(2.10) 10 Modelling the exchange of microplastics between the inner surf zone and the open sea under different beach configurations 2.2. Discretization of the problem 2.2.1. Physical system and mesh The domain of study is the nearshore, including the shoaling zone, the surf zone and the swash zone (section 1.2). As you can see in Figure 2.1 the domain is rectangular and the Cartesian coordinate system has the origin at the shoreline. The x-axis, is cross-shore directed, points seaward and the offshore boundary is set at x=Lx.The y-axis represents the longshore direction, and the lateral boundaries are located at y=0and y=Ly. The z-axis is the vertical direction and points upwards.. Figure 2.1: Physical system and frame of reference (Source: Garnier, 2006). CHAPTER 2. METHODOLOGY 11 The problem is discretized with a finite differences method. The spatial domain is discretized using constant grid sizes ∆xand ∆y, and the time is discretized with a constant step ∆t. The mesh used is known as collocated grid (all variables are computed in the center of the cells), and differs from the staggered grid used in morfo55 model (Fig. 2.2). In the staggered grid the scalar variables are usually computed in the center of the cells while other variables such as velocity or momentum are computed in the walls of the cell. Thus, it has been necessary to manipulate some morfo55 variables in order to have them located in the center of the cell. This has been the case of the depth-averaged velocities. Figure 2.2: Staggered grid To obtain these velocities in the center of the cells, a simple linear interpolation has been applied: ui,j=ui+1/2,j+ui−1/2,j 2(2.11) and, vi,j=vi,j+1/2+vi,j−1/2 2(2.12) 12 Modelling the exchange of microplastics between the inner surf zone and the open sea under different beach configurations 2.2.2. Discretization of the equation The discretization of the advection-diffusion equation (2.2) has been designed in order to achieve an approximation with the smallest possible error. The spatial derivatives (first or second order), have been approximated using central differences because they are more accurate than forward or backward differences. For the time derivative two methods have been tested. The first one is known as Euler method, it uses a first order approximation and it is normally used for simple approximations where a high accuracy is not needed. The second one, known as the Adams-Bashforth method, is a multipass method, and it gives more accuracy (while the Euler method uses only the previous time step, the AdamsBashforth method uses the two previous time steps). When using a multipass method such as Adams-Bashforth, the first iteration normally must be still done with the Euler method because in the first instant of time there are no two previous steps available. •First order spatial derivatives with central differences: ∂C(x,y,t) ∂x=Cn−1 i+1,j−Cn−1 i−1,j 2∆x(2.13) ∂C(x,y,t) ∂y=Cn−1 i,j+1−Cn−1 i,j−1 2∆y(2.14) •Second order spatial derivatives with central differences: ∂2C(x,y,t) ∂y2=Cn−1 i+1,j−2Cn−1 i,j+Cn−1 i−1,j ∆x2(2.15) ∂2C(x,y,t) ∂y2=Cn−1 i,j+1−2Cn−1 i,j+Cn−1 i,j−1 ∆y2(2.16) •Euler method: ∂C(x,y,t) ∂t=Cn i,j−Cn−1 i,j ∆t(2.17) •Adams - Bashforth method: Cn i,j−Cn−1 i,j ∆t=3 2φCn−1 i,j−1 2φCn−2 i,j(2.18) where φCn−1 i,jis the algorithm result of the previous iteration and φCn−2 i,jis the algorithm result of the two previous iterations. The governing equation discretized by Euler’s method has the following form: Cn i,j=Cn−1 i,j+∆t·(An−1+Bn−1+Cn−1), (2.19) while using the Adams method - Bashforth method, equations (2.2) reads: Cn i,j=Cn−1 i,j+∆t 2·(3(An−1+Bn−1+Cn−1)−(An−2+Bn−2+Cn−2)) (2.20) CHAPTER 2. METHODOLOGY 13 there A, B and C are as follows: A=γi+1,j−γi−1,j 2∆x Ci+1,j−Ci−1,j 2∆x+γi,j+1−γi,j−1 2∆y Ci,j+1−Ci,j−1 2∆y, (2.21) B=γi,jCi+1,j−2Ci,j+Ci−1,j ∆x2+Ci,j+1−2Ci,j+Ci,j−1 ∆y2, (2.22) C=−ui,j Ci+1,j−Ci−1,j 2∆x+vi,j Ci,j+1−Ci,j−1 2∆y(2.23) 2.2.3. Discretization of the boundary conditions The boundary conditions are different in the three types of boundaries. Note also that the governing equation can be applied at all points of the domain, including the boundaries. 2.2.3.1. Shore boundary At x=0(i=1), the shore boundary, we can use Dirichlet or mixed boundary conditions (section 2.1.3). In the former case, the value of the concentration at the boundary is imposed (Dirichlet boundary condition). The imposed value is Ci=1,j=0, which represent a situation where the plastic concentration at the dry beach is continuously cleaned up. Mixed boundary conditions have been also applied o represent the situation of an insulating boundary, i.e. without exchange of plastic between the dry beach and the surf zone. This type of boundary conditions are more complicated to impose. We start by applying the governing equation itself, where some terms appear that refer to cells that are outside the domain. Terms A, B and C of the governing equation when we use the Adam - Bashforth method are (equations 2.21 -2.23): A=γ2,j−γ0,j 2∆x C2,j−C0,j 2∆x+γ1,j+1−γ1,j−1 2∆y C1,j+1−C1,j−1 2∆y, B=γ1,jC2,j−2C1,j+C0,j ∆x2+C1,j+1−2C1,j+C1,j−1 ∆y2, C=−u1,j C2,j−C0,j 2∆x+v1,j C1,j+1−C1,j−1 2∆y The cell C0,jis known as a dummy point that does not belong to our domain. Then, we use the discretized version of the mixed boundary condition (2.9) to obtain the value of this cell. Dummypoint z}|{ C0,j=C2,j−2∆x γ1,ju1,jC1,j 14 Modelling the exchange of microplastics between the inner surf zone and the open sea under different beach configurations 2.2.3.2. Offshore boundary At x=Lx(i=Nx), offshore boundary, again we can use Dirichlet or mixed boundary conditions. The same strategies than in the shore boundary are used. In this case the terms A, B and C: A=γNx+1,j−γNx−1,j 2∆x CNx+1,j−CNx−1,j 2∆x+γNx,j+1−γNx,j−1 2∆y CNx,j+1−CNx,j−1 2∆y, B=γNx,jCNx+1,j−2CNx,j+CNx−1,j ∆x2+CNx,j+1−2CNx,j+CNx,j−1 ∆y2, C=−uNx,j CNx+1,j−CNx−1,j 2∆x+vNx,j CNx,j+1−CNx,j−1 2∆y Again, applying the mixed boundary condition we obtain the value at the dummy point outside the domain that allows us to apply the governing equation. Dummypoint z }| { CNx+1,j=CNx−1,j+2∆x γNx,juNx,jCNx,j, 2.2.3.3. Lateral boundaries When y=0(j=1), and when y=Ly(j=Ny), the lateral boundaries, periodic boundary conditions are applied. First, we must impose the governing equation in the boundary, Again, the values of C at dummy points outside the domain in the y-axis must be found by applying the lateral boundary conditions (2.10) to keep alongshore periodicity. When j = 1, A=γi+1,1−γi−1,1 2∆x Ci+1,1−Ci−1,1 2∆x+γi,2−γi,0 2∆y Ci,2−Ci,0 2∆y, B=γi,1Ci+1,1−2Ci,1+Ci−1,1 ∆x2+Ci,2−2Ci,1+Ci,0 ∆y2, C=−ui,1 Ci+1,1−Ci−1,1 2∆x+vi,1 Ci,2−Ci,0 2∆y Where we use: Ci,0=Ci,Ny−1 When j = Ny, A=γi+1,Ny−γi−1,Ny 2∆x Ci+1,Ny−Ci−1,Ny 2∆x+γi,Ny+1−γi,Ny−1 2∆y Ci,Ny+1−Ci,Ny−1 2∆y, B=γi,NyCi+1,Ny−2Ci,Ny+Ci−1,Ny ∆x2+Ci,Ny+1−2Ci,Ny+Ci,Ny−1 ∆y2, C=−ui,Ny Ci+1,Ny−Ci−1,Ny 2∆x+vi,Ny Ci,Ny+1−Ci,Ny−1 2∆y where we use Ci,Ny+1=Ci,2(2.24) CHAPTER 2. METHODOLOGY 15 2.3. Study cases 2.3.1. Analytical cases In order to verify that the chosen discretization reproduces correctly the dynamics, a few cases are initially solved that correspond to situations with an analytical solution. This analytical cases allow to test the different discretizations with the intention of validating them. Moreover, other discretizations of the time derivative, such as the Leapfrog method, can be tested. The advection-diffusion equation (2.1) in case of constant~vand γand Dirichlet boundary conditions at all the boundaries are used. Analytical solutions of this problem exist that consist of Gaussian functions in 1D or 2D. The 1D Gaussian function reads: C(x,t) = 1 √4πγ te−(x−(xo+ut))2 4γt(2.25) and the 2D Gaussian function reads: C(x,y,t) = 1 p8πγxt2γy e−(x−(xo+ut))2 4γxt+(y−(yo+vt))2 4γyt(2.26) Where (xo,yo)is the position of the initial Gaussian distribution, γis the diffusivity, γxis the diffusivity in x-axis direction and γyis the diffusivity in y-axis direction. To prove that the phenomena of diffusion and advection are faithfully described with the discretized equation, 1D and 2D Gaussian functions are used as initial conditions. Then we numerically model how they spread and migrate and compare with the analytical result. For the analytical study, the boundary conditions imposed on both axes are: C(x=0,y,t) = C(x=Lx,y,t) = 0and C(x,y=0,t) = C(x,y=Ly,t) = 0(2.27) 22 Modelling the exchange of microplastics between the inner surf zone and the open sea under different beach configurations except for the boundary. As can be seen in the zoom of Figure 3.3(b), in the boundary the analytical and numerical results are not completely equal, and this is due to the boundary conditions imposed. The condition C=0is imposed at x=0and x=Lx in the numerical case while C=0is imposed at x=−∞and x=∞in the analytical case. In the 2D case, the analytical and numerical results are also identical (Figure 3.4a,b). Thereby, it has been confirmed that the numerical schemes are able to capture the advection and he diffusion present in this analytical case using the Adams-Bashforth scheme and co≤0.25. This has been proved for both the 1D and the 2D cases. In the rest of project, the Adams-Bashforth is always used. 0 0.5 1 1.5 2 2.5 Distance 0 0.2 0.4 0.6 0.8 1 1.2 C Analytical vs Numerical solution 1D with Co = 0.5 Numerical solution Analytical solution ((a)) 2.25 2.3 2.35 2.4 2.45 2.5 Distance 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 C Analytical vs Numerical solution 1D with Co = 0.1 Numerical solution Analytical solution ((b)) Figure 3.3: Comparison of he 1D analytical and numerical solutions with Adams - Bashforth scheme for (a) co=5, and (b) co=0.1after 1 second with a velocity of 0.25 m s. Analytical solution 2D with Co = 0.1 0 0.2 0.4 0.6 0.8 1 X [m] 0 0.2 0.4 0.6 0.8 1 Y [m] 2 4 6 8 10 12 C [g/m2] ((a)) Numerical solution 2D with Co = 0.1 0 0.2 0.4 0.6 0.8 1 X [m] 0 0.2 0.4 0.6 0.8 1 Y [m] 0 2 4 6 8 10 12 C [g/m2] ((b)) Figure 3.4: Comparison between (a) analytical solution and (b) numerical solution using the Adams-Bashforth scheme with co=0.1after 0.15 second with a velocity of 0.5 m son both axes. CHAPTER 3. VALIDATION AND SENSITIVITY TO NUMERICAL METHOD 23 3.2. Morfo55 cases When the dynamics of plastic concentration is modelled using the realistic hydrodynamics obtained with morfo55 model, sensitivity analysis to the numerical parameter is performed in a similar way. Since γand ~vare now variable in space, equation (2.2) must be used, with he additional terms related with the spatial derivative of γ. This equation again has parabolic and hyperbolic parts so that the two CFL conditions (3.2) and (3.3) apply. As explained in chapter 2, morfo55 generates files where it stores the hydrodynamic data in the spatial mesh used by the morfo55 model. This mesh has specific values of ∆x,∆y, Nx,Ny,Lx,Lythat provide enough accuracy for morfo55 computations. However, this accuracy tuns out to be insufficient to model the dynamics of plastic concentration and it has been necessary to modify this mesh in order to have more cells, so that our numerical algorithm converges. The reason is related with the behaviour of Cnear the boundaries due to the imposed boundary conditions. These conditions are demanding numerically speaking, since in a small space abrupt changes occur, hence requiring a dense spatial mesh. Two situations cause these abrupt changes. In case of having mixed boundary conditions with flux equal to zero (q=0), the insulating wall condition causes the plastic to literally accumulate at the boundary, and induce huge gradients. In case of imposing Dirichlet boundary conditions (C=0) at the coast, and given that advective terms turn out to dominate and move the plastic towards the coast, plastic concentration experience a large gradient. That is why it is necessary to set a fine enough mesh so that the model is able of representing well the conditions imposed. To create more cells, we interpolate the values of the morfo55 variables of two contiguous cells. With this, we manage multiply the number of the cells along the x-axis, where these abrupt changes mainly occur. The morfo55 mesh has 47 cells on the x-axis, and 200 on the y-axis. The y-axis does not cause any problem because the lateral boundary are periodic hardly constraining the dynamics which is mainly governed by cross-shore processes conditions. The need of more accuracy occurs in the x-axis and it is necessary to increase the number of cells up to 277. The final parameter values used in the morfo55 cases (for my model) are ∆x=0.85m, ∆y=10m, Nx=277,Ny=201,Lx=237.447m, Ly=2000m. Once an adequate mesh is obtained, the sensitivity to the time step tand the conditions for numerical stabilities has been studied. The results illustrated in the Figures 3.5 and 3.6 show that numerical stability is achieved with both boundary conditions (Dirichlet BC and Mixed BC). It should be noted that the Dirichlet boundary conditions were more permissive; without the interpolation performed on the mesh, we obtained results without numerical instabilities. In contrast, mixed boundary conditions always showed instabilities without the proper mesh. 24 Modelling the exchange of microplastics between the inner surf zone and the open sea under different beach configurations ((a)) ((b)) Figure 3.5: Result of the 2D numerical solution corresponding to a morfo55 case with Dirichlet boundary condition for (a) co=0.7, and (b) co=0.6after 1 hour. ((a)) ((b)) Figure 3.6: Result of the 2D numerical solution corresponding to a morfo55 case with mixed boundary condition for (a) co=0.7, and (b) co=0.6after 1 hour. CHAPTER 4. RESULTS OF PLASTIC DYNAMICS 4.1. Input data obtained from the morfo55 model The advection-diffusion code needs the velocity and the diffusivity fields generated with the morfo55 model. These outputs of morfo55 heavily depend on the angle of incidence of the waves with respect to the shore-normal, θ. Figures 4.1,4.2 and 4.3 show these outputs for θ=0o, 5oo and 10o, respectively. It should be noted that the angle of incidence θ=20o has also been examined but is here excluded due to the great similarity with the results for the angle θ=10o. 0 50 100 150 200 250 X [m] -5 -4 -3 -2 -1 0 1 Zb [m] -1 -0.5 0 0.5 1 v [m/s] Basic state m55 alongshore-averaged H [m] Hf Mean H Mean Vel Mean Figure 4.1: Data obtained with morfo55 for θ=0o. Panel (a) shows the alongshoreaveraged sea bed level, wave height and alongshore current; panel (b) represents the bottom topographic perturbation along with the velocity field; panel (c) shows the diffusivity factor; and panel (d) displays the intensity of the velocity field where the Stokes drift has been included. 25 26 Modelling the exchange of microplastics between the inner surf zone and the open sea under different beach configurations 0 50 100 150 200 250 X [m] -5 -4 -3 -2 -1 0 1 Zb [m] 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 v [m/s] Basic state m55 alongshore-averaged H [m] Hf Mean H Mean Vel Mean ((a)) ((b)) ((c)) ((d)) Figure 4.2: Data obtained with morfo55 for θ=5o. Panel (a) shows the alongshoreaveraged sea bed level, wave height and alongshore current; panel (b) represents the bottom topographic perturbation along with the velocity field; panel (c) shows the diffusivity factor; and panel (d) displays the intensity of the velocity field where the Stokes drift has been included. CHAPTER 4. RESULTS OF PLASTIC DYNAMICS 27 0 50 100 150 200 250 X [m] -5 -4 -3 -2 -1 0 1 Zb [m] 0 0.1 0.2 0.3 0.4 0.5 v [m/s] Basic state m55 alongshore-averaged H [m] Hf Mean H Mean Vel Mean ((a)) ((b)) ((c)) ((d)) Figure 4.3: Data obtained with morfo55 for θ=10o. Panel (a) shows the alongshoreaveraged sea bed level, wave height and alongshore current; panel (b) represents the bottom topographic perturbation along with the velocity field; panel (c) shows the diffusivity factor; and panel (d) displays the intensity of the velocity field where the Stokes drift has been included. 28 Modelling the exchange of microplastics between the inner surf zone and the open sea under different beach configurations The beach profile with the longshore bar at about 75 m from the shoreline along with the wave height distribution in the cross-shore direction can be seen in panels a) of figures 4.1,4.2 and 4.3. In case of figures 4.2 and 4.3 also the longshore current profile is shown displaying two maxima, at the crest of the bar and very close to the shoreline. In panels b), c) and d), it is seen that the alongshore variability is maximum for θ=0oand it decreases by increasing the angle, being almost in-existent for θ=10o. For shore-normal waves, the crescentic shape of the bar and the rip channels are very well developed and the rip current circulation is very visible. Also, there is a strong alongshore rhythmicity in the diffusivity factor, not only on the bar but also near the shoreline (Figure 4.1). It is remarkable that the diffusivity plot bears a high resemblance with the white spots corresponding to wave breaking foam in the typical surf zone video images of crescentic bars in nature. For θ=5o, the crescentic shape is weaker and the rip channels are skewed and oriented downcurrent. The longshore component of the current dominates with only a slight meandering. The diffusivity is almost alongshore uniform but the intensity of the flow has some alongshore variability (Figure 4.2). Finally, for θ=10o, alongshore gradients are hardly seen in any of the variables (Figure 4.3). 4.2. Results for shore-normal wave incidence The dynamics of the plastic concentration is here shown for all scenarios for θ=0o. For each scenario, the 2D field of plastic concentration is displayed at the initial and the final times together with two intermediate times. The velocity field is superposed to the concentration diagrams to help understanding the behaviour. To investigate whether the plastic moves onshore, or offshore or stays in place, the total amount of plastic in some alongshore strips has been monitored over time. These strips are: A1 (from the coastline, x=0, to the breaking line, x=20 m), A2 (from the breaking line to x=80 m, which is the shallowest part of the shoaling zone), A3 (from x=80 to x=140 m, a deeper part of the shoaling zone) and A4 (from x=140 m to the offshore boundary of the integration domain). Also, the total amount of plastic in the entire domain, referred to as At, has been monitored. The scaling of the x-axis has been modified to better appreciate the variations along the cross-shore direction. CHAPTER 4. RESULTS OF PLASTIC DYNAMICS 29 ((a)) ((b)) ((c)) ((d)) 0 0.5 1 1.5 2 t [h] 0 0.2 0.4 0.6 0.8 1 1.2 C mean [g/m2] Plastic concentration per area over time A1 A2 A3 A4 At ((e)) Figure 4.4: Results for θ=0oand scenario 0. (a) numerical solution initial condition; (b) numerical solution after 30 min; (c) numerical solution after 1 hour; (d) numerical solution after 2 hours and (e) averaged plastic concentration in the pre-defined strips over time. Remember that the C scale is arbitrary. 30 Modelling the exchange of microplastics between the inner surf zone and the open sea under different beach configurations ((a)) ((b)) ((c)) ((d)) 0 0.5 1 1.5 2 t [h] 0 0.2 0.4 0.6 0.8 1 1.2 C [g/m2] Plastic concentration per area over time A1 A2 A3 A4 At ((e)) Figure 4.5: Results for θ=0oand scenario 1. (a) numerical solution initial condition; (b) numerical solution after 30 min; (c) numerical solution after 1 hour; (d) numerical solution after 2 hours and (e) averaged plastic concentration in the pre-defined strips over time. Remember that the C scale is arbitrary. CHAPTER 4. RESULTS OF PLASTIC DYNAMICS 31 ((a)) ((b)) ((c)) ((d)) 0 0.5 1 1.5 2 t [h] 0 0.2 0.4 0.6 0.8 1 1.2 C [g/m2] Plastic concentration per area over time A1 A2 A3 A4 At ((e)) Figure 4.6: Results for θ=0oand scenario 2. (a) numerical solution initial condition; (b) numerical solution after 30 min; (c) numerical solution after 1 hour; (d) numerical solution after 2 hours and (e) averaged plastic concentration in the pre-defined strips over time. Remember that the C scale is arbitrary. 38 Modelling the exchange of microplastics between the inner surf zone and the open sea under different beach configurations ((a)) ((b)) ((c)) ((d)) 0 0.5 1 1.5 2 2.5 3 t [h] 0 0.2 0.4 0.6 0.8 1 1.2 C mean [g/m2] Plastic concentration per area over time A1 A2 A3 A4 At ((e)) Figure 4.11: Results for θ=5oand scenario 3. (a) numerical solution initial condition; (b) numerical solution after 1 hour; (c) numerical solution after 2 hour; (d) numerical solution after 3 hours and (e) averaged plastic concentration in the pre-defined strips over time. Remember that the C scale is arbitrary. CHAPTER 4. RESULTS OF PLASTIC DYNAMICS 39 ((a)) ((b)) ((c)) ((d)) 0 0.5 1 1.5 2 t [h] 0 0.2 0.4 0.6 0.8 1 1.2 C mean [g/m2] Plastic concentration per area over time A1 A2 A3 A4 At ((e)) Figure 4.12: Results for θ=10oand scenario 0. (a) numerical solution initial condition; (b) numerical solution after 30 min; (c) numerical solution after 1 hour; (d) numerical solution after 2 hours and (e) averaged plastic concentration in the pre-defined strips over time. Remember that the C scale is arbitrary. 40 Modelling the exchange of microplastics between the inner surf zone and the open sea under different beach configurations ((a)) ((b)) ((c)) ((d)) ((e)) Figure 4.13: Results for θ=10oand scenario 3. (a) numerical solution initial condition; (b) numerical solution after 1 hour; (c) numerical solution after 2 hour; (d) numerical solution after 3 hours and (e) averaged plastic concentration in the pre-defined strips over time. Remember that the C scale is arbitrary. CHAPTER 5. DISCUSSION AND CONCLUSIONS 5.1. Discussion of the results Some of the present results must be regarded with care because of the lack of enough numerical accuracy. Despite having reduced these numerical errors by diminishing the spatial step ∆x(section 3.2), some situations continue to occur where part of these errors appear. For example, a clear consequence is the spurious increase of the total concentration in At (panel e) of the scenario 1 with θ=0o, Figure 4.7). These errors mainly arise when the plastic particles approach the shore where a sharp boundary layer develops. This is specially so for the mixed boundary conditions that become by this reason very demanding numerically speaking. On the other hand, the numerical inaccuracy is partially due to the onset of some numerical diffusivity. There are specific discretization methods for the advection-diffusion governing equations which are fully consistent with the conservative character of the physics (Cushman-Roisin, 2008 Finite-volume discretization, pag. 81). These discretization methods would fix the problem of numerical diffusivity that occurs and are recommended for future work on this problem. Nevertheless, although these limitations of the numerical scheme ma spoil some of our long-term model results on the plastic spreading, the transient behavior for, at least a couple of hours, is well captured. Despite these limited accuracy, the obtained results are reasonable during most of the simulated time and describe some situations that can be compared with reality. In particular, a qualitative experiment was performed by the Group of Nonlinear Fluid Dynamics on the Castelldefels beach (Southwest of Barcelona) with wave conditions relatively similar to those considered here. In this experiment, plastic bottles were throwed outside the surf zone and following their motion towards the dry beach was tracked. After 2 hours most of the bottles had arrived to the dry beach. This time scale is consistent with that obtained in the present simulations. Interestingly, a small portion of the bottles keep on being trapped by rip currents, similarly to what happen with some of our simulations. Also, the tendency of the waves to transport floating material onshore that we have found explains the amount of debris at the dry beach that is commonly seen after a storm. This indicates that the approximation taken for the Stokes drift, despite being very idealistic, captures qualitatively well the overall dynamics. However, in future studies the modelling of the Stokes drift should be improved given its importance for floating plastic dynamics. 41 42 Modelling the exchange of microplastics between the inner surf zone and the open sea under different beach configurations The wave height and period can of course affect the plastic dynamics. The present results have been done for a wave height of 0.8 m and a period of 6 s. In the case that these values were smaller. the incident wave energy would be smaller, too and the time scale of the dynamics would increase. This is due to the fact a smaller (larger) wave energy induce smaller (larger) surf zone currents (Garnier, 2006) and thereby a slower (faster) advecion of plastics. This is also reflected in the Stokes drift formula (2.6), which directly depends on the wave height. Our simulations depend strongly on the wave angle. This is because we use the final morfo55 results of a development process where the sand bar becomes crescentic or not closely depending on the angle. But in nature the situation is more complex, the morphology has some inertia or even hysteresis, so that the present morphology does not always correspond to the present wave angle. Thus, the important beach characteristics that determines plastic dynamics is the present beach morphology more than the present wave angle. CHAPTER 5. DISCUSSION AND CONCLUSIONS 43 5.2. Final conclusions In this project we have addressed the research question of how the presence of a crescentic bar together with the associated circulation affects the spreading of floating microplastics in the surf zone to both sides (dry beach and offshore). To do it, the results of a morphodynamic model called morfo55 have been used to obtain the currents and wave properties associated to the presence (or not) of a crescentic bar. Then, a Matlab code has been built to solve numerically the advection-diffusion equation for the plastic concentration. Both the initial conditions and the imposed boundary conditions have attempted to simulate situations that occur or could occur in reality. The study focuses on the crossshore transfer of plastics so that alongshore periodic boundary conditions are imposed. To test the numerical discretization of the advection-diffusion equation, an analytical solution has been initially studied. Different methods of temporal discretization were tested, from the simplest first order Euler algorithm, to multistep methods such as the Leapfrog, or the Adam-Bashforth methods. The differences and advantages of each one of them have been analyzed, and the Adams-Bashforth method has been chosen. Once we had verified that the discretized solution is almost equal to the analytical result in the test case we have used data obtained from the morfo55. The major problem has been that the results obtained from morfo55 are in a mesh that is not fine enough to simulate the plastic dynamics under some situations. It has been essential to interpolate the values obtained from morfo55 in order to have a denser mesh that could support the required numerical demand. Even so, numerical diffusivity appears due to the fact that our numerical scheme do not achieve sufficient accuracy regarding the conservation of some integrated quantities. Even so, the results obtained are reasonable during most of the simulated time and allow us to see the dynamics of microplastics in different situations. On average, the plastic tends to move onshore due to the existence of the Stokes drift (net onshore water motion on the water surface due to the waves). After a few hours, the concentration has mostly disappeared from the domain. The presence of a crescentic bar also plays an important role. The location of the rip currents, when they exist, and the corresponding alongshore gradients in diffusivity control the plastic distribution that may undergo strong alongshore gradients, with higher concentration on the current ebbs and lower at the headers. The overall effect of the crescentic bar is twofold. On the one hand the cross-shore spreading is accelerated by the presence of a crescentic bar but, on the other hand, the rip currents are able to maintain a certain residual concentration offshore. Thereby, it is possible that a small portion of the microplastic stays offshore due to the presence of a crescentic bar. 44 Modelling the exchange of microplastics between the inner surf zone and the open sea under different beach configurations BIBLIOGRAPHY Garnier R., Nonlinear modelling of surf zone morphodynamical instabilities, Barcelona 2006 Cushman-Roisin B. and Becker JM., Introduction to Geophysical Fluid Dynamics Physical and Numerical Aspects, Belgium 2008 Koelmans, A. All is not lost: deriving a top-down mass budget of plastic at sea 2017 Castelle B. and Coco G., Surf zone flushing on embayed beaches, 2013 Ribas F., Falqu´ es A., de Swart H.E., Dodd N., Garnier R., Calvete D., Understanding coastal morphodynamic patterns from depht-averaged sediment concentration, 2015 Ribas F., On the growth of nearshore sand bars as instabilit processes of equilibrium beach states, Barcelona 2003 45