Multiscale modeling and simulation of turbulent flows in porous media Open Access

A porous medium is a material that contains pores (voids). The skeletal part of the porous medium is often solid and is therefore called the solid matrix or porous matrix. The pores are often filled with fluids (gases or liquids). The transport of fluids in porous media is an important process in many applications, including CO₂ sequestration in deep saline aquifers (Benson and Cole, 2008; Huppert and Neufeld, 2014), energy extraction from geothermal reservoirs or the ambient ground (Ghoreishi-Madiseh , 2013), enhanced oil/gas recovery (Gbadamosi , 2019), transport of groundwater contaminants (Patil and Chore, 2014), and stability of snow against avalanches (Nield and Simmons, 2019).

Convection in porous media is usually laminar, since it often has a low pore-scale Reynolds number (Re_d) defined by the pore size d and the mean flow velocity u_m. However, when Re_d is of the order of 100 or higher, the flow within the pores becomes turbulent. An example of an application to which this is relevant is thermal energy storage systems that use rocks/bricks for storing thermal energy (Odenthal , 2014). The processes of heat charging and discharging in such energy storage systems are often driven by mixed (combined forced and natural) convection. Rocks or bricks have low thermal conductivity, resulting in very slow processes of charging and discharging. To overcome this disadvantage, the size of the porous element and the velocity of the fluid may be adjusted to make the flow fully turbulent. If the porous medium can be approximated by a bank of tubes, then the relationship between the Nusselt number and the Reynolds number changes from $Nu∼Red0.36$ for Re_d < 300 to $Nu∼Red0.64$ for a fully turbulent flow (Re_d > 300); see Engineering Sciences Data Unit (1984). The heat transfer is efficiently enhanced by the turbulence.

Turbulence in porous media might also occur in packed bed reactors in chemical engineering problems (Shams , 2014; Dave , 2018; and Lucci , 2017). For example, the pore-scale Reynolds number for a nuclear pebble bed reactor can reach the order of 10⁵ (Dave , 2018), which is considerably higher than the Reynolds number for laminar–turbulent transition. Likewise, flows in a catalyst packed bed can also reach the order of 10⁴ (Lucci , 2017).

Another important application of turbulence in porous media is in canopy flows. The porous canopies can be made of trees, vegetation, buildings, or wind farms (Meroney, 2007). Canopies usually have large porosity values ϕ [typically ϕ > 90%; see Ghisalberti and Nepf (2009)] and extremely high Reynolds numbers. For example, Belcher (2005) studied canopy flows in urban areas, which have pore-scale Reynolds number Re_d = 5000. Turbulence in plant canopies can have Reynolds numbers up to 10⁵ in the upper part of well-ventilated canopies and as low as 10² near the ground (Belcher , 2012).

Besides the classical turbulent flows introduced above, flow instabilities might also occur in natural convection in porous media at high Rayleigh numbers. Natural convection flows in porous media at large Rayleigh numbers are transient, chaotic, and random; these characteristics are similar to those of classical turbulence. However, transient natural convection in porous media can occur at very low pore-scale Reynolds numbers. These transient flows are called “pseudo-turbulence,” which is a term used to distinguish them from classical turbulence. “Pseudo-turbulence” of natural convection in porous media has important applications in CO₂ sequestration (Huppert and Neufeld, 2014), geothermal energy (Ghoreishi-Madiseh , 2013), and enhanced oil/gas recovery (Gbadamosi , 2019).

Driven by the need to solve industrial problems and understand physical problems in nature, there have been extensive studies of turbulent flows in porous media in the past few years. Computational fluid dynamics (CFD) is an important tool for understanding this process. According to the difference in the time and length scales of motions to be solved, CFD simulations of convection in porous media can be classified into microscopic and macroscopic simulations. In microscopic simulations, the motions smaller than the representative elementary volume (REV) should be calculated; thus, the detailed geometry of the porous elements must be considered. A pore-scale-resolved direct numerical simulation with a perfectly stirred reactor model (PSR-DNS) is a typical microscopic simulation. In a PSR-DNS, the microscopic governing equations (usually the Navier–Stokes equations) are solved directly without the introduction of additional models. A PSR-DNS method usually has higher accuracy than other simulation methods. However, only recently has it become possible to study turbulent flows in porous media with PSR-DNS methods owing to their extremely high computational cost (Chandesris , 2013; Jin , 2015; Uth , 2016; Chu , 2019; Rao and Jin, 2022; Wang , 2021; 2022; and Srikanth , 2021). Pore-scale-resolved large eddy simulation (PSR-LES) is another common microscopic simulation method. It has important applications owing to its moderate computational cost (Hutter , 2011; Shams , 2014; and Huang , 2022).

Compared with microscopic simulations, macroscopic simulations of porous medium convection are more efficient. The macroscopic equations can be obtained by time and volume averaging (or sometimes only volume averaging) of the microscopic governing equations. Macroscopic simulation is the most widely used method for calculating turbulent flows in porous media in engineering applications owing to its low computational cost (Lee and Howell, 1991; Prescott and Incropera, 1995; Lage , 1997; Kazerooni and Hannani, 2009; Kundu , 2014; and de Lemos and Braga, 2003). The problem with macroscopic simulations is that strong assumptions are often used to close the averaged momentum equation, leading to inaccuracy and uncertainty of the numerical solutions.

Despite the efforts made in recent years, microscopic simulation and macroscopic modeling of turbulent flows in porous media remain challenging tasks owing to the high complexity of the problem as well as the difficulty in validation. The purpose of this paper is to review the important work on multiscale modeling and simulation of turbulent flows carried out in the past few years. In the review, we consolidate and systemically analyze the knowledge on this topic that has been accumulated in the past years. We also discuss important unsolved problems that will need to be addressed in future work.

The structure of the paper is as follows. Following this introduction, microscopic simulation of turbulent flows in porous media is discussed in Sec. II. Then, the governing equations for macroscopic simulations are introduced in Sec. III. Section IV focuses on modeling of the macroscopic equations. Two types of turbulence models, macroscopic and microscopic turbulence models, are discussed in this section. Finally, a summary of the established knowledge and future issues is given in Sec. V.

Microscopic simulation of flows in porous media concerns simulations in which the pore-scale flows are directly resolved. DNS is the most accurate method for microscopic simulation. It is defined as a CFD simulation in which the Navier–Stokes equations are solved directly without the introduction of a turbulence model. DNS of turbulent flows in porous media requires the detailed pore-scale geometry to be resolved; thus, this DNS is also called pore-scale-resolved DNS (PSR-DNS).

A subgrid-scale (SGS) modeling term is sometimes added to the momentum equation to model small eddies to reduce the computational cost. The simulation then becomes a large eddy simulation (LES). LES also belongs to the class of microscopic simulations when the motions below the pore scale are resolved. Theoretically, the pore-scale motions can also be solved by using Reynolds-averaged Navier–Stokes (RANS) simulations (Soulaine , 2017; Ferdos and Dargahi, 2016; and Dave , 2018). Jin and Herwig (2015) simulated flows in channels with rough walls (where the walls were covered with porous elements) using eight RANS models and compared the results obtained with those from DNS. The comparison revealed that the errors in the calculated friction coefficients exceeded 20% in all-the RANS results. It was demonstrated that even k-ε and Reynolds stress models could predict incorrect trends. Therefore, it is expected that PSR-RANS results may exhibit considerable uncertainties due to model errors when used to simulate turbulent flows in porous media. PSR-DNS and PSR-LES are more suitable for understanding the physics of turbulent flows in porous media.

We focus on continuum fluids (with Knudsen number Kn < 0.01) in this review. Therefore, standard boundary conditions can be applied in the PSR-DNS and PSR-LES. No-slip boundary conditions are imposed at the surfaces of solid obstacles. Most DNS and LES studies avoid using inlet boundary conditions, primarily because it is challenging to prescribe turbulent fluctuations at the boundaries. Periodic or no-slip boundary conditions are often utilized at the boundaries of porous media. With these boundary conditions, it is possible to obtain fully developed (statistically steady) flows for statistical analysis.

To resolve the intricate pore-scale geometry, PSR-DNS and PSR-LES are often conducted using a finite volume method (FVM), as demonstrated by Jin  (2015), Gasow  (2020), and Wang  (2021, 2022), or a lattice Boltzmann method (LBM), as shown by Liu  (2021; 2023), Li  (2022), and Diao  (2023). Finite difference methods and spectral-type methods are seldom used, owing to difficulties in resolving pore-scale geometry. Among these methods, the LBM, which has emerged in recent years, has demonstrated particular potential in PSR-DNS and PSR-LES of flows in porous media owing to its excellent parallelizability. However, further efforts are needed to enhance this method for broader applications, such as calculating compressible flows and using body-fitted meshes.

Many microscopic simulation studies investigate flows in porous media with homogeneously distributed porous elements. In periodic problems, the computational domain is composed of representative elementary volumes (REVs). REVs are appropriately sized to accurately reflect the pore-scale properties of the entire volume in question (i.e., the porosity and permeability do not change if the volume is further enlarged). REVs in fundamental studies are often constructed from simple geometries; for example, REVs can be made of staggered arrangement of cylinders (Fig. 1) or spheres.

PSR-DNS and PSR-LES require the resolution of motions within both the pores and macroscopic elements, making them computationally intensive. For instance, Jin  (2015) employed 10⁷–6 × 10⁷ mesh cells to resolve each REV consisting of aligned bars, totaling 5 × 10⁷–3.8 × 10⁸ mesh cells for simulating turbulent flows in a matrix comprising six REVs. Similarly, Uth  (2016) used 4800 mesh cells to resolve each REV composed of three-dimensional spheres. Wang  (2021; 2022) utilized 8.8 × 10⁷–5.8 × 10⁸ mesh cells to simulate turbulent flows partially obstructed by a porous medium. In a DNS study of pseudo-turbulence, Gasow  (2020; 2021; and 2022) employed 3600 mesh cells to resolve a two-dimensional REV, totaling up to 7.2 × 10⁸ mesh cells.

An important motivation for microscopic simulations is to better understand the dynamics of turbulent eddies in porous media. This is crucial for predicting dissipation rates in both flow and temperature fields. According to the derivation by Jin and Herwig (2014), these dissipation rates determine both pressure drop and heat transfer rates. PSR-DNS and PSR-LES offer more accurate predictions of dissipation rates compared with RANS, since dissipation is either directly calculated (in PSR-DNS) or modeled with minimal assumptions (in PSR-LES). With a better grasp of the underlying physics, more precise modeling of turbulence transport in porous media becomes feasible.

A crucial fundamental question about the eddy dynamics in porous media is whether turbulence models should account only for turbulence structures limited in size to the pore scale or whether the size of turbulence structures could exceed the pore scale. The latter would imply the existence of macroscopic turbulence in porous media, where the turbulent eddies are larger than the size of the pores. To address this fundamental question, Jin  (2015) conducted a DNS study. Both transient and statistical results demonstrate that the length scales of turbulent structures are constrained by the pore size, thereby suggesting the pore-scale prevalence hypothesis (PSPH) when the porosity is not close to unity. Later, Uth  (2016) and Jin and Kuznetsov (2017) simulated turbulent flows in more complicated porous matrices. The DNS results further confirmed the PSPH. Figure 2 shows snapshots of turbulent structures in porous media with different pore-scale geometries. The snapshots show qualitatively that the length scales of turbulence are limited by the pore size. Rao and Jin (2022) indicate that the macroscopic turbulence may survive when the porosity is larger than a critical value, which is in the range 0.93–0.98.

Some other DNS studies have investigated the interaction between a turbulent flow in porous media and a free turbulent flow. The computational domain is partly occupied by a porous medium. For example, Motlagh and Taghizadeh (2016) studied the transverse pressure waves caused by the Kelvin–Helmholtz (K–H) instability over porous and rough walls. Kuwata and Suga (2016; 2017) studied turbulence over anisotropic porous walls. The arrays of porous elements are displayed systematically to introduce anisotropic permeability. The DNS results indicate that permeability in the streamwise (flow) and spanwise (cross-flow) directions enhances turbulence, while permeability in the wall-normal direction does not. Permeability has been shown to play a major role in determining the response of channel flow to permeable walls in DNS studies by Rosti et al. (2015, 2018), even if the permeability has a very small value. These studies also showed that making porous media anisotropic could reduce drag across a permeable wall by about 20%. The DNS results by Gómez-de-Segura and García-Mayoral (2019) reveal a linear relationship between drag reduction and the difference between the streamwise and spanwise permeabilities when the permeability is small. This linear relationship holds when the wall-normal permeability normalized by the viscous length scale $Ky+$ is above a critical value (⁠$Ky+≥0$⁠.38). Interpreting the DNS results with the transfer entropy concept, Wang  (2021; 2022) investigated the effects of permeability on the intensity, time scale, and spatial extent of surface–subsurface interactions. A neural network-based remote sensing model was used to validate the statistical results.

Flow and heat transfer in packed beds are also important topics for microscopic simulations. For example, Hutter  (2011) simulated flows in a periodic open-cell structure in the turbulent regime using LES. Pore-scale Reynolds numbers in the range between 1200 and 4500 were considered. Shams  (2014) performed an LES of flows in a randomly stacked bed composed of spherical pebbles with a Reynolds number of 9753. Since DNS is computationally very expensive at high Reynolds numbers, the flow at a high pore-scale Reynolds number can be simulated by using LES. He  (2018; 2019) performed DNS studies in a triply periodic unit cell of a face-centered cubic (FCC) lattice. The DNS results confirm the absence of macroscale turbulence structures for this configuration. In addition, the up-scaled flow statistics are well captured by the REV even if the flow in the pore is strongly anisotropic. Using both LES and DNS methods, Srikanth  (2021) studied the mechanism of a symmetry breaking phenomenon that occurs in turbulent flow in periodic porous media at a low value of porosity. The symmetry breakdown causes a deviation in the macroscale flow direction, which in turn introduces anisotropy in the macroscale Reynolds stress. Using LES, Huang  (2022) studied the effects of microscale vortices on convective heat transfer in turbulent flow within porous media. The numerical results show that both vortex shedding and secondary flow instabilities have significant influences on the Nusselt number. Local Nusselt number peaks coincided with the regions of flow mixing surrounding flow stagnation points on the solid obstacle surface. Rao and Jin (2022) simulated turbulent flows in a porous matrix characterized by two length scales using DNS. The DNS results indicate that macroscopic turbulence structures may survive when the porosity is very high. Apte  (2022) studied the clustering of inertial particles in turbulent flow through a porous unit cell. The DNS results show a large number of small volumes where cluster destruction is prominent. These small volumes are generated by the collision of particles with the wall surfaces of the porous elements.

Besides the conventional turbulent flows in porous media introduced above, it is also necessary to mention the pseudo-turbulence arising in natural convection in porous media. Figure 3 shows snapshots of the instantaneous local Reynolds number Re_K, which is based on permeability K and local velocity magnitude $u$⁠. These flows are usually characterized by very low local Reynolds numbers, and so the Forchheimer term, which characterizes the contribution of the form drag, can be neglected. However, the flow exhibits chaotic behavior at a large Rayleigh number. Conventional turbulent flows do not exhibit this behavior. Intensive DNS studies have been carried out in recent years to understand pseudo-turbulence in porous media. For example, based on a microscopic simulation study, Gasow  (2020; 2021; and 2022) found that the effects of pore-scale motions on the macroscopic instabilities of natural convection cannot be neglected. Liu  (2020) found that the scaling of the Nusselt number vs Ra is significantly affected by the ratio of the pore scale to the thickness of the thermal boundary layer. Korba and Li (2022) studied the pore-scale effects on natural convection induced by conjugate heat transfer.

A significant real-world application of pseudo-turbulence in natural convection in porous media is the CO₂ sequestration process in deep saline aquifers (Benson and Cole. 2008; Huppert and Neufeld, 2014). However, the actual CO₂ sequestration process is considerably more complex than those examined in fundamental studies. One crucial aspect to consider is the influence of multiphase flows. First, the presence of different phases affects the flow dynamics in porous media owing to surface tension effects. Additionally, the permeability of the porous medium varies when multiple phases are involved. Despite the increasing attention paid to multiphase flows in porous media in recent years, as evidenced by studies such as those by Gruener  (2012), Sun  (2016), and Li (2022), there remains a notable gap in research concerning turbulence (or pseudo-turbulence) in multiphase flows. In addition, the Darcy numbers used in the microscopic simulations are much higher than those in real applications. Calculations with lower Darcy numbers are required to provide better understanding of pseudo-turbulence in porous media.

The overbar $ā$ denotes time averaging and the angular brackets $⋅i$ denote volume averaging over the fluid region of an REV. The operation $φi$ is defined as $φi=φ−⟨φ⟩i$⁠. The different sequences of time and volume averaging result in two types of expressions; however, these are equivalent. In Eq. (4), $⟨φ̄⟩i=⟨φ⟩ī$ is the intrinsic average of the time-mean property, $⟨φ′⟩i=⟨φ⟩i′$ is the intrinsic average of the property fluctuation, $φ̄i=φī$ is the spatial deviation of the mean property, and $φ′i$ is the spatial deviation of the property fluctuation.

The modeling of the unknown terms in the macroscopic equations (5)–(8) has been extensively studied. Early studies focused on modeling the total drag, while turbulence in porous media has received little attention apart from its effect on the total drag. Some classical models for the total drag are still widely used in modern CFD simulations. We will start our discussion with these classical models, which will be introduced in Sec. IV A.

Modeling turbulence in porous media has received increasing attention in recent years. It has been reviewed in the paper by Wood et al. (2020), the books by de Lemos (2012) and Nield and Bejan (2017), and in many review chapters, such as those by Lage et al. (2002), de Lemos (2004; 2005), and Vafai (2005). In general, there are two distinct views on turbulence in porous media. The first class of models are called macroscopic turbulence models, which deal primarily with macroscopic turbulence in porous media. The second approach was originally expressed in Nield (1991) and then further developed in other studies, such as those by Nield (2001), Jin  (2015), Uth  (2016), and Jin and Kuznetsov (2017). According to the second approach, macroscopic turbulence is impossible owing to the pore-scale limitation on the size of turbulent eddies, at least in a dense porous medium. All turbulence models according to this view are microscopic turbulence models. Macroscopic turbulence models and microscopic turbulence models will be discussed in detail in Secs. IV B and IV C.

Information about turbulence transport is required in a macroscopic turbulence model. Lee and Howell (1991), Prescott and Incropera (1995), and Lage  (1997) reported representative models of this type. These macroscopic models are similar to turbulence models for clear fluid flows. The transport equations for turbulent kinetic energy and dissipation rate are often solved in these models. In addition, some turbulence models attempt to simulate both pore-scale and large-scale turbulence. Thus, macroscopic transport of turbulence is also accounted for in these models. Representative models are those developed by de Lemos and colleagues; see de Lemos and Braga (2003), de Lemos and Mesquita (2003), de Lemos and Pedras (2001), de Lemos and Rocamora (2002), de Lemos and Tofaneli (2004), Pedras and de Lemos (2000; 2001a; 2001b; 2001c; 2003), and Silva and de Lemos (2003).

To close the macroscopic equations, the tensors $R̄ijII$⁠, $R̄ijIII$⁠, and $R̄ijIV$ for the time- and volume-averaged equations need to be modeled. It should be noted that even if the flow is not turbulent, the spatial dispersion of the microscopic mean velocity $R̄ijII$ still needs to be modeled when the mean velocity gradient is not zero. However, $R̄ijII$ has received very little attention in modern turbulence models.

Macroscopic turbulence models have been widely applied to the solution of many engineering problems. For example, they have been applied to the study of flows in packed beds, in channels, in pipes, and in bundles of rods (Nakayama and Kuwahara, 2008). They have also been used in the simulation of flow in a porous matrix, which is a periodic array of square cylinders (Kazerooni and Hannani, 2009; Kundu , 2014). The macroscopic models have a common characteristic, namely, the transport of macroscopic turbulence in porous media is simulated in these models.

Besides the macroscopic models for the simulation of turbulence in porous media, there are other types of turbulence models, namely, microscopic turbulence models. In microscopic turbulence models, the effect of turbulence can be accounted for by the source terms in the momentum and energy equations. Thus, it is not necessary to calculate the transport of turbulence. These types of turbulence models originate from Nield (1991). The models developed by Nakayama and Kuwahara (1999) and Kuwahara (1998) also belong to the class of microscopic turbulence models.

Later, Jin  (2015), Uth  (2016), and Jin and Kuznetsov (2017) carried out systematic DNS studies to understand turbulence in porous media. Porous matrices composed of two-dimensional square and circular cylinders, three-dimensional spheres unbounded or bounded by walls, and a three-dimensional matrix with two characteristic pore scales have been investigated. The DNS results obtained suggest that the size of turbulent vortices is limited by the pore size, leading to the proposition of the pore-scale prevalence hypothesis (PSPH).

Natural convection in porous media is usually simulated by solving the Darcy–Oberbeck–Boussinesq (DOB) equations. These are an alternative form of the macroscopic equations (30) and (31) (Nield and Bejan, 2017). The macroscopic instabilities in natural convection are induced by a buoyancy force, which is usually modeled with the Boussinesq hypothesis. Based on the volume-averaged Navier–Stokes equations (30) and (31), Gasow  (2021) developed a two-length-scale diffusion (TSLD) model for macroscopic simulations. In contrast to the DOB model, macroscopic diffusion is accounted for in the TLSD model. Numerical studies have shown that the macroscopic instabilities of natural convection in porous media can be accurately simulated by the volume-averaged equations (Hewitt , 2012; 2013; 2014; Wen , 2015; De Paoli et al., 2016¹¹⁷; Kränzien and Jin, 2019; Pirozzoli , 2021; and Gasow , 2021; 2022).

Figure 4 shows a snapshot of the local Reynolds number Re_K obtained from a macroscopic simulation with the TLSD model. Structures made of small proto-plumes near the walls and large mega-plumes near the domain center can be observed. The macroscopic and DNS results have been compared in Gasow  (2021; 2022). However, natural convection in porous media at high Rayleigh numbers has only macroscopic instabilities. There is no microscopic turbulence, owing to the low local Reynolds number. In this respect, the macroscopic instability for natural convection in porous media is called pseudo-turbulence, which is different from conventional turbulence.

Similar to classical turbulence, pseudo-turbulence is transient, chaotic, and random. However, whether pseudo-turbulence qualifies as real turbulence remains an open question. It is still unclear whether pseudo-turbulence follows the same energy cascade pattern as classical turbulent flows, which decay with a 5/3 scaling at large Rayleigh numbers. Answering this question is crucial to determine the potential occurrence of macroscopic turbulence. Additionally, investigating whether the volume-averaged Eqs. (30) and (31) can accurately predict conventional macroscopic turbulence should be a focus of future research.

With the rapid development of high-performance computers, CFD methods are becoming an efficient tool for analyzing and predicting turbulent flows in porous media. Microscopic simulations, including pore-scale-resolved DNS and LES, have considerably improved the understanding of turbulent flows in porous media, because they provide extensive flow details, including flow details in the pores.

In the past few years, the most active microscopic simulation studies of flows in porous media have focused on the pore-scale turbulent structures, turbulence length/time scales, turbulence at the interface of a porous medium and a clear flow region, and effects of turbulence on heat/mass transfer. These studies have provided valuable knowledge for developing more efficient and accurate models for macroscopic simulations.

With time and volume averaging, or sometimes only volume averaging, the macroscopic equations for flows in porous media can be derived. It is more efficient to calculate the flows in porous media by solving the macroscopic equations; the corresponding simulations are called macroscopic simulations. Turbulence models for macroscopic simulations can be classified as macroscopic turbulence models, in which the macroscopic turbulence is simulated, or microscopic turbulence models, in which turbulence is simulated from direct solutions of the Navier–Stokes equations. DNS simulations performed by the present authors and their colleagues suggest that the pore scale limits the size of turbulent structures in porous media, unless the porosity is very large. With the new understanding of flow physics, new turbulence models for flows in porous media have emerged, including a microscopic turbulence model based on the PSPH.

Despite these advances, macroscopic modeling of turbulent flows in porous media is still a challenging task because of the complexity of the problem. Microscopic (DNS or LES) simulations are an important tool for reaching this target. Through microscopic simulations, the following fundamental questions need to be answered, which are key to developing more accurate and efficient macroscopic models:

1. How should the transport of macroscopic turbulence in porous media be modeled when it occurs? Although the DNS studies that have been conducted show no macroscopic turbulence unless the porosity is very high, there is evidence that macroscopic turbulence can occur at very large Reynolds numbers and high porosity values. In addition, it is important to know whether the pseudo-turbulence arising in natural convection in porous media is real turbulence.

2. What are the characteristics of turbulence at very high Reynolds numbers? The typical pore-scale Reynolds numbers in current DNS studies are of the order of O(10³). However, the pore-scale Reynolds numbers for nuclear pebble reactors can reach the order of O(10⁵). DNS studies of such problems are far beyond the capabilities of current high-performance computers. Whether turbulent flows at extremely high Reynolds numbers have new characteristics remains a mystery.

3. How does turbulence transport occur in inhomogeneous and anisotropic porous media? Turbulence in inhomogeneous and anisotropic porous media has different length/time scales. How motions with different length/time scales interact with each other is still unknown. A specific situation is the flow in a domain partly occupied by a porous medium and partly by a clear fluid. Modeling turbulent flow in inhomogeneous and anisotropic porous media remains a challenging task.

4. How does compressibility affect turbulent flows in porous media? Compressibility of the fluid in porous media flows has yet to be considered, partly because most previous studies have considered relatively low flow velocities. However, even if the Mach number in a porous medium is not very high, fluid compressibility still might have important effects when the temperature or pressure has strong variations. At present, very little is known about compressible flows in porous media.

5. How does turbulence or pseudo-turbulence occur, and how are the turbulence kinetic energy and thermal energy transported in multiphase flow? It is particularly important to understand how surface tension and variations in permeability caused by different phases affect turbulent flow in porous media.

6. How should fluid–structure interactions in porous media be modeled and simulated? Porous matrices may become deformable at very large Reynolds numbers. The ways in which turbulent flows interact with deformable structures should be investigated in greater depth in the future.

The authors gratefully acknowledge the support for this study by the Deutsche Forschungsgemeinschaft (DFG, Grant Nos. 408356608 and 262377115). A.V.K. acknowledges the support of the National Science Foundation (Award No. CBET-2042834) and the Alexander Humboldt Foundation Research Award.

The authors have no conflicts to disclose.

Yan Jin: Conceptualization (equal); Data curation (lead); Funding acquisition (equal); Methodology (equal); Project administration (equal); Writing – original draft (equal); Writing – review & editing (equal). Andrey V. Kuznetsov: Conceptualization (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.