Implications of Reynolds Averaging for Reactive Tracers in Turbulent Flows

1. Introduction

The ocean exhibits motions over an amazing range of scales, from the 1000 km scales of ocean gyres, to sub-centimeter scales on which viscous dissipation removes energy from the system. Modern oceanography relies on field measurements, remote sensing, and increasingly, numerical simulation of ocean processes. All of these have their own established methodology, yet all fail to provide a complete description of the ocean. However, the influence of physical processes on biogeochemical ocean dynamics is an active area of research in which numerical models play an important role as investigative and predictive tools[1,2,3,4,5]. The performance of numerical models is limited by the underlying numerical method and the model resolution. This is in turn limited by the availability of computational resources. Ocean modellers must thus make choices as to which processes to resolve, which to parameterize, and at times, even which to ignore.

The Regional Ocean Modelling System (ROMS) is a popular tool for conducting regional (on the order of 10s of kilometres) simulations due to its flexibility and stable numerical methodology. It has a topography-following coordinate system in the vertical, an implicit free-surface, and can be run as a coupled physical-biogeochemical model[6]. ROMS is a split-explicit model, meaning fast and slow dynamics are split into sub-problems. If non-dissipative time-stepping is used for both, the model becomes unstable[7,8]. Time-stepping algorithms with forward-backward feedback are used to enhance the internal stability of the model[9]. In addition to the numerical diffusivity introduced by the time–stepping method, ROMS treats horizontal and vertical diffusion differently, and has a number of options for parameterizing vertical mixing[7,10] all of which are, by nature, dissipative[11]. While not strictly necessary, the majority of regional modeling studies invoke the hydrostatic approximation, meaning some processes (e.g. density overturns) are inherently excluded from consideration. Example use cases for ROMS include the Atlantic Canada Model ($O(10$km) horizontal resolution and 30 vertical levels)[12], the Coupled-Ocean-Atmosphere-Wave-Sediment-Transport modeling system (4.7-9.0km horizontal resolution, 50 vertical layers spanning 8km)[13], and studies of storm induced circulation and hydrography changes in tropical regions (6km resolution)[14]. The diversity of studies that use ROMS indicates the importance of model flexibility and adaptability.

ROMS, and similar models, effectively solve the Reynolds Averaged Navier-Stokes (RANS) equations. Reynolds averaging is a commonly applied procedure in the study of turbulent flows [15]. It consists of decomposing a physical quantity into its ensemble mean and fluctuating parts. The procedure will be illustrated in Section 2.1.1. The appeal of RANS simulations comes from the reduced computational cost of simulating only the mean state of a system as opposed to running an ensemble of simulations. However, ensemble averaging the Navier-Stokes momentum equation introduces the Reynolds stress tensor which poses a closure problem[15]. Despite the challenge of parameterizing the Reynolds stresses, RANS is a popular method for simulating naturally occurring turbulent flows, which are ubiquitous in the oceans and atmosphere and known to cause strong mixing[16,17]. Turbulent flows are characterized by large, irregular variations in space and time, making their details unpredictable [18]. These variations often manifest in the form of vortices or filaments. For a discussion of several classes of models, as well as popular modelling techniques and closures when using RANS for turbulence modelling see [19]. For a review of the uncertainty in RANS models see [20].

Eddy viscosity models are the most common choice for parameterizing Reynolds stress terms, largely owing to their relative simplicity. These models assume that the Reynolds stress tensor is aligned with the mean strain tensor allowing the Reynolds stress term to be related to the ensemble mean velocity field[15,18]. In practice, the Reynolds stress terms in the RANS momentum equation are replaced by simply introducing additional viscosity into the system[15]. These models assume that the local gradient of the mean velocity field dictates the direction of flux. However, if the advection in a turbulent flow is strong enough, the large eddies that emerge may dominate and cause the flux to be oriented against the gradient[18], making the eddy viscosity assumption invalid. In particular, there are cases in which the eddy viscosity becomes negative and the model breaks down because anti-diffusion is mathematically ill-posed[21].

Analogously, the Reynolds fluxes that emerge from Reynolds averaging scalar transport equations are often modelled by eddy diffusivities. Eddy diffusivities have been estimated for passive tracers with satellite observations and found to be strongly space and time dependent[22]. In the upper 2000m of the oceans, horizontal eddy diffusivities on a 300km scale range between ${10}^2$m${}^2$/s and ${10}^4$m${}^2$/s[23,24,25]. The choice of expression used to calculate eddy diffusivity inevitably affects the estimate. If one were to consider the suppression of mixing by mean currents, the spatial distribution and magnitude of the eddy diffusivities changes dramatically[26] as compared to estimations that simply use mixing length theory[27]. Evidence has suggested that in stratified fluids, the ocean being a prime example, eddy diffusivity should be modelled as a decreasing function of buoyancy frequency[28].

Eddy mixing is also strongly dependent on the tracer being considered[29]. Active tracers may require counter-gradient fluxes as opposed to down gradient fluxes resulting in negative eddy diffusivities[30]. For instance, double diffusion displays counter-gradient density flux while turbulence transports density upwards[31]. Though some methods for determining eddy diffusivities and transport characteristics have shown indistinguishable performance when used for passive or active tracers[32], this cannot be the case for reactive tracers, the dynamics of which are inevitably more complicated.

The effect of turbulence on the reaction dynamics of nutrient phytoplankton zooplankton (NPZ) models has been considered from several different lenses[33,34,35]. NPZ models are almost always expressed in terms of advection-diffusion-reaction (ADR) equations. Complications arise from the impact of turbulence on the spatial distribution of a reactive tracer possibly affecting the reaction rate[35]. Further, eddy reaction terms, defined as covariances of tracer fluctuations, are introduced when the reaction function is Reynolds averaged. Eddy reactions are inevitably strongly reaction dependent and have been shown to be sensitive to seasonality, location, and spatial scales for some NPZ models[36].

In addition to being a popular tool for examining oceanic ecosystems, RANS has been used to model reactive tracers in the contexts of water disinfection[37,38,39], the spread of air pollutants[40], and combustion[41]. In this paper, we take a step back and consider a simple system in which a reactive tracer grows according to Fisher’s equation within a Rayleigh-Taylor (RT) flow. So called “toy" flows, for which velocities are specified instead of solved for, have proven to be a useful tool for conducting initial investigations into the effect of turbulence on phytoplankton[34]. Velocities are often chosen so that the flow field is composed of convective cells of multiple scales that oscillate in space [34,42]. The cells are intended to represent flows due to thermal convection and the range of scales is introduced to mimic turbulence. They are thought to be especially useful when considering under-ice convection as ice cover eliminates all wind-driven turbulent mixing[42,43]. However, this method fails to capture the cascade of energy between scales in turbulent flows nor does it account for the dissipation of turbulent kinetic energy, making it a poor representation of ocean dynamics.

We choose Rayleigh-Taylor instability-induced turbulent mixing as a non-trivial, yet well studied illustrative example. From Figure 2 in [44], it is apparent that RT flows are not cell-like and evolve into turbulent states. For a discussion of several turbulence model choices for RT instabilities, see [45,46]. RT instabilities have some memory of the spatial structure of their initial conditions when we consider large-scale effects, but the smaller-scale internal structure makes the mixing unpredictable[47]. Two dimensional simulations fail to capture the internal structure of the flow and thus cannot describe the turbulent mixing that would be seen in three dimensions[44,47,48]. For the purpose of constructing an illustrative example, however, two dimensional simulations will suffice.

The simplicity of the numerical experiment reported herein allows us to illustrate the differences between passive and reactive tracer fields evolving under turbulence. In particular, we illustrate that the evolution of a reactive tracer is poorly captured by the ensemble mean and that eddy diffusivity parameterizations are not appropriate for our simulations. The remainder of this paper is structured as follows. In Section 2 we present the governing equations and demonstrate how Reynolds averaging can be applied to ADR equations. Results from an ensemble of 50 simulations will be presented in Section 2.2 to illustrate differences in the effects of turbulent mixing on a passive tracer versus on a reactive tracer growing according to Fisher’s equation. The limitations of using RANS simulations to model reactive tracers will be discussed in Section 4.

2. Methods

2.1. Governing Equations

Consider the following set of partial differential equations

$$\begin{array}{cc}\hfill & {\overrightarrow{u}}_t+\nabla ·\left({\overrightarrow{u}}^T\overrightarrow{u}\right)=-\frac{1}{{\rho }_0}\nabla p+\frac{{\rho }_E}{{\rho }_0}\overrightarrow{g}+\nu {\nabla }^2\overrightarrow{u},\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill & {\theta }_t+\nabla ·\left(\overrightarrow{u}\theta \right)=\kappa {\nabla }^2\theta +F\left(\theta \right).\hfill \end{array}$$

(2)

We report the excess density, defined as a departure from the reference density

$$\begin{array}{c}\hfill {\rho }_E=\rho (\overrightarrow{x},t)-{\rho }_0,\end{array}$$

(3)

SPINS (Spectral Parallel Incompressible Navier–Stokes Solver)[49] will be the model used herein. The excess density is evolved through an ADR equation like that for $\theta$, but without a source term.

In equation 1, the conservation of mass, $\nabla ·\overrightarrow{u}=0$ has been used to rewrite $\overrightarrow{u}·\nabla \overrightarrow{u}$ as $\nabla ·\left({\overrightarrow{u}}^T\overrightarrow{u}\right)$. The second equation is an ADR equation describing the evolution of a scalar tracer undergoing a reaction given by the function $F\left(\theta \right)$. A reactive tracer has a nonzero reaction term, but we assume that it does not affect the properties of the fluid. Reactive tracers are an important component in biophysical models describing ecosystem dynamics which generally include components evolved by growth and interaction terms [2,3,4,5,50].

2.1.1. Reynolds Averaging

We will use overbar notation to denote Reynolds averaging. We decompose the velocity field into a slowly varying mean part, $\overrightarrow{U}$, and a rapidly fluctuating part, ${\overrightarrow{u}}^{\prime }$, so that $\overrightarrow{u}=\overrightarrow{U}+{\overrightarrow{u}}^{\prime }$ with $\overline{\overrightarrow{U}}=\overrightarrow{U}$ and $\overline{{\overrightarrow{u}}^{\prime }}=0$. Similarly for the pressure we have $p=P+p^{\prime }$. The Reynolds averaged momentum equation is

$$\begin{array}{cc}\hfill {\overrightarrow{U}}_t+\nabla ·\left({\overrightarrow{U}}^T\overrightarrow{U}\right)& =-\frac{1}{{\rho }_0}\nabla P+\nu {\nabla }^2\overrightarrow{U}-\nabla ·\overline{{\overrightarrow{u}}^{\prime T}{\overrightarrow{u}}^{\prime }}.\hfill \end{array}$$

(4)

The Reynolds stress tensor describes the effects of the fluctuations on the mean flow and is defined as $R_{ij}=\overline{u_i^{\prime }{u}_j^{\prime }}$.

The procedure can be applied to equation by decomposing the tracer field in the same way; $\theta =\Theta +{\theta }^{\prime }$. The Reynolds averaged ADR equation is

$$\begin{array}{cc}\hfill {\Theta }_t+\nabla ·\left(\overrightarrow{U}\Theta \right)& =\kappa {\nabla }^2\Theta +\overline{F(\Theta +{\theta }^{\prime })}-\nabla ·\overline{{\overrightarrow{u}}^{\prime }{\theta }^{\prime }},\hfill \end{array}$$

(5)

where $\overline{{\overrightarrow{u}}^{\prime }{\theta }^{\prime }}$ is the eddy flux representing the transport of tracer fluctuations by fluctuations in the velocity field. The treatment of the $\overline{F(\Theta +{\theta }^{\prime })}$ term will be considered in Section 2.1.2.

2.1.2. Reynolds Averaging $F\left(\theta \right)$

Fisher’s equation, $F\left(\theta \right)=\theta (1-\theta )$, will be used as an illustrative example, but the procedure presented herein can be used to treat any reaction term expressed as (or approximated by) a polynomial. We want to rewrite $\overline{{\theta }^n}$ in a way that allows us to separate the effects due to the mean tracer field, those due to fluctuations, and those due to cross-interactions of the two. By substituting $\theta =\Theta +{\theta }^{\prime }$, the reactions can be binned as follows

$$\begin{array}{cc}\hfill \overline{{\theta }^n}& =\overline{{(\Theta +{\theta }^{\prime })}^n}\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill & =\sum _{{\alpha }_1+{\alpha }_2=n;{\alpha }_1,{\alpha }_2\ge 0}\left(\genfrac{}{}{0pt}{}{n}{{\alpha }_1,{\alpha }_2}\right)\overline{{\Theta }^{{\alpha }_1}{\left({\theta }^{\prime }\right)}^{{\alpha }_2}}\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill & ={\Theta }^n+\underset{I_n}{\underbrace{{\sum }_{{\alpha }_1+{\alpha }_2=n;{\alpha }_1>0{\alpha }_2>1}\left(\genfrac{}{}{0pt}{}{n}{{\alpha }_1,{\alpha }_2}\right){\Theta }^{{\alpha }_1}\overline{{\left({\theta }^{\prime }\right)}^{{\alpha }_2}}}}+\underset{S_n}{\underbrace{\overline{{\left({\theta }^{\prime }\right)}^n}}}.\hfill \end{array}$$

(8)

The scalar quantity $I_n$ represents the interactions between the mean tracer field and the fluctuations. The correlations between the fluctuations are quantified through $S_n$. In particular, $S_2$ is the scalar variance of the tracer. Note that $I_n=0$ for $n<3$ and $S_n=0$ for $n<2$ since $\overline{{\theta }^{\prime }}=0$ and a constant function ($n=0$) does not have a fluctuating part.

The averaged reaction term can be written as

$$\begin{array}{cc}\hfill \overline{F\left(\theta \right)}& =\sum _{n=0}^N\frac{F^{\left(n\right)}\left(0\right)}{n!}{\Theta }^n+\sum _{n=3}^N\frac{F^{\left(n\right)}\left(0\right)}{n!}{I}_n+\sum _{n=2}^N\frac{F^{\left(n\right)}\left(0\right)}{n!}{S}_n\hfill \end{array}$$

(9)

and the reaction term for the fluctuation equation as

$$\begin{array}{c}\hfill F{\left(\theta \right)}^{\prime }=F\left(\theta \right)-\overline{F\left(\theta \right)}=\sum _{n=0}^N\frac{F^{\left(n\right)}\left(0\right)}{n!}({\theta }^n-{\Theta }^n)-\sum _{n=3}^N\frac{F^{\left(n\right)}\left(0\right)}{n!}{I}_n-\sum _{n=2}^N\frac{F^{\left(n\right)}\left(0\right)}{n!}{S}_n.\end{array}$$

(10)

For Fisher’s equation, the mean and fluctuation reaction terms can be expressed as follows

$$\begin{array}{cc}\hfill \overline{F\left(\theta \right)}& =\Theta (1-\Theta )-\overline{{\left({\theta }^{\prime }\right)}^2},\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill F{\left(\theta \right)}^{\prime }& ={\theta }^{\prime }(1-{\theta }^{\prime })+\overline{{\left({\theta }^{\prime }\right)}^2}-2\Theta {\theta }^{\prime }.\hfill \end{array}$$

(12)

The scalar variance, $S_2=\overline{{\left({\theta }^{\prime }\right)}^2}$, is strictly non-negative and therefore acts as a sink in the mean equation and a source in the fluctuation equation.

If the tracer in question must have non-negative values everywhere in its domain (eg. concentration, population), then its mean must also be non-negative. However, there is no such restriction on the fluctuations, the signs of which depend on the spatial distribution of the mean relative to that of each ensemble member. Fluctuations will be negative (positive) anywhere an ensemble member has less (more) tracer than the mean. At these locations, the sign of $-2\Theta {\theta }^{\prime }$ is opposite that of ${\theta }^{\prime }$ and so this term drives the fluctuations towards 0. However, it must compete with $F\left({\theta }^{\prime }\right)$ which drives ${\theta }^{\prime }\to 1$ if ${\theta }^{\prime }>0$ and ${\theta }^{\prime }\to -\infty$ if ${\theta }^{\prime }<0$ as well as $S_2$ which drives ${\theta }^{\prime }\to \infty$.

The ADR equations for the mean and fluctuations are

$$\begin{array}{cc}\hfill {\Theta }_t+\nabla ·\left(\overrightarrow{U}\Theta \right)& =\kappa {\nabla }^2\Theta +\Theta (1-\Theta )-\overline{{\left({\theta }^{\prime }\right)}^2}-\nabla ·\overline{{\overrightarrow{u}}^{\prime }{\theta }^{\prime }},\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill {\theta }_t^{\prime }+\nabla ·(\overrightarrow{U}{\theta }^{\prime }+{\overrightarrow{u}}^{\prime }\Theta )+\nabla ·({\overrightarrow{u}}^{\prime }{\theta }^{\prime }-\overline{{\overrightarrow{u}}^{\prime }{\theta }^{\prime }})& =\kappa {\nabla }^2{\theta }^{\prime }+{\theta }^{\prime }(1-{\theta }^{\prime })+\overline{{\left({\theta }^{\prime }\right)}^2}-2\Theta {\theta }^{\prime }.\hfill \end{array}$$

(14)

We note that often the fluctuation equation is multiplied by ${\theta }^{\prime }$ and averaged to produce an equation for the scalar variance[51]. However, with the reaction term, this would result in triple correlations, posing another closure problem. Since our numerical simulations provide sufficient data and resolution to treat the fluctuations directly, we consider an equation for the fluctuations.

2.2. Numerical Simulation Ensemble Setup

An ensemble of 50 simulations was performed with SPINS which allows for direct numerical simulations (DNS) of the incompressible Navier-Stokes equations. The model uses third-order accurate, variable time step Adams-Bashforth time-stepping and has spectral accuracy in space. When compared to lower order finite difference or volume schemes with the same resolution, numerical errors are much smaller[52]. Further, the numerical dispersion and diffusion errors often found in classical lower order methods are absent from spectral methods. The SPINS code has been extensively validated against linear theory (e.g., growth rate for shear instability and internal wave generation[49]), existing numerical methods (e.g., in the context of dipole–wall collisions[49]), exact solutions (e.g., for the propagation of exact internal solitary waves[53]), and experimental data (e.g., cross-boundary layer transport due to shoaling mode 2 internal waves[54]).

Simulation parameters and dimensionless numbers are given in Table 1. All simulations were run in 2D with a 2.048m wide and 0.512m tall domain. We used 4096 grid points in the horizontal and 1024 in the vertical to have 0.5mm resolution in both directions. Periodic boundary conditions were implemented in x and free slip boundary conditions were implemented in z.

To induce the formation of Rayleigh-Taylor (RT) instabilities, each ensemble member is initialized with an unstable density stratification of the form

$$\begin{array}{cc}\hfill \rho (\overrightarrow{x},t=0)& =-5\times {10}^{-4}\left(1-tanh\left(\frac{z-0.256-\epsilon \left(x\right)}{0.02}\right)\right),\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill \epsilon \left(x\right)& =2.5\times {10}^{-3}\sum _{i=1}^4{Y}_{1i}sin\left(\frac{16\pi (Y_{2i}+1)}{2.048}x+2\pi Y_{3i}\right),\hfill \end{array}$$

(16)

where $Y_{ji}\sim U(0,1)$ is a random variable chosen from a uniform distribution between 0 and 1. The distribution is seeded by the start time of each simulation. The only difference between ensemble members is the perturbation, $\epsilon$, applied to the initial density stratification. Panel (a) of Figure 1 shows an example stratification from a subdomain of one of the simulations.

Estimates of the dimensionless numbers were calculated using the domain height as the characteristic length scale and the maximum root-mean-square vertical velocity of a sample ensemble member as the velocity scale. Values are given in Table 1. The Reynolds number, $Re=UL/\nu$, represents the ratio of the inertial to viscous forces. Our simulations have $Re\approx 1.1264\times {10}^4$ indicating the flow becomes turbulent. The Schmidt number, $Sc=\nu /\kappa =7.1429$, tells us the momentum diffusivity is greater than that of the tracers. The Péclet number describes the ratio of advective to diffusive transport and is given by $Pe=Sc·Re\approx 8.0457\times {10}^4$.

In each of the simulations, the reaction term is

$$\begin{array}{c}\hfill F\left(\theta \right)=0.1·\theta (1-\theta )·tanh(\theta /0.001),\end{array}$$

(17)

where the multiplication by the hyperbolic tangent ensures that any oscillations around the unstable equilibrium $\theta =0$ are suppressed. The reaction timescale parameter $\alpha =0.1$ was chosen so that the Damköhler number, which we define as $Da=\frac{\alpha L}{U}\approx 2.3273$, is $O\left(1\right)$. The order of $Da$ was chosen to balance the reaction and advective timescales for the purposes of illustration. A fast reaction (eg. $\alpha =1$) has $\theta \to 1$ too quickly for structures created from turbulent mixing to be observed. Slow reactions (eg. $\alpha =0.01$) allow for the creation of sharp fronts which can be difficult to resolve.

The initial distributions of both the reactive and passive tracers are given by

$$\begin{array}{cc}\hfill {\theta }_0& =0.5\left(1+tanh\left(\frac{z-0.48}{0.06}\right)\right),\hfill \end{array}$$

(18)

so that the tracers lie in a band along the top of the domain as shown in panel (b) of Figure 1. Each simulation evolves both tracers with the same velocity fields and the tracers do not interact. The simulations are stopped after 40s when the system has reached a well-mixed state.

Figure 1. (a) Subdomain with sample initial density stratification perturbed by $\epsilon$ (equation 15). (b) Initial tracer concentration as a function of depth.

Figure 1. (a) Subdomain with sample initial density stratification perturbed by $\epsilon$ (equation 15). (b) Initial tracer concentration as a function of depth.

3. Results

3.1. Tracer Mean and Scalar Variance

The reactive and passive tracers start in a band along the top of the domain with concentration decreasing downwards. Initially, the system is at rest to mimic the accumulation of a tracer in calm waters. If the system started out turbulent, we would expect the tracer to be well-mixed and since turbulent mixing is inherently unpredictable, it would be challenging to draw clear comparisons between the spatial distributions of the reactive and passive tracers. The initial density stratification is centred at $z=0.256$m, about 0.16m below the bottom of the tracer band, meaning there is a delay between the start of the simulation and the time at which the RT instabilities start advecting the tracers. During this time, the tracers begin diffusing downwards and the reactive tracer grows according to the function in equation 17. The reaction both increases the concentration of the tracer and accelerates the downwards spread by enhancing Fickian diffusion. Therefore, the reactive tracer front lies lower in the domain allowing the RT plumes to reach it earlier than the passive tracer.

Figure 2 (stills taken from supplementary movies S1 and S2) shows the (a) mean and (b) scalar variance of the reactive tracer at $t=26$s. When comparing these quantities with those of the passive tracer, panels (b) and (d), it is clear that the mean concentration of the reactive tracer is nearly twice that of the passive tracer and extends further downwards. In panels (a) and (b) of Figure 3 the bulk values of the means of the reactive (blue) and passive (orange) tracers are plotted as functions of time. The mean reactive tracer concentration increases by nearly 400% while changes in the mean passive tracer concentration are $O\left({10}^{-5}\right)$ and likely due to the filter used in the numerical method. Note that in all Figures presented herein, shaded areas extend two standard deviations away from the corresponding ensemble means.

Figure 2. Mean and scalar variance of the reactive, (a) and (c), and passive, (b) and (d) tracers at $t=26$s. Panels (a) and (c) are from supplementary movie S1 and panels (b) and (d) are from S2.

Figure 2. Mean and scalar variance of the reactive, (a) and (c), and passive, (b) and (d) tracers at $t=26$s. Panels (a) and (c) are from supplementary movie S1 and panels (b) and (d) are from S2.

Figure 3. (a) Bulk mean tracer concentrations and (c) bulk scalar variances for the reactive (blue) and passive (orange) tracers. Changes in the bulk mean passive tracer concentration,(b), are $O(1e-5)$. Shaded areas extend two standard deviations away from the mean. Panels (d) and (e) depict the scalar variance integrated with respect to x for the reactive and passive tracers respectively.

Figure 3. (a) Bulk mean tracer concentrations and (c) bulk scalar variances for the reactive (blue) and passive (orange) tracers. Changes in the bulk mean passive tracer concentration,(b), are $O(1e-5)$. Shaded areas extend two standard deviations away from the mean. Panels (d) and (e) depict the scalar variance integrated with respect to x for the reactive and passive tracers respectively.

The ensemble is designed for the RT plumes in each simulation to have different spatial distributions. In each simulation, the RT bubbles (plumes of light fluid rising into heavy fluid with mushroom-shaped heads) push up against the tracer front and the tracer concentration is displaced upwards and outwards around the bubble. Meanwhile, spikes entrain the tracers into the lower parts of the domain causing the streaks observed in Figure 2. Shear causes vorticial structures to form on the undersides of the mushroom-shaped plume heads. These vortices are the first manifestation of turbulence in the flow and begin mixing the tracer that is displaced around the RT bubbles.

However, this behaviour is not obvious in the supplementary movies showing the evolution of the mean and scalar variance fields for each of the tracers. When taking the mean of an ensemble, the structures present in the flows of the individual simulations are averaged out. So, the mean tracer fields do not show any of the coherent structures that characterize turbulent mixing. Instead, the structures are captured by the fluctuations which are defined for each ensemble member as the difference between itself and the ensemble mean. The tracer fields in each of the simulations begin to differ from the mean as soon as the RT plumes start to affect their spatial distribution. Panel (c) of Figure 3 indicates a sharp increase in the scalar variance of the reactive tracer (blue) at $t=18$s as the RT bubbles reach the tracer front. It takes longer for the plumes to reach the passive tracer front.

The turbulence in each of the simulations deforms and contorts material volumes of the conserved scalar, in this case our passive tracer, and increases the magnitude of local scalar gradients. This in turn enhances molecular diffusion[51]. When this same process acts on the reactive tracer, the reaction occurring within the material volume complicates things further. For growth governed by Fisher’s equation, the concentration in each material volume increases towards the stable equilibrium value $\theta =1$ thereby increasing the rate of molecular diffusion. This doubly enhanced diffusion allows the reactive tracer to spread to areas a passive tracer may not be able to reach. Further enhancement results as the process of deformation by the flow is repeated.

Scalar variance acts as a sink in the mean ADR equation for the reactive tracer (equation13) and as a source in the fluctuation equation (equation). The fluctuations are effectively taking concentration from the mean field in addition to growing according to Fisher’s equation. In panels (a) and (c) of Figure 2 the spatial distribution of the scalar variance matches the front of the mean field. Looking at supplementary movie S1, we see that as the mean descends in the domain, its concentration is transferred to the scalar variance which evolves into a band spanning the domain horizontally. Eventually, the system becomes so well-mixed that the differences between simulations become less noticeable and thus the scalar variance decreases, Figure 3 panel (c).

Comparing the evolution of the reactive tracer with that of the passive tracer we see that in the latter case the RT instabilities quickly break the mean structure apart thereby increasing the scalar variance. However, the magnitude of the scalar variance is limited by the initial concentration of the tracer. In panels (d) and (e) of Figure 3, the reactive and passive scalar variances are integrated with respect to x and plotted as a function of depth and time. The dotted white lines indicate the inflection point of the hyperbolic tangent function describing the initial tracer distribution and the dashed line indicates the height of the mean initial density stratification. It is apparent that the reactive tracer’s scalar variance begins to increase lower in the domain as the tracer front has extended downwards. We also notice that it propagates both upwards and downwards. The upwards propagation corresponds to the transfer of concentration from the mean and to the fluctuations while the downwards propagation is uniquely due to RT induced turbulent mixing. The scalar variance of the passive tracer can only propagate downwards and is significantly smaller in magnitude.

3.2. Reynolds Stresses

The bulk absolute values of the Reynolds stresses are depicted as functions of time in panel (a) of Figure 4. The shear stress is much weaker than the other stresses (though still on the same order of magnitude) and we note $R_{22}$ varies the most between ensemble members. Panels (b) though (d) show the normalized Reynolds stresses integrated with respect to x as a function of z and time. The dashed white lines indicate the initial height of the mean density stratification. From panel (d) we see that $R_{22}$ spreads symmetrically from the location of the initial stratification which is in agreement with RT theory. The spatial distribution of $R_{11}$, panel (b), indicates that fluctuations in the horizontal velocity are symmetrical in z, greatest along the sidewalls, and take longer to increase than those in the vertical. The distribution of the shear stress is very different from those of $R_{11}$ and $R_{22}$. There is no symmetry around the initial stratification nor is there a coherent pattern.

Figure 4. (a) Bulk absolute Reynolds stresses $R_{11}$ (purple), $R_{12}$ (pink), and $R_{22}$ (green) with shaded areas extending two standard deviations away from the mean (but not crossing $R_{ij}=0$). The normalized Reynolds stresses (b) $R_{12}$, (c)$R_{11}$, and (d)$R_{22}$ integrated with respect to x. Note that panel (b) has different colorbar limits as it is the only quantity that is not positive definite.

Figure 4. (a) Bulk absolute Reynolds stresses $R_{11}$ (purple), $R_{12}$ (pink), and $R_{22}$ (green) with shaded areas extending two standard deviations away from the mean (but not crossing $R_{ij}=0$). The normalized Reynolds stresses (b) $R_{12}$, (c)$R_{11}$, and (d)$R_{22}$ integrated with respect to x. Note that panel (b) has different colorbar limits as it is the only quantity that is not positive definite.

When using an eddy viscosity parameterization for the Reynolds stresses, we have[18]

$$\begin{array}{c}\hfill \overline{u^{\prime }{u}^{\prime }}={\nu }_e^{11}{U}_x,\overline{w^{\prime }{w}^{\prime }}={\nu }_e^{22}{W}_z,\overline{u^{\prime }{w}^{\prime }}={\nu }_e^{12}(U_z+W_x)/2,\end{array}$$

(19)

and values for the eddy viscosities need to be selected. Since we know the values for the Reynolds stresses and gradients of the mean velocity fields, we can divide through by the gradients to obtain time-dependent values for the eddy viscosities. In Figure 5 the Reynolds stresses and shaded areas are scaled by the maximum values of the ensemble means. The black lines indicate the corresponding normalized right hand side of the equations in 19. Panel (d) shows the time-dependent eddy viscosities resulting from the aforementioned division. The trends captured by the Reynolds stresses are not well represented by the gradients of the mean velocity fields. In particular, while the Reynolds stresses are increasing monotonically, the gradients reach their maxima around $t=30$s and then begin to decrease. The vertical eddy viscosity, ${\nu }_e^{22}$ is the largest and the most variable. This is not entirely surprising for a RT flow as the instability manifests as plumes moving primarily in the vertical. That said, if one were to run a RANS simulation using these eddy viscosity values in their eddy viscosity parameterization, they would find the mean velocity fields to be more resistant to vertical deformations than horizontal ones, which is perhaps the opposite of what one would expect of a RT flow. In fact, in Figure 4, we saw that the vertical and horizontal Reynolds stresses have similar magnitudes (because the perturbations to the initial density stratification cause plumes to propagate diagonally) but this is clearly not captured in the gradients of the mean velocity fields.

Figure 5. Normalized bulk mean Reynolds stresses (a) $R_{11}$, (b)$R_{12}$, and (c)$R_{22}$ with shaded areas extending two standard deviations away from the mean. The corresponding normalized bulk gradient approximations (a) $U_x$, (b)$(U_z+W_z)/2$, and (c) $W_z$ are plotted with black lines. The eddy viscosity values resulting from dividing the bulk Reynolds stresses by their gradient approximations are shown in (d).

Figure 5. Normalized bulk mean Reynolds stresses (a) $R_{11}$, (b)$R_{12}$, and (c)$R_{22}$ with shaded areas extending two standard deviations away from the mean. The corresponding normalized bulk gradient approximations (a) $U_x$, (b)$(U_z+W_z)/2$, and (c) $W_z$ are plotted with black lines. The eddy viscosity values resulting from dividing the bulk Reynolds stresses by their gradient approximations are shown in (d).

3.3. Eddy Fluxes

Eddy fluxes quantify the mean transport of tracer fluctuations by fluctuations in the velocity fields. In Section 3.2 the evolution of the fluctuations in the velocity fields was discussed. We saw that fluctuations in the vertical velocity component spread symmetrically from the location of the initial density stratification while fluctuations in the horizontal one are largely found along the top and bottom boundaries. In Figure 6, we compare the bulk absolute horizontal and vertical eddy fluxes in panels (a) and (b) for the reactive (blue) and passive (orange) tracers. In Section 3.1 we saw that reactive tracer fluctuations are larger than those of the passive tracer. It thus comes as no surprise that the eddy fluxes of the reactive tracer are greater than those of the passive tracer.

Despite the magnitudes of the horizontal, $R^{11}$, and vertical, $R^{22}$, Reynolds stresses being similar, it is clear that the tracer fluctuations are more closely correlated with the vertical velocity fluctuations than the horizontal ones. Around $t=30$s, the heads of the RT plumes collide with the top and bottom boundaries of the domain. At this time, the vertical eddy fluxes reach their maxima and subsequently decrease because the tracers cannot be advected any further in the vertical. The decrease in the vertical eddy flux of the reactive tracer is more dramatic than that corresponding to the passive tracer since fluctuations in the passive tracer barely reach the bottom of the domain and thus are not significantly affected by the bottom boundary.

Figure 6. Bulk absolute value of the (a) horizontal, $\overline{u^{\prime }{\theta }^{\prime }}$, and (b) vertical, $\overline{w^{\prime }{\theta }^{\prime }}$, eddy fluxes of the reactive (blue) and passive (orange) tracers. Vertical eddy fluxes integrated with respect to x for the (c) reactive and (d) passive tracers.

Figure 6. Bulk absolute value of the (a) horizontal, $\overline{u^{\prime }{\theta }^{\prime }}$, and (b) vertical, $\overline{w^{\prime }{\theta }^{\prime }}$, eddy fluxes of the reactive (blue) and passive (orange) tracers. Vertical eddy fluxes integrated with respect to x for the (c) reactive and (d) passive tracers.

In panels (c) and (d) of Figure 6, the vertical eddy fluxes of the reactive and passive tracers are integrated with respect to x and shown as functions of depth and time. The dotted white lines indicate the inflection point of the hyperbolic tangent function describing the initial tracer distribution and the dashed lines indicate the initial height of the mean density stratification. For both tracers, the vertical eddy fluxes are negative because, as shown in the supplementary movies, the fluctuations are travelling downwards and the heights at which the eddy fluxes first become nonzero align with the locations of the tracer fronts. The reactive tracer fluctuations and fluctuations in the vertical velocity field both propagate upwards (Figure 3 and Figure 4) making the vertical eddy flux propagate upwards. This does not necessarily imply that the fluctuations are being transported upwards. Rather, it tells us that the fluctuations develop higher in the domain as they consume concentration from the front of the mean distribution. This is particularly evident in supplementary movie of the reactive tracer fields (S1). These space-time plots further support the conclusion that the decrease seen in the bulk absolute vertical eddy flux of the reactive tracer is more dramatic because fluctuations in the passive tracer barely extend to the bottom of the domain.

When running RANS simulations, eddy diffusivity parameterizations of the form

$$\begin{array}{c}\hfill \overline{u^{\prime }{\theta }^{\prime }}={\kappa }_H\frac{\partial \theta }{\partial x},\overline{w^{\prime }{\theta }^{\prime }}={\kappa }_V\frac{\partial \theta }{\partial z},\end{array}$$

(20)

are often used in place of the eddy fluxes which are unknown. In Section 3.2, we saw that the eddy viscosity parameterization does not capture the same trends as the Reynolds stresses if the eddy viscosities are assumed to be constant. This does not inspire much confidence in the performance of the analogous eddy diffusivity parameterization. In panels (a) and (b) of Figure 7 the normalized bulk absolute eddy fluxes are plotted along with the normalized right hand sides of the equations in 20. The tracer derivatives taken with respect to x capture the general behaviour of the horizontal eddy fluxes but fail to capture the details. The horizontal eddy fluxes start increasing earlier and begin to plateau at late times. The tracer derivatives taken with respect to z increase monotonically and clearly do not capture the behaviour of the vertical eddy fluxes which peak around $t=30$s and then begin to decrease. As previously discussed, the vertical eddy fluxes peak when tracer fluctuations reach the bottom of the domain. Eddy diffusivity parameterizations consider the gradient of the tracer mean which only begins to reach the bottom boundary when the simulations are stopped and the tracer is well-mixed. In panels (c) and (d), we plot the corresponding vertical and horizontal eddy diffusivities obtained from dividing through by the tracer derivatives. We immediately note that the horizontal eddy diffusivities are two orders of magnitude smaller than the vertical ones. The variability in the values brings the use of constant eddy diffusivities into question.

Figure 7. Normalized bulk absolute values of the (a) horizontal, $\overline{u^{\prime }{\theta }^{\prime }}$, and (b) vertical, $\overline{w^{\prime }{\theta }^{\prime }}$, eddy fluxes of the reactive (blue) and passive (orange) tracers. The bulk absolute values of the derivatives of the tracers with respect to (a) x and (b) z are shown in purple and yellow for the reactive and passive tracers respectively. The (c) horizontal, $k_H=\overline{u^{\prime }{\theta }^{\prime }}/{\theta }_x$, and (d) vertical, $k_V=\overline{w^{\prime }{\theta }^{\prime }}/{\theta }_z$, eddy diffusivities are plotted as functions of time.

Figure 7. Normalized bulk absolute values of the (a) horizontal, $\overline{u^{\prime }{\theta }^{\prime }}$, and (b) vertical, $\overline{w^{\prime }{\theta }^{\prime }}$, eddy fluxes of the reactive (blue) and passive (orange) tracers. The bulk absolute values of the derivatives of the tracers with respect to (a) x and (b) z are shown in purple and yellow for the reactive and passive tracers respectively. The (c) horizontal, $k_H=\overline{u^{\prime }{\theta }^{\prime }}/{\theta }_x$, and (d) vertical, $k_V=\overline{w^{\prime }{\theta }^{\prime }}/{\theta }_z$, eddy diffusivities are plotted as functions of time.

Bulk descriptions do not capture any information about the directions of the eddy fluxes or gradients. The angle between the vectors $({\Theta }_x,{\Theta }_z)$ and $(\overline{u^{\prime }{\theta }^{\prime }},\overline{w^{\prime }{\theta }^{\prime }})$ can be calculated at each grid point and indicates whether the eddy diffusivity parameterization is diffusing the tracer in the right direction. The data bins nicely into two groups, one centred around 0, and another centred around $\pm \pi$. As shown in panels (b) and (c) of Figure 8, the spread within the bins, indicated by the shaded area extending two standard deviations away from the mean of the bin, is $O\left(0.1\right)$. In panel (a), the percentage of grid points in either bin is plotted as a function of time for both the reactive and passive tracers. The percentage of points in each bin always lies within $\pm 2\%$ of 50% and, for most times, there are more points in the bin centred at 0.

Figure 8. The angle between the eddy flux and gradient vectors is calculated at each grid point. The data bins nicely into two groups centred around 0 and $\pm \pi$. In (a), the percentage of points in each bin is shown for the reactive (green, red) and passive (blue, orange) tracers. The means of the (b) 0 and (c) $\pm \pi$ bins are indicated by lines and the shaded areas extend two standard deviations away from the mean indicating the spread within each bin is reasonable. The spatial distribution of the binned points at $t=26$s and within a subdomain is shown for the (d) reactive and (c) passive tracers.

Figure 8. The angle between the eddy flux and gradient vectors is calculated at each grid point. The data bins nicely into two groups centred around 0 and $\pm \pi$. In (a), the percentage of points in each bin is shown for the reactive (green, red) and passive (blue, orange) tracers. The means of the (b) 0 and (c) $\pm \pi$ bins are indicated by lines and the shaded areas extend two standard deviations away from the mean indicating the spread within each bin is reasonable. The spatial distribution of the binned points at $t=26$s and within a subdomain is shown for the (d) reactive and (c) passive tracers.

The eddy fluxes and gradients are either almost perfectly parallel or nearly pointing in opposite directions. In the later case, a negative eddy diffusivity is needed to realign the vectors. The spatial distribution of the binned grid points at $t=26$s within a subdomain is shown for the reactive and passive tracers in panels (d) and (e) of Figure 8. In the red areas, the eddy fluxes and gradients are binned as anti-parallel and a negative eddy diffusivity is implied. The binned grid points self organize into columns, though the columns are less well-defined for the reactive tracer. The issue of negative diffusivities will be revisited in the discussion.

The plumes formed by RT instabilities primarily move vertically, though the perturbations applied to the initial density stratification do encourage some plumes to propagate diagonally. If we consider the mean initial stratification, the perturbations average out. So, the mean effect of the RT instabilities is to transport the tracers vertically until turbulent mixing becomes strong enough to induce horizontal transport. The stills in panels (d) and (e) of Figure 8 are taken just before the heads of the RT plumes reach the top of the domain after which point turbulent mixing dominates. Since eddy fluxes are an averaged quantity, they too are arranged into columns.

Meanwhile, the mean tracer concentration lies largely in a band along the top of the domain as shown in panels (a) and (b) of Figure 2. The gradient of the mean thus generally points downwards. So, while the eddy fluxes reflect properties of the flow governing the tracer transport, the gradients only reflect the spatial distribution of the mean tracer field.

The distribution of the binned data corresponding to the reactive tracer, panel (d) Figure 8, is more disorganized because the RT bubbles have already reached the mean tracer front. As previously mentioned, the increased Fickian diffusion exhibited by the reactive tracer will allow it to diffuse and then be further entrained into regions the passive tracer may not be able to reach. The reaction, then causes this process to be repeated, enriching the fluctuation field.

4. Discussion

Turbulent flows are highly unpredictable and disordered. In this study, small $O\left({10}^{-2}\right)$ perturbations to the shape of the initially unstable density stratification were enough to produce significant variability within an ensemble of 50 well-resolved simulations. RANS simulations implicitly assume that the most important dynamics are captured by the mean. We have shown that this assumption fails when considering the transport of a reactive tracer growing according to Fisher’s equation.

Under Reynolds averaging, the mean reaction equation is altered by the introduction of eddy reaction terms. The procedure in Section 2.1.2 can be applied to any polynomial reaction function. Even non-polynomial reaction functions could be subjected to the same procedure if they were approximated by Taylor polynomials. The procedure allows one to isolate the roles of the mean, fluctuations, and mean-fluctuation interactions in the Reynolds averaged reaction function. It also tells us that any terms raised to an even power will introduce sources or sinks (depending on the sign of the coefficient) into the mean ADR equation which depend only on the fluctuations. For Fisher’s equation, a sink in the mean ADR equation transfers concentration to the fluctuations. Growth, whether in the mean or fluctuations, increases diffusion allowing a reactive tracer to spread to parts of the domain that a passive tracer may not reach. These areas may include turbulent structures further transferring concentration from the mean to the fluctuations. Diffusion and turbulent mixing push the tracer concentration away from its stable equilibrium, $\theta =1$, thereby promoting further growth and creating a positive feedback loop. Numerical dispersion, while not significant in our pseudospecral simulations, has the potential to artificially increase the diffusivity of a tracer. The positive feedback would thus be intensified in an inherently dissipative model (i.e. one based on finite volumes) such as ROMS[11].

Due to the flux of concentration from the mean to the fluctuations, the mean underestimates the tracer concentrations in the individual simulations. Further, as shown by the supplementary movies, the downwards propagation of the reactive tracer is captured by the fluctuations and not the mean. In a true RANS simulation, one has no knowledge of the Reynolds stresses, eddy fluxes, or eddy reactions. In our numerical experiment, these quantities have proven to have a significant effect on the downward propagation of the reactive tracer. The bulk scalar variances (which are also the eddy reactions for Fisher’s equation) for both the passive and reactive tracers are approximately 1/4 the magnitude of the corresponding bulk mean and thus certainly non-negligible.

Given the initial distribution of the tracers (panel (b) of Figure 1), the generation of turbulent mixing with a RT instability doomed eddy diffusivity relations to fail. RT bubbles pushed up against the tracer fronts therefore causing counter-gradient fluxes and negative eddy diffusivities. As the fronts of the RT plumes hit the top and bottom of the domain, the flow transitioned from a state dominated by large scale motions to one dominated by smaller scale turbulent mixing. In the latter state, negative eddy diffusivities continue to be observed because the largest eddies are larger than the curvature of the mean vertical tracer profiles[18]. Negative eddy diffusivities have been observed in numerical studies of non-passive oceanic tracers and attributed to symmetric instabilities and other submesoscale phenomena[30] in agreement with the results presented herein. We would expect an eddy diffusivity parameterization to perform equally as bad if not worse for the same setup in three dimensions since the eddy diffusivity parameterization does not account for vortex tilting[51].

Due to the feedbacks observed between turbulent mixing and the reaction dynamics, we stress the importance of conducting studies of reactive tracers with turbulent flows. The use of convective cells to approximate turbulence is problematic primarily because the cells introduce streamlines that the tracer can only diffuse across[34,42]. Turbulent RT-induced mixing created small, irregular structures that pulled the reactive tracer into parts of the domain it otherwise could not have reached. There was a significant amount of variability in the tracer distributions of our ensemble members. With pre-specified velocity fields, it will be near impossible to capture the structures that characterize turbulent mixing and thus variability in the tracer distribution will be underestimated.

In addition to the flow used to induce turbulent mixing, the results of this study are highly specific to the reaction and Damköhler number selected. We chose $Da\sim O\left(1\right)$ to balance the effects of advection with those of the reaction. The balance causes growth and mixing processes to become coupled and eddy diffusivity parameterizations to fail[55]. Ocean tracer dynamics are complex so interactions between turbulent mixing and reaction dynamics are not as simple as the positive feedback we observed in our study[35,50,56].

When considering a slow reaction with $Da\sim O\left(0.1\right)$ or smaller, the coupling will be weaker but new challenges arise as turbulent mixing creates sharp tracer filaments. Decreasing $Da$ to $O\left(0.1\right)$ made our 0.5mm grid spacing insufficient to resolve the tracer filaments. This is problematic, given the large disparity between the scale of our simulations and any regional modeling study. Increasing the diffusivity of the tracer would smooth out the fronts thereby reducing the demand for a higher resolution. In a RANS simulation, this is done through eddy diffusivity parameterizations. Modeling studies based on ROMS generally have $O\left(1\right)$ to $O\left(10\right)$ kilometer horizontal resolution and $O\left(1\right)$ to $O\left(100\right)$ meter vertical resolution[13,14,57]. Therefore, it should be assumed that tracer distributions in these models have been over-smoothed. Further, the spacing between ROMS grid points is not uniform and the vertical resolution is often better near the surface than at the bottom (O(1m) versus O(200m) [56]). A possible improvement may be simulating near surface effects, which are weighed more heavily due to the grid spacing, while parameterizing those deeper in the domain.

Future directions for this work include running simulations with a deeper domain to allow for a longer transport period before the free slip walls begin to affect the flow. Though two dimensional RT instabilities were sufficient for this study, their mixing characteristics are well-known to differ from their three dimensional counterparts[44,47]. It would be interesting to determine whether a three dimensional RT instability affects the reactive tracer any differently, for example through organized larger scale motions. Another possible avenue for exploration may be to use the data from the ensemble generated for this study to force simulations of the mean ADR equation. These simulations would allow the effects of the eddy flux and eddy reactions to be isolated and allow for the evaluation of potential parameterizations.