Contents:
This repository contains all the code for the vertical wind mixing experiments described in Onink et al. (2022): "Empirical Lagrangian parametrization for wind-driven mixing of buoyant particulates at the ocean surface". Using the parcels (Probably A Really Computationally Efficient Lagrangian Simulator) package, a 1D model has been developed to model the vertical concentration profile of buoyant particles in the ocean surface mixed layer depending on the surface wind stress, the particle rise velocities, and the parametrization of vertical turbulent mixing.
This code was written by Victor Onink, with Charlotte Laufkoetter helping to improve the code performance.
To model stochastic particle transport, we consider two stochastic transport formulations:
- Markov - 0 (M-0): In this formulation, diffusion is implemented as a random particle displacement, and total stochastic transport is determined by this random walk together with the particle rise velocity. We use the M-0 formulation described by Grawe et al. (2012).
- Markov - 1 (M-1): In this formulation, a degree of autocorrelation set by parameter alpha is assumed in the particle velocity perturbations, and the total stochastic transport is determined by the turbulent velocity perturbation and the particle rise velocity. We follow the M-1 formulation described by Mofakham & Ahmadi (2020) In both cases, the stochastic transport is discretized using an Euler-Maruyama formulation (Maruyama, 1955)
Two approaches have been taken for the vertical diffusion profile, and in the code the diffusion parametrizations are referred as SWB and KPP diffusion
- SWB: The Surface Wave Breaking parametrization is a combination of the wind-driven surface turbulent mixing parametrization from Kukulka et al. (2012) and the work of Poulain (2020) regarding the depth dependence of the eddy viscosity.
- KPP: The K-Profile Parametrization (KPP) approach was first published by Large et al. (1994) and is commonly used in large-scale ocean models. For the 1D case, we consider a simpler approach that neglects convective fluxes presented by Boufadel et al. (2020).
The model domain is [0, 100] m, where the vertical axis z is negative upwards with z=0 indicating the ocean surface (although in the figures and the paper we follow the convention that z is positive upward with z=0 at the ocean surface). We have tested a number of boundary conditions (BCs) for the ocean surface:
- Ceiling BC: If a particle crosses the ocean surface (z=0), the particle depth is set to depth = 0.
- Reflect BC: A particle is reflected off of the ocean surface, where for a given overshoot dz of the ocean surface, the final particle position is given as |dz|.
- Mixed layer boundary condition: If the particle is within a distance L of the ocean surface, the particle particle depth is set at R x L, where R is a random number between 0 and 1. The mixing depth L is set based on the maximum depth at which a particle displacement for a timestep would still be able to cross the ocean surface.
- Reduce dt: If a particle is to cross the ocean surface, then the time step halved if the timestep is greater than a minimum timestep. If so, the particle position is recalculated using the reduced timestep. If the halved timestep is less than the minimum timestep, a reflecting boundary condition is applied.
If a particle were to cross the lower domain boundary, then we apply a reflecting boundary condition. However, since the maximum depth z=100 m is significantly deeper than the Mixed Layer Depth (MLD), particles in practice never reached the bottom of the model domain.
The particle rise velocities are assigned by the user in the src/main.py file. The corresponding spherical particle size
for a given rise velocity can be computed based on Enders et al. (2015)
In order to verify the modelled distributions, field data of vertical plastic concentration profiles have been collected from a number of sources:
- Kukulka et al. (2012): Data shared by Tobias Kukulka gave concentration depth profiles measured using trawls at a number of stations in the North Atlantic. The shared data also contained the surface wind speed and the mixed layer depth (MLD) determined from CTD using the de Boyer Montegut et al. (2004) approach.
- Pieper et al. (2020): Depth profiles measured during the PE442 cruise (August 2018) on the RV Pelagia between Terceira (the Azores) and Catania (Sicily). Unlike the other data, the concentrations were measured from Niskin bottles, and the data sampled much greater depths (up to 1150m) below the surface, but only limited data close to the surface. Also available was the wind speed at the time of the measurements and CTD profiles for almost all statons.
- Kooi et al., 2016: This is depth profiles measured for the first 5m of the water column, where the measurements were taken with a multi-depth manta trawl. Data is available for a range of wind conditions, and for each station CTD data is available for the computation of the MLD. The data is publically available here
- Amaral-Zettler (unpublished data): This is data collected during the PE448 cruise on the RV Pelagia in the South Atlantic (January 2019) using a multi-stage net for sub-surface measurements and a manta trawl for the surface measurements. Wind speed data is available, as well as CTD data for determining the MLD.
- Egger et al. (2020): This is data collected at various depths with a multi-stage neuston net in the North Pacific. Unlike the data in the paper, this data was provided by Matthias Egger without the depth correction that had been applied in their analysis.
With the exception of the Kooi et al. (2016) data, these data sets are not currently publically available and we recommend contacting the corresponding authors of the studies if one is interested in acquiring this data.
All the commands for running the model go via the main.py file within the src directory. At the beginning of the file, you can set the wind speeds w_10, the rise velocities w_rise and M-1 memory terms alpha for which you want to run the simulation (alpha values are not used in running a M-0 simulation, but alpha must still be defined). In the case of KPP diffusion, you can also set the Langmuir circulation enhancement factor theta, which is a multiplicative factor that increases the amount of mixing to account for Langmuir circulation driven mixing. All these run parameters are written as lists, since the code will run through all elements of the list to carry out the requested simulations.
The variable diffusion sets the diffusion profile that will be used in the simulation, and must be set either as 'KPP' or 'SWB'. The boundary condition is set by boundary, where the options are either 'Ceiling', 'Reflect', 'Mixed' or 'Reduce_dt' (the BC condition used for the results in the paper is typically 'Ceiling' unless otherwise stated). If you want to run a M-1 simulation, then you add '_Markov' to the diffusion string (e.g. boundary='Ceiling_Markov' would run a M-1 simulation with the Ceiling BC). Otherwise, the default is to run a M-0 simulation.
A number of run parameters (e.g. integration and model output timesteps, physical constants, particle densities, the MLD) are defined in the settings.py file. These parameters can not be changed from the main.py file.
Important: Prior to carrying out any simulation on your local computer, please update the path root_direc in settings.py to reflect where all the codes are running on your local machine.
Running the main.py file will run the simulations. First, the code will run field_data_processing(), which will check if the standardized field data files exist and print some basic statistical properties of the field data. If these files are not present, it will create these standardized files. Unless you are in possession of the original field data files, please comment out field_data_processing() prior to running any simulations.
Next, for all the wind, rise velocity and alpha values in w_10, w_rise and alpha, the function parcels_simulations(wind=wind, rise=rise, alpha=alpha_val) is called. Within parcels_simulations(), if one wishes to overwrite any prior model output that exists for the same set of run parameters, set conduct = True in utils.remove_file(conduct=True, file_name=concentration_file). Next, if an output file doesn't already exist the parcels simulation for the given parameters will be carried out, and the parcels output will be converted to a non-normalized concentration profile by analysis.depth_concentration(). Once the concentration profile is calculated, we delete the original parcels output to save storage, but if you wish to keep the original parcels file for further analysis set remove=False.
Finally, the function plotting() calls all plotting functions that are set within plotting(). The visualization functions are all contained within the visualization directory, and a brief explanation of each plotting function is given either within the plotting() function or otherwise in the function documentation of the respective visualization functions.
All these files and directories can be found within the src directory.
main.py: This is the main file used to define and run the simulationssettings.py: This is a file used to set a number of constants, dictionaries and directory paths that are used throughout the functions in the directory.parcels_simulation_functions.py: This file contains all the functions to run the parcels simulations.field_data.py: This contains all the functions to convert the original field data files into one standardized format.analysis: This directory contains all the analysis functions, split into three files:function_concentration.py: function for converting the parcels model output to concentration profiles.function_RMSE.py: all analysis functions related to calculating RMSE values between field data and concentration profilesfunction_statistics.py: analysis functions related to calculating statistics of the field data
utils: This directory contains a number of utility functions, defined as useful functions that are used at multiple points through the model code. The functions are split into three files:utils_filenames.py: These functions are used to name the parcels output, the concentration profile output and the standardized field datautils_files.py: These functions are used for a number of tasks involving files, such as checking if a file exists and saving/loading data using thepicklemodule.utils_physics.py: These functions are used to determine physical quantities, such as diffusion profiles, drag coefficients, surface wind stress, etc.
visualization: This directory contains all the files to visualize the data and create figures. Each file contains the necessary code to create one specific figure described in each function documentation, so please refer there for descriptions of what figure each function creates. The fileutils_visualization.pycontains a number of utility functions for creating the figures, such as functions for creating standardizedmatplotlib.pyplotaxes, determining linestyles and loading concentration profile data such that it can be plotted.
In order to run the model, a number of packages need to be installed. I've listed the versions that I used to run all my simulations, and in most cases its probably fine to run with newer package versions. However, this is not guaranteed to always be the case.
- parcels: v2.2.1 or higher. For the installation requirements, please consult the section Installing Parcels at oceanparcels.org. I did all my simulations with v2.2.1, and generally then the code should also work for more recent versions.
- matplotlib: v3.3.4
- progressbar: v3.37.1
- numpy: v1.19.1
- seabird: v0.11.5: This was only used to analysis CTD data to standardize the field data, so if you are not running this part of the code you don't require this package
- scipy: v1.5.2
- pandas: v1.1.3
- netCDF4: 1.5.4