# Non-intrusive, transferable model for coupled turbulent channel-porous media
flow based upon neural networks

Xu Chu  Cluster of Excellence SimTech, University of Stuttgart, Germany
Sandeep Pandey<EMAIL_ADDRESS>Institute of Thermodynamics and
Fluid Mechanics, Technische Universität Ilmenau, D-98684 Ilmenau, Germany

###### Abstract

Turbulent flow over permeable interface is omnipresent featuring complex flow
topology. In this work, a data-driven, end-to-end machine learning model has
been developed to model the turbulent flow in porous media. For the same, we
have derived a non-linear reduced order model with a deep convolution
autoencoder network. This model can reduce highly resolved spatial dimensions,
which is a prerequisite for direct numerical simulation, by 99%. A downstream
recurrent neural network has been trained to capture the temporal trend of
reduced modes, thus it is able to provide future evolution of modes. We
further evaluate the trained model’s capability on a newer dataset with a
different porosity. In such cases, fine tuning could reduce the efforts (up to
two-order of magnitude) to train a model with limited dataset (10%) and
knowledge and still show a good agreement on the mean velocity profile.
Especially, the fine-tuned model shows a better agreement in the porous domain
than the channel and interface areas indicating the topological feature is
less challenging for training than the multi-scale nature of the turbulent
flows. Leveraging the current model, we find that even quick fine-
tuning—achieving an impressive order-of-magnitude reduction in training time
by approximately $\mathcal{O}$(102)—still results in effective flow
predictions. This promising discovery encourages the fast development of a
substantial amount of data-driven models tailored for various types of porous
media. The diminished training time substantially lowers the computational
cost when dealing with changing porous topologies, making it feasible to
systematically explore interface engineering with different types of porous
media. Overall, the data-driven model shows a good agreement, especially for
the porous media which can aid the DNS and reduce the burden to resolve this
complex domain during the simulations. The fine tuning is able to reduce the
training cost significantly and maintain an acceptable accuracy when a new
flow condition comes into play.

##  I Introduction

Turbulent flow across permeable interfaces is a phenomenon that is pervasive
in both natural environments and engineering contexts, manifesting in diverse
applications ranging from sediment transport to transpiration cooling in gas
turbines. Porous media are often characterized by a complex spatial topology
(Chu  _et al._ , 2018a; Lan  _et al._ , 2021; Chu  _et al._ , 2019), and
numerous properties have been recognized to exert a substantial influence on
both mass and momentum transport at the interface between porous and free-
flowing media. Despite this understanding, the task of comprehending and
optimizing the macroscopic and microscopic characteristics of porous media
continues to present a challenge. The underlying reasons for this complexity
include the inherent difficulties associated with conducting experimental
measurements at the pore scale Terzis  _et al._ (2019) and the substantial
computational resources required for numerical simulations that adequately
resolve the entire range of relevant scales. Historical examinations of
turbulent flow over permeable substrates have investigated the impact of
variable permeability on surface flows. In scenarios where interfacial
permeability is low, the turbulent surface flow exhibits characteristics akin
to a canonical boundary layer, owing to the limited disruption of near-wall
structures. Conversely, when permeability is augmented, large-scale vortical
structures begin to manifest within the surface flow. This phenomenon has been
ascribed to Kelvin-Helmholtz (KH) type instabilities, originating from the
inflection points present in the mean velocity profile. In more recent times,
the exploration of porous media as a means of achieving drag reduction has
catalyzed a series of comprehensive studies focused on the subject of
anisotropic permeability.

Direct Numerical Simulation (DNS) provides a distinct advantage for observing
and analyzing the physics of turbulence within a constrained small spatial
domain, extending beyond traditional applications such as channel (Pandey  _et
al._ , 2020a) and pipe flows Chu  _et al._ (2016); Pandey  _et al._ (2017).
Suga  _et al._ (2020) conducted an investigation into the influence of the
anisotropic permeability tensor within porous media at a regime characterized
by higher permeability. Their findings revealed that both streamwise and
spanwise permeabilities contribute to turbulence enhancement, while vertical
permeability alone does not exert a similar influence. The effect on
turbulence is particularly pronounced in the presence of porous walls
exhibiting streamwise permeability, as this configuration promotes the
development of large-scale streamwise perturbations induced by
Kelvin–Helmholtz instability. However, despite the insights offered by DNS,
its utilization for extensive parametric studies, often required in industrial
applications, remains a costly solution. In such contexts, it is still
customary to employ either analytical models or advanced turbulence models,
which tend to provide satisfactory accuracy within established applications.
Nevertheless, these models are not without their challenges, as they can
exhibit divergent behavior, and the calibration process can prove to be both
intricate and time-consuming.

Machine learning (ML) is on surge where it has successfully tackled many
complex tasks ranging from self-driving cars Jing and Tian (2020), drug
discovery Vamathevan  _et al._ (2019), weather prediction (Chattopadhyay  _et
al._ , 2020), and more recently state-of-the-art large language model such as
GPT for natural language processing (Zhang and Li, 2021). Machine learning has
been democratized heavily in the last decade and some of the reasons include
the vast availability of open-source libraries, an active community, a huge
volume of incoming data from physical as well as numerical experiments, and an
ease of access to high-performance computing resources via cloud. In recent
years, ML has shown potential in the modeling of various fluid flow use-cases
(Pandey  _et al._ , 2022; Heyder  _et al._ , 2022; Nakamura  _et al._ , 2021;
Fang  _et al._ , 2019; Sayyari  _et al._ , 2022; Chu  _et al._ , 2018b;
Larrañaga  _et al._ , 2023; Pandey  _et al._ , 2020b; Pfeffer  _et al._ ,
2022; Heyder and Schumacher, 2021; Chang  _et al._ , 2018). Typically, big
data ML problems rely upon some kind of feature extraction process e.g.
conventional feature selection process, long-standing proper orthogonal
decomposition (POD) Wang  _et al._ (2022); Liu  _et al._ (2023a, b),
relatively new dynamic mode decomposition (DMD) Kutz  _et al._ (2016) or ML-
based methods (Storcheus  _et al._ , 2015). Feature extractors transform the
input raw data into information-preserving feature vectors with the primary
goal of reducing the dimensionality of the data. Several attempts have been
made in the past few years to combine POD or DMD with machine learning (Pandey
and Schumacher, 2020; Heyder and Schumacher, 2021; Ghazijahani  _et al._ ,
2022; Jia  _et al._ , 2022). However, these methods have certain limitations
especially when the flow advances in the turbulent flow regimes. Therefore, an
end-to-end machine learning system is being suggested where POD and DMD are
replaced by a non-linear DNN (Pandey and Schumacher, 2020; Nakamura  _et al._
, 2021; Rahman  _et al._ , 2019). A common architecture is a deep convolution
autoencoder (DAE) to extract the features for further processing and decision
making (Ryu  _et al._ , 2019; Polic  _et al._ , 2019). Turbulent flow problem
often involves the tempo-spatial data, therefore, extracted features could be
used to train a downstream network to build a data-driven dynamic model. A
common choice is recurrent neural network (RNN) which preserves the sequential
relationship in a time-series Srinivasan  _et al._ (2019). A combination of
such approaches has shown exemplary results in terms of data compression and
modeling Nakamura  _et al._ (2021).

In this work, we make an attempt to model complex flow behaviors in a porous
medium with the help of a data-driven approach. The flow is in a turbulent
regime and due to the presence of porosity, geometry is resolved to a finer
scale. This setup makes the end-to-end ML task challenging and equally
interesting. The goal is to combine the deep learning-based feature extraction
step followed by a recurrent neural network to model temporal evolution. Echo
state network was chosen as the RNN due to its proven efficacy in modeling
multi-dimensional tempo-spatial data. The trained models are restricted to the
domain of training and work poorly beyond the trained data domain which could
be a specific parameter such as the Reynolds number in fluid flow problems. A
commonly used strategy called transfer learning is typically used, where the
objective is to transfer the gained knowledge from the training of source task
to related target task with a limited dataset Zhuang  _et al._ (2020). Fine-
tuning, which can be considered a part of transfer learning, is one
methodology where an already trained and optimized model is tuned with another
similar set of data for a similar task. This could decrease the training time
significantly in addition to neural architecture search and training data
requirements Tajbakhsh  _et al._ (2016). Fine-tuning can be performed on the
entire network level where we are allowed to retrain all the layers while
initializing the parameters with the pre-trained model. Alternatively, several
layers can be frozen and retraining is allowed on only a few layers. The
former method is suitable when the network size is very large and data is
limited. While the first approach could lead to overfitting Yosinski  _et al._
(2014). Therefore, we also present the result from fine-tuning where a trained
model from a given data domain was transferred to a different data domain.

The remainder of the paper is divided into 4 sections. Section II describes
the fluid flow problem in porous media along with the underlying governing
equations and numerical method. Section III explains in detail regarding the
hierarchical data-driven model. Section IV describes the results from training
and fine-tuning. With section V, we conclude the work and give a brief
outlook.

##  II Synthetic data generation

High-fidelity simulations such as Direct Numerical Simulation (DNS) (Chu  _et
al._ , 2020; Pandey  _et al._ , 2020a; Chu  _et al._ , 2021a, b) can serve
reliable spatial and temporal resolved data for the modeling. In our DNS, the
three-dimensional incompressible Navier$-$Stokes equations are solved in a
non-dimensional form as follows:

$\frac{\partial u_{j}}{\partial x_{j}}=0$ (1) $\frac{\partial u_{i}}{\partial
t}+\frac{\partial u_{i}u_{j}}{\partial x_{j}}=-\frac{\partial p}{\partial
x_{i}}+\frac{1}{Re_{D}}\frac{\partial^{2}u_{j}}{\partial x_{i}\partial
x_{j}}+\Pi\delta_{i1}$ (2)

Here, $\Pi$ represents a constant pressure gradient in the mean-flow
direction. The governing equations are normalized by employing the half-width
of the entire simulation domain, denoted as $H$ (refer to figure 1a), and the
mean bulk velocity $U_{b}$ within the channel region specified by $y/H=[0,1]$.
The velocity components in the streamwise ($x$), wall-normal ($y$), and
spanwise ($z$) directions are henceforth expressed as $u$, $v$, and $w$,
respectively. The domain dimensions ($L_{x}/H\times L_{y}/H\times L_{z}/H$)
are fixed at $10\times 2\times 0.8\pi$ for all scenarios considered. The lower
half of the domain ($y/H=[-1,0]$) encompasses the porous media, while the
upper half ($y/H=[0,1]$) constitutes the free channel flow. The porous layer
is constructed with 50 cylindrical elements aligned in the streamwise
direction and 5 rows positioned in the wall-normal direction, as delineated in
figure 1. The distance $D$ between two adjacent cylinders is consistently set
at $D/H=0.2$. No-slip boundary conditions are enforced on the cylinders, the
upper wall, and the lower wall, whereas periodic boundary conditions are
implemented in both the streamwise and spanwise directions.

---
Figure 1: A $x-y$ cutting plane of the computational domain with intastanteous
streamwise velocity is depicted in color.

The spectral/$hp$ element solver Nektar++ is utilized to solve the equations
as referred to in Eqn. 1-2 (Cantwell  _et al._ , 2011; Chu  _et al._ , 2019;
Wang  _et al._ , 2021a, b). The geometric configuration in the $x$-$y$ plane
is discretized utilizing quadrilateral elements, with localized refinement in
the vicinity of the cylinders (refer to figure 1(a)). Localized element
expansions are implemented based on the modified Legendre basis (Karniadakis
_et al._ , 2005). Flexibility in polynomial orders is employed across the
wall-normal range through continuous Galerkin projection. Specifically, the
polynomial order within the free-flow region, defined as $y/H=[0.2,1]$, is set
at $P=6$-$7$. The proximal wall region and the upper two strata of cylinders,
where $y/H=[-0.4,0.2]$, are augmented with a superior order of $P=8$-$9$.
Conversely, an inferior order of $P=5$ is designated within the more profound
locations of the cylinder array ($y/H=[-1,-0.4]$). The spanwise orientation is
broadened utilizing a Fourier spectral method, and the 2/3 rule is implemented
to circumvent aliasing errors. The temporal progression is executed with a
second-order mixed implicit-explicit (IMEX) scheme, with the time step fixed
at $\Delta T/(H/U_{b})=5\times 10^{-4}$.

Case | $\varphi$ | $Re_{\tau}^{p}$ | $Re_{\tau}^{s}$ |  $\sqrt{K_{xx}}^{p+}$,$\sqrt{K_{yy}}^{p+}$ | $\sqrt{K_{zz}}^{p+}$ | $C_{xx}^{{p+}},C_{yy}^{{p+}}$ | $r_{c}^{{p+}}$ | $C_{f}^{p}$ | $C_{f}^{s}$
---|---|---|---|---|---|---|---|---|---
C05 | 0.5 | 336 | 180 | 4.55 | 8.86 | 0.37 | 42 | 0.0112 | 0.0084
C06 | 0.6 | 464 | 190 | 9.34 | 15.23 | 0.58 | 48 | 0.0149 | 0.0085
Table 1: Simulation parameters. The porosity of the porous medium region is
$\varphi$. The friction Reynolds numbers are $Re_{\tau}^{p}$ and
$Re_{\tau}^{s}$ for the porous and impermeable top walls respectively.
$\sqrt{K_{\alpha\alpha}}^{p+}$ and $C_{\alpha\alpha}^{{p+}}$ are the diagonal
components of the permeability tensor and Forchheimer coefficient,
respectively, in the direction of $\alpha$ ($\alpha\in\\{x,y,z\\}$), which are
normalized by wall units. $r_{c}^{{p+}}$ is the radius of the cylinders.

Two DNS cases are performed with varying porosity $\varphi=0.5,0.6,0.7$, and
$0.8$, which is defined as the ratio of the void volume to the total volume of
the porous structure. The parameters of the simulated cases are listed in
TABLE 1, where the cases are named after their respective porosity. The
superscripts $(\cdot)^{p}$ and $(\cdot)^{s}$ represent permeable wall and
smooth wall side variables, respectively. Variables with superscript + are
scaled by friction velocities $u_{\tau}$ of their respective side and
viscosity $\nu$.

It should be noted that the distance between the cylinders remains constant,
while the porosity undergoes modification through variation in the cylinder
radii. The normalized cylinder radius is found to be within the range
$r_{c}^{p+}=r_{c}u_{\tau}^{p}/\nu=42$-$52$ for all examined cases (refer to
TABLE 1), thereby implying that the effects of surface roughness are presumed
to be consistent across different scenarios. For all instances, the Reynolds
number associated with the top wall boundary layer is configured at
$Re_{\tau}^{s}=\delta^{s}u_{\tau}^{s}/\nu\approx 180$ ($\delta$ represents the
distance between the position of maximal streamwise velocity and the wall).
This configuration aids in minimizing fluctuations in the top wall boundary
layer. In the region of the upper smooth wall, the streamwise cell size spans
from $4.1\leq\Delta x^{s+}\leq 6.3$, and the spanwise cell size is confined to
below $\Delta z^{s+}=5.4$. On the side corresponding to the porous media, the
value of $\Delta z^{p+}$ is capped at 8.4, whereas $\Delta x^{p+}$ and $\Delta
y^{p+}$ are augmented by polynomial refinement of the local mesh (Cantwell
_et al._ , 2011). The cumulative number of grid points fluctuates between
88$\times 10^{6}$ (C05) and 110$\times 10^{6}$ (C06), with each cylinder
within the porous domain resolved using 80 to 120 grids along its perimeter.
The spatial resolution of the current work aligns closely with that of
preceding DNS studies Shen  _et al._ (2020); Karra  _et al._ (2022).
Additionally, the utilization of a high-order scheme in conjunction with body-
fitted mesh confers a notable advantage in resolving fine scales and
necessitates fewer grids compared to both finite volume and immersed boundary
methods to attain equivalent accuracy (Kooij  _et al._ , 2018). Furthermore, a
comparative analysis of our grid resolutions with the Kolmogorov length scale,
denoted as $\eta=(\nu/\epsilon)^{1/4}$, at the interface has been conducted.
For all cases, the findings show that $(\Delta x/\eta){y=0}\leqslant 2$,
$(\Delta y/\eta){y=0}\leqslant 1$, $(\Delta z/\eta)_{y=0}\leqslant 4.5$,
thereby validating that the current resolution is indeed sufficient

##  III Data-driven modelling

Data driven modelling (DDM) is an approach to harness the available data about
a system and establish a connection between the system state variables (input,
internal and output variables) without explicit knowledge of the physical
behaviour Solomatine  _et al._ (2008). In this work, we have the synthetic
data derived from the DNS experiments while the end model should be based upon
this data and should able to provide a model to predict the temporal evolution
of a state variable. Figure 2 depicts our end-to-end 2-step ML framework. To
have an extendable framework, we chose to first extract the features with the
help of a deep encoder and feed the extracted features to a second-level
recurrent neural network (RNN). This architecture is inspired from our earlier
work Pandey  _et al._ (2022).

Figure 2:  Two-step machine learning architecture for end-to-end ROM model for
porous media. In the first step, DAE is trained as a feature extractor cum
decoder while the second step trains the extracted feature in an
autoregressive manner.

###  III.1 Feature extraction with autoencoder

The deep convolution autoencoder (DAE) operates within the framework of
unsupervised machine learning, meaning that the model’s training does not
depend on labeled data. The learning objective is to distill efficient
representations or features within a low-dimensional hyper-plane, with the aim
of reproducing the original data while minimizing reconstruction error. In
more formal terms, an autoencoder encodes input data $\bm{I}$ into a compact
form $\bm{c}=\textnormal{encode}(\bm{I})$, where $\bm{c}\in\mathbb{R}^{l}$. It
then decodes this compact form back into a reconstructed representation
$\hat{\bm{I}}=\textnormal{decode}(\bm{c})$. The difference between the
reconstructed and original data, given by $\hat{\bm{I}}-\bm{I}$, constitutes
the reconstruction error.

In the specific implementation discussed, 2-dimensional data were sourced from
DNS, as delineated in II. The architecture of the autoencoder is deep,
comprising several convolution layers. Each layer functions as a local feature
extractor and shares its weights across the network. Following each
convolution layer is a pooling layer—specifically, a max pooling layer in this
instance—which diminishes the dimensionality of the data and emphasizes non-
repetitive key features. This contributes to the network’s unique translation
invariance property. The encoder segment of the network concludes with a
bottleneck layer, which harbors the latent space or extracted features. These
features are subsequently channeled into a decoder section, which attempts to
reconstruct the original DNS data. The training process initiates with random
weight initialization across the network and iteratively refines these weights
according to the loss function, employing backpropagation for optimization.

Consider three-dimensional input data $\bm{I}\in\mathbb{R}^{N_{x}\times
N_{y}\times C_{\textrm{in}}}$ with $C_{\textrm{in}}$ because the number of
input channels. Here, $C_{\textrm{in}}=1$ as we load the streamwise velocity
field $u_{x}$ only into the network and thus $\bm{I}\in\mathbb{R}^{N_{x}\times
N_{y}}$. This input is convoluted with a kernel $k_{m}$ as shown in Eqn. 3.

$\textrm{Conv}(m,\bm{I})=\psi\left(b_{m}+\sum_{i=1}^{C_{\textrm{in}}}k_{m}\ast\bm{I}_{i}\right)\,\,\,\textrm{for}\,m\in
M$ (3)

In the encoder, multiple convolutional layers are typically followed by
pooling layers. A window of configurable size and sliding is used with a
configurable step across the input. This process combines the input within
each window to produce a single output element using a chosen aggregation
function. The common max-pooling layer’s aggregation function retains only the
maximum element within each window. The pooling layer’s step size is usually
chosen to reduce the input’s dimensionality, denoted as $\bm{I}$.

Within the encoder segment of a DAE, multiple convolutional layers are
commonly succeeded by pooling layers. These layers operate by deploying a
configurable sliding window across the input, aggregating the information
within each window into a single output element via a specified aggregation
function. A frequent choice for this function is the retention of the maximum
element within each window, as implemented in a max-pooling layer. This
approach serves to diminish the dimensionality of the input, represented by
$\bm{I}$, typically through the strategic selection of the pooling layer’s
step size.

The cumulative effect of the encoding sequence is to map the input data into a
latent space $\mathbb{R}^{l}$, where $l\leq N_{x}\times N_{y}$, and $N_{x}$
and $N_{y}$ denote the dimensions of the initial input. Following the encoder,
the decoder of the DAE primarily mirrors the encoder’s design, with a critical
distinction: it substitutes max-pooling layers with upsampling layers to
expand the data back to its original domain size. This expansion process,
known as unpooling, essentially doubles the data’s dimensions. In the context
of the work discussed, nearest-neighbor linear interpolation is employed to
facilitate the upsampling. The processing of intricate turbulence data
necessitates a sophisticated, deep architecture composed of multiple
convolutional layers, augmented by corresponding upsampling layers. Such an
arrangement ensures the attainment of a faithful representation, minimizing
information loss. In the comprehensive system, subsequent to training all
relevant networks, the DAE’s decoder is applied to extract a compressed form,
denoted by $\hat{\bm{c}}_{t}$, from a specific turbulence snapshot
$\bm{I}_{t}$. This condensed snapshot is then input into an RNN, which
predicts the subsequent reduced representation $\hat{\bm{c}}_{t+1}$. The DAE’s
decoder translates $\hat{\bm{c}}_{t+1}$ back into the input domain, recreating
the anticipated turbulence snapshot. Through the integrated use of deep
convolutional autoencoders and recurrent neural networks, the system attains
efficient and insightful compression of complex turbulence data. This approach
enhances the capacity for predictive modeling and nuanced analysis within the
realm of fluid mechanics.

###  III.2 Echo State Network (ESN)

Within the framework of the present study, ESN serves as a downstream
recurrent neural network, functioning to predict the temporal evolution of the
features extracted by the encoder segment of a DAE in an autoregressive
fashion. The forecasted values are reintroduced to the decoder segment of the
DAE, facilitating the recreation of the target feature. During the operational
phase, the focus is exclusively on the trained ESN and the decoder, while the
encoder is disregarded.

The structural construct of the reservoir is represented as an $N\times N$
adjacency matrix, denoted by $\textbf{W}_{\eta}^{(r)}$. The initialization of
this reservoir matrix is predicated on a specific vector of hyperparameters,
symbolized by $\eta$, and is applied to translate the ESN’s input into a
concealed representation known as the reservoir state, amassing information
from previous inputs. The reservoir dimension, referred to as $N$, is usually
selected to substantially exceed the latent space’s dimension, following the
criterion $N\geq l$. Additionally, this reservoir is dynamically updated with
every time increment, receiving new input. Specifically, at each discrete time
$t$, the input vector $\bm{c}_{t}$ modulates the computation of an updated
reservoir state $\bm{r}_{t}\in\mathbb{R}^{N}$.

The evolution of the reservoir state is defined by:

$\bm{r}t=(1-\alpha)\bm{r}{t-1}+\alpha\textrm{tanh}\left(\textbf{W}^{(\textrm{in})}\bm{c}t+\textbf{W}^{(\textrm{r})}\eta\bm{r}{t-1}\right)$
(4)

where $\alpha$ denotes a leakage rate, dictating the integration of previous
states and current inputs. The random matrices $\textbf{W}^{(\textrm{in})}$
and $\textbf{W}^{(\textrm{r})}\eta$ are initialized at the commencement of
training and are retained as constants. The ESN’s output at the time step $t$
is deduced by:

$\hat{\bm{c}}_{t+1}=\textbf{W}^{(\textrm{out})}\bm{r}_{t}$ (5)

Contrasting with convolutional networks, in the ESN case, only the final
output layer requires training, obviating the need for backpropagation across
several epochs. This implies that the elements of the output matrix
$\textbf{W}^{(\textrm{out})}$ must be meticulously optimized to minimize a
cost function, denoted by $C$, which quantifies the discrepancy between the
training data and the ESN’s output.

The cost function and its minimization are formally represented as:

$C[\textbf{W}^{(\textrm{out})}]=\frac{1}{n_{\textrm{tr}}}\sum_{t=1}^{n_{\textrm{tr}}}(\textbf{W}^{(\textrm{out})}\bm{r}_{t}-\bm{c}_{t})^{2}+\beta\sum_{i=1}^{l}\|W^{(\textrm{out})}_{i}\|_{2}^{2}$
(6)

Here, $W^{(\textrm{out})}_{i}$ represents the $i$-th row of
$\textbf{W}^{(\textrm{out})}$, and $||_{2}^{2}$ signifies the $L_{2}$ norm.
The regression problem is resolved by:

$\textbf{W}^{(\textrm{out})}\ast=\textbf{Y}\textbf{S}^{\textrm{T}}(\textbf{S}\textbf{S}^{\textrm{T}}+\beta\textbf{I}\textbf{d})^{-1}.$
(7)

The ESN’s hyperparameters include the reservoir size $N$, the node density
$D$, the spectral radius of the reservoir
$\rho(\textbf{W}^{(\textrm{r})}\eta)$, the leakage rate $\alpha$, and the
Tikhonov regularization parameter $\beta$ in the cost function $C$.
Collectively, these are captured in the vector
$\bm{\eta}=(N,D,\rho,\alpha,\beta)$. This training methodology is inherently
expeditious relative to other forms of RNNs. Following successful training and
optimization of hyperparameters, the ESN operates as an autonomous dynamical
system. In this mode, a forecasted compressed snapshot $\bm{c}_{t+1}$ can be
iteratively used as the succeeding ESN input to predict $\bm{c}_{t+2}$ and
continue thereafter. This mode of operation is recognized as the closed-loop
scenario.

##  IV Results and Discussion

###  IV.1 Data driver flow model

As mentioned earlier, we first trained a DAE network to extract the feature in
lower dimensions. For the same, input and output layers have the same dataset
i.e. the snapshots from the DNS consist of $u$-velocity component with a
dimension of $1501\times 109\times 1$. The neural network has multiple
convolution layers combined with a max pooling layer while the decoder section
has a convolution layer and up-sampling layers, giving a converging-diverging
shape. The features or latent modes are obtained at the end of the encoder
section and it has a shape of $24\times 2\times 34$, thereby giving a
dimensionality reduction of 99%. This could increase further on the expense of
reconstruction loss. TABLE 2 illustrates the various layers and their shape.
After fixing the DAE architecture, various hyperparameters were obtained by
using Bayesian optimization. This particular approach is part of AutoML
paradigm He  _et al._ (2021), and it assists in achieving an optimized set of
parameters in a relatively smaller number of iterations when compared to the
grid or random search Turner  _et al._ (2021). It is worth mentioning that
data was scaled between 0-1 before the training. For the entire process, we
used 2000 snapshots, 60% are used for training, 20% for validation while
remaining 20% for blind testing. The model has around 1.9 Million trainable
parameters and the model reaches a loss of $O(10^{-4})$ for both training and
validation set in 50 epochs.

Encoder |  | Decoder
---|---|---
Layer | Output size |  | Layer | Output size
Encoder input | $1501\times 109\times 1$ |  | Decoder input | $24\times 2\times 34$
2D Conv-E1 | $1501\times 109\times 256$ |  | 2D Conv-D1 | $24\times 2\times 34$
Max Pool-E1 | $751\times 55\times 256$ |  | 2D Upsamp-D1 | $48\times 4\times 34$
2D Conv-E2 | $751\times 55\times 256$ |  | 2D Conv-D2 | $48\times 4\times 32$
Max Pool-E2 | $376\times 28\times 256$ |  | 2D Upsamp-D2 | $96\times 8\times 32$
2D Conv-E3 | $376\times 28\times 128$ |  | 2D Conv-D3 | $96\times 8\times 64$
Max Pool-E3 | $188\times 14\times 128$ |  | 2D Upsamp-D3 | $192\times 16\times 64$
2D Conv-E4 | $188\times 14\times 64$ |  | 2D Conv-D4 | $192\times 16\times 128$
Max Pool-E4 | $94\times 7\times 64$ |  | 2D Upsamp-D4 | $384\times 32\times 128$
2D Conv-E5 | $94\times 7\times 32$ |  | 2D Conv-D5 | $384\times 32\times 256$
Max Pool-E5 | $47\times 4\times 32$ |  | 2D Upsamp-D5 | $768\times 64\times 256$
2D Conv-E6 | $47\times 4\times 34$ |  | 2D Conv-D6 | $768\times 64\times 256$
Max Pool-E6 | $24\times 2\times 34$ |  | 2D Upsamp-D6 | $1536\times 128\times 256$
|  |  | Output with 2D Conv | $1536\times 128\times 1$
|  |  | Output with Cropping | $1501\times 109\times 1$
Table 2: DAE architecture used at first level for feature extraction. The
table summarizes the encoder and decoder architectures (Conv = convolution,
Max Pool = max pooling, Upsamp = upsampling). Symbol E2 denotes for example
encoder hidden layer No. 2.

Once DAE is trained, we extract the features for the entire dataset using its
encoder section. These encoded data which have 99% reduced dimension were
further used to train the second level of RNN i.e. ESN. Unlike DAE, training
an ESN is relatively easy and computationally cheaper because of
straightforward design choices. The optimized network of ESN has 3726
reservoirs, 0.93 as spectral radius, 0.94 as leaking rate and 0.16 as
reservoir density. For the best results, the look back history period is 4
timesteps i.e., the model requires 4 timesteps data as an initialization then
the system can predict auto-regressively. As a first step of validation, we
visualize the low-dimensional features from DAE along with the generated data
from ESN on the blind validation set. We call the method generative because
the model works auto-repressively once we provide the initial data. Figure 3
shows such a time series for 2 of the modes out of 1632. It is impressive that
the model maintains a value of 0 without any minor fluctuations for $\phi_{1}$
while the series fluctuates with the right phase and amplitude for
$\phi_{10}$. Figure 4 further illustrates the probability density function
(PDF) of these 2 modes for the DAE and ESN. It reaffirms the generative
capability of ESN where it is able to match the overall distribution with GT
(i.e., DAE for ESN).

Figure 3:  A comparison of extracted features from DAE which acts as ground
truth for the training, and prediction from ESN for 2 of the modes. Figure 4:
Comparison between DAE and ESN modes in terms of PDF for 2 of the modes

As a next step of validation, a comparison in high-dimensional space where
flow physics needs to be captured is performed. Figure 5 displays contour
plots of the original DNS data in panels (a) and (d), reconstructed data from
the features space of DAE in panels (b) and (e), and the decoded value from
the prediction of ESN is shown in panels (c) and (f). One can observe a good
qualitative agreement in terms of capturing the flow fields. Both regions i.e.
main channel flow and porous flow can be predicted. Both DAE and ESN were
trained from the DNS, and are not able to fully capture the strong turbulent
field in the channel. While the flow in the porous media section is much
better because it falls within the laminar region mostly. The obvious reason
is the limited capability of DAE to capture the entire flow field because it
discards 99% of input data and uses only 1% of the features, similar to
conventional reduced-order modeling approaches such as POD. The derived model
with ESN shows an excellent result while matching the field compared to its GT
i.e. DAE. It suggests that the entire model can further be improved if a
state-of-the-art feature extractor such as a transformer model can be used
instead of DAE. It is worth mentioning that DAE surpasses conventional as well
as some advanced feature extractors e.g., variational autoencoder for this
use-case, which was observed in our experimentation. Often, averages and
fluctuations are critical values for the analysis in industrial applications.
Therefore, we further validate the results with mean and fluctuation profiles
over the cross-section of the channel. Figure 6 shows these profiles and the
mean velocity profile overlaps with the DNS for both DAE and ESN. The
fluctuation of velocity is slightly deviated in the main channel while they
remain zero in the porous section due to dominating laminar flow. The
distinction made by the model is also great where it can predict the interface
where flow is transiting from laminar to turbulent.

Figure 5:  Contours of $u_{x}$ by (a,d) DNS (GT), (b,e) DAE (feature
extractor) and (c,f) ESN (final end-to-end model) at 2 different timestamps.
Figure 6:  Comparison among ground truth (DNS), reduced order model (DAE) and
temporal model for DAE modes using ESN. (a): Average fields (b) Fluctuations

For a comprehensive analysis, we include findings from Fourier spectra in our
comparison. Figure 7 contrasts the streamwise spectra $E_{uu}$ near both the
smooth wall at $y=19.5$ and the permeable wall at $y=10.2$. Both the DAE and
ESN models perform well in capturing the energy at larger scales ($k_{x}<1$).
However, neither DAE nor ESN successfully models the abundant turbulent
kinetic energy found in the inertial subrange, particularly at $y=19.5$. This
limitation is anticipated, given that the DAE model achieves a 99%
dimensionality reduction.

Figure 7: Streamwise spectra of velocity fluctuation at different locations,
(a) close to the smooth wall, (b) close to the permeable wall

###  IV.2 Fine tuning to a various porosity

The trained model has shown a good agreement in terms of quantitative as well
as qualitative characteristics of the flow. However, the trained model is
limited to a given data domain that was used for the training which means that
the model will not predict accurate results if we use different parametric
conditions. This is one of the biggest downsides of such a data-driven model.
Considering the current setup, different porous media with various porosities
(Wang  _et al._ , 2021b, 2022) and topologies (Wang  _et al._ , 2021a) are
often compared together, whereas the channel remains unchanged. Variations in
porous topology, porosity, and permeability can lead to distinct influences on
turbulence modulation, which in turn could impact functionalities like drag
reduction, noise absorption, and heat transfer enhancement. We evaluate the
use of fine-tuning as a method to extend the versatility of these models to
broader data sets. Fine-tuning involves adapting a pre-trained model to new,
limited data. Given that the DAE is computationally intensive within this
workflow, we are seeking a substantial reduction of the training time by using
fine-tuning. The DAE model was tested using a distinct data set for this
assessment. Here, we used case C06 which has higher level of porosity
$\varphi=0.6$. The other parameters such as turbulent Reynolds number on both
wall sides $Re_{\tau}^{p},Re_{\tau}^{s}$, permeability
$\sqrt{K_{\alpha\alpha}}^{p+}$ and Forchheimer coefficients
$C_{\alpha\alpha}^{{p+}}$ are also quite different between C05 and C06 , as
shown in TABLE 1.

Case Name | Description | Training time | Mean Square Error
---|---|---|---
DNS | Simulation results which acts as GT for case C06 | - | -
case A | Direct prediction using model from case C05 without any training | 0 | 0.081
case B | Model from case C05 with only 1 trainable layer for case C06 |  $\mathcal{O}$(0.5 $\times$ 101) | 0.064
case C | Model from case C05 retrained with 10% data points for case C06 |  $\mathcal{O}$(101) | 0.049
case D | Model from case C05 retrained with all data points for case C06 |  $\mathcal{O}$(102) | 0.030
case E | Model trained from scratch similar to case C05 |  $\mathcal{O}$(103) | 0.025
Table 3: Various cases considered for fine-tuning.

TABLE 3 shows the cases used for the fine-tuning study. It compares the
training time and the mean square error of these cases. The case numbering is
based on an increasing order of the training time. Case A uses the trained
model with C05 without any change. Therefore, its training time is 0 while the
error is also the highest. Case E is trained from scratch, which is the same
as the procedure in the last section. It is featured with the longest training
time and the lowest mean square error showing the high accuracy and high
training cost. Cases A and E represent two extreme conditions, whereas cases
B, C and D fall between the two extremes. Case D is based on the same model as
the original model from C05, however, it is retrained with all data points for
case C06. Its training cost $\mathcal{O}$(102) is one order of magnitude lower
than case A whereas the mean square error is only slightly higher which
indicates a worthy trade-off. Case C is fundamentally close to case D however
only retrained with 10$\%$ data for case C06. Hence, its training time
$\mathcal{O}$(101) is again one order of magnitude lower than case D but with
a significant increase of the mean square error. Case B is inherited with the
model from case C05 with only 1 trainable layer for case C06. The training
time is approximately half of that of case C and shows a significant increase
in the error. In general, both case C and case D are valuable approaches with
a reasonable trade-off of accuracy and training cost. Case D offers a high
accuracy and still a significant reduction of training cost, whereas case C
can be valuable when the training resource is extremely limited.

Figure 8:  Contours of (a, b, c, d, e, f): Instantaneous velocity, (g, h, i,
j,k,l): temporal average velocity. Labels are as per Tab. 3

Figure 8 compares the instantaneous velocity and time-averaged velocity
obtained from the ground truth (DNS) with all the fine-tuning cases from TABLE
3. In the DNS result (figure 8(a)), turbulent structures on multi-scales can
be resolved by fine resolutions. On the other hand, low-rank modeled data by
DAE can capture the dominant structures such as the blow and suction events,
for instance, the low-velocity region at the interface ($x\approx 20,y\approx
10$). The fine details are missing due to the truncation of the original data.
Compared to the model case A without additional training (figure 8(c)), cases
B,C,D,E certainly show more details, which corresponds to the accuracy. From
cases A to E, an increase in fine motions is observed. On the average results,
all the mean velocity contours show similar behavior with the maxima in the
channel center and the minima down to the porous domain. It is somehow
difficult to tell the quality of modeling.

Moreover, figure 9 shows the quantification of the time-averaged streamwise
velocity profile and its fluctuation. In the channel region, all the models
are able to reproduce the DNS results even in case A without any training on
the DNS data. It is because the channel domain is topologically the same even
though the porosity in the connected domain is different. A clear difference
is seen in the porous domain as well as in the vicinity of the interface. Case
A in red color is not able to follow the DNS data of the mean velocity as well
as the velocity fluctuation. On the other hand, cases B, C and case D,
empowered by the fine-tuning method, show an excellent modeling of the flow
inside the porous domain regarding the mean velocity profile. This is quite
impressive considering case B has only one trainable layer and with more than
two orders of magnitude lower training cost. This suggests that a small amount
of data is able to fine-tune the model to model and predict the mean profile
in the entire domain. The prediction of velocity fluctuation in figure 9(b)
seems much more challenging. Case B shows only a slight improvement but is
still close to that of case A. Even though, case B can model the mean velocity
profile with no clear discrepancy. Case C (10% new training data and
$\mathcal{O}$(102) reduction of the training time) shows a significant
improvement in both the channel and porous region. Its performance is quite
close to the case D and E with much higher training cost. Therefore, case C is
considered as a balanced case with efficiency and accuracy. On the other hand,
a further cost reduction of from case C to case B is not quite significant
where the accuracy loss is noticeable.

Figure 9:  Comparison among ground truth (DNS), reduced order model (DAE)
models. (a): Average fields (b) Fluctuations

For a more in-depth examination of the fine-tuning results in turbulent
channel flows, Figure 10 shows the streamwise spectra $E_{uu}$ for all fine-
tuned cases against DNS results at $y=19.5$ (close to the top wall) and
$y=10.2$ (close to the porous wall). Across all fine-tuned models, the
predictions for large-scale motion at $k_{x}<1$ are highly accurate. However,
deviations become apparent at higher streamwise wavenumbers. Although
differences in time-averaged velocity fluctuations are noticeable (as seen in
figure 9), they are less discernible in frequency space when plotted on a log-
log scale. Among the cases, cases D and E stand out for their superior
performance in the inertial subrange at both positions.

---
Figure 10: Streamwise spectra $E_{uu}$ close to the smooth wall (a) and
permeable wall (b)

To further investigate the performance of fine-tuned models, figure 11
displays the instantaneous velocity field $u_{x}$ to the porous media domain.
Case A fails to accurately capture the flow pattern within the porous domain
due to the missing input of the new training data. Starting from case B, fine-
tuning is able to significantly increase the modeling of the flow inside
porous media even though the training data is quite limited. Increasing the
training effort from case B to case E can progressively enhance the fidelity.
Even though the flow in porous media is significantly influenced by dispersion
caused by the porous structure, there are still instances of non-uniform
velocity around $y\approx 8$. This is indicative of a transitional regime
between laminar and turbulent flow, due to its proximity to the interface.
This phenomenon is better captured by Cases D and E.

---
Figure 11: A detailed look of instantaneous velocity $u_{x}$ to the porous
media domain

The interfacial behavior is of great interest since it enables bi-directional
interactions between channel flow and porous media (Wang  _et al._ , 2021b,
2022). On the other hand, it is featured with high complexity involving two
different types of geometries as well as a mixing state of turbulent,
transition and laminar flows. Figure 12 depicts the instantaneous velocity
fluctuation at the interface from DNS and various reduced-order models. In the
boundary layer region ($x>10$), DNS shows a rich scale from the large-scale
structures to the fine motions, whereas ROM can only reproduce the bigger
ones. This is reasonable considering that ROM has to ignore these details to
deliver a low-rank model. Among them, case A shows a filtered modeling of the
velocity fluctuation field, whereas more details are retained by other cases.
Below the interface, the strength of fluctuation is much weaker. But case C, D
and E can still deliver the model in a certain way. Case A has received no
training data from this DNS, therefore, the porous topology is unclear and the
prediction is far from satisfactory.

---
Figure 12: A detailed look of instantaneous velocity fluctuation
$u_{x}^{\prime}$ to the channel flow-porous media interface

##  V Conclusions and outlook

Flows characterized by turbulence over permeable interfaces are ubiquitous and
present intricate flow patterns. Addressing the challenges posed by the
complex geometry and mixed state of laminar, transitional and turbulent
regimes, our study introduces a data-driven, end-to-end machine learning
framework specifically designed to model turbulent flows within a channel flow
coupled with porous media. To achieve this, we’ve established a non-linear
reduced-order model powered by a deep convolution autoencoder network.
Impressively, this model can reduce the spatial dimensions required for highly
resolved simulations—integral for direct numerical simulation—by 99%. To
capture the time-varying characteristics of the reduced modes, we employed a
downstream recurrent neural network. This strategic integration allows the
model to accurately predict the future trajectories of these modes.
Furthermore, to test the adaptability and robustness of our model, we
evaluated its performance on a dataset that presents a different porosity from
the original training set. Through fine-tuning, our system has demonstrated
the potential to significantly cut down on computational resources, reducing
training efforts by up to two orders of magnitude when working with a limited
dataset (constituting merely 10% of the full training data). Even with this
reduced data input, the results achieved were commendable, particularly with
respect to the mean velocity profile. The efficient fine-tuning of our current
model significantly reduces computational time, making it much easier to study
the effect of different types of porous media for drag reduction, enhanced
heat transfer and noise absorption. This efficiency is particularly
advantageous for systematic investigations into interface engineering,
allowing for more focused and quicker research cycles. Of particular note is
the observation that the fine-tuned model excels within the porous domain in
comparison to the channel and interface regions. This suggests that the
inherent topological features of the porous media are more amenable to
training compared to the challenging multi-scale attributes typical of
turbulent flows.

###### Acknowledgements.

The work is supported by the funding by Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) project SFB1313 (project No. 327154368) and under
Germany’s Excellence Strategy-EXC2075-390740016. In addition, all authors
gratefully acknowledge the access to the high performance computing facility
Hawk at HLRS, Stuttgart.

## References

  * Chu  _et al._ (2018a) X. Chu, B. Weigand, and V. Vaikuntanathan, “Flow turbulence topology in regular porous media: From macroscopic to microscopic scale with direct numerical simulation,” Physics of Fluids 30, 065102 (2018a).
  * Lan  _et al._ (2021) Yu-Hsiang Lan, Paul Fischer, Elia Merzari, and Misun Min, “All-hex meshing strategies for densely packed spheres,” arXiv preprint arXiv:2106.00196 (2021).
  * Chu  _et al._ (2019) X. Chu, G. Yang, S. Pandey, and B. Weigand, “Direct numerical simulation of convective heat transfer in porous media,” International Journal of Heat and Mass Transfer 133, 11–20 (2019).
  * Terzis  _et al._ (2019) A. Terzis, I. Zarikos, K. Weishaupt, G. Yang, X. Chu, R. Helmig, and B. Weigand, “Microscopic velocity field measurements inside a regular porous medium adjacent to a low reynolds number channel flow,” Physics of Fluids 31, 042001 (2019).
  * Pandey  _et al._ (2020a) S. Pandey, X. Chu, B. Weigand, E. Laurien, and J. Schumacher, “Relaminarized and recovered turbulence under nonuniform body forces,” Physical Review Fluids 5, 104604 (2020a).
  * Chu  _et al._ (2016) Xu Chu, Eckart Laurien, and Sandeep Pandey, “Direct numerical simulation of heated pipe flow with strong property variation,” in  _High Performance Computing in Science and Engineering´ 16_ (Springer, 2016) pp. 473–486.
  * Pandey  _et al._ (2017) Sandeep Pandey, Eckart Laurien, and Xu Chu, “A modified convective heat transfer model for heated pipe flow of supercritical carbon dioxide,” International Journal of Thermal Sciences 117, 227–238 (2017).
  * Suga  _et al._ (2020) K. Suga, Y. Okazaki, and Y. Kuwata, “Characteristics of turbulent square duct flows over porous media,” Journal of Fluid Mechanics 884 (2020).
  * Jing and Tian (2020) Longlong Jing and Yingli Tian, “Self-supervised visual feature learning with deep neural networks: A survey,” IEEE Trans. Pattern Anal. Mach. Intell (2020).
  * Vamathevan  _et al._ (2019) Jessica Vamathevan, Dominic Clark, Paul Czodrowski, Ian Dunham, Edgardo Ferran, George Lee, Bin Li, Anant Madabhushi, Parantu Shah, Michaela Spitzer, _et al._ , “Applications of machine learning in drug discovery and development,” Nat. Rev. Drug Discov. 18, 463–477 (2019).
  * Chattopadhyay  _et al._ (2020) Ashesh Chattopadhyay, Pedram Hassanzadeh, and Saba Pasha, “Predicting clustered weather patterns: A test case for applications of convolutional neural networks to spatio-temporal climate data,” Sci. Rep. 10, 1–13 (2020).
  * Zhang and Li (2021) Min Zhang and Juntao Li, “A commentary of gpt-3 in mit technology review 2021,” Fundamental Research 1, 831–833 (2021).
  * Pandey  _et al._ (2022) Sandeep Pandey, Philipp Teutsch, Patrick Mäder, and Jörg Schumacher, “Direct data-driven forecast of local turbulent heat flux in rayleigh–bénard convection,” Physics of Fluids 34, 045106 (2022).
  * Heyder  _et al._ (2022) Florian Heyder, Juan Pedro Mellado, and Jörg Schumacher, “Generalizability of reservoir computing for flux-driven two-dimensional convection,” Physical Review E 106, 055303 (2022).
  * Nakamura  _et al._ (2021) Taichi Nakamura, Kai Fukami, Kazuto Hasegawa, Yusuke Nabae, and Koji Fukagata, “Convolutional neural network and long short-term memory based reduced order surrogate for minimal turbulent channel flow,” Physics of Fluids 33, 025116 (2021).
  * Fang  _et al._ (2019) R. Fang, D. Sondak, P. Protopapas, and S. Succi, “Neural network models for the anisotropic Reynolds stress tensor in turbulent channel flow,” J. Turb. , 1–19 (2019).
  * Sayyari  _et al._ (2022) Mohammad Javad Sayyari, Jinyul Hwang, and Kyung Chun Kim, “Unsupervised deep learning of spatial organizations of coherent structures in a turbulent channel flow,” Physics of Fluids 34, 115138 (2022).
  * Chu  _et al._ (2018b) Xu Chu, Wanli Chang, Sandeep Pandey, Jiayu Luo, Bernhard Weigand, and Eckart Laurien, “A computationally light data-driven approach for heat transfer and hydraulic characteristics modeling of supercritical fluids: From dns to dnn,” International Journal of Heat and Mass Transfer 123, 629–636 (2018b).
  * Larrañaga  _et al._ (2023) A Larrañaga, SL Brunton, J Martínez, S Chapela, and J Porteiro, “Data-driven prediction of the performance of enhanced surfaces from an extensive cfd-generated parametric search space,” Machine Learning: Science and Technology 4, 025012 (2023).
  * Pandey  _et al._ (2020b) Sandeep Pandey, Jörg Schumacher, and Katepalli R Sreenivasan, “A perspective on machine learning in turbulent flows,” Journal of Turbulence 21, 567–584 (2020b).
  * Pfeffer  _et al._ (2022) Philipp Pfeffer, Florian Heyder, and Jörg Schumacher, “Hybrid quantum-classical reservoir computing of thermal convection flow,” Physical Review Research 4, 033176 (2022).
  * Heyder and Schumacher (2021) Florian Heyder and Jörg Schumacher, “Echo state network for two-dimensional turbulent moist rayleigh-bénard convection,” Physical Review E 103, 053107 (2021).
  * Chang  _et al._ (2018) W. Chang, X. Chu, A. Fareed, S. Pandey, J. Luo, B. Weigand, and E. Laurien, “Heat transfer prediction of supercritical water with artificial neural networks,” Applied Thermal Engineering 131, 815–824 (2018).
  * Wang  _et al._ (2022) Wenkang Wang, Adrián Lozano-Durán, Rainer Helmig, and Xu Chu, “Spatial and spectral characteristics of information flux between turbulent boundary layers and porous media,” Journal of Fluid Mechanics 949, A16 (2022).
  * Liu  _et al._ (2023a) Yanchao Liu, Wenkang Wang, Guang Yang, Hassan Nemati, and Xu Chu, “The interfacial modes and modal causality in a dispersed bubbly turbulent flow,” Physics of Fluids 35 (2023a).
  * Liu  _et al._ (2023b) Y. Liu, X. Chu, W. Wang, and B. Weigand, “Large-eddy simulation, convective instability, and modal causality of coaxial supersonic air–water jets considering a swirl effect,” Physics of Fluids 35 (2023b).
  * Kutz  _et al._ (2016) J Nathan Kutz, Steven L Brunton, Bingni W Brunton, and Joshua L Proctor,  _Dynamic mode decomposition: data-driven modeling of complex systems_ (SIAM, 2016).
  * Storcheus  _et al._ (2015) Dmitry Storcheus, Afshin Rostamizadeh, and Sanjiv Kumar, “A survey of modern questions and challenges in feature extraction,” in  _Feature Extraction: Modern Questions and Challenges_ (PMLR, 2015) pp. 1–18.
  * Pandey and Schumacher (2020) Sandeep Pandey and Jörg Schumacher, “Reservoir computing model of two-dimensional turbulent convection,” Physical Review Fluids 5, 113506 (2020).
  * Ghazijahani  _et al._ (2022) Mohammad Sharifi Ghazijahani, Florian Heyder, Jörg Schumacher, and Christian Cierpka, “On the benefits and limitations of echo state networks for turbulent flow prediction,” Measurement Science and Technology 34, 014002 (2022).
  * Jia  _et al._ (2022) Xuyi Jia, Chunna Li, Wen Ji, and Chunlin Gong, “A hybrid reduced-order model combing deep learning for unsteady flow,” Physics of Fluids 34, 097112 (2022).
  * Rahman  _et al._ (2019) Sk M Rahman, Suraj Pawar, Omer San, Adil Rasheed, and Traian Iliescu, “Nonintrusive reduced order modeling framework for quasigeostrophic turbulence,” Physical Review E 100, 053306 (2019).
  * Ryu  _et al._ (2019) Seunghyoung Ryu, Hyungeun Choi, Hyoseop Lee, and Hongseok Kim, “Convolutional autoencoder based feature extraction and clustering for customer load analysis,” IEEE Transactions on Power Systems 35, 1048–1060 (2019).
  * Polic  _et al._ (2019) Marsela Polic, Ivona Krajacic, Nathan Lepora, and Matko Orsag, “Convolutional autoencoder for feature extraction in tactile sensing,” IEEE Robotics and Automation Letters 4, 3671–3678 (2019).
  * Srinivasan  _et al._ (2019) P. A. Srinivasan, L. Guastoni, H. Azizpour, P. Schlatter, and R. Vinuesa, “Predictions of turbulent shear flows using deep neural networks,” Phys. Rev. Fluids 4, 054603 (2019).
  * Zhuang  _et al._ (2020) Fuzhen Zhuang, Zhiyuan Qi, Keyu Duan, Dongbo Xi, Yongchun Zhu, Hengshu Zhu, Hui Xiong, and Qing He, “A comprehensive survey on transfer learning,” Proceedings of the IEEE 109, 43–76 (2020).
  * Tajbakhsh  _et al._ (2016) Nima Tajbakhsh, Jae Y Shin, Suryakanth R Gurudu, R Todd Hurst, Christopher B Kendall, Michael B Gotway, and Jianming Liang, “Convolutional neural networks for medical image analysis: Full training or fine tuning?” IEEE transactions on medical imaging 35, 1299–1312 (2016).
  * Yosinski  _et al._ (2014) Jason Yosinski, Jeff Clune, Yoshua Bengio, and Hod Lipson, “How transferable are features in deep neural networks?” Advances in neural information processing systems 27 (2014).
  * Chu  _et al._ (2020) X. Chu, Y. Wu, U. Rist, and B. Weigand, “Instability and transition in an elementary porous medium,” Physical Review Fluids 5, 044304 (2020).
  * Chu  _et al._ (2021a) X. Chu, W. Wang, G. Yang, A. Terzis, R. Helmig, and B. Weigand, “Transport of turbulence across permeable interface in a turbulent channel flow: interface-resolved direct numerical simulation,” Transport in Porous Media 136, 165–189 (2021a).
  * Chu  _et al._ (2021b) Xu Chu, Johannes Müller, and Bernhard Weigand, “Interface-resolved direct numerical simulation of turbulent flow over porous media,” in  _High Performance Computing in Science and Engineering ’19_ (Springer International Publishing, 2021) pp. 343–354.
  * Cantwell  _et al._ (2011) C. D. Cantwell, S. J. Sherwin, R. M. Kirby, and P. HJ. Kelly, “From h to p efficiently: Strategy selection for operator evaluation on hexahedral and tetrahedral elements,” Computers & Fluids 43, 23–28 (2011).
  * Wang  _et al._ (2021a) Wenkang Wang, Guang Yang, Cenk Evrim, Alexandros Terzis, Rainer Helmig, and Xu Chu, “An assessment of turbulence transportation near regular and random permeable interfaces,” Physics of Fluids 33, 115103 (2021a).
  * Wang  _et al._ (2021b) W. Wang, X. Chu, A. Lozano-Durán, R. Helmig, and B. Weigand, “Information transfer between turbulent boundary layer and porous media,” Journal of Fluid Mechanics 920, A21 (2021b).
  * Karniadakis  _et al._ (2005) George Em Karniadakis, George Karniadakis, and Spencer Sherwin,  _Spectral/hp element methods for computational fluid dynamics_ (Oxford University Press on Demand, 2005).
  * Shen  _et al._ (2020) G. Shen, J. Yuan, and M. S. Phanikumar, “Direct numerical simulations of turbulence and hyporheic mixing near sediment water interfaces,” Journal of Fluid Mechanics 892, A20 (2020).
  * Karra  _et al._ (2022) Shashank K Karra, Sourabh V Apte, Xiaoliang He, and Timothy D Scheibe, “Pore-resolved simulations of turbulent boundary layer flow over permeable and impermeable sediment beds,” arXiv preprint arXiv:2204.13875 (2022).
  * Kooij  _et al._ (2018) Gijs L Kooij, Mikhail A Botchev, Edo MA Frederix, Bernard J Geurts, Susanne Horn, Detlef Lohse, Erwin P van der Poel, Olga Shishkina, Richard JAM Stevens, and Roberto Verzicco, “Comparison of computational codes for direct numerical simulations of turbulent rayleigh–bénard convection,” Computers & Fluids 166, 1–8 (2018).
  * Solomatine  _et al._ (2008) Dimitri Solomatine, Linda M See, and RJ Abrahart, “Data-driven modelling: concepts, approaches and experiences,” Practical hydroinformatics: Computational intelligence and technological developments in water applications , 17–30 (2008).
  * He  _et al._ (2021) Xin He, Kaiyong Zhao, and Xiaowen Chu, “Automl: A survey of the state-of-the-art,” Knowledge-Based Systems 212, 106622 (2021).
  * Turner  _et al._ (2021) Ryan Turner, David Eriksson, Michael McCourt, Juha Kiili, Eero Laaksonen, Zhen Xu, and Isabelle Guyon, “Bayesian optimization is superior to random search for machine learning hyperparameter tuning: Analysis of the black-box optimization challenge 2020,” in  _NeurIPS 2020 Competition and Demonstration Track_ (PMLR, 2021) pp. 3–26.


* Garraway (2011) B. M. Garraway, “The Dicke model in quantum optics: Dicke model revisited,” Phil. Trans. R. Soc. A 369, 1137–1155 (2011).
  * Ruggenthaler  _et al._ (2014) M. Ruggenthaler, J. Flick, C. Pellegrini, H. Appel, I. V. Tokatly, and A. Rubio, “Quantum-electrodynamical density-functional theory: Bridging quantum optics and electronic-structure theory,” Phys. Rev. A 90, 012508 (2014).
  * Ruggenthaler (2015) M. Ruggenthaler, “Ground-State Quantum-Electrodynamical Density-Functional Theory,” arXiv:1509.01417 [quant-ph] (2015).
  * Tokatly (2013) I. V. Tokatly, “Time-dependent density functional theory for many-electron systems interacting with cavity photons,” Phys. Rev. Lett. 110, 233001 (2013).
  * Schäfer  _et al._ (2021) C. Schäfer, F. Buchholz, M. Penz, M. Ruggenthaler, and A. Rubio, “Making ab initio QED functional(s): Non-perturbative and photon-free effective frameworks for strong light-matter coupling,” arXiv:2106.07507 [quant-ph] (2021).
  * Aschroft and Mermin (1976) N. W. Aschroft and N. Mermin,  _Solid State Physics_ (Harcourt College Publishers, 1976).
  * Callaway (1964) J. Callaway,  _Energy Band Theory_ (Academic Press Inc., 1964).
  * Kohmoto (1985) M. Kohmoto, “Topological invariant and the quantization of the Hall conductance,” Annals of Physics 160, 343 – 354 (1985).
  * Cahill and Glauber (1969) K. E. Cahill and R. J. Glauber, “Ordered expansions in boson amplitude operators,” Phys. Rev. 177, 1857–1881 (1969).
  * Harper (1955) P. G. Harper, “Single band motion of conduction electrons in a uniform magnetic field,” Proc. Phys. Soc. A 68, 874–878 (1955).
  * Rokaj  _et al._ (2022) Vasil Rokaj, Michael Ruggenthaler, Florian G. Eich, and Angel Rubio, “Free electron gas in cavity quantum electrodynamics,” Phys. Rev. Research 4, 013012 (2022).
  * Eckhardt  _et al._ (2021) C. J. Eckhardt, G. Passetti, M. Othman, C. Karrasch, F. Cavaliere, M. A. Sentef, and D. M. Kennes, “Quantum floquet engineering with an exactly solvable tight-binding chain in a cavity,” arXiv:2107.12236 [cond-mat.str-el] (2021).
  * Claro and Wannier (1979) F. H. Claro and G. H. Wannier, “Magnetic subband structure of electrons in hexagonal lattices,” Phys. Rev. B 19, 6068–6074 (1979).
  * Du  _et al._ (2018) Liang Du, Qi Chen, Aaron D. Barr, Ariel R. Barr, and Gregory A. Fiete, “Floquet hofstadter butterfly on the kagome and triangular lattices,” Phys. Rev. B 98, 245145 (2018).
  * Zhao  _et al._ (2020) Ming Zhao, Qi Chen, Xue-Dong Tian, and Liang Du, “Floquet engineering the hofstadter butterfly in the square lattice and its effective hamiltonian,” arXiv:2007.01071 [cond-mat.mes-hall] (2020).
  * Kooi  _et al._ (2018) S. H. Kooi, A. Quelle, W. Beugeling, and C. Morais Smith, “Genesis of the Floquet Hofstadter butterfly,” Phys. Rev. B 98, 115124 (2018).
  * Sentef  _et al._ (2020) Michael A. Sentef, Jiajun Li, Fabian Künzel, and Martin Eckstein, “Quantum to classical crossover of floquet engineering in correlated quantum systems,” Phys. Rev. Research 2, 033033 (2020).
  * Schäfer  _et al._ (2018) C. Schäfer, M. Ruggenthaler, and A. Rubio, “Ab initio nonrelativistic quantum electrodynamics: Bridging quantum chemistry and quantum optics from weak to strong coupling,” Physical Review A 98, 043801 (2018).
  * Lloyd-Hughes  _et al._ (2021) J. Lloyd-Hughes, P. M. Oppeneer, T. Pereira dos Santos, A. Schleife, S. Meng, M. A. Sentef, M. Ruggenthaler, A. Rubio, I. Radu, M. Murnane, X. Shi, H. Kapteyn, B. Stadtmüller, K. M. Dani, F. H. da Jornada, E. Prinz, M. Aeschlimann, R. L. Milot, M. Burdanova, J. Boland, T. Cocker, and F. Hegmann, “The 2021 ultrafast spectroscopic probes of condensed matter roadmap,” J. Phys.: Condens. Matter 33, 353001 (2021).
  * Hone  _et al._ (1997) Daniel W. Hone, Roland Ketzmerick, and Walter Kohn, “Time-dependent floquet theory and absence of an adiabatic limit,” Phys. Rev. A 56, 4045–4054 (1997).
  * Schäfer  _et al._ (2019) C. Schäfer, M. Ruggenthaler, V. Rokaj, and A. Rubio, “Relevance of the quadratic diamagnetic and self-polarization terms in cavity quantum electrodynamics,” arXiv:1911.08427 [quant-ph] (2019).
  * Luttinger (1951) J. M. Luttinger, “The effect of a magnetic field on electrons in a periodic potential,” Phys. Rev. 84, 814–817 (1951).
  * Pfannkuche and Gerhardts (1992) Daniela Pfannkuche and Rolf R. Gerhardts, “Theory of magnetotransport in two-dimensional electron systems subjected to weak two-dimensional superlattice potentials,” Phys. Rev. B 46, 12606–12626 (1992).
  * Langbein (1969) D. Langbein, “The tight-binding and the nearly-free-electron approach to lattice electrons in external magnetic fields,” Phys. Rev. 180, 633–648 (1969).
  * Rauh (1975) A. Rauh, “On the broadening of Landau levels in crystals,” phys. stat. sol. (b) 69, K9–K13 (1975).
  * Petschel and Geisel (1993) G. Petschel and T. Geisel, “Bloch electrons in magnetic fields: Classical chaos and Hofstadter’s butterfly,” Phys. Rev. Lett. 71, 239–242 (1993).
  * Wannier (1978) G. H. Wannier, “A result not dependent on rationality for Bloch electrons in a magnetic field,” physica status solidi (b) 88, 757–765 (1978).
  * Brown (1964) E. Brown, “Bloch electrons in a uniform magnetic field,” Phys. Rev. 133, A1038–A1044 (1964).
  * Zak (1964a) J. Zak, “Magnetic translation group,” Phys. Rev. 134, A1602–A1606 (1964a).
  * Zak (1964b) J. Zak, “Magnetic translation group. II. Irreducible representations,” Phys. Rev. 134, A1607–A1611 (1964b).
  * Griffiths (1995) D. J. Griffiths,  _Introduction to Quantum Mechanics_ (Prentice Hall, 1995).
  * Tong (2016) D. Tong, “Lectures on the quantum Hall effect,” arXiv:1606.06687 [hep-th] (2016).
  * Peierls (1955) R. E. Peierls,  _Principles of the Theory of Solids_ (Oxford University Press, 1955).
  * Barrier  _et al._ (2020) J. Barrier, P. Kumaravadivel, R. Krishna Kumar, L. A. Ponomarenko, N. Xin, M. Holwill, C. Mullan, M. Kim, R. V. Gorbachev, M. D. Thompson, J. R. Prance, K. Taniguchi, K. Watanabe, I. V. Grigorieva, K. S. Novoselov, A. Mischenko, V. I. Fal’ko, A. K. Geim, and A. I. Berdyugin, “Long-range ballistic transport of Brown-Zak fermions in graphene superlattices,” Nature 11, 5756 (2020).
  * Ikeda (2018) K. Ikeda, “Hofstadter’s butterfly and Langlands duality,” J. Math. Phys. 59, 061704 (2018).
  * Lamb and Retherford (1947) Willis E. Lamb and Robert C. Retherford, “Fine structure of the hydrogen atom by a microwave method,” Phys. Rev. 72, 241–243 (1947).