2.1. The CNN method

We utilise a CNN (Lecun et al. Reference Lecun, Bottou, Bengio and Haffner1998), which is widely used to process grid-like data to obtain their spatial features, to model the mapping $\mathcal {F}$ described above for flow reconstruction purposes. The overall architecture of the NN is sketched in figure 1(a). The model consists of two parts: an encoder that maps the surface quantities into a representation with reduced dimensions and a decoder that reconstructs the three-dimensional flow field from the reduced-dimension representation. The input fed into the encoder includes the surface quantities $\boldsymbol {E}$ sampled on a grid of size $256 \times 128$. With $\boldsymbol {E}$ consisting of four variables, the total dimensions of the input values are $256 \times 128 \times 4$. The encoding process, which outputs a compressed representation of dimensionality $32 \times 16\times 24$, extracts the most important surface features for the subsequent reconstructions, which may alleviate the problem of the model overfitting insignificant or non-generalisable features (Verleysen & François Reference Verleysen and François2005; Ying Reference Ying2019). The decoder reconstructs velocities with three components from the output of the encoder onto a three-dimensional grid of size $128\times 64 \times 96$. The total dimensions of the output data are $128\times 64 \times 96 \times 3$. The input and output of the CNN are based on the same time instant, and each instant is considered independently; therefore, no temporal information is used for the reconstruction process.

Figure 1. Building blocks of the network for flow reconstruction. (a) Overview of the reconstruction process. (b) Structure of the residual block in the encoder. Here, SE refers to the squeeze-and-excitation layer (Hu, Shen & Sun Reference Hu, Shen and Sun2018).

The detailed architectures of the CNN layers are listed in table 1. The basic building block of a CNN is the convolutional layer, which is defined as

(2.3)\begin{equation} X^{out}_n = b_n + \sum_{k=1}^{K^{in}} W_n(k)*X^{in}_k,\quad n=1,\ldots,K^{out}, \end{equation}

where $X^{in}_k$ denotes the $k$th channel of an input consisting of $K^{in}$ channels, $X^{out}_n$ denotes the $n$th channel of the $K^{out}$-channel output, $b$ denotes the bias, $W$ is the convolution kernel and ‘$*$’ denotes the convolution operator. Strictly speaking, ‘$*$’ calculates the cross-correlation between the kernel and the input, but in the context of a CNN, we simply refer to it as convolution. We note that, in (2.3), the input may consist of multiple channels. For the first layer that takes the surface quantities as inputs, each channel corresponds to one surface variable. Each convolutional layer can apply multiple kernels to the input to extract different features and as a result, its output may also contain multiple channels. Each channel of the output is thus called a feature map.

Table 1. Detailed architecture of the CNN with its parameters, including the input size, kernel size and stride of each block (if applicable). The input size is written as $N_x\times N_y \times N_c$ for two-dimensional data for the encoder or $N_x \times N_y \times N_z \times N_c$ for three-dimensional data for the decoder, where $N_x$, $N_y$ and $N_z$ denote the grid dimensions in the $x$-, $y$- and $z$-directions, respectively, and $N_c$ denotes the number of channels. Note that the output of each block is the input of the next block. The output size of the last block is $128\times 64 \times 96 \times 3$.

These feature maps are often fed into a nonlinear activation function, which enables the CNN to learn complex and nonlinear relationships. We have compared the performance of several types of nonlinear activation functions, including the rectified linear unit function, the sigmoid function and the hyperbolic tangent function. The hyperbolic tangent function, which performs better, is selected in the present work. We note that Murata et al. (Reference Murata, Fukami and Fukagata2020) also found the hyperbolic tangent function performs well for the mode decomposition of flow fields.

In the encoder, we use strided convolution layers and blur pooling layers to reduce the dimensionality of feature maps. In a strided convolution, the kernel window slides over the input by more than one pixel (grid point) at a time, and as a result of skipping pixels, the input is downsampled. A blur pooling layer can be considered a special strided convolution operation with fixed weights. The blur kernel we use here is the outer product of $[0.25,0.5,0.25]$ (Zhang Reference Zhang2019), which applies low-pass filtering and downsampling to the input at the same time. The blur pooling can improve the shift invariance of the CNN (Zhang Reference Zhang2019), which is discussed further below. In the encoder part, we also utilise a design called residual blocks, which is empirically known to ease the training processes of NNs (He et al. Reference He, Zhang, Ren and Sun2016). The structure of a residual block is illustrated in figure 1(b). The output of the block is the sum of the original input and the output of several processing layers. Following Vahdat & Kautz (Reference Vahdat and Kautz2020), the residual block in this work consists of two rounds of activation–convolution layers, followed by a squeeze-and-excitation layer that enables the network to learn crucial features more effectively (Hu et al. Reference Hu, Shen and Sun2018). In our experiments, we find that the use of residual blocks in the encoder part can reduce the number of epochs needed to train the network. The residual blocks can also be implemented for the decoder as in Vahdat & Kautz (Reference Vahdat and Kautz2020),but we do not observe any significant differences in the training or the performance of the network.

The reconstruction process of the decoder uses a transposed convolution operation, which is essentially the adjoint operation of convolution. To interpret the relationships between a convolution and a transposed convolution, we first note that a convolution operation can be computed as a matrix-vector product, with the matrix defined by the kernel weights and the vector being the flattened input. By transposing the matrix, one obtains the adjoint operation, i.e. the transposed convolution, associated with the initial convolution operation. It should also be noted that the dimensions of the input and output of a transposed convolution are the reverse of those of the corresponding convolution. In the decoder, strided transposed convolutions have the opposite effect relative to that of the strided convolutions in the encoder, i.e. they are used to upsample the feature maps.

To ensure that the CNN can better represent physical correlations between the surface and subsurface flows, the CNN architecture is designed to maintain good translation invariance, i.e. the ability to produce equivalent reconstructions when the surface input is shifted in space. Although the convolution operations are shift invariant by themselves, the downsampling by strided convolutions can produce different results when inputs are shifted (Azulay & Weiss Reference Azulay and Weiss2019; Zhang Reference Zhang2019). From the perspective of the Fourier analysis, downsampling may introduce aliasing errors at high wavenumber modes of the feature maps, which may be amplified by subsequent layers, and as a result, break the translation invariance of the network. Therefore, following Zhang (Reference Zhang2019), we employ the blur pooling layers, which perform antialiasing while downsampling to reduce aliasing errors. The blur pooling layers act as the first two downsampling layers in the encoder (see table 1). As reported in Appendix A, the present network can produce more consistent reconstructions for shifted inputs than a network that only uses strided convolutions for downsampling. The blur pooling is not applied to the last downsampling operation in the encoder because we find that it barely improves the translation invariance but decreases the reconstruction accuracy as the blurring can result in loss of information (Zhang Reference Zhang2019). We also apply periodic paddings in the horizontal directions to reduce the translation-induced errors owing to the edge effect. Furthermore, we apply random translations in the horizontal directions to the training data, known as data augmentation, which is detailed in § 2.3.2, to help improve the translation invariance (Kauderer-Abrams Reference Kauderer-Abrams2017). We note that complete translation invariance is difficult to achieve because nonlinear activation functions can also introduce aliasing effects (Azulay & Weiss Reference Azulay and Weiss2019). Nevertheless, our tests indicate that the reconstructions are consistent for shifted inputs and therefore the translation-induced errors should be small.

The choice of model parameters, such as the kernel sizes and strides, is often a balance between performance and computational cost. Most existing CNNs use small kernels, with sizes between $3$ and $7$ (see e.g. Simonyan, Vedaldi & Zisserman Reference Simonyan, Vedaldi and Zisserman2014; He et al. Reference He, Zhang, Ren and Sun2016; Lee & You Reference Lee and You2019; Tan & Le Reference Tan and Le2019), for computational efficiency. In the present work, we use kernels of sizes varying from $3$ and $5$. For transposed convolutions, the kernel size is chosen to be divisible by the stride to reduce the checkerboard artefacts (Odena, Dumoulin & Olah Reference Odena, Dumoulin and Olah2016). For both strided convolutions and blur pooling layers, we use a stride of $2$, i.e. a downsampling ratio of 2, which is a common choice when processing images (Tan & Le Reference Tan and Le2019) and fluid fields (Lee & You Reference Lee and You2019; Park & Choi Reference Park and Choi2020). Similarly, in the decoder, the spatial dimensions are in general expanded by a ratio of $2$, except for one layer where the features maps are expanded by a ratio of $3$ in the vertical direction owing to the fact that the final dimension $96$ has a factor of $3$. Some transposed convolution layers only expand one spatial dimension to accommodate the different expansion ratios needed in different dimensions to obtain the output grid. Although there are many combinations for how the dimension expansions are ordered in the decoder, we find that the network performance is not particularly sensitive to the ordering.

The performance of the CNN is also affected by the number of channels (feature maps) in each layer. A network with more layers can in theory provide a greater capacity to describe more types of features. As presented in § 1 of the supplementary material available at https://doi.org/10.1017/jfm.2023.154, a network with $50\,\%$ more channels has almost the same performance as the one presented in table 1, indicating that the present network has enough number of feature maps for expressing the surface–subsurface mapping. A special layer worth more consideration is the output of the encoder or the input of the decoder, which is the bottleneck of the entire network and encodes the latent features that are useful for reconstructions. The effect of the dimension of this bottleneck layer on the reconstruction performance is studied by varying the number of channels in this layer. It is found that the present size with $24$ layers can retain most surface features necessary for the network to achieve an optimal performance. More detailed comparisons are provided in § 1 of the supplementary material. We have also considered some variants of the network architecture to investigate whether a deeper network can improve the reconstruction performance. These modified networks show no significant performance differences from the current model, indicating that the present network already has enough expressivity to describe the mapping for subsurface flow reconstructions. Detailed comparisons are reported in § 2 of the supplementary material.

The entire network is trained to minimise the following loss function:

(2.4)\begin{equation} J = \frac{1}{3 L_z}\sum_{i=1}^3 {\int_z \frac{\displaystyle\iint_{x,y}{\left| \tilde{u}_i - u_i \right|}^2 \mathrm{d}\kern0.06em x\,\mathrm{d}y}{\displaystyle\iint_{x,y} u_i^2 \,\mathrm{d}\kern0.06em x,\mathrm{d}y} \mathrm{d}z}, \end{equation}

which first computes the mean squared errors normalised by the plane-averaged Reynolds normal stresses for each component and for each horizontal planes, and then averages these relative mean squared errors. This loss function is defined to consider the inhomogeneity of the open-channel flow in the vertical direction and the anisotropy among the three velocity components. Compared with the mean squared error of the velocity fluctuations (equivalent to the error of the total turbulent kinetic energy), the above loss function gives higher weights to the less energetic directions and less energetic flow regions.

The incompressibility constraint is embedded into the network to produce a reconstruction that satisfies mass conservation. At the end of the decoder, the reconstructed velocity is determined by

(2.5)\begin{equation} \tilde{\boldsymbol{u}} = \boldsymbol{\nabla} \times \tilde{\boldsymbol{A}}, \end{equation}

where $\tilde {\boldsymbol {A}}$ is a vector potential estimated by the network. The solenoidality of $\tilde {\boldsymbol {u}}$ is automatically satisfied following the vector identity $\boldsymbol {\nabla } \boldsymbol {\cdot } (\boldsymbol {\nabla }\times \tilde {\boldsymbol {A}}) = 0$. It should be noted that the vector potential of a velocity field is not unique, which can be seen from $\boldsymbol {\nabla }\times \tilde {\boldsymbol {A}} = \boldsymbol {\nabla }\times (\tilde {\boldsymbol {A}} + \boldsymbol {\nabla }\psi )$, where $\psi$ is an arbitrary scalar function. Considering that $\tilde {\boldsymbol {A}}$ is just an intermediate quantity and it is the velocity field we are interested in, we do not add any constraints on $\boldsymbol {A}$ to fix the choice of $\psi$. In other words, we let the network architecture and the optimisation process freely estimate the vector potential. We find that predicting the velocity field through its vector potential and predicting the velocity directly yields no discernible differences in the reconstruction accuracy (see detailed results in § 4 of the supplementary material), indicating that the imposed incompressibility constraint does not negatively impact the CNN's reconstruction capability.

The spatial derivatives in (2.5) are computed using the second-order central difference scheme, which can be easily expressed as a series of convolution operations with fixed kernel coefficients. We have also tested the Fourier spectral method to calculate the derivatives in the $x$- and $y$-directions, which barely affect the reconstruction accuracy but significantly increases the training time. This indicates that a higher-order scheme does not provide additional reconstruction accuracy, which we believe is because most reconstructable structures are low-wavenumber structures and the second-order central scheme is adequately accurate for resolving these structures.

The total number of trainable parameters in the designed network is 279 551. The network is trained using the adaptive moment estimation with decoupled weight decay optimiser (Loshchilov & Hutter Reference Loshchilov and Hutter2019) with early stopping to avoid overfitting. The training process is performed on an NVIDIA V100 GPU.

2.2. The LSE method

For comparisons with the CNN model, we also consider a classic approach based on the LSE method (Adrian & Moin Reference Adrian and Moin1988), which is commonly used to extract turbulent coherent structures from turbulent flows (Christensen & Adrian Reference Christensen and Adrian2001; Xuan, Deng & Shen Reference Xuan, Deng and Shen2019). It has also been applied to reconstruct the velocity on an off-wall plane in a turbulent channel flow from wall measurements (Suzuki & Hasegawa Reference Suzuki and Hasegawa2017; Encinar & Jiménez Reference Encinar and Jiménez2019; Guastoni et al. Reference Guastoni, Encinar, Schlatter, Azizpour and Vinuesa2020). Wang et al. (Reference Wang, Gao, Wang, Pan and Wang2021) demonstrated that this method can be used to transform a three-dimensional vorticity field into a velocity field.

The LSE method, as the name suggests, estimates a linear mapping from the observations of some known random variables to unknown random variables. In the present problem, the known random variables are the surface quantities $\boldsymbol {E}$ and the unknowns are the estimated $\tilde {\boldsymbol {u}}$. The reconstruction of $\tilde {\boldsymbol {u}}$ by a linear relationship can be expressed as

(2.6)\begin{equation} \tilde{u}_i(\boldsymbol{x}) = \mathcal{L}_{ij}(\boldsymbol{x}; \boldsymbol{x}_s) E_j(\boldsymbol{x}_s),\quad i=1,2,3;\enspace j=1,2,3,4, \end{equation}

where $\mathcal {L}_{ij}$ (written as $\mathcal {L}$ below for conciseness) are linear operators. In other words, (2.6) treats the mapping $\mathcal {F}$ in (2.1) as a linear function of $\boldsymbol {E}$. Note that two distinct coordinates are contained in the above equation: $\boldsymbol {x}$ and $\boldsymbol {x}_s$. The coordinate $\boldsymbol {x}$ denotes the position to be reconstructed and $\mathcal {L}$ is a function of $\boldsymbol {x}$ because the mapping between $\boldsymbol {E}$ and $\tilde {\boldsymbol {u}}$ should depend on the location to be predicted. The coordinate $\boldsymbol {x}_s$ denotes the horizontal coordinates, i.e. $\boldsymbol {x}_s=(x,y)$, which signifies that the observed surface quantity $\boldsymbol {E}(\boldsymbol {x}_s)$ is a planar field. Here, we denote the surface plane as $S$, i.e. $\boldsymbol {x}_s \in S$. For each $\boldsymbol {x}$, $\mathcal {L}$ is a field operator on $S$ that transforms $\boldsymbol {E}(\boldsymbol {x}_s)$ into the velocity $\boldsymbol {\tilde {u}}(x)$. The LSE method states that $\mathcal {L}$ can be determined by

(2.7)\begin{equation} \left\langle{E_m(\boldsymbol{r}') E_j(\boldsymbol{x}_s)}\right\rangle \mathcal{L}_{ij}(\boldsymbol{x};\boldsymbol{x}_s) = \left\langle{u_i(\boldsymbol{x}) E_m(\boldsymbol{r}')}\right\rangle,\quad \forall\ r'\in S;\enspace m=1,2,3,4 ,\end{equation}

where $\langle {\cdot }\rangle$ is defined as averaging over all the ensembles (instants) in the dataset. The above equation is obtained by minimising the mean squared errors between $\tilde {u}_i(\boldsymbol {x})$ and ${u}_i(\boldsymbol {x})$ over all the ensembles (instants) in the dataset; i.e. the operator $\mathcal {L}$ defined above is the optimal estimation of the mapping $\mathcal {F}$ in the least-squares sense. Note that, on the left-hand side of (2.7), $\langle {E_m(\boldsymbol {r}') E_j(\boldsymbol {x}_s)}\rangle$ should be considered a planar field function of $\boldsymbol {x}_s$ when multiplied by the operator $\mathcal {L}$. By varying $\boldsymbol {r}'$ and $m$, one can obtain a linear system to solve for $\mathcal {L}$. Similar to the CNN model in § 2.1, the LSE model only describes the spatial correlations between the subsurface and surface quantities.

As pointed out by Wang et al. (Reference Wang, Gao, Wang, Pan and Wang2021), the operator $\mathcal {L}$ in (2.6) can be alternatively written in an integral form as

(2.8)\begin{equation} \tilde{u}_i(\boldsymbol{x}) = \iint_S {l_{ij}(\boldsymbol{x}; \boldsymbol{x}_s) E_j(\boldsymbol{x}_s)}\,\mathrm{d}\kern0.7pt \boldsymbol{x}_s. \end{equation}

Utilising the horizontal homogeneity of the open-channel turbulent flow, it can be easily shown that $l_{ij}$ is a function that is dependent on $\boldsymbol {x}-\boldsymbol {x}_s$ and, therefore, (2.8) can be written as

(2.9)\begin{equation} \tilde{u}_i(\boldsymbol{x}) = \iint_S {l_{ij} (\boldsymbol{x}-\boldsymbol{x}_s) E_j(\boldsymbol{x}_s)}\,\mathrm{d}\kern0.7pt \boldsymbol{x}_s. \end{equation}

Therefore, the LSE of the subsurface velocity field can be obtained by a convolution of the kernel $l_{ij}$ and the surface quantities $E_j$ over the horizontal plane $S$. As a result, the vertical coordinate $z$ is decoupled, and (2.9) can be evaluated independently for each $x$–$y$ plane.

For efficiently computing $l_{ij}$, we apply the convolution theorem to (2.9) and transform it into pointwise multiplications in the Fourier wavenumber space. Similarly, (2.7) can also be expressed in the convolutional form. This yields

(2.10)$$\begin{gather} \hat{\tilde{u}}_i = \hat{l}_{ij}\hat{E}_j, \end{gather}$$

(2.11)$$\begin{gather}\left\langle{\hat{E}_m^* \widehat{E_j}}\right\rangle \hat{l}_{ij} = \left\langle{\hat{E}_m^* \hat{u}_i}\right\rangle, \end{gather}$$

where $\hat {f}$ denotes the Fourier coefficients of a quantity $f$ and $\hat {f}^*$ denotes the complex conjugate of $\hat {f}$. The above equations can be solved independently for each wavenumber in the horizontal directions, each vertical location and each velocity component.

2.3. Dataset descriptions

The reconstruction models proposed above are trained and tested with datasets obtained from the DNS of a turbulent open-channel flow. Figure 2(a) shows the configuration of the turbulent open-channel flow, which is doubly periodic in the streamwise ($x$) and spanwise ($y$) directions. In the vertical ($z$) direction, the bottom wall satisfies the no-slip condition, and the top surface satisfies the free-surface kinematic and dynamic boundary conditions. The utilised simulation parameters are listed in table 2. The domain size is $L_x\times L_y \times L_z = 4{\rm \pi} h \times 2{\rm \pi} h \times h$, with $L_z$ being the mean depth of the channel. The flow is driven by a constant mean pressure gradient $G_p$ in the streamwise direction. The friction Reynolds number is defined as ${Re}_\tau =u_\tau h/\nu =180$, where $u_\tau =\sqrt {G_p/\rho }$ is the friction velocity, $h$ is the mean depth of the channel, $\rho$ is the density and $\nu$ is the kinematic viscosity. Flows with different Froude numbers, which are defined based on the friction velocity as $Fr_\tau = u_\tau /\sqrt {gh}$ with g being the gravitational acceleration, are considered. The domain size $4{\rm \pi} h \times 2{\rm \pi} h \times h$ is sufficiently large to obtain domain-independent one-point turbulence statistics in open-channel flows at ${Re}_\tau =180$ (Wang, Park & Richter Reference Wang, Park and Richter2020). The autocorrelations studied by Pinelli et al. (Reference Pinelli, Herlina, Wissink and Uhlmann2022) indicate that a spanwise domain size of $6h$ is sufficiently large for the spanwise scales to be decorrelated; in the streamwise direction, a much larger domain length, $L_x>24h$, is needed for complete decorrelation. Owing to the higher computational cost of the deformable free-surface simulations compared with the rigid-lid simulations, here we use a smaller domain length $L_x=4{\rm \pi} h$. The normalised autocorrelation of the streamwise velocity fluctuations, $R_u(x', z) = \langle u(x') u(x+ x') \rangle /\langle u^2\rangle$, where $\langle {\cdot }\rangle$ denotes the plane average, is smaller than $0.065$ at the largest separation distance $L_x/2$ across the depth of the channel, which is small enough to allow most streamwise structures to be captured, consistent with the numerical result of Pinelli et al. (Reference Pinelli, Herlina, Wissink and Uhlmann2022).

Figure 2. (a) Set-up of the turbulent open-channel flow. The contours illustrate the instantaneous streamwise velocity $u$ of one snapshot from the simulation. (b) Sketch of the boundary-fitted curvilinear coordinate system. For illustration purposes, the domain is stretched, and the grid resolution is lowered in the plots.

Table 2. The simulation parameters of the turbulent open-channel flows. The superscript ${}^+$ denotes the quantities normalised by wall units $\nu /u_\tau$.

2.3.1. Simulation method and parameters

To simulate free-surface motions, we discretise the governing equations on a time-dependent boundary-fitted curvilinear coordinate system, $\boldsymbol {\xi } = (\xi _1, \xi _2, \xi _3)$, which is defined as (Xuan & Shen Reference Xuan and Shen2019)

(2.12a–c)\begin{equation} \xi_1 = x,\quad \xi_2 = y, \quad \xi_3 = \frac{z}{\eta(x,y,t)+h}h. \end{equation}

The above transformation maps the vertical coordinate $z\in [0, \eta +h]$ to the computational coordinate $\xi _3\in [0, h]$. The mass and momentum equations are written using the curvilinear coordinates as

(2.13a)$$\begin{gather} \frac{{\partial} (J^{{-}1} u_i)}{{\partial} t} - \frac{{\partial} (J^{{-}1} U^j_g u_i)}{{\partial} \xi^j} + \frac{{\partial} (J^{{-}1} U^j u_i)}{{\partial} \xi^j} ={-}\frac{{\partial}}{{\partial} \xi^j}\left(J^{{-}1} \frac{{\partial} \xi^j}{{\partial} x_i} p\right) + \frac{1}{\mbox{{Re}}} \frac{{\partial}}{{\partial} \xi^j}\left({J^{{-}1} g^{ij} \frac{{\partial} u_i}{{\partial} \xi^j}}\right), \end{gather}$$

(2.13b)$$\begin{gather}\frac{{\partial} U^j}{{\partial} \xi^j} = 0, \end{gather}$$

where $J=\det ( {\partial } \xi ^i / {\partial } x_j )$ is the Jacobian determinant of the transformation; $U^j=u_k ( {\partial } \xi ^j / {\partial } x_k)$ is the contravariant velocity; $U^j_g= \dot {x}_k({ {\partial } \xi ^j}/{ {\partial } x_k})$ is the contravariant velocity of the grid, where $\dot {x}_k(\xi _j)$ is the moving velocity of $\xi _j$ in the Cartesian coordinates; $g^{ij}=( {\partial } \xi ^i / {\partial } x_k)( {\partial } \xi ^j / {\partial } x_k)$ is the metric tensor of the transformation. The above curvilinear coordinate system, which is illustrated in figure 2(b), enables the discretised system to more accurately resolve the surface boundary layers.

The solver utilises a Fourier pseudo-spectral method for discretisations in the horizontal directions, $\xi _1$ and $\xi _2$, and a second-order finite-difference scheme in the vertical direction $\xi _3$. The temporal evolution of (2.13) is coupled with the evolution of the free surface $\eta (x, y, t)$; the latter is determined by the free-surface kinematic boundary condition, written as

(2.14)\begin{equation} \eta_t = w - u \eta_x - v \eta_y. \end{equation}

The above equation is integrated using a two-stage Runge–Kutta scheme. In each stage after the $\eta$ is updated, (2.13) is solved to obtain the velocity in the domain bounded by the updated water surface. This simulation method has been validated with various canonical wave flows (Yang & Shen Reference Yang and Shen2011; Xuan & Shen Reference Xuan and Shen2019) and has been applied to several studies on turbulent free-surface flows, such as research on the interaction of isotropic turbulence with a free surface (Guo & Shen Reference Guo and Shen2010) and turbulence–wave interaction (Xuan, Deng & Shen Reference Xuan, Deng and Shen2020; Xuan & Shen Reference Xuan and Shen2022). More details of the numerical schemes and validations are described in Xuan & Shen (Reference Xuan and Shen2019).

The horizontal grid, following the Fourier collocation points, is uniform. Near the bottom of the channel and the free surface, the grid is clustered with a minimum vertical spacing $\varDelta ^+_{min}=0.068$, where the superscript ${}^+$ denotes the length normalised by wall units $\nu /u_\tau$; the maximum vertical spacing at the centre of the channel is $\Delta z^{+}_{max}=2.7$. The above grid resolutions, as listed in the last four columns in table 2, are comparable to those used in the DNS of turbulent rigid-lid open-channel flows (Yao, Chen & Hussain Reference Yao, Chen and Hussain2022) and turbulent free-surface open-channel flows (Yoshimura & Fujita Reference Yoshimura and Fujita2020). Figure 3 compares the free-surface fluctuation magnitudes $\eta_{rms}$, where rms denotes the root-mean-sqaure value, in the present simulations with the literature (Yokojima & Nakayama Reference Yokojima and Nakayama2002; Yoshimura & Fujita Reference Yoshimura and Fujita2020). We find that the present results agree well with their data. The profiles of the mean velocity and Reynolds stresses are also in agreement with the results reported by Yoshimura & Fujita (Reference Yoshimura and Fujita2020) (not plotted).

Figure 3. Comparison of the root-mean-square (rms) free-surface fluctuations $\eta_{rms}$ with the literature: the present DNS results ($\blacksquare$) and the DNS results of Yoshimura & Fujita (Reference Yoshimura and Fujita2020) ($\bullet$) and Yokojima & Nakayama (Reference Yokojima and Nakayama2002) ($\blacktriangle$) are compared. The dashed line (– – –) is an approximated parameterisation of $\eta _{rms}/h$ proposed by Yoshimura & Fujita (Reference Yoshimura and Fujita2020); i.e. $\eta _{rms}/h= 4.5\times 10^{-3}{Fr}^2$ where $Fr=\bar {U}/\sqrt {gh}$ is defined based on the bulk mean velocity $\bar {U}$ and, in the present study, $\bar {U}=15.6 u_\tau$.

2.3.2. Dataset preprocessing and training

For each case, the dataset consists of a series of instantaneous flow snapshots stored at constant intervals for a period of time. The intervals and the total number of samples collected, denoted by $\Delta t u_\tau /h$ and $N$, respectively, are listed in table 2. We note that the simulations are performed on a computational grid of size $N_x \times N_y \times N_z=256\times 256 \times 128$. The simulation data are interpolated onto the input and output grids, which are uniform grids with coarser resolutions, as listed in table 1, to reduce the computational cost of the training process and the consumption of GPU memory. Although the reduced resolutions remove some fine scale structures from the training data compared with the original DNS data, we find that an increased resolution has a negligible effect on the reconstruction (see details in § 3 of the supplementary material), indicating that the present resolution is adequate for resolving the reconstructed structures. We apply random translations to the flow fields in the horizontal periodic directions, i.e. in the $x$ and $y$ directions, and obtain a total of $N_i$ samples, as listed in table 3. This process, known as data augmentation, which artificially increases the amount of data via simple transformations, can help reduce overfitting and increase the generalisation performance of NNs, especially when it is difficult to obtain big data (Shorten & Khoshgoftaar Reference Shorten and Khoshgoftaar2019). Here, owing to the high computational cost of free-surface flow simulations, we adopt this technique to expand the datasets. As discussed in § 2.1, this process can also improve the translation invariance of the CNN (Kauderer-Abrams Reference Kauderer-Abrams2017). The surface elevation in the input is normalised by $\eta _{rms}$ of each case; the surface velocities are normalised by the root-mean-square value of all three velocity components, ${\langle (u_s^2 + v_s^2 + w_s^2)/3 \rangle }^{1/2}$; the subsurface velocities are normalised by the root-mean-square magnitude of velocity fluctuations in the subsurface flow field, ${\langle (u'^2+v'^2+w'^2)/3 \rangle }^{1/2}$. This normalisation uniformly scales three velocity components because it is necessary to obtain a divergence-free velocity.

Table 3. The sampling parameters and the sizes of the datasets, including: the non-dimensionalised sampling interval $\Delta t u_\tau /h$, the total number of snapshots sampled from simulations $N$, the total number of snapshots after augmentation $N_i$, the number of snapshots in the training set $N_{train}$, the number of snapshots in the validation set $N_{{val}}$ and the number of snapshots in the test set $N_{test}$.

Each dataset is split into training, validation and test sets with $75\,\%$, $10\,\%$ and $15\,\%$ of the total snapshots, respectively. The numbers of snapshots in each set are listed in table 3. The training set is used to train the network and the validation set is used to monitor the performance of the network during the training process. The training procedure is stopped when the loss function (2.4) computed over the validation set no longer improves, an indication that the model is overfitting (see § 5 of the supplementary material for learning curves). The test set is used to evaluate the network. It should be noted that the training, validation and test sets are obtained from different time segments of the simulations such that the flow snapshots in the test set do not have strong correlations with those in the training set.

The LSE also uses the same training set as in CNN to determine the operator $\mathcal {L}$, i.e. the LSE also uses the interpolated simulation data for reconstructions. However, random translations are not applied to the dataset for the LSE reconstructions because the LSE method is inherently translation-invariant.

3.1. Instantaneous flow features

We first inspect the instantaneous flow fields to qualitatively assess the performance of the CNN and LSE methods. In this section, unless specified otherwise, we use one snapshot in the test set for the case with ${Fr}_\tau =0.08$ as an example to compare the reconstructed flow field and the ground truth from DNS. It should also be noted that hereafter, unless specified otherwise, the DNS results as the ground truth refer to the velocity fields on the uniform reconstruction grid as described in § 2.3.2.

Figure 4 compares the subsurface streamwise velocity fluctuations estimated by the CNN and LSE methods with the DNS results. Near the surface at $z/h=0.9$, the DNS result (figure 4a) shows that the regions with negative $u'$ values, i.e. the regions with low-speed streamwise velocities, are characterised by patchy shapes. The high-speed regions distributed around the low-speed patches can appear in narrow-band-like patterns. The patchy low-speed and streaky high-speed structures of the near-surface streamwise motions are consistent with the open-channel flow simulation results obtained by Yamamoto & Kunugi (Reference Yamamoto and Kunugi2011). Despite some minor differences, both the CNN and LSE methods provide reasonably accurate reconstructions. In the high-speed regions, the LSE method (figure 4c) predicts the amplitudes of the velocity fluctuations slightly more accurately than the CNN method (figure 4b). Moreover, the contours of the LSE reconstruction, especially in high-speed regions, have sharper edges, indicating that the LSE method can predict higher gradients than the CNN method. This result suggests that the LSE approach can capture more small-scale motions, thereby yielding slightly better reconstruction performance in the high-speed regions with narrow-band patterns. On the other hand, the CNN method can more accurately predict the low-speed patches, whereas the LSE method tends to underestimate the fluctuations in the low-speed regions. Overall, both methods perform similarly well in the near-surface regions.

Figure 4. Comparisons among the instantaneous streamwise velocity fluctuations $u'$ (a,d,g) obtained from the DNS results and the fields reconstructed by (b,e,h) the CNN and (c, f,i) the LSE methods. The $x$–$y$ planes at (a–c) $z/h=0.9$, (d–f) $z/h=0.6$ and (g–i) $z/h=0.3$ are plotted for the case of ${Fr}_\tau =0.08$.

Away from the free surface, close to the centre of the channel ($z/h=0.6$) (figure 4d), we can observe increases in the streamwise scales of $u'$ compared with those of the near-surface region ($z/h=0.9$). The high-speed and low-speed regions start to appear as relatively long streamwise streaks and alternate in the spanwise direction. At this height, the reconstructions provided by both the CNN and LSE methods are not as accurate as their reconstructions near the surface. However, the CNN method (figure 4e) provides a better reconstruction than the LSE method (figure 4f) in that the former better captures the streaky pattern of $u'$. The structures predicted by the LSE method generally appear shorter than those in the ground truth. Furthermore, although the reconstructed velocities provided by both methods exhibit underpredicted magnitudes, the velocity magnitudes predicted by the CNN method are clearly higher than those predicted by the LSE method and follow the DNS results more closely.

Farther away from the surface at $z/h=0.3$, the streamwise streaks become more pronounced, and the intensities of the fluctuations become stronger (figure 4g); these are the typical features of turbulent flows affected by strong shear near the wall. The CNN method (figure 4h), despite its underestimated velocity magnitudes and missing small-scale structures, can still predict distributions for the large-scale streaks that are similar to the DNS result. In contrast, the streaky structures in the reconstruction provided by LSE (figure 4i) are far less discernible, and the magnitudes of the velocity fluctuations are also weaker.

The results obtained for the vertical and spanwise velocity fluctuations are presented in Appendix B. The performance difference between the CNN and LSE methods is similar to what we have observed for the streamwise velocity component.

Based on the above observations, we find that the CNN method can reconstruct the characteristic turbulent coherent structures away from the surface, e.g. the streamwise streaks near the wall, to some extent. To further assess the capabilities of the reconstruction methods regarding instantaneous flow structures, we extract vortices from the DNS data and the reconstruction data using the $\lambda _2$-criterion (Jeong & Hussain Reference Jeong and Hussain1995). It should be noted that both the CNN and LSE models are optimised to minimise the errors in the velocities and do not guarantee similarity in the vortical structures. Considering the imperfect velocity reconstruction results observed above, we are interested in whether the two models can reconstruct physically reasonable vortical structures.

As shown in figures 5 and 6, compared with the DNS results, the reconstructed flow fields contain far fewer vortical structures. Specifically, most small-scale vortices and near-wall vortices are not captured by the reconstructions. This is another indication that the two reconstruction methods have difficulties in capturing small-scale structures. However, we note that the vortical structures recovered by the CNN method extend farther away from the surface than the vortices recovered by the LSE method. Most vortices in the LSE reconstruction are located near the surface. This difference indicates that the CNN method is more capable than the LSE method of reconstructing the flow structures away from the free surface. Moreover, the vortices that span the entire channel, which can be seen more clearly in the side view plotted in figure 6(b), are inclined towards the flow direction ($+x$-direction), similar to what can be observed in the DNS result (figure 6a). The inclination angles of these vortices, approximately $45^\circ$, are consistent with the features of the rolled-up vortices that originate from the hairpin vortices in the shear-dominated near-wall region (Christensen & Adrian Reference Christensen and Adrian2001). This result suggests that the vortical structures reconstructed by the CNN method are indeed turbulent coherent structures that would be expected in a turbulent open-channel flow.

Figure 5. Comparisons of the instantaneous vortex structures among the (a) DNS, (b) CNN reconstruction and (c) LSE reconstruction. The vortex structures are educed with the criterion $\lambda _2<0$ (Jeong & Hussain Reference Jeong and Hussain1995). The isosurface with $0.9\,\%$ of the minimum $\lambda _2$ value is plotted and is coloured by $\omega _z$.

Figure 6. Side views of the inclined vortices in the (a) DNS, (b) CNN reconstruction and (c) LSE reconstruction, as plotted in figures 5(a), 5(b) and 5(c), respectively. The vortex structures are educed by the $\lambda _2$-criterion and the isosurfaces are coloured by $\omega _z$.

It is evident from our inspection of the instantaneous velocity fields and vortical structures that the CNN method has a notable advantage over the traditional LSE method in reconstructing subsurface flows, especially away from the free surface. We find that the CNN method captures the characteristic features of the near-wall turbulence more accurately than the LSE method.

3.2. Normalised mean squared errors of reconstruction

To quantify the performance of the CNN and LSE methods, we compute a normalised mean squared error defined as

(3.1)\begin{equation} e_i(z) = \frac{\left\langle \displaystyle\iint_{x,y} {\left| \tilde{u}_i - u_i \right|}^2 \mathrm{d}\kern0.06em x\,\mathrm{d}y \right\rangle}{\left\langle \displaystyle\iint_{x,y} u_i^2 \,\mathrm{d}\kern0.06em x\,\mathrm{d}y \right\rangle}, \end{equation}

which measures the reconstruction error of each velocity component on each horizontal plane. We note that the performance evaluation is based on the test set of the corresponding case, i.e. the averaging $\langle {\cdot }\rangle$ in (3.1) is conducted on the snapshots in the test set. The reconstruction errors of the CNN and LSE models for different cases are plotted in figure 7. We find that the reconstruction performances for cases with different $Fr$ are comparable and, therefore, here, we focus on the differences between the CNN and LSE reconstructions.

Figure 7. Normalised mean squared errors (3.1) of the CNN (——) and LSE (– – –) reconstructions at $Fr_\tau =0.08$ ($\circ$, blue), $Fr_\tau =0.03$ ($\square$, orange) and $Fr_\tau =0.01$ ($\triangle$, green) for the (a) streamwise, (b) spanwise and (c) vertical velocity fluctuations between the DNS and reconstructed fields.

As shown in figure 7, the overall reconstruction accuracy of the CNN reconstructions is significantly higher than the accuracy of the LSE method. Near the free surface, both the CNN and LSE methods have low reconstruction errors, indicating that the reconstruction of the near-surface flow is relatively more accurate. The LSE even slightly outperforms the CNN in predicting the vertical velocity fluctuations. Farther away from the free surface, the errors of the LSE increase more rapidly than the errors of the CNN. As a result, the advantage of the CNN method becomes more pronounced near the centre of the channel. For example, at $z/h=0.6$, the reconstruction of $u'$ provided by the LSE method has a normalised mean squared error of $e_{u'}=0.68$, which is $45\,\%$ larger than the CNN reconstruction result of $e_{u'}=0.47$. This result is also consistent with our observation of the instantaneous flow fields (figure 4d, f).

Both the CNN and LSE reconstructions become less accurate away from the free surface, suggesting that the flow field away from the free surface is physically more difficult to predict. This behaviour is because the flow structures away from the surface have fewer direct effects on the free-surface motions and, therefore, become increasingly decorrelated with the free-surface motions. Such an accuracy decrease when reconstructing flows away from where the measurements are taken was also observed by Guastoni et al. (Reference Guastoni, Güemes, Ianiro, Discetti, Schlatter, Azizpour and Vinuesa2021) and Wang et al. (Reference Wang, Wang and Zaki2022) for wall-bounded flows when these authors used the CNN method and adjoint-variational method, respectively, to predict turbulent flows based on wall measurements. However, it should be noted that the mean squared error defined in (3.1) is an aggregated measure of the reconstruction accuracy for the different structures in a flow. As discovered in § 3.1 above, the CNN reconstructions of large-scale structures near the bottom wall still have a decent similarity with the ground truth despite that the small-scale motions are missing. In other words, the loss of the reconstruction accuracy varies for structures at different scales, and, therefore, the scale-specific reconstruction performance is investigated in § 3.3 below.

3.3. Spectral properties of reconstruction

In this section, we further evaluate the CNN and LSE reconstruction performance with respect to spatial scales based on the Fourier spectrum of the reconstructed velocities. Here, the coefficients of the Fourier modes of the velocity $u_i$, which are denoted by $\hat {u}_i$, are computed with respect to the streamwise wavenumber $k_x$ and spanwise wavenumber $k_y$, thereby representing the turbulence fluctuations at different streamwise scales and spanwise scales, respectively. A two-dimensional spectrum with respect to ($k_x$, $k_y$) can also be computed. But to conduct easier comparisons, here, we are interested only in the one-dimensional spectrum.

We first assess the reconstruction accuracy in terms of the amplitude of the spectrum by computing the relative Fourier amplitude, $\hat {r}_i = |\hat {\tilde {u}}_i|/|\hat {u}_i|$, which is the ratio of the amplitude of the reconstructed mode $|\hat {\tilde {u}}_i|$ to the ground truth $|\hat {u}_i|$. Because the Fourier amplitude spectrum is the square root of the energy spectrum, the relative amplitude spectrum can also be viewed as a measure of how much turbulent kinetic energy at each scale is recovered by the reconstructions.

Figure 8 plots the relative amplitude spectrum $\hat {r}_i$ for several representative vertical locations of the case with $Fr_\tau =0.08$. In general, for each scale, $\hat {r}_i$ of either the CNN or LSE method decreases with increasing depth, which is consistent with the conclusion drawn above from the mean squared errors. Regarding the scale-wise performance, we observe that $\hat {r}_i$ at small wavenumbers are generally higher than $\hat {r}_i$ at high wavenumbers, indicating that large-scale motions can be predicted more accurately than small-scale motions.

Figure 8. Relative amplitudes of the Fourier coefficients of reconstructed velocity fluctuations compared with the ground truth: (a,d) $\hat {r}_u$, (b,e) $\hat {r}_v$ and (c, f) $\hat {r}_w$, for the streamwise, spanwise and vertical velocity fluctuations, respectively. The spectra are evaluated at $z/h=0.9$ (——, blue), $z/h=0.6$ (– – –, orange) and $z/h=0.3$ (—${\cdot }$—, green) for the CNN reconstructions ($\circ$) and LSE reconstructions ($\square$). The first row (a–c) and the second row (d–f) plot the spectra with respect to the streamwise wavenumber $k_x$ and the spanwise wavenumber $k_y$, respectively.

Near the surface at $z/h=0.9$, the CNN method outperforms the LSE method in predicting the streamwise and spanwise velocity fluctuations at small wavenumbers, yet the LSE still has good accuracy and has a slightly better performance than CNN for vertical velocity fluctuations. At larger wavenumbers, approximately $k_x h>10$ and $k_y h>10$, the performance of the CNN method drops rapidly and becomes worse than that of the LSE method. This result indicates that, although the presented CNN model can accurately predict near-surface flow motions over a wide range of scales, it has difficulty reconstructing small-scale motions. This result is also consistent with our observation regarding the instantaneous streamwise velocity fields (figure 4a–c) where the CNN method more accurately predicts the velocities in the low-speed regions, which appear more frequently in patchy shapes that have large scales, while the LSE method can better reconstruct the narrow-band high-speed regions that are more skewed towards small scales.

Away from the surface at $z/h=0.6$ and $z/h=0.3$, the CNN method outperforms the LSE method at almost all scales, and its performance lead is even more pronounced for large-scale motions. For example, at $z/h=0.6$, the relative amplitudes of the CNN reconstructions at $k_x h < 3$ and $k_y h < 3$ for the streamwise, spanwise and vertical velocity components are at least $38\,\%$, $24\,\%$ and $19\,\%$ higher than the relative amplitudes of the LSE reconstructions, respectively.

Figure 8(d) shows an interesting phenomenon concerning the spanwise structures of the streamwise velocity fluctuations. The relative amplitude $\hat {r}_u$ for the CNN reconstructions exhibits a peak around $k_y h \approx 3$, indicating that the CNN-reconstructed streamwise velocity around this spanwise scale agrees with the DNS result better than at other scales. Moreover, the relative amplitude around $k_y h \approx 3$, being greater than $0.9$ at $z/h=0.6$ and still as high as $0.6$ at $z/h=0.3$, decreases only slowly with the increasing depth (not plotted for other examined depths). This behaviour indicates that the corresponding spanwise structures in the streamwise velocity fluctuations can be reconstructed by the CNN with good accuracy across the channel depth. This result quantitatively supports our observation that the CNN method can accurately capture the large-scale streaky pattern exhibited by the streamwise velocity fluctuations that alternate in signs in the spanwise direction (see figure 4e,h). By comparison, the performance of the LSE method at this scale is not as good, which is consistent with the lack of streaky patterns in the reconstructed streamwise velocity fields shown in figure 4( f,i).

Next, we measure the phase errors of the reconstructions. Denoting $\theta$ as the phase angle of a Fourier mode, a Fourier coefficient can be expressed as $\hat {u}_i = |\hat {u}_i|\,{\rm e}^{{\rm i}\theta }$. The phase difference between the reconstruction $\hat {\tilde {u}}_i$ and the ground truth $\hat {u}_i$ can then be calculated as $\Delta \theta =\tilde {\theta } - \theta =\mathrm {arg}(\hat {\tilde {u}}_i/\hat {u}_i) = \mathrm {arg}({\rm e}^{{\rm i}\Delta \theta })$, where $\mathrm {arg}({\cdot })$ denotes the angle of a complex value. A large phase difference means that the predicted structure is shifted by a large distance relative to the true structure. To assess the phase errors for all snapshots in the dataset, we calculate a mean phase difference between the reconstructions and the DNS results defined as

(3.2)\begin{equation} \varTheta_i = \mathrm{arg}\left(\frac{\left\langle |\hat{u}_i| \,{\rm e}^{{\rm i}\Delta\theta} \right\rangle}{\left\langle |\hat{u}_i| \right\rangle}\right). \end{equation}

In the above equation, the phase difference is weighted by the amplitude of the fluctuations such that the phase errors of stronger modes have more weights in the mean value.

The results of the mean phase errors are plotted in figure 9. It should be noted that because both the CNN and LSE methods cannot reconstruct small-scale motions with consistent accuracy, we focus only on the phase errors at low wavenumbers. It is found that in the wavenumber range $k_x h < 6$ and $k_y h< 6$, the phase errors of the CNN reconstructions are within $7.8^{\circ }$ at $z/h=0.9$ and $z/h=0.6$, and are within $16.9^{\circ }$ at $z/h=0.3$, indicating that the CNN reconstructions have low phase errors. Furthermore, when we compare the full error $|\hat {u}_i - \hat {\tilde {u}}_i|$ that accounts for both the amplitude and phase difference with the amplitude difference $|\hat {u}_i (1-\hat {r}_i)|$ at $z/h=0.3$ in the above wavenumber range $k_x h < 6$ and $k_y h< 6$, we find that the amplitude error is close to the full error, indicating that the amplitude error is dominant. In other words, although the CNN model under-estimates the fluctuating amplitudes of these large-scale structures, it can predict the phases of these fluctuation modes with relatively good accuracy. By comparison, the phase errors of the LSE reconstructions are significantly larger in the above wavenumber range, reaching $31^{\circ }$ at $z/h=0.6$ and exceeding $63^{\circ }$ at $z/h=0.3$. We also note that the LSE produces considerably larger phase errors in $u'$ at $k_y h \approx 3$, corresponding to the streaky patterns, than the CNN.

Figure 9. Mean phase errors of the reconstructed velocity fluctuations (3.2): (a,d) $\hat {\varTheta }_u$, (b,e) $\hat {\varTheta }_v$ and (c, f) $\hat {\varTheta }_w$ for the streamwise, spanwise and vertical velocity fluctuations, respectively. The phase errors are evaluated at $z/h=0.9$ (——, blue), $z/h=0.6$ (– – –, orange) and $z/h=0.3$ (—${\cdot }$—, green) for the CNN reconstructions ($\circ$) and LSE reconstructions ($\square$). The first row (a–c) and the second row (d–f) plot the errors with respect to the streamwise wavenumber $k_x$ and the spanwise wavenumber $k_y$, respectively.

The above quantitative assessments further confirm that the CNN method is more accurate than the LSE method in predicting subsurface flows from free-surface measurements. For turbulence motions across a wide range of scales, the fluctuating amplitudes and phases predicted by the CNN method are more accurate than the reconstructions provided by the LSE method. In particular, we find that the low-wavenumber spanwise structures, which correspond to streamwise streaks characterised by positive–negative $u'$ values alternating in the spanwise direction, are effectively captured by the CNN method.

The structures captured by the reconstructions also depict how subsurface turbulence influences free-surface motions. Because the reconstructions are essentially mappings from the surface quantities to the subsurface velocities (2.1), the structures that can be captured must influence the surface motions, either by directly impacting the surface or by correlating with the near-surface structures that have manifestations on the surface. This enables us to gain some understanding of the interactions between subsurface turbulence and the free surface. The CNN can capture large-scale streamwise streaks better than other near-wall structures, indicating that these streaks have strong relations with free-surface motions. Considering that the near-wall streaks are closely related to hairpin vortices that originate from the near-wall shear layer, we can infer that the correlations between the streamwise streaks and the free-surface motions are due to the evolving hairpin vortices that impact the free surface; this is consistent with the findings by Nagaosa & Handler (Reference Nagaosa and Handler2003) and Sanjou & Nezu (Reference Sanjou and Nezu2011) that the hairpin vortices can roll up from the wall and rise to the free surface. The relations between the rolled-up hairpin vortices and the surface motions are supported by the presence of tilted vortices in the flow field reconstructed by the CNN method (figure 6). Therefore, we see that the CNN reconstruction method can serve as a promising tool for revealing the dynamics of free-surface turbulence. In the next section, how the free-surface dynamics affects the reconstruction models is discussed further.