2 The parent DNS dataset

We model a stratified mixing layer by assuming initial velocity and density distributions that have a hyperbolic tangent form, as

(2.1a,b ) $$\begin{eqnarray}\overline{u}(z,0)=U_{0}\tanh \left(\frac{z}{d}\right),\quad \overline{\unicode[STIX]{x1D70C}}(z,0)=\unicode[STIX]{x1D70C}_{0}\left[1-\tanh \left(\frac{z}{\unicode[STIX]{x1D6FF}}\right)\right],\end{eqnarray}$$

in the Boussinesq approximation such that $\unicode[STIX]{x1D70C}_{0}\ll \unicode[STIX]{x1D70C}_{r}$ (note that here density represents departures from a hydrostatic state associated with  $\unicode[STIX]{x1D70C}_{r}$ ). Also, $U_{0}$ and $\unicode[STIX]{x1D70C}_{0}$ denote, respectively, half the total velocity and density jumps across the shear layer (with a total depth  $2d$ ) and the density layer (with a total depth of  $2\unicode[STIX]{x1D6FF}$ ). As a result of this canonical setting, the dimensionless Boussinesq equations are governed by four important non-dimensional parameters, namely the (initial) Reynolds number $Re$ ; the bulk Richardson number $Ri_{b}$ ; the Prandtl number $Pr$ ; and the initial scale ratio  $R$ , defined together as

(2.2a-d ) $$\begin{eqnarray}Re=\frac{U_{0}d}{\unicode[STIX]{x1D708}},\quad Ri_{b}=\frac{g\unicode[STIX]{x1D70C}_{0}d}{\unicode[STIX]{x1D70C}_{r}U_{0}^{2}},\quad Pr=\frac{\unicode[STIX]{x1D708}}{\unicode[STIX]{x1D705}},\quad R=\frac{d}{\unicode[STIX]{x1D6FF}},\end{eqnarray}$$

in which $\unicode[STIX]{x1D708}$ is the kinematic viscosity, $\unicode[STIX]{x1D705}$ is the molecular diffusivity and $g$ is the gravitational acceleration. Table 1 lists all the DNS analyses from which data will be employed for training and validation of the proposed artificial neural networks. These simulations have been thoroughly analysed and discussed previously in a number of recent publications on KHI (Salehipour & Peltier Reference Salehipour and Peltier2015; Salehipour, Peltier & Mashayek Reference Salehipour, Peltier and Mashayek2015; Salehipour et al. Reference Salehipour, Peltier, Whalen and MacKinnon2016b ) and HWI (Salehipour et al. Reference Salehipour, Caulfield and Peltier2016a , Reference Salehipour, Peltier and Caulfield2018). For details of each simulation, interested readers are referred to the relevant papers.

Table 1.

The collection of initial parameters (as defined in (2.2)) employed for conducting DNS experiments associated with either KHI or HWI. $n_{s}$ indicates the number of saved snapshots for each individual simulation. The split between training and validation sets is also highlighted.

For the supervised machine learning application to be discussed herein, we have further subdivided these datasets into training and validation sets with an approximate 80 %–20 % ratio, as indicated in table 1. Both these subsets include examples of flow evolution due to KHI and HWI. We have intentionally chosen the validation dataset to include all DNS cases with extreme values for their initial parameters, which are well beyond the range of similar parameters employed for training purposes. This enables us to investigate the extent to which our trained model is generalizable and thus robust. Note that our training dataset had a very limited number of HWI examples (compared to KHI) and that these examples are also at much smaller values of $Ri_{b}$ and  $R$ .

3 Preprocessing of DNS data

The result of each three-dimensional DNS experiment associated with the evolution of either KHI or HWI is comprised of $n_{s}$ snapshots in time, where each saved snapshot represents three-dimensional fields of flow quantities, namely the density $\unicode[STIX]{x1D70C}$ and velocity fields $\boldsymbol{u}=(u,v,w)$ (table 1 lists $n_{s}$ for each simulation). The intensity of turbulent activity may be represented by the pointwise dissipation rate of total kinematic energy, $\unicode[STIX]{x1D716}(\boldsymbol{x},t)$ , defined as

(3.1) $$\begin{eqnarray}\unicode[STIX]{x1D716}(\boldsymbol{x},t)=2\unicode[STIX]{x1D708}\unicode[STIX]{x1D634}_{ij}\unicode[STIX]{x1D634}_{ij},\end{eqnarray}$$

in which $\unicode[STIX]{x1D634}_{ij}=(\unicode[STIX]{x2202}u_{i}/\unicode[STIX]{x2202}x_{j}+\unicode[STIX]{x2202}u_{j}/\unicode[STIX]{x2202}x_{i})/2$ is the total strain rate tensor. We may also reduce the above three-dimensional fields into a one-dimensional profile by performing horizontal averaging (to be denoted here by an overbar). Thus the horizontally averaged dissipation rate of total kinematic energy, $\overline{\unicode[STIX]{x1D716}}(z,t)$ , and the mean flow density, $\overline{\unicode[STIX]{x1D70C}}(z,t)$ , are defined as

(3.2a,b ) $$\begin{eqnarray}\overline{\unicode[STIX]{x1D716}}(z,t)=\frac{1}{L_{x}L_{y}}\displaystyle \int \unicode[STIX]{x1D716}(\boldsymbol{x},t)\,\text{d}x\,\text{d}y,\quad \overline{\unicode[STIX]{x1D70C}}(z,t)=\frac{1}{L_{x}L_{y}}\displaystyle \int \unicode[STIX]{x1D70C}(\boldsymbol{x},t)\,\text{d}x\,\text{d}y,\end{eqnarray}$$

where $L_{x}$ and $L_{y}$ denote the size of the computational domain in the streamwise and spanwise directions.

The (generally) time-dependent mixing efficiency, $\mathscr{E}$ , may be computed precisely by invoking the concept of irreversible diapycnal mixing (originally introduced by Winters et al. Reference Winters, Lombard, Riley and D’Asaro1995), which relies on a special kind of reduction operator, namely a three-dimensional sorting of the density field into a notional state that is strictly stably stratified (Peltier & Caulfield Reference Peltier and Caulfield2003), and is defined as

(3.3) $$\begin{eqnarray}\mathscr{E}(t)=\frac{\mathscr{M}(t)}{\mathscr{M}(t)+\langle \overline{\unicode[STIX]{x1D716}}(z,t)\rangle },\end{eqnarray}$$

where $\langle \,\rangle$ denotes vertical averaging. For a precise definition of $\mathscr{M}(t)$ refer to (2.18) of Salehipour et al. (Reference Salehipour, Caulfield and Peltier2016a ) and the cited discussions therein. The required parallel implementation of the sorting procedure is described in Salehipour et al. (Reference Salehipour, Peltier and Mashayek2015) (see, for example, their figure 1). Such an elaborate technique for calculating $\mathscr{E}$ is only viable in numerical simulations such as those employed in this work because in practice oceanographers only measure one-dimensional profiles in depth and are therefore unable to perform the same analysis. Indeed, there is a similar subtlety in defining the ‘background’ buoyancy frequency, $N^{2}(z,t)$ , as described in Salehipour & Peltier (Reference Salehipour and Peltier2015) (see their discussion leading to (2.23)) and more recently in Arthur et al. (Reference Arthur, Venayagamoorthy, Koseff and Fringer2017). In order to distinguish between irreversible mixing and reversible stirring, $N^{2}(z,t)$ must be defined based on the same notional state obtained by the three-dimensional sorting procedure. For consistency with common practice in oceanography, in this paper we may define $N^{2}$ using the mean flow density introduced in (3.2) such that $N^{2}(z,t)=-(g/\unicode[STIX]{x1D70C}_{r})\,\text{d}\overline{\unicode[STIX]{x1D70C}}/\text{d}z$ .

We seek a mapping between the instantaneous vertical profiles of $\overline{\unicode[STIX]{x1D716}}(z,t_{0})$ and $N^{2}(z,t_{0})$ (i.e. at a given time  $t_{0}$ ) and the precisely computed values of mixing efficiency, $\mathscr{E}(t_{0})$ . Once the network is trained, this mapping would essentially reveal a reduction operator that is conceivably very different from a straightforward vertical averaging; one which also incorporates the structural pattern and length scales that implicitly exist and are thus ‘hidden’ in these profiles. Thus the inputs to our artificial neural network are tuples of $({\mathcal{X}}_{1},{\mathcal{X}}_{2})$ defined respectively as

(3.4a,b ) $$\begin{eqnarray}{\mathcal{X}}_{1}(z,t_{0})\equiv \frac{\overline{\unicode[STIX]{x1D716}}(z,t_{0})}{\unicode[STIX]{x1D705}\displaystyle \int N^{2}(z,t_{0})\,\text{d}z},\quad {\mathcal{X}}_{2}(z,t_{0})\equiv \frac{N^{2}(z,t_{0})}{\displaystyle \int N^{2}(z,t_{0})\,\text{d}z}.\end{eqnarray}$$

Furthermore, the true ‘labels’ in our supervised learning setting are the instantaneous values of mixing efficiency, namely

(3.5) $$\begin{eqnarray}{\mathcal{Y}}(t_{0})\equiv \mathscr{E}(t_{0}).\end{eqnarray}$$

It is important to highlight that $({\mathcal{X}}_{1},{\mathcal{X}}_{2})$ appear in a normalized form to render $\overline{\unicode[STIX]{x1D716}}(z,t_{0})$ and $N^{2}(z,t_{0})$ comparable in terms of their dimensionality and physical relevance, and furthermore to extend the applicability of the trained network to oceanographic profiles. In (3.4), ${\mathcal{X}}_{1}$ represents the vertical profile of kinetic energy dissipation rate relative to the molecular diffusion rate in the absence of mean flow shear. The vertical profiles are assumed to have a fixed length of $i=512$ points in which the ‘dead’ regions of the simulation (near top and bottom boundaries) have been excluded by focusing on the largest segments of the profiles where $|\text{d}\overline{\unicode[STIX]{x1D70C}}(z)/\text{d}z|\geqslant 10^{-3}$ . This approach is analogous to identifying ‘patches’ of turbulence from DNS calculations (Smyth, Moum & Caldwell Reference Smyth, Moum and Caldwell2001).

4 Deep convolutional neural networks

Deep learning methods involve a multilayer stacking of simple modules that perform linear or nonlinear input–output mappings whose weights and biases are subject to ‘training’ through an optimization procedure (LeCun, Bengio & Hinton Reference LeCun, Bengio and Hinton2015). These techniques became widely popularized after Krizhevsky, Sutskever & Hinton (Reference Krizhevsky, Sutskever and Hinton2012), from University of Toronto, employed a ‘deep convolutional neural network’ to classify a dataset of 1.2 million images and won first place in the 2012 ImageNet competition. A convolutional neural network (CNN) is a special type of neural network architecture that relies on the convolution operator in lieu of general matrix multiplication in at least one layer of its configuration (Goodfellow, Bengio & Courville Reference Goodfellow, Bengio and Courville2016). In two dimensions, this operator may be defined as

(4.1) $$\begin{eqnarray}\widetilde{\boldsymbol{{\mathcal{X}}}}(i,j)=(\boldsymbol{{\mathcal{X}}}\ast \mathscr{K})(i,j)=\mathop{\sum }_{m}\mathop{\sum }_{n}\boldsymbol{{\mathcal{X}}}(m,n)\mathscr{K}(i-m,j-n),\end{eqnarray}$$

where the input field $\boldsymbol{{\mathcal{X}}}\in \mathbb{R}^{i\times j}$ and the convolution kernel $\mathscr{K}\in \mathbb{R}^{m\times n}$ is represented by characteristic filter lengths of size $m\leqslant i$ and $n\leqslant j$ . A convolved (i.e. filtered) field is constructed by traversing the kernel $\mathscr{K}$ over the dimensions of  $\boldsymbol{{\mathcal{X}}}$ . To keep the filtered field, $\widetilde{\boldsymbol{{\mathcal{X}}}}$ , the same size as $\boldsymbol{{\mathcal{X}}}$ , zero-padding is often employed.

Figure 2.

Illustration of the type of convolutional neural network investigated in this paper with increasing number of stacked layers. Refer to the text for definition of the various layers. For any given example data point, the tuples of $\boldsymbol{{\mathcal{X}}}=({\mathcal{X}}_{1},{\mathcal{X}}_{2})$ are of size ( $512\times 2$ ), representing vertical profiles of normalized $\unicode[STIX]{x1D716}(z)$ and $N^{2}(z)$ as defined in (3.4).

Figure 2 illustrates the schematic configuration of the selected neural network architectures to be employed in this paper, which may consist of one to seven convolution layers labelled as CNN1 to CNN7. Each configuration receives the input $\boldsymbol{{\mathcal{X}}}=({\mathcal{X}}_{1},{\mathcal{X}}_{2})\in \mathbb{R}^{512\times 2}$ , as defined in (3.4), and passes it to a ‘batch normalization’ layer (Ioffe & Szegedy Reference Ioffe and Szegedy2015) that, for any given ‘batch’ of the training dataset (with size  $n_{b}$ ), normalizes ${\mathcal{X}}_{1}$ and ${\mathcal{X}}_{2}$ individually by subtracting the batch mean and dividing by its variance. The batch-normalized data are then subsequently fed into a series of convolution layers each having 64 filters with a kernel $\mathscr{K}\in \mathbb{R}^{4\times 2}$ . Each convolution kernel undergoes a nonlinear activation function of type $f(x)=\max (0,x)$ , also known as a Rectified Linear Unit (ReLU). The output of each convolution layer is followed by an ‘average pooling’ operator which effectively reduces the size of the profile by half through averaging any two adjacent data in the profile (i.e. averaging window of size $(2\times 1)$ ). The reduced outputs are then reshaped appropriately to be fed into a dense (or fully connected) layer with 64 neurons that also employs ReLU as its nonlinearity function. To avoid overfitting and to improve the model predictions, we regularize the network by a method known as ‘dropout’ (Hinton et al. Reference Hinton, Srivastava, Krizhevsky, Sutskever and Salakhutdinov2012), which randomly turns off 50 % of the neurons, thereby eliminating their contribution in the ‘backpropagation’ procedure (a method to apply the chain rule to derive gradients of the loss function with respect to trainable parameters in the network) of the optimization process. As a result the network is forced to learn robust features that emerge more frequently in the random subsets during training. Finally we use a single neuron to represent the network output, $\hat{{\mathcal{Y}}}$ . We have chosen a sigmoid activation function for the output layer because efficiency values must be between 0 and 1. We have used the Adam optimizer (Kingma & Ba Reference Kingma and Ba2014) for performing stochastic gradient descent to minimize the loss function, defined as the mean squared error $\sum _{n_{b}}({\mathcal{Y}}-\hat{{\mathcal{Y}}})^{2}/n_{b}$ , where the batch size is set as $n_{b}=100$ .

Notice that the number of trainable parameters, $n_{p}$ , in the network decreases from CNN1 to CNN6 and increases slightly from CNN6 to CNN7 (see figure 3 a). Taking CNN1 for instance, the network has copious neurons that link the first convolution layer (after pooling) to a fully connected dense layer, leading to $n_{p}>2\times 10^{6}$ . As the network becomes deeper, increasingly more structure is built into the network due to the locality of the convolution operators that may be contrasted to the global connectivity of the dense layers. By definition (4.1), the convolution layer shares its kernel parameters ( $m\times n$ weights for a two-dimensional kernel) across a given ‘tiling’ of its input field. Notice that for deeper convolution layers whose inputs are derivative of the previous pooling operator with $q$ filters, there are $m\times n\times q$ trainable weights and one trainable bias that are shared by each tiling of the convolution layer. For instance in CNN2, $n_{p}\sim 64\times (4\times 2+1)_{conv1}+64\times (4\times 2\times 64+1)_{conv2}+(128\times 2\times 64\times 64)_{dense}$ . Refer to Goodfellow et al. (Reference Goodfellow, Bengio and Courville2016, Chap. 9) for further details on parameter sharing in convolutional networks. Figure 3(a) also evaluates the effect of increasing the depth of the network in so far as the validation data is concerned. Clearly CNN6 outperforms others, which might be explained by the observed saturation of network training capacity also shown in this figure.

Figure 3.

(a) Comparing the predictions of CNN1–CNN7 on the validation set. (b) The results of CNN6 (this study) are compared quantitatively with other methods in terms of mean squared error for various sets of validation data, labelled (a–h) as per figure 5.

Obviously there are many parameters (or hyper-parameters) that we have assumed to be fixed within the above networks. Moreover, there are many other types of deep neural networks (DNN) (Goodfellow et al. Reference Goodfellow, Bengio and Courville2016) that could be exploited – an alternatively good candidate being the recurrent neural network, which enables handling input sequences of vertical profiles with varying dimensions. We only note in passing that we have also investigated a deep feed-forward neural network (that consists of deep stacking of dense layers only) in this work, but the CNN results were substantially more accurate. We wish to emphasize that our focus in this paper has not been to find the optimal configuration (or architecture) for producing the least possible error on the validation set. This paper rather intends to make the first step in introducing the idea of employing deep learning methods for the purpose of parameterizing subgrid scale processes using DNS datasets. We have open-sourced our code and postprocessed dataset in the hope of encouraging the community to further enhance such a data-driven approach to parameterization. Indeed, we believe the ultimate success of these efforts will rely upon a cohesive community-driven collaboration. Fortunately, these data-driven ideas are being embraced most recently in the broader context of Earth system modelling (Schneider et al. Reference Schneider, Lan, Stuart and Teixeira2017).

6 Summary

Using properly normalized vertical structures of $N^{2}(z)$ and $\overline{\unicode[STIX]{x1D716}}(z)$ , we have proposed a data-driven approach based on deep convolutional neural networks (CNN) that can accurately predict the value of mixing efficiency for the entire life cycle of KHI and HWI (i.e. two ‘atoms’ of turbulence in stratified flows) beyond the range of initial conditions that have been employed for training the network. The large overturns that are convectively unstable are no longer ignored in such an approach. We have also shown that the results of the CNN model for KHI and HWI are more reliable and accurate than those based on the Osborn–Cox method.

Deep neural networks have a compositional hierarchy in which low-level features are composed to form higher-level features (for example, for image recognition, the first layers of a CNN detect basic abstract features such as edges, then deeper layers combine edges to form motifs, and subsequent layers assemble parts from motifs). We believe the proposed CNN model has similarly discovered such an ‘abstract’ level of stratified turbulence with characteristics that are so universal that even with a small portion of data associated with HWI, the generic behaviour of its induced mixing efficiency can be predicted robustly for wildly different initial conditions.

What makes such a data-driven approach especially appealing is its capability to become increasingly more accurate, robust and generic. This is foreseen to be achieved by (i) experimenting with many other types of DNN architectures, (ii) tuning the hyper-parameters (of which there are many) and, perhaps most importantly, (iii) further enriching the training dataset by adding additional examples of KHI and HWI, as they become available, or perhaps more excitingly by including more ‘atoms’ of ocean turbulence such as those induced by, for example, double-diffusion, Taylor and Rayleigh–Taylor instabilities. Another exciting future direction would involve using observed profiles (either from laboratory or real environments) to estimate mixing efficiency based on the proposed model, especially due to the relative inaccuracies of the Osborn–Cox method.