acceptance probability across the batch is greater than $0.65$, we set
$\hat{\epsilon}_{n}=1.05\hat{\epsilon}_{n}$, and
$\epsilon_{\text{shared}}=1.02\epsilon_{\text{shared}}$. If the average
acceptance probability across a batch is lower than $0.65$, we set
$\hat{\epsilon}_{n}=\hat{\epsilon}_{n}/1.05$, and
$\epsilon_{\text{shared}}=\epsilon_{\text{shared}}/1.02$. Adapting shared
parameters across the AIS forward pass violates Markov property, as the
transitions late in the MCMC chain will have a weak dependency on the earlier
transitions. However, the effect of this is minor as the step size changes are
relatively small for each run. For evaluation the AIS parameters are frozen,
so it respects the Markov property.

In FAB with prioritized buffer, we use a total of $L=8$ gradient update steps
per AIS sampling step. We initialise the buffer with $65,536$ samples from the
initialized flow-AIS combination and use a maximum buffer length of $512,000$.
We do not use any clipping for $w_{\text{correction}}$.

CRAFT specific details: For CRAFT we run our experiments with two different
setups, which we simply refer to as CRAFT (config 1) and CRAFT (config 2). The
first uses a configuration that is similar to the that of FAB in terms of the
flow architecture, and MCMC. This CRAFT model uses 5 temperatures, and uses an
auto-regressive affine flow. The second setup uses the configuration provided
in the CRAFT implementation at
https://github.com/deepmind/annealed_flow_transport, which uses neural spline
flows, with 11 temperatures (10 flow/MCMC steps) and HMC containing 10 leap-
frog steps. We use the default HMC implementation and configuration for CRAFT,
which uses a fixed step sizes of 0.3 for the first half of the intermediate
distributions and 0.2 for the rest. Both models use 3 flow layers per
temperature.

Table 7: Number of flow and target evaluations during training for each method on the Many Well Problem. For SNF and CRAFT we report the number of model forward passes - each of which contains multiple flow transport and MCMC steps.  | Number of flow/model evaluations | Number of target evaluations
---|---|---
Flow w/ ML | $10^{10}$ | $10^{7}$
Flow w/ $D_{\alpha=2}$ | $10^{10}$ | $10^{10}$
Flow w/ KLD | $10^{10}$ | $10^{10}$
Flow w/ RBD | $10^{10}$ | $10^{10}$
SNF w/ KLD | $10^{10}$ | $2.5\times 10^{11}$
CRAFT (config 1) | $2.5\times 10^{9}$ | $10^{10}$
CRAFT (config 2) | $10^{9}$ | $10^{10}$
_FAB w/o buffer_ | $10^{10}$ | $10^{10}$
_FAB w/ buffer_ | $10^{10}$ | $7.2\times 10^{9}$

###  E.4 Further Results

Figure 14 and Figure 15 show contour plots of 2D marginals from the 32
dimensional Many Well target. The contours are for pairs of variables in the
first four elements of $\mathbf{x}$ belonging to different copies of the
Double Well distribution. Each plot is obtained by scanning two variables
while the other ones are kept fixed to zero. These figures also shows samples
from each analyzed method. FAB based methods and the flow trained by maximum
likelihood place samples at each of the modes in the contour plots, while the
other methods fail to do so.

Figure 14: Target contours and model samples for 2D marginals over the first
four elements of $\mathbf{x}$ in the 32 dimensional Many Well Problem. The
plots are for pairs of variables belonging to different copies of the Double
Well distribution. For $D_{\alpha=2}(p\|q_{\theta})$ minimization with samples
from the flow, we plot results at iteration 56 of training as the final model
samples were outside of the plotting regions due to training instabilities.
Figure 15: Target contours and model samples for 2D marginals over the first
four elements of $\mathbf{x}$ in the 32 dimensional Many Well Problem for the
Resampling Base Distribution (RBD) flow model trained with the KL divergence.

##  Appendix F Alanine dipeptide experiments

###  F.1 Setup

Coordinate transformation Boltzmann generators usually do not operate on
Cartesian coordinates. In particular, Noé et al. (2019) introduced a
coordinate transformation whereby a subset of the coordinates are mapped to
internal coordinates, i.e., bond lengths, bond angles, and dihedral angles,
see also Appendix G.1 and Figure 11 in (Stimper et al., 2022). The internal
coordinates are normalized and the respective means and standard deviations
for these coordinates are computed on samples from the target distribution
generated with MD. For the remaining Cartesian coordinates, principal
component analysis is applied to the samples and the six coordinates with the
lowest variance are eliminated. The rationale behind this is that the
Boltzmann distribution is invariant in six degrees of freedom, i.e., three of
translation and three of rotation, and consequently, the corresponding
unnecessary coordinates should be removed. However, the mapping of vectors
onto a fixed set of principal components is generally neither invariant to
translations nor to rotations, and, therefore, the transformed coordinates do
not satisfy these invariances. When training Boltzmann generators with
samples, this is not a problem since the flow will learn to generate molecular
configurations for a specific rotation or translation, but when only using the
target distribution to train the flow, the model will spend some of its
capacity to sample different translational and rotational states, which is
unnecessary since they can easily be sampled independently. This will harm
performance.

Instead, we transform all Cartesian coordinates to internal coordinates, which
is a representation invariant to translations and rotations. Since we do not
want to use MD samples for our model, we use the position with the minimum
energy instead as shift and fix values for the scale when normalizing the
coordinates. The former can be easily estimated with gradient descent using
less than 100 steps. As scale parameters, we used $0.005\text{\,}\mathrm{nm}$
for the bond lengths, $0.15\text{\,}\mathrm{rad}$ for the bond angles, and
$0.2\text{\,}\mathrm{rad}$ for the dihedral angles. Coordinates which are
treated as circular are not scaled.

Model architecture We use Neural Spline Flows with rational quadratic splines
having 8 bins each. The parameter mapping is done through coupling (Durkan et
al., 2019). Dihedral angles which can freely rotate, e.g., because it is not a
double bond, are treated as periodic coordinates (Rezende et al., 2020). For
these coordinates, we use a uniform base distribution, while we pick a
Gaussian for the other ones. The flow has 12 layers and the parameter maps are
residual networks with one residual block, while the two linear layers in each
block have 256 hidden units. The flow layers were initialized in a way that
they correspond to the identity map.

One model uses a RBD, which has a residual network with two blocks having 512
hidden units per layer as acceptance function. The truncation parameter is set
to the common value $T=100$.

The models trained with FAB do AIS with 8 intermediate distributions, which
are linear interpolations between the flow and target log-densities (Neal,
2001), where the latter one is unnormalized. We use HMC with 4 Leapfrog steps
as the MCMC operator in AIS. The same procedure is used when we use AIS in the
other baseline models, see Table 10, Figure 20, Figure 21, and Figure 24. The
HMC parameters are initialized and tuned using the same procedure as the Many
Well problem, see Appendix E.3. The SNF model does additionally 10 Metropolis-
Hastings steps every two flow layers. Since this renders sampling from this
model already expensive, we do not do AIS with this model.

Dataset Since the energy surface of alanine dipeptide in an implicit solvent
has several modes of different sizes with large energy barriers between them,
see Figure 19, a very long MD simulation would be required to obtain samples
that represent the target distribution well. To get around this problem, which
is well known in computational physics and chemistry, we do a replica exchange
MD simulation (Mori & Okamoto, 2010), which is a parallel tempering technique
(Earl & Deem, 2005). We use 21 replicas starting at a temperature of
$300\text{\,}\mathrm{K}$ and increasing the temperature by an increment of
$50\text{\,}\mathrm{K}$. The replicas are exchanged every 200 iterations and
use the state at each multiple of 1000 time steps as samples. To reduce the
time it takes to generate the data, we run many of these simulations in
parallel with different seeds. Since the initial condition is always the same,
i.e., the position with minimum energy as it is usually done, we let the
system equilibrate for $2\times 10^{5}$ iterations and run the simulation
subsequently for $2\times 10^{6}$ iterations.

We split the data into 1) a training set, which consists of $10^{6}$ samples
and is only used to train the baseline flow model with ML; 2) a validation set
consisting of $10^{6}$ samples as well, which is used to find a suitable set
of hyperparameters for our experiments; and 3) a test set with $10^{7}$
samples, which is used to evaluate all models.

To generate the training data alone we had to evaluate the target distribution
and its gradients $2.3\times 10^{10}$ times, which is what we report as cost
in terms of target evaluations in Table 8.

Filtering chiral forms As mentioned in Section 4.2, alanine dipeptide is a
chiral molecule, i.e., it can occur in two different forms that are mirror
images of each other, see Figure 16. They cannot be easily converted into each
other as this would involve breaking existing and forming new bonds. In
nature, we find almost exclusively, while the D-form typically only exists in
synthetically created compounds. Hence, whenever alanine dipeptide is
considered in the literature, it is almost always as the L-form (Wu et al.,
2020; Campbell et al., 2021; Stimper et al., 2022; Dibak et al., 2022; Köhler
et al., 2022). Therefore, we aim to train our model on only this form as well.
However, since the energy of the molecule does not change when creating a
mirror image of it, models trained to approximate its Boltzmann distribution
will a priori learn to generate both forms.

Figure 16: Visualization of alanine dipeptide in its two chiral forms. In
nature, we see almost exclusively the L-form and, likewise, we aim to only
generate samples of this form.

The two forms can be separated by using the following procedure. The two
chiral forms differ by the positioning of the neighboring atoms at a chiral
center, i.e., the center carbon atom. Hence, the difference of the dihedral
angles of those atoms with respect to a fixed reference will differ relative
to each other, i.e., their difference will change. Hence, we compute this
difference and check whether it is close to a reference configuration for
which we know that it corresponds to the L-form. As reference configuration,
we use the position with minimum energy, which we already determined for the
coordinate transformation.

We use this procedure to filter the configurations generated by the flow model
during training and included only the samples that correspond to the L-form
when computing the loss. Thereby, the model learns to only generate this
chiral form.

In Appendix F.2 we will investigate a model trained on both chiral forms and
compare it to one that was only trained on the L-form.

Training All models were trained using the Adam optimizer (Kingma & Ba, 2015)
with a batch size of 1024. A learning rate of $5\times 10^{-4}$ is initially
linearly warmed up over 1000 iterations and decayed with a cosine annealing
schedule over the course of training. We use a weight decay of $10^{-5}$ and
clip gradients at a value of $10^{3}$. When training the models with the
prioritized replay buffer, we ensured a minimum buffer length of 64 batches
and started replacing the oldest samples once its length reached 512 batches.

Evaluation To evaluate the models, we draw $10^{7}$ samples from the models
with and without the use of AIS. Since there were some outliers of the
importance weights due to flow numerics, we took the $10^{3}$ highest weights
and clipped them to the lowest value in this set to compute the ESS (Koblents
& Míguez, 2015; Dibak et al., 2022). This corresponds to a fraction of
$10^{-4}$, or $0.01\%$, of the weights. For the flow trained with the
$\alpha$-divergence with $\alpha=2$, the resulting ESS is close to $10^{-4}$,
or $0.01\%$, indicating that the true ESS is even lower. We estimated the
Ramachandran plots, i.e., made a histogram of the dihedral angles $\phi$ and
$\psi$, see Figure 3a, with $100\times 100$ bins, and used them to compute the
KL divergence between the test samples and the samples from the model. We
repeated this with the reweighted samples, whereby we also used the clipped
weights. We computed the log-likelihood on the test set with the models, but
we did not do so for the SNF, as it only computes importance weights and does
not directly estimate the density.

Computational cost The two main contributors to the computational expenses
necessary to train flows approximating Boltzmann distributions are the number
of evaluations of the flow and the target. Typically, we need both the value
and the gradient and, hence, we regard obtaining them as one operation.
Moreover, the flows that we use take the same time for sampling and for pure
likelihood computation, which is why we regard them as the same operation as
well. The flow with RBD is an exception, as sampling from it is more expensive
due to learned rejection sampling being used in the base distribution. The
number of flow and target evaluations for each model and training procedure
are listed in Table 8.

In general, we trained all models using an equal number of flow evaluations,
with the exception of SNF. SNF requires a large number of target evaluations
due to the sampling layers. Because of this, we reduced the number of flow
evaluations done in total by this method.

Table 8: Number of flow and target evaluations needed to train the models. For the flow being trained with ML on MD samples, we report the number of target evaluations that are needed to generate the training dataset with MD. Our methods are marked in _italic_. | Number of flow evaluations | Number of target evaluations
---|---|---
Flow w/ ML | $2.5\times 10^{8}$ | $2.3\times 10^{10}$
Flow w/ $D_{\alpha=2}$ | $2.5\times 10^{8}$ | $2.5\times 10^{8}$
Flow w/ KLD | $2.5\times 10^{8}$ | $2.5\times 10^{8}$
RBD w/ KLD | $2.5\times 10^{8}$ | $2.5\times 10^{8}$
SNF w/ KLD | $6.0\times 10^{7}$ | $3.6\times 10^{9}$
_FAB w/o buffer_ | $2.5\times 10^{8}$ | $2.5\times 10^{8}$
_FAB w/ buffer_ | $2.5\times 10^{8}$ | $2.0\times 10^{8}$

Computational resources and runtime To generate the MD dataset, we ran the
replica exchange MD simulations on servers with an Intel Xeon IceLake-SP 8360Y
processors having 72 cores and 256 GB RAM. We used a total of 100 nodes which
needed roughly 15.7h each adding up to about 113 kCPUh.

The flow models were trained on servers with an NVIDIA A100 GPU and an Intel
Xeon IceLake-SP 8360Y processor with 18 cores and 128 GB RAM. The training
time for each model is listed in Table 9. In total, we invested around 1.02
kGPUh in the experiments. Although training the flow with ML is faster than
training with FAB, note that generating the data for ML training with MD
requires an additional 9.4 kCPUh, which would take 131h when executed on one
server.

Table 9: Runtime for training the models on the same server type as specified in the text. Our methods are marked in _italic_. | Runtime
---|---
Flow w/ ML | 13.8h
Flow w/ $D_{\alpha=2}$ | 20.0h
Flow w/ KLD | 20.0h
RBD w/ KLD | 82.5h
SNF w/ KLD | 170h
_FAB w/o buffer_ | 18.8h
_FAB w/ buffer_ | 15.7h

###  F.2 Further results

Figure 17: Marginal distribution of the dihedral angle $\phi$ for selected
models. The visualized data is the same as in Figure 3b, but here a normal
scale instead of a log scale is used for the density.

Model trained on both chiral forms To demonstrate the importance of filtering
for the L-form during training, we trained a model on both chiral forms using
FAB with a replay buffer with the same setting as in the other experiments. We
drew $10^{7}$ samples from the model and found exactly 50% of them correspond
to the L- and 50% to the D-form. As can be seen in Figure 18, the marginal
distributions of the dihedral angles $\psi$ and $\phi$ for the two forms are
mirror images of each other, while the flow model, generating both forms, is a
mixture of the two.

(a)

(b)

Figure 18: Marginal distribution of the dihedral angles $\psi$ and $\phi$ of a
model which has been trained on both chiral forms of the alanine dipeptide,
i.e., the L- and D-form. We plot the density obtained with all samples from
the model in blue and separate the samples in the two forms, yielding the pink
and green curves. The distributions of the two forms are mirror images of each
other.

The log-likelihood on the test set is $210.80$, which is $0.74$ less than the
corresponding model only trained on the L-form. This is close to
$\log(2)\approx 0.69$, i.e. the density is roughly by a factor of two lower,
confirming once more that the flow density is a mixture of the density of the
two forms.

Model performance with AIS As mentioned in the previous section, we do AIS
with all our trained models except the SNF, which already involves sampling
layers. We adopt the same AIS setting used for training the flow models with
FAB, i.e., we use 8 intermediate distributions given by linear interpolations
between the flow and the target distributions and sample from them with HMC
performing 4 Leapfrog steps. For comparison, we provide the performance when
using the untrained base distribution of our flow as proposal for AIS. The
results are shown in Table 10. When comparing tables 2 and 10, we observe that
AIS improves performance for those models that approximate the target
distribution at least fairly well. Again, the flow trained with FAB with a
replay buffer outperforms the baselines.

Table 10: ESS, and the KL divergence of the Ramachandran plots with and without reweighting (RW) for flow models when sampling from them with AIS. The results are averages over 3 runs and the standard error is given as uncertainty. Our methods are marked in _italic_ and the highest ESS or log-likelihood or lowest KL divergence values are emphasized in bold. | ESS (%) | KLD | KLD w/ RW
---|---|---|---
Base untrained | $0.013\hskip 1.42271pt\pm\hskip 1.42271pt0.000$ | $1.96\hskip 1.42271pt\pm\hskip 1.42271pt0.00$ | $8.5\hskip 1.42271pt\pm\hskip 1.42271pt0.0$
Flow w/ ML | $11.5\hskip 1.42271pt\pm\hskip 1.42271pt0.5$ | $\bm{(4.92}\hskip 1.42271pt\pm\hskip 1.42271pt\bm{2.13)\times 10^{-3}}$ | $(1.88\hskip 1.42271pt\pm\hskip 1.42271pt0.75)\times 10^{-2}$
Flow w/ $D_{\alpha=2}$ | $0.012\hskip 1.42271pt\pm\hskip 1.42271pt0.000$ | $2.52\hskip 1.42271pt\pm\hskip 1.42271pt0.06$ | $11.3\hskip 1.42271pt\pm\hskip 1.42271pt0.2$
Flow w/ KLD | $80\hskip 1.42271pt\pm\hskip 1.42271pt8$ | $2.99\hskip 1.42271pt\pm\hskip 1.42271pt0.19$ | $2.95\hskip 1.42271pt\pm\hskip 1.42271pt0.19$
RBD w/ KLD | $61\hskip 1.42271pt\pm\hskip 1.42271pt23$ | $2.84\hskip 1.42271pt\pm\hskip 1.42271pt0.05$ | $2.81\hskip 1.42271pt\pm\hskip 1.42271pt0.04$
_FAB w/o buffer_ | $70.6\hskip 1.42271pt\pm\hskip 1.42271pt0.7$ | $(2.65\hskip 1.42271pt\pm\hskip 1.42271pt0.12)\times 10^{-2}$ | $(2.45\hskip 1.42271pt\pm\hskip 1.42271pt0.91)\times 10^{-2}$
_FAB w/ buffer_ | $\bm{96.7}\hskip 1.42271pt\pm\hskip 1.42271pt\bm{0.2}$ | $\bm{(2.53}\hskip 1.42271pt\pm\hskip 1.42271pt\bm{0.41)\times 10^{-3}}$ | $\bm{(2.17}\hskip 1.42271pt\pm\hskip 1.42271pt\bm{0.15)\times 10^{-3}}$

FAB with varying values of $\bm{\alpha}$ In Appendix B.1, we derived a variant
of FAB having an arbitrary value for the $\alpha$ parameter of the
$\alpha$-divergence. Here, we want to investigate how the performance of
models trained with FAB changes as we vary $\alpha$. Therefore, we leave the
setup the same as used in Section 4.2 and only changed the $\alpha$ parameter
of FAB. The results when using a replay buffer are given in Table 11 and
without it they are reported in Table 12. We see that $\alpha=2$ outperforms
the other values in almost all performance metrics, no matter whether a replay
buffer is used or not. This justifies our theoretical arguments for picking
$\alpha=2$ empirically.

Table 11: ESS, log-likelihood on the test set, and KL divergence (KLD) of Ramachandran plots with and without reweighting (RW) for each method. Here, all models are trained with FAB using a replay buffer, but we varied the value of $\alpha$ being used, see Appendix B.1. The results are averaged over 3 runs and the standard error is given as uncertainty. Best results are emphasized in bold. $\alpha$ | ESS (%) | $\operatorname{E}_{p(\mathbf{x})}\left[\log q(\mathbf{x})\right]$ | KLD | KLD w/ RW
---|---|---|---|---
$0.25$ | $0.021\hskip 1.42271pt\pm\hskip 1.42271pt0.000$ | $-546\hskip 1.42271pt\pm\hskip 1.42271pt14$ | $5.36\hskip 1.42271pt\pm\hskip 1.42271pt0.72$ | $10.9\hskip 1.42271pt\pm\hskip 1.42271pt1.4$
$0.5$ | $0.023\hskip 1.42271pt\pm\hskip 1.42271pt0.001$ | $-606\hskip 1.42271pt\pm\hskip 1.42271pt20$ | $4.84\hskip 1.42271pt\pm\hskip 1.42271pt0.95$ | $10.4\hskip 1.42271pt\pm\hskip 1.42271pt1.9$
$1$ | $0.027\hskip 1.42271pt\pm\hskip 1.42271pt0.000$ | $60.1\hskip 1.42271pt\pm\hskip 1.42271pt9.0$ | $2.79\hskip 1.42271pt\pm\hskip 1.42271pt0.81$ | $6.73\hskip 1.42271pt\pm\hskip 1.42271pt1.76$
$1.5$ | $89.4\hskip 1.42271pt\pm\hskip 1.42271pt0.3$ | $211.49\hskip 1.42271pt\pm\hskip 1.42271pt0.02$ | $(1.44\hskip 1.42271pt\pm\hskip 1.42271pt0.90)\times 10^{-2}$ | $(1.52\hskip 1.42271pt\pm\hskip 1.42271pt0.93)\times 10^{-2}$
$2$ | $\bm{92.8}\hskip 1.42271pt\pm\hskip 1.42271pt\bm{0.1}$ | $\bm{211.54}\hskip 1.42271pt\pm\hskip 1.42271pt\bm{0.00}$ | $\bm{(3.42}\hskip 1.42271pt\pm\hskip 1.42271pt\bm{0.45)\times 10^{-3}}$ | $\bm{(2.51}\hskip 1.42271pt\pm\hskip 1.42271pt\bm{0.39)\times 10^{-3}}$
$3$ | $\bm{93.9}\hskip 1.42271pt\pm\hskip 1.42271pt\bm{1.1}$ | $211.49\hskip 1.42271pt\pm\hskip 1.42271pt0.03$ | $(2.58\hskip 1.42271pt\pm\hskip 1.42271pt0.96)\times 10^{-2}$ | $(2.19\hskip 1.42271pt\pm\hskip 1.42271pt0.82)\times 10^{-2}$
Table 12: ESS, log-likelihood on the test set, and KL divergence (KLD) of Ramachandran plots with and without reweighting (RW) for each method. Here, all models are trained with FAB without a replay buffer, but we varied the value of $\alpha$ being used, see Appendix B.1. The results are averaged over 3 runs and the standard error is given as uncertainty. Best results are emphasized in bold. $\alpha$ | ESS (%) | $\operatorname{E}_{p(\mathbf{x})}\left[\log q(\mathbf{x})\right]$ | KLD | KLD w/ RW
---|---|---|---|---
$0.25$ | $1.7\hskip 1.42271pt\pm\hskip 1.42271pt0.3$ | $198.71\hskip 1.42271pt\pm\hskip 1.42271pt0.09$ | $3.94\hskip 1.42271pt\pm\hskip 1.42271pt0.40$ | $8.45\hskip 1.42271pt\pm\hskip 1.42271pt0.20$
$0.5$ | $9.4\hskip 1.42271pt\pm\hskip 1.42271pt4.4$ | $-192\hskip 1.42271pt\pm\hskip 1.42271pt326$ | $7.70\hskip 1.42271pt\pm\hskip 1.42271pt5.59$ | $7.99\hskip 1.42271pt\pm\hskip 1.42271pt5.49$
$1$ | $15.8\hskip 1.42271pt\pm\hskip 1.42271pt4.9$ | $210.16\hskip 1.42271pt\pm\hskip 1.42271pt0.40$ | $(1.09\hskip 1.42271pt\pm\hskip 1.42271pt0.19)\times 10^{-1}$ | $(5.60\hskip 1.42271pt\pm\hskip 1.42271pt0.75)\times 10^{-2}$
$1.5$ | $34.8\hskip 1.42271pt\pm\hskip 1.42271pt3.6$ | $210.87\hskip 1.42271pt\pm\hskip 1.42271pt0.07$ | $(7.72\hskip 1.42271pt\pm\hskip 1.42271pt0.72)\times 10^{-2}$ | $(4.04\hskip 1.42271pt\pm\hskip 1.42271pt0.04)\times 10^{-2}$
$2$ | $\bm{52.2}\hskip 1.42271pt\pm\hskip 1.42271pt\bm{1.3}$ | $\bm{211.13}\hskip 1.42271pt\pm\hskip 1.42271pt\bm{0.03}$ | $\bm{(6.28}\hskip 1.42271pt\pm\hskip 1.42271pt\bm{0.33)\times 10^{-2}}$ | $\bm{(2.66}\hskip 1.42271pt\pm\hskip 1.42271pt\bm{0.90)\times 10^{-2}}$
$3$ | $24.8\hskip 1.42271pt\pm\hskip 1.42271pt4.6$ | $210.48\hskip 1.42271pt\pm\hskip 1.42271pt0.14$ | $(1.31\hskip 1.42271pt\pm\hskip 1.42271pt0.05)\times 10^{-1}$ | $(3.88\hskip 1.42271pt\pm\hskip 1.42271pt0.03)\times 10^{-2}$

Ramachandran plots Figure 19 shows the Ramachandran plot of the test set and
Figure 20, Figure 21, Figure 22, Figure 23, Figure 24, Figure 25, and Figure
26 show the Ramachandran plots of all the models we trained for the first run,
including the samples drawn from them via AIS.

Figure 19: Ramachandran plot of the test data.

(a) Flow

(b) Flow, reweighted

(c) AIS

(d) AIS, reweighted

Figure 20: Ramachandran plots of a flow trained with the
$\alpha=2$-divergence.

(a) Flow

(b) Flow, reweighted

(c) AIS

(d) AIS, reweighted

Figure 21: Ramachandran plots of a flow trained with the KL divergence.

(a) Flow

(b) Flow, reweighted

(c) AIS

(d) AIS, reweighted

Figure 22: Ramachandran plots of a flow with a resampled base distribution
trained with the KL divergence.

(a) SNF

(b) SNF, reweighted

Figure 23: Ramachandran plots of a SNF trained with the KL divergence. Since
the SNF already has layers which do sampling, we did not do AIS with it.

(a) Flow

(b) Flow, reweighted

(c) AIS

(d) AIS, reweighted

Figure 24: Ramachandran plots of the flows trained with ML.

(a) Flow

(b) Flow, reweighted

(c) AIS

(d) AIS, reweighted

Figure 25: Ramachandran plots of the flows trained with FAB without the use of
a replay buffer.

(a) Flow

(b) Flow, reweighted

(c) AIS

(d) AIS, reweighted

Figure 26: Ramachandran plots of the flows trained with FAB with the use of a
replay buffer.


We investigate the mechanisms by which a particle bed modifies the oscillatory boundary layer. We perform Euler-Lagrange simulations of an OBL over a particle bed at Reynolds numbers $\Rey_{\delta}= 200$, 400, and 800, density ratio $\rho_p/\rho_f = 2.65$, maximum Shields numbers ranging from $1.42 \times 10^{-2}$ to $7.40 \times 10^{-2}$, Galileo number $\mathrm{Ga} = 51.9$, Stokes number ranging from 3.99 to 15.97, and Kuelegan-Carpenter number ranging from $1.7 \times 10^3$ to $6.8 \times 10^3$, and showed large modulation to the underlying boundary layer. We show two mechanisms of modification to the oscillatory boundary layer: (I) the permeability of the bed allows fluid flow to penetrate into the bed, introducing a slip velocity at the bed-fluid interface, an inflection point and the expansion of the boundary layer, and (II) for $\Rey_{\delta} = 400$, and 800, particle motion creates an evolving bedform. For $\Rey_{\delta} = 200$, particles are nearly motionless and the permeability of the bed drives the alteration of the velocity profile. The initial bed form drives velocity fluctuations not seen for flow over a smooth, impermeable wall. For Reynolds 400 and 800 particle momentum is associated with the peak velocity fluctuations, indicating that the particle feedback force contributes to velocity fluctuations. We show that the coefficient of friction, computed using the bed shear stress, is decreased by nearly half for Reynolds 200 and 400, and nearly an eighth for Reynolds 800.

keyword 1, keyword 2, keyword 3

MSC Codes   (Optional) Please enter your MSC Codes here

§ INTRODUCTION

In shallow areas of the ocean, the seafloor maybe subject to large oscillating pressure gradients and strong shear forces. This causes sediment to become suspended and transported to new locations, where they are deposited as the shear force oscillates. A model flow often used to investigate this process is the oscillatory boundary layer (OBL) problem. Stokes, 1855 was amongst the first to address this problem, specifically  in the limit where viscous effects dominate and where the bottom surface is represented as a smooth flat wall. Under these assumptions, Stokes, 1855 derived analytical solutions which show the establishment of a boundary layer with characteristic thickness $\delta =\sqrt{2\nu/\omega}$, where $\nu$ is the fluid kinematic viscosity and $\omega$ the angular frequency of the oscillations. Due to the assumption of dominating viscous effects, these solutions apply only in the limit of very small Reynolds numbers $Re_\delta=U_0 \delta/\nu$, where $U_0$ is the velocity amplitude of the oscillations. Later, many researchers investigated the dynamics of oscillatory boundary layers over smooth and rough walls at larger Reynolds numbers, including when the Reynolds number is sufficiently high for turbulence to emerge [Akhavan et al., 1991, Carstensen et al., 2010, Carstensen et al., 2012, Costamagna et al., 2003, Fytanidis et al., 2021, Ghodke & Apte, 2016, Ghodke & Apte, 2018, Hino et al., 1976, Mazzuoli & Vittori, 2019, Mazzuoli et al., 2020, Mazzuoli et al., 2016, Ozdemir et al., 2014, Pedocchi et al., 2011, Salon et al., 2007, Sarpkaya, 1993, Vittori & Verzicco, 1998, Vittori et al., 2020]. However, it is unclear whether these results are applicable to seafloors. Unlike the previously studied configurations with impermeable and fixed smooth or rough walls, seafloors are made of sediment particles that together form a porous bed. Depending on the details of the flow over it, the bed may be static, with or without bedforms, and may even evolve dynamically as sediment particles saltate or get suspended by the flow  [Finn et al., 2016]. In this manuscript, we investigate the extent to which the OBL is altered by a bottom sediment bed at different Reynolds number.

There has been significant effort devoted to the characterization of the boundary layer that develops over a smooth or rough wall under the action of an oscillatory forcing.  Depending on the Reynolds number  $Re_{\delta}$, different regimes have been identified as detailed in [Vittori & Verzicco, 1998, Pedocchi et al., 2011, Akhavan et al., 1991, Ozdemir et al., 2014, Fytanidis et al., 2021]. To summarize, an OBL developing over an impermeable wall exhibits four different regimes. In experiments and simulations with wavy or rough botom walls, with Reynolds number $Re_{\delta} \leq 85$ the flow is laminar throughout the oscillation cycle [Blondeaux & Seminara, 1979, Vittori & Verzicco, 1998, Akhavan et al., 1991], and is well described the analytical solutions of Stokes, 1855. For larger Reynolds numbers up to $\sim 500$ with rough or wavy bottom walls, [Hino et al., 1976, Jensen et al., 1989] or $\sim 700$ [Blennerhassett & Bassom, 2002], the disturbed laminar regime is observed. The latter is characterized by the appearance of small amplitude perturbations superimposed upon the Stokes flow [Vittori & Verzicco, 1998]. Fytanidis et al., 2021 found that the limiting Reynolds number for this regime is dependent on background disturbances, such as initialization parameters. With increasing Reynolds number, the flow enters an intermittent turbulent regime characterized by sudden eruptions of turbulence during the decelerating portion of the oscillatory period before relaminarizing again. The exact value for the transition to intermittent turbulence is dependent upon wall imperfections, background disturbances, and initialization conditions in numerical studies [Vittori & Verzicco, 1998].
For Reynolds numbers greater than $\sim3460$, Jensen et al., 1989 show that the OBL presents sustained velocity fluctuations and a logarithmic-layer for at least 90 percent of the cycle. This criterion defines the turbulent regime.

A number of studies have considered the effects of bottom interfaces upon the structure of the OBL. Jensen et al., 1989 experimentally studied smooth beds, rough sandpaper beds, and sand beds. They found that increased roughness relative to the wave height lead to a thicker boundary layer and larger logarithmic layer. Vittori & Verzicco, 1998 investigated the transition to turbulence by simulating an OBL near an imperfect wall. The wall imperfections were modeled through waviness and roughness elements. These imperfections created instabilities that lead to transition. Xiong et al., 2020 showed that a temporary roughness element can trigger the formation of a vortex tube which causes the transition to turbulence. Additional studies of roughness effects can be found in [Ghodke & Apte, 2016, Ghodke & Apte, 2018, Mazzuoli & Vittori, 2019]. In the case of freely moving particles, Mazzuoli et al., 2016 considered a chain of particles resting on a horizontal plane. The particles are shown to form structured chains when set in motion by the OBL. If the particle bed has some waviness then there is also vortex shedding caused by the bed dunes. Finn & Li, 2016 proposed a regime map based on the Shields and Galileo number for sediment-turbulence interaction. The Shields number is a nondimensional shear stress, while the Galileo number is the ratio of gravitational forces to viscous forces upon a particle. The Shields number is dependent upon the particle parameters and the flow parameters, while the Galileo number only depends on particle parameters. This leads to a range of potential bedforms with different waviness that may emerge dependent upon the underlying fluid flow. Mazzuoli et al., 2020 performed interface-resolved direct numerical simulations of an OBL over a particle bed. They found that particle transport varies between the accelerating and decelerating portion of the flow when the Shields number is small, and is virtually the same when the Shields number is large. The limitation of interface resolved simulations due to high resolution requirements prevents the capture of several bedform periods, which may modulate turbulence. This may further alter particle transport.

To study the OBL, several models can be applied to the study of sediment transport. Namely, Euler-Euler (EE) methods, Euler-Lagrange (EL) methods, and Particle-Resolved DirectNumerical Simulations (PR-DNS). EE methods consider the particle phase as a continuum. Lee et al., 2016 developed a 3D EE method using phased based averaging. They showed good predictive capability in sheet flow and scouring simulations. Chiodi et al., 2014 derived an EE model for dense viscous flows and turbulent suspensions. They used this method to investigate sediment transport and captured the transition from bedload and suspended load. Maurin et al., 2016 compared  an EE method with an EL in simulations of turbulent bedload transport and showed agreement between the two methods over a range of Shields numbers. Kasbaoui et al., 2019 compared EE simulations with EL simulations for sedimenting particles. They showed that so long as particle clustering is weak, the two methods agree well. Additionally, EE methods are computationally efficient, but are difficult to extract information about particle momentum from. EL formulations provide a high fidelity method for simulating particle-laden flows. EL formulations track individual particles on a separate Lagrangian grid, while the fluid phase is solved on a Eulerian mesh. Capecelatro & Desjardins, 2013 used EL formulations previously to simulate dense fluidized beds, and Arolla & Desjardins, 2015 used EL methods to reproduce bed-form regimes in slurries in pipe flows. PR-DNS captures detailed information about the particle boundary layer. They have a very high degree of fidelity as they do not rely upon models for the forces on particles. However, this comes at a high computational cost, which restricts their application to smaller domains. Compared to PR-DNS, the EL formulation provides a computationally efficient method for investigating larger fluid domains, which allows meso-scales of the flow to be captured. Importantly, this allows bedforms to emerge as a result of the driving flow, with several periods of the bedform.

In this paper, we study the effect of a cohesionless particle bed on the dynamics of an OBL from the laminar regime to the onset of turbulence. The structure of the manuscript is as follows. In section <ref>, we describe the equations that dictate the dynamics of the flow and sediment particles. In section <ref>, we provide information on the computations and parameters used in this study. In order to understand the role played by the particle bed, we analyze the dynamics of an OBL over a flat plate in section <ref>, i.e., without particles, as a baseline for comparison. In the following section <ref>, we analyze the simulations with a particle-bed  and highlight the role of suspended particles and bedforms in the modification of the flow regimes. Finally, we give concluding remarks in section <ref>.

§ GOVERNING EQUATIONS

We use the volume-filtering approach of Anderson & Jackson, 1967 and Euler-Lagrange methodology of Capecelatro & Desjardins, 2013  to describe the dynamics of the sediment-laden flow. The carrier phase is an incompressible fluid with density $\rho_f$ and viscosity $\mu_{f}$. The volume-filtered Navier-Stokes equations read

\begin{eqnarray}
\frac{\partial}{\partial t} (\alpha_{f}  \rho_{f})+\nabla \cdot (\alpha_{f} \rho_{f} \boldsymbol{u}_{f}) &=& 0 \label{eq:continuity},\\
\frac{\partial}{\partial t}(\alpha_f {\rho_f} \boldsymbol{u}_{f})+\nabla \cdot (\alpha_f {\rho_f} \boldsymbol{u}_{f} \boldsymbol{u}_{f}) &=& \nabla \cdot {\boldsymbol{\tau}}+\alpha_f {\rho_f} \boldsymbol{g} - \boldsymbol{F}_{p}+ \nabla \cdot \boldsymbol{R}_{\mu}+ \boldsymbol{A}\label{eq:momentum},
\end{eqnarray}

where $\alpha_f$ is the fluid volume fraction, $\boldsymbol{u}_{f}$ is the volume-filtered fluid velocity, ${\boldsymbol{\tau}}=-p\boldsymbol{I}+\mu [\nabla \boldsymbol{u}_{f}+\nabla \boldsymbol{u}_{f}^{T}-\frac{2}{3}(\nabla \cdot \boldsymbol{u}_{f})\boldsymbol{I}]$ is the resolved fluid stress tensor [Capecelatro & Desjardins, 2013], $\boldsymbol{g}$ is the gravitational acceleration, and $\boldsymbol{F}_{p}$ is the momentum exchange between the particles and the fluid. The tensor $\boldsymbol{R}_{\mu}$ represents the so-called residual viscous stress tensor. This term arises from filtering the point-wise stress tensor. It includes sub-filter scale terms which require closure. This term is believed to be responsible for the apparent enhanced viscosity observed in viscous fluids containing suspended solid particles. For this reason Capecelatro & Desjardins, 2013 proposed a closure using an effective viscosity, which when combined with the effective viscosity model of Gibilaro et al., 2007 leads to an expression for the residual viscous stress tensor

\begin{equation}
\boldsymbol{R}_{\mu} = \mu_{f}(\alpha_{f}^{-2.8}-1)[\nabla \boldsymbol{u}_{f}+\nabla \boldsymbol{u}_{f}^{T}-\frac{2}{3}(\nabla \cdot \boldsymbol{u}_{f})\boldsymbol{I}].
\end{equation}

In order to study the dynamics of an oscillatory boundary layer, we drive the flow using the last term in equation (<ref>), which  expresses as

\begin{equation}
\boldsymbol{A} = \alpha_f \rho_f U_0 \omega \cos(\omega t)\boldsymbol{e}_x,
\end{equation}

which represents a harmonic pressure gradient forcing with angular frequency $\omega$ and amplitude $U_0$. Here, $x$ is the coordinate in the streamwise direction along the unitary vector $\boldsymbol{e}_x$, $y$ is the coordinate in the wall normal direction, $z$ is the coordinate in the cross-stream direction.

The momentum exchange term $\boldsymbol{F}_p$ contains contributions due to resolved stresses and drag force. It reads

\begin{equation}
\boldsymbol{F}_{p} = \alpha_{p} \nabla \cdot \boldsymbol{\tau}|_{p} + f_d \rho_{p} \alpha_{p} \frac{\boldsymbol{u}_{f}|_{p}-\boldsymbol{u}_{p}}{\tau_{p}}
\end{equation}

where $\alpha_p=1-\alpha_f$ is the particle volume fraction, $\rho_p$ is the particle density, $\tau_p=\rho_p d_p^2/(18\mu)$ is the response time of a particle with diameter $d_p$, and $f_d$ represents an inertial drag correction. In this study we, use the correction proposed by Tenneti et al., 2011 and derived from particle-resolved direct numerical simulations. This correction accounts for inertial and volume fraction effects.

The particles are described in the Lagrangian frame. Following Maxey & Riley, 1983, the equations of motion of a particle “$i$” are given by

\begin{eqnarray}
\frac{d \boldsymbol{x}^i_p}{d t}(t) &=& \boldsymbol{u}^i_p(t)\label{eq:part_pos}\\
m_p \frac{d \boldsymbol{u}^i_{p}}{d t}(t) &=& V_{p} \nabla \cdot \tau + m_p f_d\frac{\boldsymbol{u}_{f}(\boldsymbol{x}^i_p,t)-\boldsymbol{u}^i_{p}}{\tau_{p}}  + \boldsymbol{f}_{p}^\mathrm{col} + m_p\boldsymbol{g}\label{eq:part_vel}
\end{eqnarray}

where $\boldsymbol{x}^{p}$, $\boldsymbol{u}^{p}$, $m_p$, and $V_p$ are the particle position, velocity, mass, and volume, respectively. The term $\boldsymbol{f}_p^{col}$ represents the collisional force exerted on the particle due to particle-particle and particle-wall collisions. These collision are modeled using the soft sphere model, as detailed in Capecelatro & Desjardins, 2013.

Schematic of the configuration with a bottom sediment bed. The latter is generated in precursor runs where the particles are seeded towards the middle of the domain and allowed to settle on the bottom boundary.

Note that the governing equations (<ref>) and (<ref>) for the fluid phase are solved in both simulations with particles and without. In the latter case, $\alpha_f=1$ throughout the domain, which recovers the standard Navier-Stokes equations.

We calculate the particle volume fraction and particle velocity from the Langrangian quantities using

\begin{eqnarray}
\alpha_{p}(\boldsymbol{x},t) &=& \sum_{i=1}^{N} V_{p} g(||\boldsymbol{x}-\boldsymbol{x}^i_{p}||) \\
\alpha_{p} \boldsymbol{u}_{p}(\boldsymbol{x},t) &=& \sum_{i=1}^N \boldsymbol{u}^i_{p}(t) V_{p} g(||\boldsymbol{x}-\boldsymbol{x}^i_{p}||)
\end{eqnarray}

where $V_{p} = \pi d^{3}_{p}/6$ is the particle volume. In these equations, $g$ represents a Gaussian filter with width $\delta_{f} = 7 d_{p}$. Additional details on the computation of these terms, as well as on the solver, can be found in [Capecelatro & Desjardins, 2013].

§ NUMERICAL EXPERIMENTS

§.§ Configuration

Case     $\Rey_{\delta}$     $\rho_p/\rho_{f}$     $\Theta_\mathrm{max}$     $\mathrm{Ga}$     $\mathrm{St}_{\delta}$     $\mathrm{KC}$

1     200      –       –         –       –      –

2     400      –       –         –       –      –

3     800      –       –         –       –      –

4     200      2.6     1.42$\times 10^{-2}$     51.9     3.99      1.7$\times 10^{3}$

5     400      2.6     2.85$\times 10^{-2}$     51.9     7.99      3.4$\times 10^{3}$

6     800      2.6     7.40$\times 10^{-2}$     51.9     15.97     6.8$\times 10^{3}$

Summary of the non-dimensional parameters for the presents runs of an OBL over a smooth wall (runs 1-3) and over a sediment bed (runs 4-6).

We consider a total of six cases where the Reynolds number is progressively increased from $\Rey_\delta=200$ to $\Rey_\delta=800$. A summary of the relevant non-dimensional parameters for each run is listed in table <ref>. In order to provide a baseline for comparisons with the sediment-laden cases, we first investigate the dynamics of an oscillatory boundary layer over a smooth wall, where we vary the forcing amplitude $U_0$ to yield $\Rey_\delta=200$, 400, and 800. The first two cases are in the laminar regime, so the flow field is expected obey the Stokes solutions [Carstensen et al., 2010]. The case at $\Rey_\delta=800$ is in the intermittent turbulent regime.  Each of these cases is then compared to a companion run at the same Reynolds number with a particle bed in order to elucidate the impact of sediment motion, bedforms, and bed permeability on the flow statistics.

The presence of particles introduces additional dimensionless parameters. These are:

* the density ratio $\rho_p/\rho_f$,

* the Shields number $\Theta_\mathrm{max}=\tau_{b,\mathrm{max}}/((\rho_{p}-\rho_{f})g d_{p})$, where $\tau_{b,\mathrm{max}}$ is the maximum bed shear stress observed in simulations,

* the Galileo number $\mathrm{Ga}=d_{p}\sqrt{(\rho_{p}/\rho_{f}-1)gd_p}/\nu$,

* the Stokes number $\mathrm{St}_{\delta}=\tau_{p}U_{0}/\delta$, where $\tau_{p}$ is the particle response time,

* and the Keulegan-Carpenter number, $KC=U_{0}T/d_{p}$.

The values for each case is shown in table <ref>. The work of Finn et al., 2016 suggests regimes of particle transport determined by the combination of the Shields and Galileo numbers. Based on their work and the combination of the present parameters, case 4 falls into the “no motion regime”. Cases 5 and 6 fall in the gravitational settling regime. We expect particle motion in both of these cases, with case 6 having higher suspended sediment concentration. Keulegan-Carpenter number is large in all cases, meaning that particle drag will dominate inertial forces associated with the flow oscillations.

Figure <ref> shows a schematic of the computational domain used in these simulations. In the single-phase flow simulations, the domain and discretization follow the description in Salon et al., 2007. The dimensions are $50\delta$ in the stream-wise direction, $40\delta$ in the wall normal direction, and $25\delta$ in the span-wise direction. For $\Rey_{\delta} = 200, 400$ the discretization is $64$ points in the stream-wise direction, $256$ in the wall normal direction, and $64$ in the spanwise direction. For $\Rey_{\delta} = 800$ we use a discretization of $256$ points in the streamwise direction, $256$ in the wall-normal direction, and $128$ in the spanwise direction. In order to capture the flow structures near the wall, the mesh is stretched in the wall-normal direction using a hyperbolic tangent function such that the minimum grid spacing is $\Delta y_{min}/\delta = 0.03$. In simulations with a particle bed, we use a uniform discretization $672\times158\times 67$ such that the mesh spacing $\Delta x$ is equal to half a particle diameter. The bedheight is kept at $18.59 \delta$, which gives a depth of 25 particles. This gives a total number of particles $N = 3.1\times10^{5}$. The particle restitution coefficient is held constant at 0.9, as is the particle friction coefficient at 0.1.

The simulation timestep differs from single phase to particle-laden cases. In the single phase cases, we use a period to timestep ratio of $T/\Delta t = 5.6\times10^{3}$. The same ratio in the particle-laden case is $T/\Delta t = 2.24\times10^5$. The timestep requirement in the particle-laden case is significantly more restricted. This restriction is driven by the requirement that the bottom layer of the particles must support the weight of the particle bed above them. For more details, see [Capecelatro & Desjardins, 2013].

The protocol to initialize the simulations and gathering statistics is as follows. For the single-phase cases at $\Rey_\delta=200$ and $\Rey_\delta=400$, we initialize the simulations with a quiescent flow and simulate 10 periods of the flow. We discard the first two periods, as they have transient effects from the initialization, and perform phase averaging on the remaining eight periods. For the single-phase case at $\Rey_\delta=800$, we find that intermittent turbulence does not emerge within a reasonable simulation time when initializing from a quiescent flow. For this reason, we initialize the flow from an already turbulent flow generated at $\Rey_\delta = 1790$, that we have previously validated with the experimental data of Jensen et al., 1989. This provides sufficient initial disturbances for the flow to trigger intermittent turbulence at $\Rey_\delta = 800$. We run the simulation for several periods until the flow reaches a stationary state, which requires 2 periods for cases 1-2 and 4-6. Case 3 requires 4 periods to reach stationary state because it is initialized from a highly turbulent flow field. This is validated by confirming that fluid velocity statistics no longer vary with additional periods. After this, we collect data over the following 8 periods for the purpose of phase-averaging. For cases with a sediment-bed, we perform precursor simulations to generate a realistic poured-bed as described in <ref>. Then, we carry out simulations initialized  from quiescent flow. Monitoring changes to the flow statistics from period to period shows that transient effects due to the initialization are gone with 4 periods. We then, simulate for additional 8 periods to collect statistics.

§.§ Bed formation and bed-fluid interface

The particle bed is initialized by letting particles settle onto the bottom wall. (a) This procedure results in a volume fraction profile that is consistent with that of a poured bed. (b) The isosurface $\alpha_p=0.2$ represents a good indicator of the location of the bed-fluid interface.

To form the sediment bed in cases 4-6, we perform precursor simulations that serve to generate a realistic bed volume fraction that matches the volume fraction of a poured bed, about $63\%$  [Scott & Kilgour, 1969]. In these runs, the oscillatory forcing is turned off and the particles are initially uniformly distributed towards the middle of the domain at an average volume fraction of 40% and with small random velocity fluctuations. We then integrate the governing equations (<ref>), (<ref>), (<ref>), and (<ref>) until the particles fully settle down. Particle-particle collisions and fluid-mediated particle-particle interactions lead to the formation of the poured bed in figure <ref>.

Figure <ref> shows the average particle volume fraction $\langle \alpha_p\rangle_{xz}$ profile as a function of the wall normal distance. Note that here and onward, the notation $\langle \cdot\rangle_{xz}$ refers to ensemble and spatial averaging over the streamwise ($x$) and spanwise ($z$) directions. As anticipated, the volume fraction within the bed matches the random poured packing  [Scott & Kilgour, 1969]. It smoothly transitions to zero away from the bed. Further, we conduct the simulations with particle beds that are sufficiently deep to ensure that  the interaction between the particle bed and the turbulent flow above is captured without interference from the bottom boundary. In the present study, the sediment bed is thick by about 25 particle diameters, which corresponds to about $\sim 8.30 \delta$.

At this point, we must address the way we define the bed-fluid interface. We follow the approach of Kidanemariam & Uhlmann, 2014, where we define the bed-fluid interface using an isosurface of the particle volume fraction $\alpha_{p}=\alpha_{p,b}<0.63$. However, it is important to recognize that the choice of the isosurface $\alpha_{p,b}$ demarcating the bed-fluid interface is somewhat arbitrary since the computation of the volume fraction field $\alpha_{p}$ depends on numerical choices. For example, the shape and size of the filter kernel used to compute $\alpha_{p}$ control the width of the transition region in figure <ref>. With the filtering described in <ref>, the isosurface $\alpha_{p,b}=0.2$ provides a good indicator of the approximate location of the bed-fluid interface. We determine this by verifying that this surface lies right on top of the particles as shown in figure <ref>.

§ STRUCTURE OF AN OSCILLATORY BOUNDARY LAYER OVER A SMOOTH WALL

We now turn our attention to cases 1-3, where the bottom boundary is an impermeable smooth wall rather than a sediment bed. The reasons for this are twofold. First, analytical solutions exist for the laminar regime, which allows the validation of the computational approach. Second, these runs serve as benchmark to elucidate the changes to the flow in presence of a sediment bed.

$\omega t=0^{\circ}$
$\omega t=30^{\circ}$
$\omega t=60^{\circ}$
$\omega t=90^{\circ}$
$\omega t=120^{\circ}$
$\omega t=150^{\circ}$
Normalized spanwsie vorticity fields at $\Rey_{\delta}=400$ for a smooth, impermeable wall, at phases 0, 30, 60, 90, 120 and 150 degrees. The vorticity is arranged in laminae at all phases.

Normalized stream-wise velocity for the case of an OBL over an impermeable smooth wall at $\Rey_\delta=400$  during (a) the positive velocity half cycle, and (b) the negative velocity half cyle. Solid black lines with symbols correspond to the Stokes solution, while dashed red lines with symbols correspond to simulation data. The symbols indicate the phase: ($\square$,$\omega t = 0,180$), ($\circ$,$\omega t = 45,225$), ($\triangle$,$\omega t = 90,270$), ($\diamond$,$\omega t = 135,315$). The strong agreement between the simulated data and the Stokes solution indicates that the flow is fully laminar in this case.

Figure <ref> shows the normalized spanwise vorticity field at phases $\omega t = 0, 60, 90, 120,$ and $180^{\circ}$ for the case at Reynolds 400. The solution at Reynolds 200 shows similar vorticity structure to the Reynolds 400 and, thus, is not included here. The vorticity in these low Reynolds number cases, is organized into sheets in the near-wall region. This indicates laminar flow, and so the flow at Reynolds 200 and 400 should obey the Stokes solution,

\begin{eqnarray}
u_{f,x}/U_0 &=& \cos(\omega t)-e^{-y/\delta}\cos({\omega t - y/\delta}),\\
u_{f,y}/U_0	&=& 0.
\end{eqnarray}

To verify this, we compare the Stokes solution to vertical profiles of the phase-averaged fluid velocity from the simulation at $\Rey_\delta=400$ in figure <ref>. The agreement between the simulated data and the analytical solution is excellent, showing that the flow is indeed laminar at these Reynolds numbers. This also validates the computational approach for an impermeable smooth wall oscillatory boundary layer at low Reynolds numbers.

Scaled coefficient of friction over one period at $Re_{\delta} = 200$ and 400, for a smooth, impermeable wall. The black line correspond to the Stokes solution. Symbols correspond to the numerical solution at ($\square$,$\Rey_{\delta} = 200$) and ($\triangle$,$\Rey_{\delta} = 400$). All cases collapse onto the Stokes solution.

As an additional comparison, we compute the coefficient of friction $C_f$ defined as,

\begin{equation}
C_f=\frac{|\tau_w|}{(1/2)\rho_{f} U_0^{2}},
\end{equation}

where $\tau_w$ is the wall shear stress. Figure <ref> shows the variation of the coefficient of friction scaled by $\Rey_{\delta}$. This cancels the Reynolds number dependence of the coefficient of friction. The scaled coefficient is plotted over a period for the cases of an OBL over an impermeable smooth wall at $\Rey_\delta = 200$ and 400, alongside the Stokes solution. Here too, the agreement between numerical and Stokes solution is excellent which further demonstrates that the OBL at these Reynolds number is fully laminar.

$\omega t=0^{\circ}$
$\omega t=30^{\circ}$
$\omega t=60^{\circ}$
$\omega t=90^{\circ}$
$\omega t=120^{\circ}$
$\omega t=150^{\circ}$
Normalized spanwise vorticity field in DNS of an OBL over an impermeable smooth wall at $\Rey_{\delta}=800$, for a smooth, impermeable wall. The eruption of velocity fluctuations during the decelerating portion of the cycle ($120^{\circ}$ and $150^{\circ}$) indicates that this flow is in the intermittent turbulent regime.

Unlike the lower Reynolds number cases, figure <ref> shows significant vorticity throughout the cycle for the case at $\Rey_{\delta}=800$. Of particular interest is the range of scales seen at phase $120^{\circ}$, which is in the decelerating portion of the cycle. This eruption of velocity fluctuations, followed by partial relaminarization,  is characteristic of the intermittent turbulence regime. Similar observations were made by Jensen et al., 1989, Vittori & Verzicco, 1998, and Salon et al., 2007.

$\omega t = 0^{\circ}$

$\omega t = 30^{\circ}$

$\omega t = 60^{\circ}$

$\omega t = 90^{\circ}$

$\omega t = 120^{\circ}$

$\omega t = 150^{\circ}$

Wall scaled mean velocity profiles for $Re_{\delta} = 800$, for a smooth, impermeable wall. No logarithmic layer is observed.

The vorticity field for Reynolds 800 indicates the possibility of turbulence. To investigate this further, we plot the spatially averaged streamwise velocity in semilog, shown in figure <ref>. The velocity is normalized using wall units, i.e., the friction velocity $u_{\tau} = \sqrt{\tau_{w}/\rho_{f}}$ as velocity scale and $\nu/u_{\tau}$ as length scale. The profiles do not show a logarithmic layer. The smooth wall baseline flow does not exhibit turbulence.

§ OBL OVER A COHESIONLESS PARTICLE BED

Having characterized the dynamics in the case of a bottom smooth flat wall, we now analyze the changes that occur when the oscillatory boundary layer develops over a cohesionless particle bed.

§.§ Case at $\Rey_\delta=200$

Span-wise vorticity fields at $\Rey_{\delta}=200$, at phases 0, 30, 60, 90, 120, and 150 degrees. Small ripples in the bedform cause repeated shedding vortices to emerge.

While in the case of an OBL over a smooth impermeable wall at $\Rey_\delta=200$ the flow is devoid of any fluctuations, the presence of a sediment bed leads to
notable vortex shedding. These dynamics are visible in figure <ref> depicting the spanwise vorticity in a wall normal plane at different phases of the cycle. To highlight the bedform, figure <ref> also shows the volume fraction contour $\alpha_p=\alpha_{p,b}=0.2$ that demarcates the sediment bed-fluid interface.  The latter displays small waviness which results in some vortex shedding. These structures are small departures from the single phase case, showing a shift in the presence of a bedform, towards the disturbed laminar regime.

$\omega t = 0^{\circ}$

$\omega t = 30^{\circ}$

$\omega t = 60^{\circ}$

$\omega t = 90^{\circ}$

$\omega t = 120^{\circ}$

$\omega t = 150^{\circ}$

Profiles of the mean streamwise velocity at $\Rey_\delta=200$. The profiles are given at phases 0, 30, 60, 90, 120, and 150. Red lines with symbols correspond to simulations with a particle bed, black lines correspond to smooth wall simulations.

$\omega t = 0^{\circ}$

$\omega t = 30^{\circ}$

$\omega t = 60^{\circ}$

$\omega t = 90^{\circ}$

$\omega t = 120^{\circ}$

$\omega t = 150^{\circ}$

Profiles of the rms fluctuations of the streamwise velocity at $\Rey_\delta=200$, in presence of a particle bed. The profiles are given at phases 0, 30, 60, 90, 120, and 150. Even at Reynolds 200 velocity fluctuations are non-zero in the presence of a particle bed, whereas fluctuations are zero throughout the water columns for a smooth wall.

The mean velocity profile is altered by the cohesionless particle bed. Figure <ref> shows the average velocity profiles for Reynolds 200. To compare the two-phase velocity profile to the analytical solutions, the vertical distance is shifted by the average bed height at initialization, that is, the average $y$ location for volume fraction 0.2. In these cases, the bed does not evolve significantly and so the average bed height does not meaningfully change throughout the simulation. Reynolds 200 is close to the analytical solution, but with some important deviations. Specifically, the boundary layer is thicker, and there is an inflection point near the bed interface. Additionally, the velocity does not become zero immediately within the bed, showing flow intrusion within the bed. The flow intrusion leads to a thicker boundary layer in the presence of a bed. Finally, the maximum velocities within the boundary layer is decreased.

The presence of a particle bed leads to greater velocity fluctuations than in the smooth wall cases. Figure <ref> shows the streamwise velocity fluctuations. The rms fluctuations for a smooth wall at Reynolds 200 is zero throughout the water column. However, the presence of a particle bed leads to non-zero fluctuations. For Reynolds 200, there are peaks near 10 percent of the forcing amplitude at $1 d_{p}$ above the bed. The fluctuations show multiple peaks for all phases except 150. This indicates repeated structures at fixed heights above the bed.

The particle bed leads to a different condition at the fluid-bed interface as compared to a smooth wall. In the smooth wall case, no-slip applies at the wall, while the particle bed is porous, which leads to a slip velocity at the fluid-bed interface. This causes the bed shear stress to be reduced compared to the smooth wall case. We define the bed shear stress as the shear stress conditioned on an isosurface corresponding to the bed interface $\alpha_{p} = \alpha_{p,b}$.

\begin{equation}
\tau_{b}= \vert \vert \langle \boldsymbol{n}\cdot\boldsymbol{\tau'} \rvert_{\alpha_{p}=\alpha_{p,b}} \rangle \rvert \rvert
\end {equation}

where $\boldsymbol{n}$ is the normal vector on the isosurface $\alpha_{p}=\alpha_{p,b}$, and $\boldsymbol{\tau'}=\mu [\nabla \boldsymbol{u} +\nabla \boldsymbol{u}^T- (2/3)(\nabla\cdot\boldsymbol{u})\boldsymbol{I}]+\boldsymbol{R_\mu}$ is the deviatoric stress tensor. With the closure of Gibilaro et al., 2007, this tensor reads

\begin{equation}
\boldsymbol{\tau'}=\mu \alpha_f^{-2.8} \left(\nabla \boldsymbol{u} +\nabla \boldsymbol{u}^T- (2/3)(\nabla\cdot\boldsymbol{u})\boldsymbol{I}\right)
\end{equation}

The deviatoric stress is computed throughout the domain, and then interpolated to the bed surface. We spatially average the stress over the bed, and then take the magnitude. We define the coefficient of friction at the bed interface as

\begin{equation}
C_{f} = \frac{\tau_{b}}{(1/2)\alpha_{f} \rho_f U_{0}^{2}}
\end {equation}

The coefficient of friction scaled by $\Rey_{\delta}$ is plotted in figure <ref> as a function of time for Reynolds 200. The amplitude of the coefficient of friction is reduced by 40 percent. The phase is shifted by $30^{\circ}$. The analytical phase shift for the coefficient of friction in the laminar regime is $45^{\circ}$, showing a departure even at low Reynolds number form the analytical behavior. The reduction in coefficient of friction is caused by the change in boundary condition near the wall, which leads to an inflection point in the streamwise velocity.

Coefficient of friction for an OBL over a particle bed at $\Rey_{\delta}=200$. Red symbols corresponds to the particle bed case, blue symbols correspond to the smooth wall case. The coefficient of friction is reduced in amplitude by nearly half in the presence of a particle bed. The phase is shifted by $30^{\circ}$.

§.§ Case at $\Rey_\delta=400$

Span-wise vorticity fields at $\Rey_{\delta}=400$, at phases 0, 30, 60, 90, 120, and 150 degrees. Increasing Reynolds number leads to greater vortex shedding, with small non-repeating vortices emerging at Reynolds 400.

The disturbances caused by the small bedform increase significantly at $\Rey_\delta=400$.  Figure <ref>  shows the instantaneous spanwise vorticity at $\Rey_\delta=400$. We note significantly greater vortex shedding, particularly at phases $90^{\circ}$ and $120^{\circ}$. The vortex structures observed at those phases show a chaotic behavior. The larger structures spin off and break down into smaller scale structures. However, the range of scales is limited compared to what may be expected for a fully turbulent flow. Still, the presence of vortex shedding shows a process through which turbulence can emerge at much lower Reynolds numbers than over an impermeable smooth wall, as large vortices shed by small ripples in the particle bed break into smaller flow structures. The vorticity does not penetrate more than one Stokes thickness. Flow intrusion leads to particle transport by creating a seepage flow. Where the flow intrusion leaves the bedform, the flow exerts drag on particles, which causes particle transport, as described in Jewel et al., 2019. The outflow from the bed may set some particles into motion. The bedform, demarked by a black line in figure <ref>, changes with the phase. This shows that particles are transported by the flow.
$\omega t = 0^{\circ}$

$\omega t = 30^{\circ}$

$\omega t = 60^{\circ}$

$\omega t = 90^{\circ}$

$\omega t = 120^{\circ}$

$\omega t = 150^{\circ}$

Profiles of the mean streamwise velocity at $\Rey_\delta=400$, in the presence of a particle bed. The profiles are given at phases 0, 30, 60, 90, 120, and 150. Red lines with symbols correspond to simulations with a particle bed, black lines correspond to smooth wall simulations.

$\omega t = 0^{\circ}$

$\omega t = 30^{\circ}$

$\omega t = 60^{\circ}$

$\omega t = 90^{\circ}$

$\omega t = 120^{\circ}$

$\omega t = 150^{\circ}$

Profiles of the rms fluctuations of the streamwise velocity at $\Rey_\delta=400$, in presence of a particle bed. The profiles are given at phases 0, 30, 60, 90, 120, and 150. Velocity fluctuations show a thick single peak near the bed, indicating randomness in the fluctuations.

Figure <ref>  shows the average velocity profiles for Reynolds 400. While the smooth wall case showed good agreement with the laminar solution for $\Rey_{\delta} = 400$, the OBL over a particle bed the profile is significantly altered. The velocity profile shows flow penetration within the bed, to a depth of up to $2.5 \delta$. This leads to an inflection point within the boundary layer, and a reduction in the gradient of the streamwise velocity. For all phases, the boundary layer is significantly thicker in the presence of the bed. The maximum velocity in the boundary layer is decreased in the presence of a particle bed for phases 0 and 30 degrees.

The presence of a particle bed leads to greater velocity fluctuations than in the smooth wall cases. Figure <ref> shows the streamwise velocity fluctuations. The rms fluctuations for a smooth wall at Reynolds 400 are zero throughout the water column. However, the presence of a particle bed leads to non-zero fluctuations. Reynolds 400 shows a peak at 20 percent of the forcing amplitude at phase 90 degrees, located at $0.25 d_p$ above the bed. This indicates that the larger fluctuations are generated at the fluid-bed interface, possibly by saltating particles. Additionally, we observe single peaks near the bed interface for all phases except $\omega t = 30^{\circ}$. The thickened single peak indicates randomness in the fluctuations.

To investigate the particle motion we show the particle momentum in figure <ref>. The particle momentum is negligible $5.9 d_p$ below the bed interface. For all phases, the particle momentum begins to increase near $1.5 d_p$ below the bed interface.  Phase $0^{\circ}$ shows a momentum peak at $1.4 d_p$ below the bed interface, and decays to zero by $2.0 d_p$ above the bed surface. The momentum peak increases through phases 30, 60, and $90^{\circ}$, and the reagion of particle transport thickens. Phase $90^{\circ}$ shows a momentum peak at $0.46 d_p$ below the bed interface, and decays to zero by $5.0 d_p$. This indicates a small layer of of particles transported close to the bed surface. Phase $120^{\circ}$ shows a decreased momentum maximum, and the region of particle transport is thinned. This indicates that particles are slowing and being redeposited in the bed. Phase $150$ shows low particle momentum throughout the column, and a reversal in the direction of the momentum. The particles have been redeposited in the bed and are nearly at rest. The mobile particles at certain phases are responsible for the fluid velocity fluctuations observed in <ref>.

Particle momentum as a function of vertical direction at $\Rey_{\delta} = 400$. ($\omega t=0$,),($\omega t=30$,),($\omega t=60$,),($\omega t=90$,),($\omega t=120$,),($\omega t=150$,). Significant particle momentum is seen near the bed interface, indicating particle motion at around the bed interface.

$\omega t = 0^{\circ}$

$\omega t = 30^{\circ}$

$\omega t = 60^{\circ}$

$\omega t = 90^{\circ}$

$\omega t = 120^{\circ}$

$\omega t = 150^{\circ}$

Top-down view of particle bed at $\Rey_{\delta} = 400$ and phases 0, 30, 60, 90, 120, 150 degrees. The particles are colored by their normalized stream-wise velocity. Particles are set into motion by the fluid flow, leading to particles rolling over the bed surface.

Figure <ref> shows a top-down view of the particle bed where the particles are colored by their normalized velocities. Beginning at phase 0 degrees, particles are primarily stationary within the bed. By phase 60 degrees the fluid velocity increases and some particles are set into motion, rolling over the bed surface. By phase 90 degrees, a greater number of particles are set in motion, with some rolling over the bed, and some skipping off the bed surface. After this the fluid velocity decreases, as does the particle velocity, until most particles rest in the bed at phase 150 degrees.

$\omega t = 0^{\circ}$

$\omega t = 30^{\circ}$

$\omega t = 60^{\circ}$

$\omega t = 90^{\circ}$

$\omega t = 120^{\circ}$

$\omega t = 150^{\circ}$

Reynolds 400 bedform height deviations. Small ripples rise and fall below the average bed height.

At Reynolds 400 the bedform evolves due to particle transport. We use isosurfaces of the fluid-bed interface, $\alpha_{p,b} = 0.2$, shown in figure <ref>, to visualize the bedform. The bedform is colored by the deviation in local bed height from the average bed height at the corresponding phase. At phase $\omega t = 0^{\circ}$ to $30^{\circ}$, when the majority of particles are at rest in the bed, we see a relatively flat bed, close to the average bed height throughout the domain. As particles are set into motion small dunes emerge at phase $\omega t = 60^{\circ}$. These dunes grow in height, to a maximum of one Stokes thickness above or below the average bed height at phase $\omega t = 90^{\circ}$. After this the particles begin to be redeposited in the bed, and the dunes are reduced. This continues through phase $\omega t = 150^{\circ}$, which shows only small deviations from the average bed height, and is relatively flat. The presence of hills and dunes at the bed interface may contribute to the increased fluctuations observed at the fluid-bed interface.
The coefficient of friction scaled by $\Rey_{\delta}$ is plotted in figure <ref> as a function of time for 400. Compared to the particle free case, Reynolds 400 shows a reduction in coefficient of friction of about 40 percent with a particle bed. The phase shift is $30^{\circ}$. The analytical phase shift for the coefficient of friction in the laminar regime is $45^{\circ}$, showing a departure even at low Reynolds number from the analytical behavior. As with the case at $\Rey_{\delta} = 200$, the reduction in coefficient of friction is caused by the change in boundary condition at the fluid-bed interface.

Coefficient of friction for an OBL over a particle bed at $\Rey_{\delta}=400$. Red symbols corresponds to the particle bed case, blue symbols correspond to the smooth wall case. The coefficient of friction is reduced in amplitude nearly half in the presence of a particle bed. The phase is shifted by $30^{\circ}$.

§.§ Case at $Re_{\delta}= 800$

$\omega t = 0^{\circ}$

$\omega t = 30^{\circ}$

$\omega t = 60^{\circ}$

$\omega t = 90^{\circ}$

$\omega t = 120^{\circ}$

$\omega t = 150^{\circ}$

Top-down view of the particle bed at $\Rey_{\delta} = 800$ and phases 0, 30, 60, 90, 120, 150 degrees. Particles are colored by their normalized stream-wise velocity. Particles are set into motion and lifted into the flow before being redeposited in the bed, driving bedform evolution.

At Reynolds 800 the fluid flow is able to set a greater number of particles into motion compared to Reynolds 400. Figure <ref> shows a top-down view of the particle bed where the particles are colored by their normalized velocities. The particles first saltate across the bed surface, and then are suspended. Later, as fluctuations decay, the particles are deposited in the bed. The particles are at rest in the bed at phase $0$. As the velocity of the flow increases a layer of particles are set in motion, which can be seen at phase $30^{\circ}$. Most of the particles in motion at this point saltate over the top of the bed. As the velocity continues to increase to phase $60^{\circ}$ more particles are set in motion. At this phase a large number of the particles are lifted into the flow and transported through the domain. This continues through phase $120^{\circ}$, where the fluid velocity has begun to decrease. As the velocity decreases to phase $150^{\circ}$ less particles are in motion, most having been redeposited in the bed. Saltation is once again the dominant mode of particle transport. This cycle continues in the opposite direction in the following half cycle, which is not visualized here. The process of particle transport leads to the formation of ripples in the bed for some phases, while at others, such as phase 90 a transport layer emerges at the top of the bed.

Particle momentum is shown in figure <ref> at phases 0, 30, 60, 90, 120, and $150^{\circ}$. For phase $\omega t = 0$ there is small levels of motion from $-7.5 \delta$ to $-2 \delta$ deep into the bed. The particle momentum then increases to a maximum at $-0.5 \delta$. Afterwards, the particle momentum decreases to zero far from the bed. The particle momentum maximum increases through phases 30, 60 and $90^{\circ}$, and the region of particle transport thickens. The phase $\omega t = 90$ shows negligible particle momentum within the bed up to $-3 \delta$ from the average bed height. The momentum peaks at $-0.5 \delta$ for phase $90^{\circ}$, and then  decreases to zero near $4 \delta$ above the initial bed height. This shows a large number of particles moving with significant velocity at phase 90. The maximum particle momentum decreases for phase 120, and the region of particle transport is thinned. This indicates particle deposition. Phase 150 shows a further decrease in the momentum maximum, to a similar level as in phase 0. The maximum is higher in the flow, however, with small particle transport above the fluid-bed interface. Compared to Reynolds 400, this shows continued particle transport, rather than particles being largely at rest in the bed. Particles are lifted out of the bed and transported in the flow, within a region of six stokes thickness.

Particle momentum as a function of vertical direction at $\Rey_{\delta} = 800$. ($\omega t=0$,),($\omega t=30$,),($\omega t=60$,),($\omega t=90$,),($\omega t=120$,),($\omega t=150$,). Particle momentum shows small particle transport near the bed surface for phases 0 and $180^{\circ}$, while phases 90 and $270^{\circ}$ show high particle transport in a layer of $6 \delta$.

$\omega t = 0^{\circ}$

$\omega t = 30^{\circ}$

$\omega t = 60^{\circ}$

$\omega t = 90^{\circ}$

$\omega t = 120^{\circ}$

$\omega t = 150^{\circ}$

Reynolds 800 bedform height deviations at phases $0$, $30$, $60$, $90$, $120$, $150$. Large dunes rise and fall below the average bed height.

Figure <ref> shows snapshots of the bed interface. We note that the bed interface forms into dunes that are much larger than those seen at Reynolds 200-400. The dunes expand and roll through the domain. Large numbers of particles are set into motion at certain phases. Those particles are transported through the flow before being deposited into dunes. This is different from the Reynolds 200 case, where the particles only move small amounts and the bedform itself is static. This creates a lower boundary that greatly varies throughout the simulation, and is different from both the smooth wall case and the low Reynolds number cases. That is, the bed interface becomes much more wavy, and at certain phases much of the bed interface is sliding across the bed. At other phases there are particles lifted out of the and carried above the interface itself.

Span-wise vorticity at $\Rey_{\delta}$ and phases 0, 30, 60, 90, 120, and 150 degrees. The bedform shifts into dunes at various phases. The shedding vortices create a large range of scales. The eddies penetrate the bed interface.

At Reynolds 800, the dunes in the bed-form cause greatly increased vortex shedding which leads to stronger turbulence sustained over a larger portion of the period. Figure <ref> shows the vorticity at Reynolds 800, overlaid with the isocontour $\alpha_p=\alpha_{p,b}$ denoting the location of the bed-fluid interface. The field drastically deviates from the single phase simulations, with both large scale and small scale structures. The bed generates vorticity through the small dunes as seen in lower Reynolds numbers. Whereas at low Reynolds numbers the dunes caused large shedding vortices to form in the region just above the bed, Reynolds 800 shows large vortical structures throughout the water column. In the near bed region, at phases $60$, $90$, and $120^{\circ}$, there are small scale structures.  While the smooth wall case also saw significant fluctuations and a range of scales, the fluctuations were not as significant throughout the water column as they are here. The vorticity penetrates deeper into the bed than in the Reynolds 100-400 cases. This indicates greater flow intrusion, large enough at this Reynolds number to cause significant particle transport, as described in Jewel et al., 2019. The penetration is most notable at phases $60$, $90$, and $120^{\circ}$. These correspond to the phases with greatest particle transport, as seen above.

Coefficient of friction for an OBL over a particle bed at $\Rey_{\delta}=800$. Red symbols represent particle-laden flow, blue symbols represent single phase flow. The particle bed significantly reduces coefficient of friction.

The coefficient of friction is shown in figure <ref>, scaled by $\Rey_{\delta}$. The turbulent eruption can be clearly seen in the single phase case as the dramatic increase in coefficient of friction at phase $75^{\circ}$. In the case of a particle bed there is a clear eruption of turbulence which occurs earlier in the cycle compared to the single phase case, at phase $45^{\circ}$. This represents a phase shift in the turbulent portion of the cycle of $30^{\circ}$. The eruption occurs during the accelerating portion of the cycle, rather than in the decelerating portion as in the smooth wall case. Additionally, where the presence of turbulence leads to an increase in the coefficient of friction compared to the laminar cases in the smooth wall configuration, in the particle bed case the coefficient of friction is decreased in the turbulent case compared to the lower Reynolds number cases. This is due to the flow penetration and transport layer, which reduces the rate the fluid velocity goes to zero within the bed.

$\omega t = 0^{\circ}$

$\omega t = 30^{\circ}$

$\omega t = 60^{\circ}$

$\omega t = 90^{\circ}$

$\omega t = 120^{\circ}$

$\omega t = 150^{\circ}$

Streamwise velocity profile for Reynolds 800. There is no logarithmic layer present.

The stream-wise velocity is plotted in the figure <ref>, in semi-logarithmic format. As before, the $y$ direction is re-centered by the average bed height for that phase. In all six cases there is clear flow penetration shown by the nonzero velocity at the bed surface. Indeed, for $\omega t = 60$, $90$, and $120^{\circ}$ the wall scaled streamwise velocity is 40 at the bed interface. While we don't observe turbulence at any of the phases shown here, the high flow penetration and

$\omega t = 0^{\circ}$

$\omega t = 30^{\circ}$

$\omega t = 60^{\circ}$

$\omega t = 90^{\circ}$

$\omega t = 120^{\circ}$

$\omega t = 150^{\circ}$

Profiles of the mean streamwise velocity at $\Rey_\delta=800$. The profiles are given at phases 0, 30, 60, 90, 120, and 150. Red lines with symbols correspond to simulations with a particle bed, black lines correspond to smooth wall simulations.

$\omega t = 0^{\circ}$

$\omega t = 30^{\circ}$

$\omega t = 60^{\circ}$

$\omega t = 90^{\circ}$

$\omega t = 120^{\circ}$

$\omega t = 150^{\circ}$

Profiles of the rms fluctuations of the streamwise velocity at $\Rey_\delta=800$. The profiles are given at phases 0, 30, 60, 90, 120, and 150. Velocity fluctuations show a thick single peak near the bed, indicating randomness in the fluctuations.

We also show the velocity profiles scaled using the outer scaling, $U_0$ and $\delta$, seen in <ref>. This highlights the thickened boundary layer, with a much lower gradient in the bed case than in the smooth wall case. There is an inflection point in the fluid velocity profile. Additionally, we note flow penetration, up to $5 \delta$ below the bed interface. The permeability of the particle bed, as well as the particle suspension, represents a very different condition as compared to the no slip and no penetration conditions present for a smooth wall.

As in the low Reynolds cases, the presence of a particle bed leads to greater velocity fluctuations than in the smooth wall cases. Figure <ref> shows the stream-wise velocity fluctuations for Reynolds 800. Phase $\omega t = 90^{\circ}$ shows a thick single peak, at the bed interface. The thickness indicates randomness in the fluctuations, rather than the periodic fluctuations observed at $\Rey_{\delta} = 200$. When a particle bed is present, the peak is thicker than in the smooth wall case, which indicates that the random fluctuations are present higher in the water column. At phase $\omega t = 0^{\circ}$ shows a thick region of fluctuations that extends throughout the water column. This is compared to the smooth wall case, where a similar region of fluctuations extends $12 \delta$ away from the wall. While the velocity profiles above show that the flow is no longer turbulent at these phase $\omega t = 0^{\circ}$ fluctuations remain throughout the water column. The high fluctuations and thicker region of fluctuations indicate a greater level of turbulence in the presence of a particle bed as compared to the smooth wall case.

§ CONCLUSIONS

We investigate the mechanisms by which a particle bed modifies the oscillatory boundary layer. Using Euler-Lagrange simulations of an OBL over a particle bed at Reynolds numbers $\Rey_{\delta}= 200$, 400, and 800, density ratio $\rho_p/\rho_f = 2.65$, maximum Shields numbers ranging from $1.42 \times 10^{-2}$ to $7.40 \times 10^{-2}$, Galileo number $\mathrm{Ga} = 51.9$, Stokes number ranging from 3.99 to 15.97, and Kuelegan-Carpenter number ranging from $1.7 \times 10^3$ to $6.8 \times 10^3$, and showed large modulation to the underlying boundary layer. Additionally, we performed simulations with a smooth wall to provide a benchmark of the OBL for comparison with the OBL over a particle bed. We observed two mechanisms of modulation to the boundary layer: (I) the permeability of the bedform allows fluid flow to penetrate into the bed, causing an expansion of the boundary layer, and (II) resuspension of particles at high Reynolds number leads to significant fluctuations in the fluid phase, as a result of the feedback force exerted by the particle.

For the case at $\Rey_{\delta} = 200$, there is  significant flow modulation, despite a stationary bed. In comparison to the case of an OBL over a smooth, impermeable wall, we observe significant velocity fluctuations, changes to the velocity profile, and reduction in coefficient of friction. These changes are due to the permeability of the bed and the rough wavy fluid-bed interface. The bed permeability allows flow penetration, which creates a slip velocity at the bed-fluid interface, an inflection point in the velocity profile, and an expansion of the boundary layer. The alteration of the velocity profile drives a reduction in the coefficient of friction by nearly half. This, in turn leads to reduced Shields numbers compared to predictions from cases with an impermeable wall, be it smooth, rough, or wavy.

The case at $\Rey_{\delta} = 400$ has a greater extent of flow penetration, and a thicker boundary layer. The coefficient of friction is reduced by nearly half. In this case, the particle transport is cyclic. The particles are at rest when the fluid velocity is minimum. As the latter increases, the particles are set into motion and the saltate over the particle bed. This leads to an evolving bedform and the emergence of dunes. These dunes contribute to vortex shedding, in contrast to the smooth wall case, where no vortex shedding is observed. The particle transport also contributes to the increased fluid velocity fluctuations through the particle feedback force, as evidenced by the correspondence of particle momentum and fluid velocity fluctuation peaks.

For $\Rey_{\delta} = 800$, we see more particles set in motion over a greater portion of the period and higher above the bed surface. The maximum particle momentum is near 3 times as at $\Rey_{\delta} = 400$. Particle transport is still cyclic. The particles are at rest when the fluid velocity is minimum. As the fluid velocity increases, particles are set into motion and saltate over the bed, as at $\Rey_{\delta} = 400$. This leads to dune formation as before, but then gives way to particle resuspension. The particles are transported in a layer of the flow above the particle bed. As the fluid velocity decreases, particles slow and are redeposited in the bed. This creates a boundary layer that is much thicker than the smooth wall case, an inflection point in the velocity profile and a lower slope in the velocity profile in the boundary layer. The alteration to the fluid velocity profile causes a drastic reduction in the coefficient of friction, by nearly an eighth. Instead of the increased fluctuations accompanying an increase in the coefficient of friction, as in the smooth, impermeable wall case, the particle bed case shows a reduction in the coefficient of friction. This is an important departure from the known OBL behavior over an impermeable wall for the intermittent turbulent regime. The maximum fluctuations are associated with maximum particle transport. The region of particle transport is associated with the region of velocity fluctuations, showing that the particle feedback force contributes to the increased fluctuations.

These results show that the permeability of the particle bed introduces additional physics which are not well captured by impermeable wall treatments, even if the wall has some degree of roughness and waviness. Predictions of sediment transport may be inaccurate if the bed permeability is ignored, as inclusion of permeability leads to lower coefficient of friction, and thereby the Shields number.

[Acknowledgements]Computing resources through ACCESS award EES230041 and Research Computing at Arizona State University are gratefully acknowledged.

[Funding]This research received no specific grant from any funding agency, commercial or not-for-profit sectors.

[Declaration of interests]The authors report no conflict of interest.

[Author ORCIDs]
M. H. Kasbaoui, https://orcid.org/0000-0002-3102-0624;

natexlab#1#1#1#1 #1#1 #1#1#1#1#1#1
#1#1#1#1 #1#1 #1#1

[Akhavan et al., 1991]
Akhavan, R., Kamm, R. D. & Shapiro, A. H.
1991a  An investigation of transition to
turbulence in bounded oscillatory Stokes flows Part 1.
Experiments.  Journal of Fluid Mechanics  225,

[Akhavan et al., 1991]
Akhavan, R., Kamm, R. D. & Shapiro, A. H.
1991b  An investigation of transition to
turbulence in bounded oscillatory Stokes flows Part 2. Numerical
simulations.  Journal of Fluid Mechanics  225,  423–444.

[Anderson & Jackson, 1967]
Anderson, T. B. & Jackson, R. 1967  Fluid
Mechanical Description of Fluidized Beds. Equations of
Motion.  Industrial & Engineering Chemistry Fundamentals
6 (4),  527–539.

[Arolla & Desjardins, 2015]
Arolla, S. K. & Desjardins, O. 2015  Transport
modeling of sedimenting particles in a turbulent pipe flow using
Euler–Lagrange large eddy simulation.  International Journal of
Multiphase Flow  75,  1–11.

[Blennerhassett & Bassom, 2002]
Blennerhassett, P. J. & Bassom, A. P. 2002  The
linear stability of flat Stokes layers.  Journal of Fluid Mechanics
464,  393–410.

[Blondeaux & Seminara, 1979]
Blondeaux, P. & Seminara, G. 1979  Transizione
incipiente al fondo di un'onda di gravità.  Atti della Accademia
Nazionale dei Lincei. Classe di Scienze Fisiche, Matematiche e Naturali.
Rendiconti  67 (6),  408–417.

[Capecelatro & Desjardins, 2013]
Capecelatro, J. & Desjardins, O. 2013  An
Euler–Lagrange strategy for simulating particle-laden flows.
Journal of Computational Physics  238,  1–31.

[Carstensen et al., 2010]
Carstensen, S., Sumer, B. M. & Fredsøe, J. 2010
Coherent structures in wave boundary layers. Part 1. Oscillatory
motion.  Journal of Fluid Mechanics  646,  169–206.

[Carstensen et al., 2012]
Carstensen, S., Sumer, B. M. & Fredsøe, J. 2012
A note on turbulent spots over a rough bed in wave boundary layers.
Physics of Fluids  24 (11),  115104.

[Chiodi et al., 2014]
Chiodi, Filippo, Claudin, Philippe & Andreotti, Bruno
2014  A two-phase flow model of sediment transport: Transition from
bedload to suspended load.  Journal of Fluid Mechanics  755,
561–581,  arXiv: 1401.0807.

[Costamagna et al., 2003]
Costamagna, P., Vittori, G. & Blondeaux, P. 2003
Coherent structures in oscillatory boundary layers.  Journal of
Fluid Mechanics  474,  1–33.

[Finn & Li, 2016]
Finn, J. R. & Li, M. 2016  Regimes of
sediment-turbulence interaction and guidelines for simulating the multiphase
bottom boundary layer.  International Journal of Multiphase Flow
85,  278–283.

[Finn et al., 2016]
Finn, J. R., Li, M. & Apte, S. V. 2016  Particle
based modelling and simulation of natural sand dynamics in the wave bottom
boundary layer.  Journal of Fluid Mechanics  796,

[Fytanidis et al., 2021]
Fytanidis, D. K., García, M. H. & Fischer, P. F.
2021  Mean flow structure and velocity–bed shear stress maxima
phase difference in smooth wall, transitionally turbulent oscillatory
boundary layers: Direct numerical simulations.  Journal of Fluid
Mechanics  928,  A33.

[Ghodke & Apte, 2016]
Ghodke, C. D. & Apte, S. V. 2016  DNS study of
particle-bed-turbulence interactions in an oscillatory wall-bounded flow.
Journal of Fluid Mechanics  792,  232–251.

[Ghodke & Apte, 2018]
Ghodke, C. D. & Apte, S. V. 2018  Roughness effects
on the second-order turbulence statistics in oscillatory flows.
Computers & Fluids  162,  160–170.

[Gibilaro et al., 2007]
Gibilaro, L. G., Gallucci, K., Di Felice, R. &
Pagliai, P. 2007  On the apparent viscosity of a fluidized
bed.  Chemical Engineering Science  62 (1),  294–300.

[Hino et al., 1976]
Hino, M., Sawamoto, M. & Takasu, S. 1976
Experiments on transition to turbulence in an oscillatory pipe flow.
Journal of Fluid Mechanics  75 (2),  193–207.

[Jensen et al., 1989]
Jensen, B. L., Sumer, B. M. & Fredsøe, J. 1989
Turbulent oscillatory boundary layers at high Reynolds numbers.
Journal of Fluid Mechanics  206,  265–297.

[Jewel et al., 2019]
Jewel, A., Fujisawa, K. & Murakami, A. 2019
Effect of seepage flow on incipient motion of sand particles in a bed
subjected to surface flow.  Journal of Hydrology  579,

[Kasbaoui et al., 2019]
Kasbaoui, M. Houssem, Koch, Donald L. & Desjardins,
Olivier 2019  Clustering in Euler–Euler and
Euler–Lagrange simulations of unbounded homogeneous particle-laden
shear.  Journal of Fluid Mechanics  859,  174–203.

[Kidanemariam & Uhlmann, 2014]
Kidanemariam, A. G. & Uhlmann, M. 2014
Interface-resolved direct numerical simulation of the erosion of a
sediment bed sheared by laminar channel flow.  International Journal of
Multiphase Flow  67,  174–188.

[Lee et al., 2016]
Lee, Cheng-Hsien, Low, Ying Min & Chiew, Yee-Meng
2016  Multi-dimensional rheology-based two-phase model for sediment
transport and applications to sheet flow and pipeline scour.  Physics of
Fluids  28 (5),  053305.

[Maurin et al., 2016]
Maurin, Raphael, Chauchat, Julien & Frey, Philippe
2016  Dense granular flow rheology in turbulent bedload transport.
Journal of Fluid Mechanics  804,  490–512,  arXiv:

[Maxey & Riley, 1983]
Maxey, Martin R. & Riley, James J. 1983  Equation of
motion for a small rigid sphere in a nonuniform flow.  The Physics of
Fluids  26 (4),  883–889.

[Mazzuoli et al., 2020]
Mazzuoli, M., Blondeaux, P., Vittori, G., Uhlmann, M.,
Simeonov, J. & Calantoni, J. 2020  Interface-resolved
direct numerical simulations of sediment transport in a turbulent oscillatory
boundary layer.  Journal of Fluid Mechanics  885.

[Mazzuoli et al., 2016]
Mazzuoli, M., Kidanemariam, A. G., Blondeaux, P.,
Vittori, G. & Uhlmann, M. 2016  On the formation of
sediment chains in an oscillatory boundary layer.  Journal of Fluid
Mechanics  789,  461–480.

[Mazzuoli & Vittori, 2019]
Mazzuoli, M. & Vittori, G. 2019  Turbulent spots in
an oscillatory flow over a rough wall.  European Journal of Mechanics -
B/Fluids  78,  161–168.

[Ozdemir et al., 2014]
Ozdemir, C., Hsu, T. J. & Balachandar, S. 2014
Direct numerical simulations of transition and turbulence in
smooth-walled Stokes boundary layer.  Physics of Fluids  26.

[Pedocchi et al., 2011]
Pedocchi, F., Cantero, M. I. & García, M. H.
2011  Turbulent kinetic energy balance of an oscillatory boundary
layer in the transition to the fully turbulent regime.  Journal of
Turbulence  12,  N32.

[Salon et al., 2007]
Salon, S., Armenio, V. & Crise, A. 2007  A
numerical investigation of the Stokes boundary layer in the turbulent
regime.  Journal of Fluid Mechanics  570,  253–296.

[Sarpkaya, 1993]
Sarpkaya, T. 1993  Coherent structures in oscillatory
boundary layers.  Journal of Fluid Mechanics  253,

[Scott & Kilgour, 1969]
Scott, G. D. & Kilgour, D. M. 1969  The density of
random close packing of spheres.  Journal of Physics D: Applied Physics
2 (6),  863–866.

[Stokes, 1855]
Stokes, G. G. 1855  On the effects of internal friction of
fluids on the motion of pendulums.  Transactions of the Cambridge
Philosophical Society  9.

[Tenneti et al., 2011]
Tenneti, S., Garg, R. & Subramaniam, S. 2011
Drag law for monodisperse gas–solid systems using particle-resolved
direct numerical simulation of flow past fixed assemblies of spheres.
International Journal of Multiphase Flow  37 (9),

[Vittori et al., 2020]
Vittori, G., Blondeaux, P., Mazzuoli, M., Simeonov, J.
& Calantoni, J. 2020  Sediment transport under oscillatory
flows.  International Journal of Multiphase Flow  133,

[Vittori & Verzicco, 1998]
Vittori, G. & Verzicco, R. 1998  Direct simulation of
transition in an oscillatory boundary layer.  Journal of Fluid
Mechanics  371,  207–232.

[Xiong et al., 2020]
Xiong, Chengwang, Qi, Xiang, Gao, Ankang, Xu, Hui,
Ren, Chengjiao & Cheng, Liang 2020  The bypass
transition mechanism of the Stokes boundary layer in the intermittently
turbulent regime.  Journal of Fluid Mechanics  896,  A4.