1. Introduction
Convection in fluid–porous systems is a universally observed phenomenon, with applications arising in technological, geophysical and astrophysical settings. In technological and industrial applications, fluid–porous convection is prevalent in heat sinks and cooling technologies found in laptops and computers (Yu, Lee & Yook Reference Yu, Lee and Yook2010, Reference Yu, Lee and Yook2011; Al-Zamily Reference Al-Zamily2017), as well as in the solidification of alloys (Le Bars & Worster Reference Le Bars and Worster2006a,Reference Le Bars and Worsterb). Geophysical examples can be seen in the coupled fluid–porous flow of subglacial or dry salt lakes (Hirata, Goyeau & Gobin Reference Hirata, Goyeau and Gobin2012; Couston Reference Couston2021; Lasser, Ernst & Goehring Reference Lasser, Ernst and Goehring2021), in carbon dioxide sequestration or the flow of oil in underground reservoirs (Allen Reference Allen1984; Ewing Reference Ewing1997; Huppert & Neufeld Reference Huppert and Neufeld2014) and in contaminant transport in subsoil water reservoirs (Curran & Allen Reference Curran and Allen1990; Allen & Khosravani Reference Allen and Khosravani1992). Other geophysical applications that feature natural convection include plate tectonics (Zhang & Libchaber Reference Zhang and Libchaber2000; Mac Huang et al. Reference Mac Huang, Zhong, Zhang and Mertz2018), cave ventilation (Khazmutdinova et al. Reference Khazmutdinova, Nof, Tremaine, Ye and Moore2019) and morphological formation from solute-laden flows (Wykes et al. Reference Wykes, Mac Huang, Hajjar and Ristroph2018; Mac Huang et al. Reference Mac Huang, Tong, Shelley and Ristroph2020; Mac Huang & Moore Reference Mac Huang and Moore2022). The phenomenon of convection also extends far beyond Earth and into astrophysical applications, such as in moons of Saturn and Jupiter (Choblet et al. Reference Choblet, Tobie, Sotin, Běhounková, Čadek, Postberg and Souček2017; Le Reun & Hewitt Reference Le Reun and Hewitt2020; Vilella et al. Reference Vilella, Choblet, Tsao and Deschamps2020).
Although flow in fluid–porous systems has been a staple of the research community since Saffman's and Jones’s work in the 1970s (Saffman Reference Saffman1971; Jones Reference Jones1973), convection in these systems still poses unique and timely questions. Recent years have seen a resurgence of research on fluid–porous convection from a variety of viewpoints, including conducting nonlinear stability analysis, exploring bifurcations and developing stable numerical schemes for solutions (McCurdy, Moore & Wang Reference McCurdy, Moore and Wang2019; Chen et al. Reference Chen, Han, Wang and Zhang2020; Han, Wang & Wang Reference Han, Wang and Wang2020; Chen et al. Reference Chen, Han, Wang and Zhang2022; Wang & Wu Reference Wang and Wu2021). One recurring theme observed in many of these works is the contrast between deep convection, in which convection cells occupy the entirety of the coupled domain, and shallow convection, in which cells only circulate in the free-flow region.
The arrangement that we consider is a saturated porous layer lying beneath a clear fluid region that is free from obstructions, illustrated with the schematic in figure 1. This entire coupled system is heated from below. Each domain has a governing set of equations for the fluid flow – Navier–Stokes for the free-flow region and Darcy's equations for the porous medium – with a set of conditions imposed along the interface between the two. If the temperature difference between the lower and upper plates rises above a critical threshold, then the conductive state becomes unstable and gives way to natural convection. If the temperature difference is sufficiently small, the conductive state remains stable.
Figure 1. Schematic of the domain 
$\varOmega = \{(x,y)\in \mathbb {R}^2 \times z\in (-d_m,d_f)\}$, comprising a free-flow region 
$\varOmega _f$ and a porous medium 
$\varOmega _m$. The two subdomains meet at an interface 
$\varGamma _i$. The upper and lower boundaries are impermeable and held at constant temperatures 
$T_U$ and 
$T_L$, respectively, with 
$T_L>T_U$. We assume periodicity of the velocity and temperature in the horizontal direction(s) as well.
In previous work, McCurdy et al. (Reference McCurdy, Moore and Wang2019) conducted linear and nonlinear stability analyses of the superposed system. With certain parameter regimes, the bimodal marginal stability curves suggested that a small change in the depth ratio of the two regions could trigger a drastic change in the convection patterns. The interested reader can look ahead to figure 2(a,b) to see the stark contrast between the flow profiles for depth ratios of 
$\hat {d} = 0.18$ and 
$\hat {d} = 0.19$, respectively. Indeed, slightly altering the depth ratio induces a qualitative shift in flow behaviour, as originally observed by Chen & Chen (Reference Chen and Chen1988). This drastic change spurred McCurdy et al. (Reference McCurdy, Moore and Wang2019) to develop a simple, coarse-grained model to narrow down the parameter ranges for which the transition occurs. This simple model, which neglected any coupling between the free-flow and porous regions, provided promising results and laid the groundwork for the current study.
Figure 2. Marginally stable flow configurations and temperature profiles (colour) for (a) 
$\hat {d}=0.18$ (deep convection) and (b) 
$\hat {d}=0.19$ (shallow convection), with 
$\sqrt {Da} = 5.0\times 10^{-3},\ Pr_{m} =0.7,\ \epsilon _T=0.7,\ \alpha =1.0$ fixed. (c) Marginal stability curves while varying 
$\hat {d}$ from 
$\hat {d}=0.15$ to 
$\hat {d}=0.22$ by increments of 0.05 with the respective critical Rayleigh numbers 
$Ra_m^*$ shown in red. As 
$\hat {d}$ crosses the critical depth ratio of 
$\hat {d}^*=0.181$ for this parameter regime, the most unstable wavenumber jumps from 
$a_m\approx 2.1$ to 
$a_m\approx 14.5$. This signifies a sudden transition from deep to shallow convection (at the onset of convection) as 
$\hat {d}$ increases from 
$\hat {d}<\hat {d}^* \rightarrow \hat {d} > \hat {d}^*$.
Here, we improve upon the model to better account for the flow conditions at the interface between the fluid and the porous medium. The analysis results in a simple, explicit formula for the critical depth ratio at which shallow convection transitions to deep convection. We expect the new formula to be relevant for geophysical applications, such as predicting the penetration of tracers into groundwater, and industrial applications, for instance controlling heat dissipation by appropriately choosing the depth and/or porosity of a heat sink. To test the new model, we conduct numerical simulations of the fully nonlinear, coupled Navier–Stokes–Darcy–Boussinesq system. Computing numerical solutions to this system presents several challenges, for example accurately representing the sharp transitions in physical properties (e.g. density, conductivity, diffusivity) across interfaces and achieving stable time-stepping in face of the nonlinearities present in the Navier–Stokes equations. As the numerical scheme is detailed, we explain how we address each of these challenges. Ultimately, the fully nonlinear numerical simulations provide a more comprehensive picture of fluid–porous convection, revealing novel flow configurations not easily predicted by stability analysis or the coarse-grained model.
A few recent works have focused on convection in superposed fluid–porous layers from a numerical perspective. Zhang, Shan & Hou (Reference Zhang, Shan and Hou2020) used a finite element method (FEM) to study the stationary Navier–Stokes–Darcy–Boussinesq system and investigated the well-posedness of their finite element approximation. Other works have numerically examined convection in related systems (Tatsuo et al. Reference Tatsuo, Toru, Mitsuhiro, Yuji and Hiroyuki1986; Valencia-Lopez & Ochoa-Tapia Reference Valencia-Lopez and Ochoa-Tapia2001; Al-Zamily Reference Al-Zamily2017; Le Reun & Hewitt Reference Le Reun and Hewitt2021), albeit with different governing equations (e.g. the Brinkman system instead of the Darcy system, or the Stokes equations instead of Navier–Stokes). Some works have explored schemes to couple the Cahn–Hilliard equations with fluid–porous flow (Chen et al. Reference Chen, Han, Wang and Zhang2020, Reference Chen, Han, Wang and Zhang2021), while others have taken an analytical approach to convection in coupled layers. For example, Han et al. (Reference Han, Wang and Wang2020) recently examined transitions in the same Navier–Stokes–Darcy–Boussinesq system considered here, although through a different lens. Their main focus was the transition from a conductive to a convective state as the Raleigh number increases. The work of Han et al. (Reference Han, Wang and Wang2020) rigorously showed that transitions exist between different convective profiles, like deep and shallow convection, and noted how the transitions behaved – as continuous transitions or jump transitions – around their critical Rayleigh numbers. In this work, we examine a similar kind of transition as we develop a model to predict the parameter regimes where the switch between deep and shallow convection takes place. This kind of transition, dubbed a ‘dynamic transition’ by Han et al. (Reference Han, Wang and Wang2020), is a bifurcation of the system as one moves through the parameter space and has been observed in a number of papers (Chen & Chen Reference Chen and Chen1988, Reference Chen and Chen1989, Reference Chen and Chen1992; McKay Reference McKay1998; Straughan Reference Straughan2002; Hirata et al. Reference Hirata, Goyeau, Gobin, Carr and Cotta2007; Yin, Fu & Tan Reference Yin, Fu and Tan2013; McCurdy et al. Reference McCurdy, Moore and Wang2019). Our model improves the prediction of our previous model by utilizing an open boundary condition at the interface when calculating the critical Rayleigh numbers in our theory. Additionally, we note a second kind of transition in this paper. This transition, which we refer to as a ‘dynamic shift’, is associated with the time evolution of the system where the conductive state transitions through a metastable shallow-convection state en route to its steady state of deep convection. Our numerical simulations shed light on this new flow configuration: shallow-to-deep convection.
The article is organized as follows. In § 2, we present the system of equations, interface conditions and the non-dimensionalized system. Then, the transition theory – one of the two main contributions of this work – is introduced in § 3 along with results showing its efficacy. Next, we detail our numerical scheme using a FEM in § 4, and we note how the treatment of interfacial and nonlinear terms allows us to write the system as a set of linear, sequentially decoupled equations. Results are shown in § 5 along with a discussion of how our three complementary approaches – stability analysis, fully nonlinear numerical simulations and the coarse-grained model – agree to provide a more complete picture of convection in fluid–porous systems. Our results showcase the second main contribution of this paper: the novel convection pattern of shallow-to-deep convection.
2. The coupled Navier–Stokes–Darcy system
In this section, we present the governing equations, detail the boundary and interface conditions and introduce the non-dimensional system.
2.1. Governing equations
In the free-flow zone, we use the same system as that studied by McCurdy et al. (Reference McCurdy, Moore and Wang2019) and Han et al. (Reference Han, Wang and Wang2020) – the incompressible Navier–Stokes equations with constant viscosity and the Boussinesq approximation, coupled with the advection–diffusion equations for heat:
\begin{equation} \left.\begin{gathered} \rho_0\left( \dfrac{\partial \boldsymbol{u_f} }{\partial t_f} + \left(\boldsymbol{u_f}\ \boldsymbol{\cdot}\ \boldsymbol{\nabla}\right) \boldsymbol{u_f} \right) = \boldsymbol{\nabla}\ \boldsymbol{\cdot}\ {{\boldsymbol{\mathsf{T}}}}\left(\boldsymbol{u_f} , p_f \right) - g\rho_0 \left[1-\beta\left(T_f-T_0 \right) \right]\boldsymbol{k}, \\ \boldsymbol{\nabla}\ \boldsymbol{\cdot}\ \boldsymbol{u_f} =0, \\ \dfrac{\partial T_f}{\partial t_f} + \boldsymbol{u_f}\ \boldsymbol{\cdot}\ \boldsymbol{\nabla} T_f = \frac{\kappa_f}{ \left(\rho_0c_p\right)_f} \boldsymbol{\nabla}^2T_f, \end{gathered}\right\} \end{equation}
where 
$\boldsymbol {u_f} = (u_f, v_f, w_f)$, 
$p_f$ and 
$T_f$ are the free-flow velocity, pressure and temperature, respectively, with 
$g$, 
$\rho _0$, 
$\beta$ and 
$T_0$ as acceleration due to gravity, the reference density of the fluid, the coefficient of thermal expansion and the temperature of the conductive state at the interface, respectively. The stress tensor and rate of strain tensor are defined as 
${{\boldsymbol{\mathsf{T}}}}(\boldsymbol {u_f} , p_f)=2\mu _0{{\boldsymbol{\mathsf{D}}}}(\boldsymbol {u_f} )-p_f{{\boldsymbol{\mathsf{I}}}}$ and 
${{\boldsymbol{\mathsf{D}}}}(\boldsymbol {u_f} )=\frac {1}{2}(\boldsymbol {\nabla }\boldsymbol {u_f} + \boldsymbol {\nabla } \boldsymbol {u_f} ^{\textrm {T}})$, respectively, with 
$\mu _0$ as dynamic viscosity and 
$\boldsymbol {k}$ as the upward-pointing unit normal. Additionally, 
$\kappa _f$, 
$c_p$ and 
$\lambda _f = \kappa _f/(\rho _0c_p)_f$ are the thermal conductivity of the fluid, specific heat capacity of the fluid and thermal diffusivity of the fluid, respectively.
For fluid flow in the porous medium, we assume the medium has a small porosity, as is generally applicable to geophysical systems (Bear Reference Bear1972; Nield & Bejan Reference Nield and Bejan2017). We therefore employ the Darcy–Boussinesq system with the advection–diffusion equation for heat:
\begin{equation} \left.\begin{gathered} \frac{\rho_0}{\chi}\dfrac{\partial \boldsymbol{u_m} }{\partial t_m} + \frac{\mu_0}{\varPi} \boldsymbol{u_m} ={-}\boldsymbol{\nabla} p_m - g\rho_0 \left[1-\beta\left(T_m-T_L \right) \right] \boldsymbol{k},\\ \boldsymbol{\nabla}\ \boldsymbol{\cdot}\ \boldsymbol{u_m} =0,\\ \frac{\left(\rho_0c_p\right)_m}{\left(\rho_0c_p\right)_f}\dfrac{\partial T_m}{\partial t_m} + \boldsymbol{u_m}\ \boldsymbol{\cdot}\ \boldsymbol{\nabla} T_m =\frac{\kappa_m}{ \left(\rho_0c_p\right)_f} \boldsymbol{\nabla}^2T_m, \end{gathered}\right\} \end{equation}
where 
$\boldsymbol {u_m} = (u_m, v_m, w_m),$ 
$p_m$ and 
$T_m$ are the velocity, pressure and temperature in the porous medium, respectively, 
$\chi$ and 
$\varPi$ are the porosity and permeability, 
$\lambda _m = \kappa _m/(\rho _0 c_p)_f$ is the thermal diffusivity of the medium and 
$T_L$ is the temperature at the lower boundary of the domain. We assume the medium to be homogeneous and isotropic so that the permeability 
$\varPi$ is constant and scalar-valued. The thermal conductivity 
$\kappa _m$ and specific heat capacity 
$(\rho _0 c_p)_m$ of the porous medium are defined as averages of the fluid and solid components. While many studies neglect the time derivative 
$\partial _t \boldsymbol {u_m}$ in the first equation of (2.2), retaining this term is useful for energy analysis of the system (see McCurdy et al. Reference McCurdy, Moore and Wang2019).
2.2. Boundary and interface conditions
The domain, shown in figure 1, consists of flat, horizontal, non-penetrable plates at the top and bottom with a non-deforming interface between the two regions, 
$\varOmega _f=\{(x,y,z)\in \mathbb {R}^2 \times (0,d_f)\}$ for the free flow and 
$\varOmega _m=\{(x,y,z)\in \mathbb {R}^2 \times (-d_m,0)\}$ for the porous medium. The temperature is held constant at the top and bottom plates. The flow satisfies a free-slip condition at the top and an impermeable condition at the bottom:
\begin{equation}\left.\begin{gathered} T_f = T_U,\quad \boldsymbol{u_f}\boldsymbol{\cdot}\ \boldsymbol{n} =\frac{{\partial} \boldsymbol{u_f}{_\tau}}{{\partial} \boldsymbol{n}}=0\quad\mbox{on}\ z=d_f,\\ T_m = T_L, \quad \boldsymbol{u_m}\boldsymbol{\cdot}\ \boldsymbol{n} =0\quad \mbox{on}\ z={-}d_m, \end{gathered}\right\}\end{equation}
where 
${\boldsymbol {u_f} }_\tau =(u_f, v_f)$ denotes the tangential (horizontal) components of the velocity at the top of the domain with 
$\boldsymbol {n}$ as the unit normal vector.
At the interface 
$\varGamma _i$ 
$(z=0)$, we require continuity of temperature, heat flux and the normal component of velocity:
\begin{equation} \left.\begin{gathered} T_f = T_m, \\ \kappa_f \boldsymbol{\nabla} T_f\ \boldsymbol{\cdot}\ \boldsymbol{n} = \kappa_m \boldsymbol{\nabla} T_m\ \boldsymbol{\cdot}\ \boldsymbol{n},\\ \boldsymbol{u_f}\ \boldsymbol{\cdot}\ \boldsymbol{n} = \boldsymbol{u_m}\ \boldsymbol{\cdot}\ \boldsymbol{n}. \end{gathered}\right\} \end{equation}For the last two conditions, we use the Beavers–Joseph–Saffman–Jones (BJSJ) condition (Saffman Reference Saffman1971) and the Lions interface condition to specify the tangential and normal stresses, respectively. The BJSJ condition, also known as the Navier-slip condition, is
\begin{equation} -\boldsymbol{\tau}\ \boldsymbol{\cdot}\ {{\boldsymbol{\mathsf{T}}}} \left(\boldsymbol{u_f} ,p_f \right)\boldsymbol{n} = \frac{\mu_0 \alpha}{\sqrt{\varPi}}\boldsymbol{\tau}\ \boldsymbol{\cdot}\ \boldsymbol{u_f} , \end{equation}
where 
$\alpha$ is an empirically determined coefficient and 
$\boldsymbol {\tau }$ denotes the unit tangent vectors. Lastly, the Lions interface condition is
\begin{equation} -\boldsymbol{n}\ \boldsymbol{\cdot}\ {{\boldsymbol{\mathsf{T}}}}\left(\boldsymbol{u_f} ,p_f \right)\boldsymbol{n} + \frac{\rho_0}{2}\left|\boldsymbol{u_f} \right|^2 = p_m. \end{equation}
The inclusion of the 
$({\rho _0}/{2})|\boldsymbol {u_f} |^2$ term in this interface condition is essential in conducting the nonlinear stability analysis (Chidyagwai & Rivière Reference Chidyagwai and Rivière2009; Çeşmelioğlu & Rivière Reference Çeşmelioğlu and Rivière2008, Reference Çeşmelioğlu and Rivière2009; Discacciati & Quarteroni Reference Discacciati and Quarteroni2009; Girault & Rivière Reference Girault and Rivière2009; McCurdy et al. Reference McCurdy, Moore and Wang2019).
2.3. System of perturbed variables
Instead of working with the physical variables, we consider the perturbed variables; that is, we consider the deviation of the velocity, temperature and pressure profiles from their conductive steady states. The conductive state, denoted with an overhead bar, is a stationary fluid and a piecewise-linear temperature:
\begin{equation} \left.\begin{gathered} \bar{\boldsymbol{u}}_{\boldsymbol{f}} = \bar{\boldsymbol{u}}_{\boldsymbol{m}} = 0, \\ \bar{T}_f = T_0 + z\frac{T_U - T_0 }{d_f},\\ \bar{T}_m = T_0 + z\frac{T_0 - T_L }{d_m}. \end{gathered}\right\} \end{equation}
Here, 
$T_0$ represents the interface temperature of the conductive solution
\begin{equation} T_0 = \frac{\kappa_m d_f T_L + \kappa_f d_m T_U}{\kappa_m d_f + \kappa_f d_m}. \end{equation}
If 
$T_U>T_L$, the conductive state is stable, but if 
$T_L>T_U$, buoyancy can destabilize the system. Throughout this work, we only consider the case where 
$T_L>T_U$. Additionally, we choose 
$\bar {p}_f$ and 
$\bar {p}_m$ to satisfy
\begin{equation} \left.\begin{gathered} \boldsymbol{\nabla} \bar{p}_f ={-}g\rho_0\left(1-\beta\left(\bar{T}_f-T_0\right)\right)\boldsymbol{k}, \\ \boldsymbol{\nabla} \bar{p}_m ={-}g\rho_0\left(1-\beta\left(\bar{T}_m-T_L\right)\right)\boldsymbol{k}. \end{gathered}\right\} \end{equation}
With the perturbation variables 
$\boldsymbol {v}_j$, 
$\theta _j$ and 
${\rm \pi} _j$ for 
$j=\{f,m\}$ regions, we perturb the steady-state solutions:
\begin{equation} \left.\begin{gathered} \boldsymbol{u_f} =\bar{\boldsymbol{u}}_{\boldsymbol{f}} + \boldsymbol{v_f} ,\quad \boldsymbol{u_m} = \bar{\boldsymbol{u}}_{\boldsymbol{m}} +\boldsymbol{v_m} ,\\ T_f = \bar{T}_f +\theta_f,\quad T_m = \bar{T}_m+\theta_m, \\ p_f = \bar{p}_f +{\rm \pi}_f,\quad p_m= \bar{p}_m + {\rm \pi}_m. \end{gathered}\right\} \end{equation}Substituting (2.10) into the original system produces a system for the perturbed variables:
\begin{equation} \left.\begin{gathered} \rho_0\left(\dfrac{\partial \boldsymbol{v_f} }{\partial t_f} + \left(\boldsymbol{v_f}\ \boldsymbol{\cdot}\ \boldsymbol{\nabla}\right) \boldsymbol{v_f} \right) = \boldsymbol{\nabla}\ \boldsymbol{\cdot}\ {{\boldsymbol{\mathsf{T}}}}\left(\boldsymbol{v_f} , {\rm \pi}_f \right) + g\rho_0 \beta\theta_f\boldsymbol{k}, \\ \boldsymbol{\nabla}\ \boldsymbol{\cdot}\ \boldsymbol{v_f} =0, \\ \dfrac{\partial \theta_f}{\partial t_f} + \boldsymbol{v_f}\ \boldsymbol{\cdot}\ \boldsymbol{\nabla} \theta_f =\lambda_f \boldsymbol{\nabla}^2\theta_f -\boldsymbol{v_f}\ \boldsymbol{\cdot}\ \boldsymbol{k} \left(\frac{T_U-T_0}{d_f}\right) \end{gathered}\right\} \end{equation}
for 
$\varOmega _f$;
\begin{equation} \left.\begin{gathered} \frac{\rho_0}{\chi}\dfrac{\partial \boldsymbol{v_m} }{\partial t_m} + \frac{\mu_0}{\varPi} \boldsymbol{v_m} ={-}\boldsymbol{\nabla} {\rm \pi}_m + g\rho_0 \beta\theta_m\boldsymbol{k},\\ \boldsymbol{\nabla}\ \boldsymbol{\cdot}\ \boldsymbol{v_m} =0\\ \frac{\left(\rho_0c_p\right)_m}{\left(\rho_0c_p\right)_f}\dfrac{\partial \theta_m}{\partial t_m} + \boldsymbol{v_m}\ \boldsymbol{\cdot}\ \boldsymbol{\nabla} \theta_m =\lambda_m \boldsymbol{\nabla}^2\theta_m -\boldsymbol{v_m}\ \boldsymbol{\cdot}\ \boldsymbol{k} \left(\frac{T_0-T_L}{d_m}\right) \end{gathered}\right\} \end{equation}
for 
$\varOmega _m$;
\begin{equation} \left.\begin{gathered} \theta_f = \theta_m, \\ \kappa_f \boldsymbol{\nabla} \theta_f\ \boldsymbol{\cdot}\ \boldsymbol{n} = \kappa_m \boldsymbol{\nabla} \theta_m\ \boldsymbol{\cdot}\ \boldsymbol{n},\\ \boldsymbol{v_f}\ \boldsymbol{\cdot}\ \boldsymbol{n} = \boldsymbol{v_m}\ \boldsymbol{\cdot}\ \boldsymbol{n}, \\ -\boldsymbol{\tau}\ \boldsymbol{\cdot}\ {{\boldsymbol{\mathsf{T}}}} \left(\boldsymbol{v_f} ,{\rm \pi}_f \right)\boldsymbol{n} = \frac{\mu_0 \alpha}{\sqrt{\varPi}}\boldsymbol{\tau}\ \boldsymbol{\cdot}\ \boldsymbol{v_f} ,\\ -\boldsymbol{n}\ \boldsymbol{\cdot}\ {{\boldsymbol{\mathsf{T}}}}\left(\boldsymbol{v_f} ,{\rm \pi}_f \right)\boldsymbol{n} + \frac{\rho_0}{2}\left|\boldsymbol{v_f} \right|^2 = {\rm \pi}_m \end{gathered}\right\} \end{equation}
for the interface conditions on 
$\varGamma _i$; and
\begin{equation} \left.\begin{gathered} \theta_f = 0,\quad \boldsymbol{v_f}\ \boldsymbol{\cdot}\ \boldsymbol{n}= \frac{{\partial} \boldsymbol{v_f}{ _\tau}}{{\partial} \boldsymbol{n}}=0 \quad \mbox{on}\ z=d_f,\\ \theta_m = 0, \quad \boldsymbol{v_m}\ \boldsymbol{\cdot}\ \boldsymbol{n} =0\quad \mbox{on}\ z={-}d_m \end{gathered}\right\} \end{equation}for the boundary conditions.
2.4. Non-dimensional system
As in previous work, we non-dimensionalize the system using the porous values as a reference (Chen & Chen Reference Chen and Chen1988; Straughan Reference Straughan2002; McCurdy et al. Reference McCurdy, Moore and Wang2019), with non-dimensional variables denoted by tildes:
\begin{equation} \boldsymbol{v}_{\boldsymbol{j}} = \tilde{\boldsymbol{v}}_{\boldsymbol{j}} \frac{\nu}{d_m},\quad \boldsymbol{x}_{\boldsymbol{j}} = \tilde{\boldsymbol{x}}_{\boldsymbol{j}} d_m,\quad t_j = \tilde{t}\frac{d_m^2}{\lambda_m},\quad \theta_j = \tilde{\theta}_j \frac{\left(T_0-T_L\right)\nu}{\lambda_m},\quad {\rm \pi}_j = \tilde{\rm \pi}_j \frac{\rho_0 \nu^2}{d_m^2}, \end{equation}
for 
$j=\{f,m\}$, where 
$\nu = \mu _0/\rho _0$ is the kinematic viscosity. We also note that the stress tensor has been altered slightly from when it was first introduced; the non-dimensional stress tensor is defined as 
$\tilde {{{\boldsymbol{\mathsf{T}}}}}(\boldsymbol {v_f} , {\rm \pi}_f)=2{{\boldsymbol{\mathsf{D}}}}(\boldsymbol {v_f} )-{\rm \pi} _f{{\boldsymbol{\mathsf{I}}}}$.
Substituting the non-dimensional variables into (2.11)–(2.14) results in the governing equations (after dropping the tildes) for 
$\varOmega _f = \{(x,y,z) \in \mathbb {R}^2\times (0,\hat {d})\}$, where 
$\hat {d}$ is the ratio of the free-flow depth to that of the porous medium:
\begin{equation} \left.\begin{gathered} \frac{1}{Pr_{m}}\dfrac{\partial \boldsymbol{v_f} }{\partial t} + \left(\boldsymbol{v_f}\ \boldsymbol{\cdot}\ \boldsymbol{\nabla} \right)\boldsymbol{v_f} = \boldsymbol{\nabla}\ \boldsymbol{\cdot}\ \tilde{{{\boldsymbol{\mathsf{T}}}}}\left(\boldsymbol{v_f} , {\rm \pi}_f \right) + \frac{Ra_m}{Da}\theta_f \boldsymbol{k}, \\ \boldsymbol{\nabla}\ \boldsymbol{\cdot}\ \boldsymbol{v_f} =0, \\ \dfrac{\partial \theta_f}{\partial t} + Pr_{m} \boldsymbol{v_f}\ \boldsymbol{\cdot}\ \boldsymbol{\nabla} \theta_f = \epsilon_T \boldsymbol{\nabla}^2 \theta_f + \frac{1}{\epsilon_T}\boldsymbol{v_f}\ \boldsymbol{\cdot}\ \boldsymbol{k}, \end{gathered}\right\} \end{equation}
for 
$\varOmega _m = \{(x,y,z) \in \mathbb {R}^2\times (-1,0)\}$:
\begin{equation} \left.\begin{gathered} \frac{Da}{Pr_{m} \chi }\dfrac{\partial \boldsymbol{v_m} }{\partial t} + \boldsymbol{v_m} ={-}Da\boldsymbol{\nabla} {\rm \pi}_m +Ra_m\theta_m \boldsymbol{k},\\ \boldsymbol{\nabla}\ \boldsymbol{\cdot}\ \boldsymbol{v_m} =0,\\ \varrho \dfrac{\partial \theta_m}{\partial t} + Pr_{m}\boldsymbol{v_m}\ \boldsymbol{\cdot}\ \boldsymbol{\nabla} \theta_m = \boldsymbol{\nabla}^2 \theta_m + \boldsymbol{v_m}\ \boldsymbol{\cdot}\ \boldsymbol{k}, \end{gathered}\right\} \end{equation}
and at 
$\varGamma _i = \{(x,y,z) \in \mathbb {R}^2\times (z=0)\}$:
$$\begin{gather} \theta_f = \theta_m, \end{gather}$$
$$\begin{gather}\epsilon_T \boldsymbol{\nabla} \theta_f\ \boldsymbol{\cdot}\ \boldsymbol{n} =\boldsymbol{\nabla} \theta_m\ \boldsymbol{\cdot}\ \boldsymbol{n}, \end{gather}$$
$$\begin{gather}\boldsymbol{v_f}\ \boldsymbol{\cdot}\ \boldsymbol{n} = \boldsymbol{v_m}\ \boldsymbol{\cdot}\ \boldsymbol{n} \end{gather}$$
$$\begin{gather}-\boldsymbol{\tau}\ \boldsymbol{\cdot}\ \tilde{{{\boldsymbol{\mathsf{T}}}}} \left(\boldsymbol{v_f} ,{\rm \pi}_f \right)\boldsymbol{n} = \frac{\alpha}{\sqrt{Da}}\left(\boldsymbol{\tau}\ \boldsymbol{\cdot}\ \boldsymbol{v_f} \right), \end{gather}$$
$$\begin{gather}-\boldsymbol{n}\ \boldsymbol{\cdot}\ \tilde{{{\boldsymbol{\mathsf{T}}}}} \left(\boldsymbol{v_f} ,{\rm \pi}_f \right)\boldsymbol{n} + \frac{1}{2}|\boldsymbol{v_f} |^2 ={\rm \pi}_m. \end{gather}$$Boundary conditions at the top and bottom of the domain are
\begin{equation} \left.\begin{gathered} \theta_f = 0,\quad \boldsymbol{v_f}\ \boldsymbol{\cdot}\ \boldsymbol{n} = \frac{{\partial} \boldsymbol{v_f}{ _\tau}}{{\partial} \boldsymbol{n}}=0 \quad \mbox{on}\ z=\hat{d},\\ \theta_m = 0,\quad \boldsymbol{v_m}\ \boldsymbol{\cdot}\ \boldsymbol{n} =0\quad \mbox{on}\ z={-}1. \end{gathered}\right\} \end{equation}
The non-dimensional numbers 
$\hat {d}$, 
$Pr_{m}$, 
$Da$, 
$\epsilon _T$ and 
$\varrho$ are defined by
\begin{equation} \hat{d} = \frac{d_f}{d_m},\quad Pr_{m} = \frac{\nu}{\lambda_m},\quad Da = \frac{\varPi}{d_m^2},\quad \epsilon_T = \frac{\lambda_f}{\lambda_m},\quad \varrho = \frac{\left(\rho_0c_p\right)_m}{\left(\rho_0c_p\right)_f}, \end{equation}and the Rayleigh numbers of both regions are
\begin{equation} Ra_m = \frac{g\beta \left(T_L-T_0 \right) Da d_m^3}{\nu \lambda_m},\quad Ra_{f} = \frac{\hat{d}^4}{Da\epsilon_T^2}Ra_m. \end{equation}While the Rayleigh number of the porous medium is used throughout this work, the corresponding Rayleigh number of the free-flow region can be determined via the relationship above.
The depth ratio of the two layers 
$\hat {d}$ plays a central role in our study. Other dimensionless parameters, like 
$Pr_{m},\epsilon _T,\varrho$, represent values inherent to the fluid and/or porous medium. In industrial applications, these cannot easily be altered, or it may be impractical to do so. However, the depth ratio can more readily be changed; for example, by adding or removing fluid from the free-flow layer, the depth ratio varies. The coarse-grained model presented in the next section details how changing the depth ratio can significantly affect convection profiles.
The second parameter of interest for geophysical applications is the Darcy number 
$Da$ – especially the small-Darcy regime since many materials have small permeability values. Relatively impervious materials, like limestone or granite, can have permeabilities between 
$10^{-18}$ and 
$10^{-15}\,\textrm {m}^2$, while more porous media, like sorted gravel or sand, can have permeabilities in the range 
$10^{-10}\lesssim \varPi \lesssim 10^{-7}\,\textrm {m}^2$ (Bear Reference Bear1972; Nield & Bejan Reference Nield and Bejan2017). The resulting Darcy numbers in experiments and/or simulations are small (typically 
$10^{-10}\leq Da \leq 10^{-5}$), highlighting the necessity to study the relevant small-Darcy regime in more depth, as done in Lyu & Wang (Reference Lyu and Wang2021).
3. Coarse-grained model for the transition between deep and shallow convection
McCurdy et al. (Reference McCurdy, Moore and Wang2019) conducted linear and nonlinear stability analyses of the Navier–Stokes–Darcy–Boussinesq system to determine the threshold Rayleigh number needed for the conductive state to become unstable at a given wavenumber, 
$a_m$. For a fixed depth ratio, figure 2(c) shows the collection of these points in the 
$a_m\unicode{x2013}a_m$ plane, called the marginal stability curve, which delineates regions of stability from instability. The critical Rayleigh number 
$Ra_m^*$ is the global minimum of this curve, representing the lowest 
$Ra_m$ value that produces a flow instability. The corresponding critical wavenumber 
$a_m^*$ dictates the flow configuration at the onset of convection: small 
$a_m^*$ (i.e. large wavelength) indicates deep convection cells that occupy the entire coupled domain as seen in figure 2(a), while large 
$a_m^*$ (i.e. small wavelength) indicates narrow convection cells that occupy the fluid region only, as seen in figure 2(b). The bimodal nature of the marginal stability curve can create an abrupt shift between these two flow configurations as the depth ratio changes.
The marginal stability curves shown in figure 2(c) are calculated for a range of 
$\hat {d}$ values using the linear stability analysis outlined by McCurdy et al. (Reference McCurdy, Moore and Wang2019). The global minimum of each curve corresponds to 
$Ra_m^*$ and is indicated by a red dot. The red arrow illustrates the path that connects sequential values of 
$Ra_m^*$ as 
$\hat {d}$ increases from 0.15 to 0.22. As seen in the figure, a jump in 
$a_m^*$ occurs as 
$\hat {d}$ changes from 0.18 to 0.19, indicating an abrupt transition from deep convection to shallow convection.
This observation spurred McCurdy et al. (Reference McCurdy, Moore and Wang2019) to develop a simple model to predict the critical depth ratio 
$\hat {d}^*$ that distinguishes shallow convection from deep convection. A brief outline of this model is as follows. First, due to the dimensionless definitions, the depth ratio can be expressed as the following combination of the other parameters:
\begin{equation} \hat{d} = \left[\frac{Ra_f Da \epsilon_T^2}{Ra_m}\right]^{1/4}. \end{equation}
Of particular importance is the appearance of the two Raleigh numbers, 
$Ra_{f}$ and 
$Ra_m$. If the Rayleigh number exceeds its critical value in the fluid region, 
$Ra_{f} > Ra_{f}^*$, but remains below the critical value in the porous medium, 
$Ra_m < Ra_m^*$, then only shallow convection cells would be expected to form. On the other hand, if the Rayleigh number exceeds the critical value in both regions, then convection cells can penetrate deeper into the porous domain. Hence, the transition between shallow and deep convection should occur when both Rayleigh numbers, 
$Ra_{f}$ and 
$Ra_m$, become critical simultaneously. In reality, the critical values, 
$Ra_{f}^*$ and 
$Ra_m^*$, are coupled to one another through the flow details at the interface. However, McCurdy et al. (Reference McCurdy, Moore and Wang2019) showed that useful estimates could be obtained by treating the fluid and porous regions as uncoupled for the sake of calculating 
$Ra_{f}^*$ and 
$Ra_m^*$, and then inputting these values into (3.1) to estimate the critical depth ratio 
$\hat {d}^*$ that defines the transition between shallow and deep convection:
\begin{equation} \hat{d}^* = \left[\frac{Ra_{f}^*}{Ra_m^*}Da\epsilon_T^2\right]^{1/4}. \end{equation}
In particular, McCurdy et al. (Reference McCurdy, Moore and Wang2019) assumed no-slip and no-penetration conditions along the boundaries of the fluid domain and no-penetration conditions along the boundaries of the porous medium. These assumptions yield 
$Ra_{f}^*=1707.8$ and 
$Ra_m^* = 4{\rm \pi} ^2\approx 39.5$, which gives
\begin{equation} \hat{d}^* = \left[\frac{1707.8}{39.5}Da\epsilon_T^2\right]^{1/4}. \end{equation}This formula was shown to predict the actual critical depth ratio to within a relative error of 13 %–17 % for the cases tested by McCurdy et al. (Reference McCurdy, Moore and Wang2019). This level of accuracy, while not exceptional, is promising considering that the formula neglects any kind of coupling between the two regions.
Here, we further refine this model by introducing asymptotically weak coupling at the interface (see e.g. Moore & Shelley Reference Moore and Shelley2012). As demonstrated below, the new model yields improved estimates for 
$\hat {d}^*$, with relative errors of the order of 0 %–4 % for sufficiently small 
$Da$. We first asymptotically expand both velocity fields, i.e. inside the fluid and the medium regions, in small 
$Da:=\varepsilon ^2$:
\begin{equation} \boldsymbol{v_j} ^{\varepsilon} = \boldsymbol{v_j} ^{(0)} + \varepsilon \boldsymbol{v_j} ^{(1)} + \varepsilon^2 \boldsymbol{v_j} ^{(2)} +\cdots,\quad \textrm{for }j\in\{f,m\}. \end{equation}
If 
$Da=0$, the Darcy system (2.17) is degenerate, giving a porous-medium velocity that vanishes throughout the domain, 
$\boldsymbol {v_m} ^{(0)} \equiv 0$. Therefore, interface condition (2.18c) gives 
$\boldsymbol {v_f} ^{(0)}\ \boldsymbol {{\cdot }}\ \boldsymbol {n}=0$ at leading order. Meanwhile, the leading-order equation of the BJSJ condition (2.18d) gives 
$\boldsymbol {v_f} ^{(0)}\ \boldsymbol {{\cdot }}\ \boldsymbol {\tau }=0$. Thus, both components of the fluid velocity 
$\boldsymbol {v_f} ^{(0)}$ vanish at the interface to leading order in 
$\varepsilon$, and so the no-slip and no-penetration interface conditions assumed by McCurdy et al. (Reference McCurdy, Moore and Wang2019) are valid to leading order in the fluid domain. We therefore continue to use the corresponding value 
$Ra_{f}^*=1707.8$ in the new model.
Approaching the interface from the porous-medium side, though, the condition for an impenetrable boundary, 
$\boldsymbol {v_m}\ \boldsymbol {{\cdot }}\ \boldsymbol {n}=0$, is not recovered as the leading-order non-trivial dynamics in small Darcy number. In the improved model, we instead view the top of the porous domain as an open boundary, along which the pressure is uniform (Nield & Bejan Reference Nield and Bejan2017). Due to Darcy's law (2.17), the condition of constant pressure along a horizontal interface implies that the tangential velocity vanishes 
$\boldsymbol {v_m} \times \boldsymbol {n}=0$, which produces a critical value of 
$Ra_m^*= 27.1$ at the porous wavenumber of 
$a_{m,1}^*=2.3$ for the uncoupled porous medium (Tyvand Reference Tyvand2002; Nield & Bejan Reference Nield and Bejan2017). Although we have been so far unable to rigorously justify the use of this condition from first principles, we observe that, in practice, it yields significantly improved estimates for 
$Ra_m^*$ and 
$a_{m,1}^*$ as 
$Da \to 0$. Table 1 demonstrates this idea by showing the true critical values 
$Ra_m^*$ for the coupled system and their critical wavenumbers as calculated numerically by the linear stability analysis, for a range of Darcy numbers and three values of 
$\epsilon _T$. The table shows that as 
$Da \to 0$, 
$Ra_m^*$ is well approximated by 27.1 in each case, with relative errors of the order of 1 %–2 %.
Table 1. Critical Rayleigh numbers 
$Ra_{m,c}$ and their critical wavenumbers 
$a_{m,1}^*$ at respective 
$\hat {d}^*$ values with 
$\epsilon _T=\{0.5, 0.7,1.0\}$ and Darcy numbers 
$Da\to 0$.
With these new critical Rayleigh numbers, we obtain the more accurate coarse-grained model for predicting the critical depth ratio:
\begin{equation} \hat{d}^* = \left[\frac{1707.8}{27.1}Da \epsilon_T^2 \right]^{1/4}. \end{equation}As we show below, the results produced with this model become increasingly accurate in the small-Darcy limit, which is reasonable given that our intuition in choosing boundary conditions came from the small-Darcy limit.
Figure 3 shows data for three different 
$\epsilon _T$ values while varying 
$Da$, since our formula is a function of these two variables. As 
$Da\to 0,$ the critical depth ratio 
$\hat {d}^*$ goes to zero as well. We plot the critical depth ratios found with two methods – the first predicted from the heuristic theory, compared with the second obtained from the marginal stability results from McCurdy et al. (Reference McCurdy, Moore and Wang2019). To show that our model predicts 
$\hat {d}^*$ going to zero at the same rate as the actual values found using the linear stability analysis 
$\hat {d}^*_{LSA}$, we plot the data 
$Da$ versus 
$\hat {d}^*$ with a log–log plot along with a reference triangle to illustrate the slope of 
$1/4$.
Figure 3. Critical depth ratios for various 
$\epsilon _T$ (ratio of thermal diffusivities) values, 
$\epsilon _T=0.5, 0.7, 1.0$. The solid lines represent the predicted 
$\hat {d}^*$ values from our theory (3.5), and the circles are the 
$\hat {d}^*$ values calculated from the marginal stability curves determined by McCurdy et al. (Reference McCurdy, Moore and Wang2019).
To quantify the error of our predicted 
$\hat {d}^*$ values, we define the relative error between the predicted 
$\hat {d}^*$ values from (3.5), 
$\hat {d}^*_{Thry}$, and values found using the linear stability analysis, 
$\hat {d}^*_{LSA}$, with
\begin{equation} e_{rel} =\frac{\left|\hat{d}^*_{LSA}-\hat{d}^*_{Thry} \right|}{\hat{d}^*_{LSA}}. \end{equation}
The relative errors are noted in table 2. For 
$Da=1.0\times 10^{-4}$, the model has about 10 % relative error, while the relative error drops to less than 2 % when 
$Da=1.0\times 10^{-8}$. The worst relative error with our new formula outperforms the best relative error of the formula presented in McCurdy et al. (Reference McCurdy, Moore and Wang2019).
Table 2. Relative errors – calculated with (3.6) – between predicted 
$\hat {d}^*$ values and values found using the linear stability analysis with 
$\epsilon _T=\{0.5, 0.7,1.0\}$ and Darcy numbers 
$Da\to 0$.
This coarse-grained model allows us to accurately predict the depth ratio that triggers a shift between deep and shallow convection. In applications with technology, being able to harness convection in this manner could have considerable impact with the design of heat sinks. By selecting physical properties of the porous medium, a suitable depth ratio could be chosen to obtain the desired flow configuration to prevent overheating and enhance circulation. Rather than conduct costly experiments to determine a depth ratio that yields the intended results, our formula provides a reference for an appropriate range of 
$\hat {d}$ values needed for deep or shallow convection.
This bifurcation from shallow to deep convection is highlighted in several works, using linear stability or numerical simulations to show the deep versus shallow convection profiles. Finding the critical depth ratio with either of these methods – linear stability and numerical simulations – can be computationally expensive and time-consuming. Our model allows us to accurately predict a range of depth ratios where this shift in convection occurs, each of which are in good agreement with the critical depth ratios found by previous researchers. For example, with the same governing equations we consider, Chen & Chen (Reference Chen and Chen1988) showed via linear stability that a transition occurred between depth ratios of 
$\hat {d} = 0.12$ and 
$\hat {d}=0.13$ using the fixed parameters 
$\sqrt {Da}=0.003, \epsilon _T=0.7$. Substituting these values into (3.5) allows our model to predict that a transition occurs at a depth ratio of
\begin{equation} \hat{d}^* = \left[\frac{1707.8}{27.1}\left(0.003\right)^2 \left(0.7\right)^2 \right]^{1/4}\approx 0.1291, \end{equation}
agreeing well with the work conducted by Chen & Chen (Reference Chen and Chen1988). The recent work by Han et al. (Reference Han, Wang and Wang2020) presented an example where the convection profiles shifted from shallow to deep convection as the depth ratio was altered from 
$\hat {d}=0.16$ to 
$\hat {d}=0.17$, found via numerics and their centre manifold theory (see figure 5.2 of Han et al. (Reference Han, Wang and Wang2020)). With the values of 
$Da=25\times 10^{-6}$ and 
$\epsilon _T=0.7$ from this example in their work, we predict a critical depth ratio of 
$\hat {d}^* \approx 0.1667$, which falls in the range determined by Han et al. The formula we developed greatly reduces the computational time needed to determine where the shallow to deep transition occurs. Using our coarse-grained model bypasses the necessity of searching through large parameter spaces while solving stability problems or conducting numerical simulations, as it quickly narrows the parameter regime where the transition occurs, especially for small Darcy numbers. In the Appendix, we present data showing good agreement between true critical depth ratios and those predicted from our model. Additionally, we show that our formula (3.5) can be used to predict how altering the ratio of thermal diffusivities 
$\epsilon _T$ or the Darcy number 
$Da$ can trigger a shift in convection.
Although each of the papers noted above use two-dimensional simulations or experiments in their work, our model can also be used in three-dimensional settings as well. Since our theory comes from the stability of the uncoupled systems under the assumption of infinite horizontal plates, both the two- and three-dimensional problems reduce to one spatial dimension – the vertical direction. As such, our model can be applied to both two- and three-dimensional systems, albeit with some error since the assumption of infinite horizontal plates cannot be satisfied in physical settings or numerical simulations.
The main limitation of our coarse-grained model is that the results are confined to the onset of convection, where Raleigh numbers are close to their critical counterparts. To develop a more comprehensive understanding of this phenomenon over a more broad parameter regime, we turn to numerical simulations of the full, nonlinear system.
4. Numerical scheme
The stability analyses and coarse-grained model can predict flow configurations near the onset of convection, but numerical simulations are needed to examine dynamics outside of this regime. In this section, we present a numerical scheme to simulate convection in the superposed fluid–porous medium system using a FEM. The main advantage to this scheme is that, by lagging the nonlinear terms and the terms associated with interface conditions, we produce a set of linear, sequentially decoupled equations. A similar scheme was studied by Chen et al. (Reference Chen, Han, Wang and Zhang2020) with the Cahn–Hilliard–Navier–Stokes–Darcy–Boussinesq system – or the same system we study with the inclusion of a phase-field function using the Cahn–Hilliard equations. Using a numerical method related to that of Chen et al. (Reference Chen, Han, Wang and Zhang2020) is beneficial since the method is unconditionally long-time stable. We outline the scheme below, and note the time-lagged terms as they are shown. Next, we detail how the variational forms are obtained, using the following notation for vector-valued functions 
$\boldsymbol {f}$ and 
$\boldsymbol {g}$ and matrix-valued functions 
${{\boldsymbol{\mathsf{A}}}}$ and 
${{\boldsymbol{\mathsf{B}}}}$:
\begin{equation} \left(\boldsymbol{f} , \boldsymbol{g} \right)_j = \int_{\varOmega_j}\boldsymbol{f}\ \boldsymbol{\cdot}\ \boldsymbol{g} \,{\rm d}\varOmega_j,\quad \langle{{\boldsymbol{\mathsf{A}}}}, {{\boldsymbol{\mathsf{B}}}} \rangle_j = \int_{\varOmega_j}{{\boldsymbol{\mathsf{A}}}}\,\colon{{\boldsymbol{\mathsf{B}}}}\,{\rm d}\varOmega_j,\quad \|\boldsymbol{f} \|_j^2=(\boldsymbol{f} ,\boldsymbol{f} )_j,\quad |\boldsymbol{f} |^2 = \boldsymbol{f}\ \boldsymbol{\cdot}\ \boldsymbol{f} , \end{equation}
for domains 
$j \in \{f,m\}$ where the 
$f$ and 
$m$ subscripts correspond to the fluid and porous medium regions, respectively.
We introduce the following finite element spaces:
(i)
$W_f = \{\boldsymbol {w}_{\boldsymbol {f}} \in [H^1(\varOmega _f )]^2: \boldsymbol {w}_{\boldsymbol {f}}\ \boldsymbol {{\cdot }}\ \boldsymbol {n} =0 \textrm { at top }+\textrm{ periodic on left and right} \}$,(ii)
$Q_f = \{ q_f \in L^2(\varOmega _f ): \int _{\varOmega _f} q_f\,\textrm {d}\kern0.06em \boldsymbol {x}=0 \}=L^2_0( \varOmega _f)$,(iii)
$Q_m = \{ q_m \in L^2(\varOmega _m ): \int _{\varOmega _m} q_m\,\textrm {d}\kern0.06em \boldsymbol {x}=0 \}=L^2_0( \varOmega _m)$,(iv)
$\varPsi = \{ \psi \in H^1(\varOmega ): \psi =0 \textrm { at top and bottom }+\textrm{ periodic on left and right} \}$.
The space 
$\varPsi$ spans the entire domain and is reserved for the temperature field over 
$\varOmega$, since the advection–diffusion equation for heat can be written as one equation over the entirety of the domain. Instead of solving two problems for 
$\theta _f$ and 
$\theta _m,$ we solve only one for
\begin{equation}
\theta = \left\{\begin{array}{@{}ll@{}} \theta_f & \textrm{for }
\boldsymbol{x}\in\varOmega_f, \\
\theta_m & \textrm{for }\boldsymbol{x}\in\varOmega_m. \end{array}\right.
\end{equation}With superscripts denoting the time iteration, we present the variational problems corresponding to the system (2.16)–(2.19).
First, given 
$(\boldsymbol {v_f} ^{(n)},\boldsymbol {v_f} ^{(n-1)}, \theta _m^{(n)})$, we find the perturbed pressure in the porous medium 
${\rm \pi} _m^{(n+1)}\in Q_m$ with
\begin{align} &
Da\left(\boldsymbol{\nabla} {\rm \pi}_m^{(n+1)},
\boldsymbol{\nabla} q_m\right)_m -Ra_m
\left(\theta_m^{(n)}\boldsymbol{k}, \boldsymbol{\nabla} q_m\right)_m \nonumber\\
&\quad +\int_{\varGamma_i}\left[\frac{Da}{Pr_{m} \chi }
\frac{\boldsymbol{v_f} ^{(n)}-\boldsymbol{v_f}
^{(n-1)}}{{\rm \Delta} t} + \boldsymbol{v_f} ^{(n)}\right]\
\boldsymbol{\cdot}\ \boldsymbol{n} q_m \,{\rm
d}\varGamma_i=0\quad \forall\ q_m\in Q_m.
\end{align}
To begin decoupling the problems, we use the previous 
$\boldsymbol {v_f} ^{(n)}$ values in the integral along the interface. Our treatment of the interface term, along with the method in which we solve for the Darcy velocity, differs from the method used in Chen et al. (Reference Chen, Han, Wang and Zhang2020). We still observe experimentally that our method is stable, with the time-step restriction of 
${\rm \Delta} t \sim O(Da)$ or smaller. Having solved for 
${\rm \pi} _m^{(n+1)}$, we can use the previous Darcy velocity 
$\boldsymbol {v_m} ^{(n)}$ to find the updated velocity 
$\boldsymbol {v_m} ^{(n+1)}$. The Darcy equation is
\begin{equation} \frac{Da}{Pr_{m} \chi}\frac{\boldsymbol{v_m} ^{(n+1)} - \boldsymbol{v_m} ^{(n)} }{{\rm \Delta} t} + \boldsymbol{v_m} ^{(n+1)} ={-}Da\boldsymbol{\nabla} {\rm \pi}_m^{(n+1)} +Ra_m \theta_m^{(n)}\boldsymbol{k}, \end{equation}
which we can solve for 
$\boldsymbol {v_m} ^{(n+1)}$ with
\begin{equation} \boldsymbol{v_m} ^{(n+1)} =\left[{-}Da\boldsymbol{\nabla} {\rm \pi}_m^{(n+1)} + Ra_m\theta_m^{(n)}\boldsymbol{k} + \frac{Da}{Pr_{m}\chi {\rm \Delta} t} \boldsymbol{v_m} ^{(n)}\right]\left(\frac{Pr_{m} \chi{\rm \Delta} t}{Da +Pr_{m} \chi{\rm \Delta} t} \right). \end{equation} Next, using 
$(\boldsymbol {v_f} ^{(n)},{\rm \pi} _m^{(n+1)},\theta _f^{(n)})$ – where the previously-found 
${\rm \pi} _m^{(n+1)}$ term is used in an integral along the interface – we find 
$(\boldsymbol {v_f} ^{(n+1)},{\rm \pi} _f^{(n+1)})\in W_f\times Q_f$ by solving
\begin{align} & \left(\frac{1}{Pr_{m}}\frac{\boldsymbol{v_f} ^{(n+1)}-\boldsymbol{v_f} ^{(n)}}{{\rm \Delta} t},\boldsymbol{w}_{\boldsymbol{f}} \right)_f+ {{\boldsymbol{\mathsf{B}}}}_f\left(\boldsymbol{v_f} ^{(n)},\boldsymbol{v_f} ^{(n+1)},\boldsymbol{w}_{\boldsymbol{f}} \right)+2 \left\langle{{\boldsymbol{\mathsf{D}}}}\left(\boldsymbol{v_f} ^{(n+1)} \right), {{\boldsymbol{\mathsf{D}}}}\left(\boldsymbol{w}_{\boldsymbol{f}} \right) \right\rangle_f \nonumber\\ &\quad - \left({\rm \pi}_f^{(n+1)},\boldsymbol{\nabla}\ \boldsymbol{\cdot}\ \boldsymbol{w}_{\boldsymbol{f}} \right)_f + \left(\boldsymbol{\nabla}\ \boldsymbol{\cdot}\ \boldsymbol{v_f} ^{(n+1)}, q_f \right)_f \nonumber\\ &\quad +\frac{Ra_m}{Da}\left(\theta_f^{(n)}, \boldsymbol{w}_{\boldsymbol{f}}\ \boldsymbol{\cdot}\ \boldsymbol{k} \right)_f + \int_{\varGamma_i}{\rm \pi}_m^{(n+1)} \left(\boldsymbol{w}_{\boldsymbol{f}}\ \boldsymbol{\cdot}\ \boldsymbol{n} \right) {\rm d}\varGamma_i \nonumber\\ &\quad + \int_{\varGamma_i}\frac{\alpha}{\sqrt{Da}} \left(\boldsymbol{v_f} ^{(n+1)}\ \boldsymbol{\cdot}\ \boldsymbol{\tau} \right) \left(\boldsymbol{w}_{\boldsymbol{f}}\ \boldsymbol{\cdot}\ \boldsymbol{\tau} \right) {\rm d}\varGamma_i= 0 \quad \forall \boldsymbol{w}_{\boldsymbol{f}}\in W_f,\enspace q_f\in Q_f, \end{align}
where the trilinear term 
${{\boldsymbol{\mathsf{B}}}}_f$ is defined with
\begin{align} 2{{\boldsymbol{\mathsf{B}}}}_f\left(\boldsymbol{u},\boldsymbol{v}, \boldsymbol{w} \right) &= \int_{\varOmega_f}\left(\boldsymbol{u}\ \boldsymbol{\cdot}\ \boldsymbol{\nabla} \boldsymbol{v} \right)\boldsymbol{w} - \left(\boldsymbol{u}\ \boldsymbol{\cdot}\ \boldsymbol{\nabla} \boldsymbol{w} \right)\boldsymbol{v}\,{\rm d}\varOmega_f \nonumber\\ &\quad +\int_{\varGamma_i}\left(\boldsymbol{u}\ \boldsymbol{\cdot}\ \boldsymbol{w} \right)\left(\boldsymbol{v}\ \boldsymbol{\cdot}\ \boldsymbol{n} \right) - \left(\boldsymbol{u}\ \boldsymbol{\cdot}\ \boldsymbol{v} \right)\left(\boldsymbol{w}\ \boldsymbol{\cdot}\ \boldsymbol{n} \right){\rm d}\varGamma_i. \end{align}
The first integral of the trilinear term 
${{\boldsymbol{\mathsf{B}}}}_f$ was first used by Temam in the late 1960s (Temam Reference Temam1968, Reference Temam1969a,Reference Temamb) and has remained a critical tool in approximating solutions to the Navier–Stokes equations since then. The second integral is included to deal with the non-homogeneous boundary value in our application. The skew-symmetric form of the trilinear term 
${{\boldsymbol{\mathsf{B}}}}(\boldsymbol {v_f} ^{(n)},\boldsymbol {v_f} ^{(n+1)},\boldsymbol {w}_{\boldsymbol {f}} )_f$ allows for partial time-lagging of the nonlinear term of Navier–Stokes, which linearizes the problem while still conserving the energy of the system. Treating the nonlinear term in this manner has been used in a number of recent works (Chen et al. Reference Chen, Han, Wang and Zhang2020, Reference Chen, Han, Wang and Zhang2021), albeit for an additional reason as well. With the coupled Navier–Stokes–Darcy equations, use of this trilinear form allows for cancellation of the 
$\frac {1}{2}|\boldsymbol {v_f} |^2$ term from the Lions interface condition (2.18e) with integration by parts on the 
$(\boldsymbol {v}\ \boldsymbol {{\cdot }}\ \boldsymbol {\nabla })\boldsymbol {v}$ term of Navier–Stokes.
With the velocity of both subdomains known, we write them as a single, updated velocity field:
\begin{equation}
\boldsymbol{v}^{(n+1)} = \left\{\begin{array}{@{}ll@{}}
\boldsymbol{v_f} ^{(n+1)} & \textrm{for }
\boldsymbol{x}\in\varOmega_f, \\
\boldsymbol{v_m} ^{(n+1)} & \textrm{for }\boldsymbol{x}\in\varOmega_m.
\end{array}\right.
\end{equation}
With 
$(\boldsymbol {v}^{(n+1)}, \theta ^{(n)})$, we obtain 
$\theta ^{(n+1)}\in \varPsi$ by solving
\begin{align} & \left(\delta_1\frac{\theta^{(n+1)}-\theta^{(n)}}{{\rm \Delta} t}, \psi \right)_\varOmega + Pr_{m}{{\boldsymbol{\mathsf{B}}}}\left(\boldsymbol{v}^{(n+1)},\theta^{(n+1)},\psi \right) \nonumber\\ &\quad + \left(\delta_2 \boldsymbol{\nabla} \theta^{(n+1)},\boldsymbol{\nabla} \psi \right)_\varOmega - \left(\delta_3 \boldsymbol{v}^{(n+1)}\ \boldsymbol{\cdot}\ \boldsymbol{k}, \psi \right)_\varOmega =0\quad \forall\ \psi\in\varPsi, \end{align}
where the coefficients 
$\delta _i$ are defined by
\begin{equation}
\delta_1=\begin{cases} 1
& \textrm{for } \boldsymbol{x}\in\varOmega_f, \\ \varrho &
\textrm{for }\boldsymbol{x}\in\varOmega_m, \end{cases}\quad
\delta_2=\begin{cases} \epsilon_T & \textrm{for }
\boldsymbol{x}\in\varOmega_f, \\ 1 & \textrm{for }
\boldsymbol{x}\in\varOmega_m, \end{cases} \quad
\delta_3=\begin{cases} 1/\epsilon_T & \textrm{for }
\boldsymbol{x}\in\varOmega_f, \\ 1 & \textrm{for
}\boldsymbol{x}\in\varOmega_m,
\end{cases}\end{equation}
and the trilinear term 
${{\boldsymbol{\mathsf{B}}}}$ is defined with
\begin{equation} 2{{\boldsymbol{\mathsf{B}}}}\left(\boldsymbol{u}, v, w \right) = \int_{\varOmega}\boldsymbol{u}\ \boldsymbol{\cdot}\ \left(w\boldsymbol{\nabla} v\right) - \boldsymbol{u}\ \boldsymbol{\cdot}\ \left(v\boldsymbol{\nabla} w\right){\rm d}\varOmega. \end{equation}
These coefficients are discontinuous over the interface since they are related to physical properties of each region, and are chosen to enforce the interface conditions – continuity of the temperature and heat flux. With the updated velocity field over the entire domain, we are also able to solve for the streamlines 
$\phi$:
\begin{equation} \left( \boldsymbol{\nabla}\phi^{(n+1)}, \boldsymbol{\nabla} \varphi\right)_\varOmega - \left(\boldsymbol{\nabla} \times \boldsymbol{v}^{(n+1)}, \varphi \right)_\varOmega = 0\quad \forall \varphi \in \varPhi,\end{equation}
where 
$\varPhi =\{\varphi \in H^1(\varOmega ): \varphi =0 \textrm { on top and bottom }+\textrm{ periodic on left and right}\}$.
Each simulation begins by perturbing the conductive steady state. That is, we apply a seeded, random perturbation (of the order of 
$10^{-6}$) to a stationary fluid and a piecewise-linear temperature profile. The simulations begin with
\begin{equation} \left.\begin{gathered} \boldsymbol{v_f} ^{(0)} = \boldsymbol{v_m} ^{(0)} = \boldsymbol{0} + \boldsymbol{\varepsilon}_1(\boldsymbol{x}), \\ \theta_f^{(0)}=\theta_m^{(0)} = 0 + \varepsilon_2(\boldsymbol{x}), \end{gathered}\right\} \end{equation}
where 
$\boldsymbol {\varepsilon }_1(\boldsymbol {x})$ and 
$\varepsilon _2(\boldsymbol {x})$ are small, seeded, random perturbations. These perturbations – the components of the vector-valued function 
$\boldsymbol {\varepsilon }_1(\boldsymbol {x})$ and the scalar 
$\varepsilon _2(\boldsymbol {x})$ – take the form
\begin{equation} A\sum_{k,\ell={-}N}^{N}\frac{1}{\sqrt{\rm \pi}}\exp\left({-\frac{1}{2}\left(k^2+\ell^2\right)}\right) \sin\left(\frac{k{\rm \pi} x}{L_x} + \alpha_{k,\ell}\right) \sin\left(\frac{\ell{\rm \pi} (y-\hat{d})}{1+\hat{d}} + \beta_{k,\ell}\right), \end{equation}
where 
$\alpha _{k,\ell },\ \beta _{k,\ell }$ are randomized phases from a uniform distribution on 
$[0,2{\rm \pi} )$ and 
$A$ sets the overall magnitude of the perturbation.
For numerical stability, the time step is heavily restricted by the requirement 
${\rm \Delta} t \sim O(Da)$. Our method is first-order in time and second-order in space. The time integrator could be swapped out for a high-order method. However, in practice, we find that the spatial error dominates errors from time discretization, essentially rendering higher-order time integrators unnecessary. Each simulation conducted in this work has a length scale in 
$x$ of 
$2.1$ units, which, through experimentation, we determined to be the smallest domain able to support a pair of counter-rotating, deep convection cells. The length scale in 
$y$ has height 
$1+\hat {d},$ with 
$1$ unit for the porous height and 
$\hat {d}$ for the fluid region. For the spatial discretization, we use a uniform grid with 
$32 \times 32$ triangular elements over a unit area in the non-dimensional domain. Future work is being conducted on using a non-uniform mesh for this system to investigate the associated errors with this discretization. All of the numerical simulations are run using FreeFem++ (Hecht Reference Hecht2012).
5. Results and discussion
From stability analyses and our theory, there are two types of convection we expect – deep and shallow convection. In figure 4, these two flow configurations are shown at their steady states with streamlines shown as contours and the temperature profiles in colour. The system with deep convection exhibits cells which extend throughout the entirety of the domain, evidenced by the wave-like temperature profile and streamlines circulating throughout the whole domain. In contrast, with shallow convection, almost all of the velocity and temperature deviations from the conductive state are confined to the fluid region alone, leaving the porous medium largely unchanged with the exception of a small region immediately below the interface. For these two cases, the systems begin with random initial data and go directly to their preferred convection states.
Figure 4. Shallow and deep convection with temperature (colour) and streamlines (contour) at their respective steady states, with (a) 
$Ra_m=10,\ \hat {d}=0.35$ at 
$t=3.0$ and (b) 
$Ra_m=30,\ \hat {d}=0.20$ at 
$t=3.0$, respectively. Fixed parameters: 
$Da=1.0\times 10^{-4},\ Pr_{m} =0.7,\ \epsilon _T=0.7,\ \alpha =1.0,\ {\rm \Delta} t = 2.5\times 10^{-4}$.
With the case shown in figure 4(a), the fluid evolves from random, disorganized data to a steady state with shallow convection. Conversely, for the simulation in figure 4(b), the system goes from the perturbed conductive state to a stable steady state of deep convection. That is, in these two respective cases, the shallow convection case's cells originate and stay in the fluid region, while cells from the deep convection case occupy the whole domain for the duration of the simulations. This phenomenon is expected; with these two types of convection, instabilities from one group of wavenumbers outperform the other by a significant margin. For example, with deep convection, the growth from small wavenumbers of the perturbation dominates any contribution from the larger wavenumbers, resulting in large cells. Alternatively, with shallow convection, the growth from small wavenumbers pales in comparison to that of the larger wavenumbers.
There is another form of convection that the stability analyses cannot predict though: shallow-to-deep convection, as shown in figure 5. This kind of convection takes place in a parameter regime where both the small and large wavenumbers are unstable, and the interaction between their instabilities gives rise to this hybrid category of convection. In these situations, the system goes from random initial data to a metastable shallow convection state, followed by the collapse and consolidation of cells in the fluid region. The joined cells gain enough strength to penetrate into the porous medium, forming larger-scale cells which occupy the fluid and porous regions alike, and the system finally steadies out at the preferred stable configuration of deep convection. For this choice of parameters, this route – from random initial data to a metastable shallow state followed by a dynamic shift to deep convection – is the route taken from perturbing the unstable conductive state. While we observe shallow-to-deep convection, we never observe a deep-to-shallow dynamic shift; this leads us to speculate that when both small and large wavenumbers are unstable, deep convection is the preferred steady state of this system.
Figure 5. Shallow-to-deep convection with 
$Ra_m=20,\ \hat {d}=0.30$ with temperature (colour) and streamlines (contour). Fixed parameters: 
$Da=1.0\times 10^{-4},\ Pr_{m} =0.7,\ \epsilon _T=0.7,\ \alpha =1.0,\ {\rm \Delta} t = 2.5\times 10^{-4}$.
Although shallow-to-deep convection and deep convection ultimately both achieve a steady state with cells that extend throughout the domain, they are distinct flow progressions. Their most notable differences concern the routes taken to their steady state and the efficiency of heat transfer; shallow-to-deep convection has metastable shallow cells before forming deep convection cells for its steady state, and these flow profiles can exhibit higher Nusselt values compared with their counterparts with deep convection.
To measure convection, we define two quantities. The first is a mathematical energy that allows us to find the 
$Ra_m$ value where 
${\textrm {d} E}/{\textrm {d} t}= 0$, distinguishing between the regions of stability and instability – with 
${\textrm {d} E}/{\textrm {d} t}< 0$ and 
${\textrm {d} E}/{\textrm {d} t}> 0$, respectively – as done in McCurdy et al. (Reference McCurdy, Moore and Wang2019). This uses the volume-averaged norms of the perturbed velocity and temperature profiles, measuring how these profiles compare with their conductive steady states:
\begin{equation} 2E(t) = \frac{1}{Pr_{m}}\frac{1}{V} \|\boldsymbol{v} \|_{\varOmega}^2 +\frac{1}{V} \|\theta\|_{\varOmega}^2, \end{equation}
with 
$V$ as the volume. Qualitatively, the onset of convection is noted by a jump in the energy profile. Second, we use the physically motivated Nusselt number 
$Nu(t)$ as the ratio of vertical convective and conductive fluxes with
\begin{equation} \left.\begin{gathered} Nu(t) = \frac{J_{cnv } + J_{cnd}}{J_{cnd}} \quad \textrm{with } J_{cnv} = \left(\boldsymbol{u}\ \boldsymbol{\cdot}\ \boldsymbol{k},\ \left(\theta + \bar{T}\right)\right)_{\varOmega}\\ \textrm{and}\quad J_{cnd}={-}\left(\delta_2\boldsymbol{\nabla} \left(\theta + \bar{T}\right)\ \boldsymbol{\cdot}\ \boldsymbol{k},1\right)_{\varOmega}, \end{gathered}\right\} \end{equation}
where 
$\bar {T}$ is the conductive temperature and 
$\delta _2$ is defined in (4.10a–c), which allows us to write the conductive fluxes of each region as a single term.
In figure 6, we plot the energy and Nusselt profiles for the three different qualitative states presented in figures 4 and 5. For each curve, the initial jump indicates the time of noticeable formation of convection cells in their respective cases with their steady states achieved as the profiles level out. The deep and shallow cases (shown in red and blue, respectively) show the formation of cells, immediately followed by a steady state. The phenomenon of shallow convection cells forming faster than deep cells is expected. Shallow cells have nothing obstructing their formation, while deep convection is hindered by the porous medium. As a result, it takes longer for the fluid to gain enough velocity in the porous region for noticeable deep convection cells to form. More precisely, the characteristic time scale for the porous-medium flow given in (2.15a–e) is 
${d_m^2}/{\lambda _m}$, which has been normalized to unity in our analysis. The flow time scale for the fluid region, meanwhile, is 
${d_f^2}/{\lambda _f}$, which is shorter by a factor of 
$\hat {d}^2/\epsilon _T$.
These time scales can be used to better understand the shallow-to-deep convection case, which is shown by the yellow curves in figure 6. The first jump in 
$E$ and 
$Nu$ occurs as the convection cells form in the free-flow zone, followed sometime later by a secondary jump whereby the cells coalesce into the larger, deep convection cells and push into the porous medium. As noted above, the time scale for porous-medium flow has been normalized to unity and the time scale for free-fluid flow is shorter by a factor of 
$\hat {d}^2/\epsilon _T$ or about 0.09 in the case shown. These values are consistent with a quick initial onset free-flow convection, followed by a slower shift to deep convection.
The jumps in 
$Nu$ with the shallow-to-deep case demonstrate that these configurations have the potential to transfer heat more efficiently than their counterparts with deep or shallow convection. Further investigations into this hybrid flow pattern could determine the parameter regimes which produce shallow-to-deep convection to maximize heat flow in various applications, like heat sink technology. For example, in comparing the deep convection and shallow-to-deep convection cases with the 
$E(t)$ and 
$Nu(t)$ profiles, we see both cases eventually steady out to approximately the same energy; that is, they have comparable deviations from their conductive states. However, the case with shallow-to-deep convection has a larger Nusselt number than its deep-convection counterpart, indicating a greater efficiency in transferring heat. We speculate that shallow-to-deep convection is more efficient in its heat transfer than deep convection due to the route each takes to its steady state. For shallow-to-deep convection, as cells form and gain speed in the free-flow region, there is a spike in the 
$Nu$ profile. As the shallow cells consolidate and penetrate into the medium, a second spike in the Nusselt value occurs.
The parameters for the simulations with deep and shallow convection are chosen so that a single group of wavenumbers was contributing to the perturbations. These calculations come from looking at the marginal stability curves that identify the 
$Ra_m$ values for the onset of convection at a given depth ratio. For larger Rayleigh numbers, where the small and large wavenumbers both contribute to perturbation growth, the stability analyses cannot predict the preferred flow configuration. To provide a more holistic picture of flow behaviour, we use marginal stability analyses and numerical simulations to identify regions in a 
$\hat {d}$–
$Ra_m$ phase space where each configuration occurs.
Marginal stability analyses are able to predict the 
$Ra_m$ values for the onset of convection at a given depth ratio (shown with the solid black lines of figure 7). There is a small region above the critical Rayleigh numbers where a single group of wavenumbers is unstable, resulting in either deep or shallow convection, and the upper bound is denoted in the phase-space diagrams with the dashed lines. For example, with 
$\hat {d} = 0.2$ and 
$Ra_m=25$, only the small wavenumbers are unstable, resulting in deep convection. However, once the Rayleigh number is increased past the dashed line, say to 
$Ra_m=50$ or 
$Ra_m=80$, the larger wavenumbers are unstable as well. In regions where both groups of wavenumbers are unstable, linear stability analysis and our coarse-grained model are not able to accurately predict convection patterns. In practice, we find the deep or shallow convection regimes extend above the predicted regions from the stability analyses, as evidenced with the green and red triangles extending into the blue region where both groups of wavenumbers are unstable. The distance in which the red triangles persist into the blue region tells us that for these Rayleigh numbers, the small wavenumbers (which produce the deep convection cells) are only slightly unstable, and not enough to overpower the instabilities coming from the large wavenumbers which produce the shallow cells. For Rayleigh numbers well above these regions in the phase space, where larger and smaller wavenumbers can interact and our analytical tools are ineffective, shallow-to-deep convection takes place.
Figure 7. Phase space with different types of convection noted. Solid black lines denote the 
$Ra_m$ values for marginal stability, and dashed lines mark the 
$Ra_m$ values where both large and small wavenumbers become unstable. Fixed parameters: 
$Da=1.0\times 10^{-4},\ Pr_{m} =0.7,\ \epsilon _T=0.7, \alpha =1.0,\ {\rm \Delta} t = 2.5\times 10^{-4}$.
The phase space shown in figure 7 ties our theory, the stability analyses and numerical simulations together nicely. Marginal stability predicts the Rayleigh numbers needed for the onset of convection, and the wavenumbers associated with the unstable modes ultimately dictate the flow configuration for Rayleigh numbers near their corresponding critical value. This is represented in the phase diagram by the region between the solid and dashed lines where the marginal stability curves can effectively determine the behaviour of convection. The shift in the unstable wavenumbers occurs where the solid and dashed lines come together, noting the shift in convection, that can accurately be determined by our theory. The numerical simulations agree with results from the stability analyses and theory for Rayleigh numbers around the onset of convection. Further, the numerics allow for a more comprehensive understanding of superposed fluid–porous convection with higher Rayleigh numbers.
6. Conclusions
Based on the framework developed in McCurdy et al. (Reference McCurdy, Moore and Wang2019), we proposed a coarse-grained model to predict the critical depth ratio needed for the transition from deep to shallow convection, and we observed the formula is increasingly accurate in the limit of 
$Da\rightarrow 0$. The previous stability analyses and the new transition theory, when viewed in tandem with these novel numerical simulations, provide a more complete view into the phenomenon of convection in superposed fluid–porous media systems.
The coarse-grained model presented was formed under the assumption that heuristically, the shift between deep convection and shallow convection takes place when the Rayleigh numbers of the fluid and porous regions are equivalent in some sense. In the physically relevant small-Darcy-number regime, we deduced appropriate boundary conditions for the uncoupled regions to determine their critical Rayleigh numbers, and as a result, we saw that the formula became more accurate in the small-Darcy limit of 
$Da\rightarrow 0$.
To explore the extensive parameter regime outside the limitations of our model, we outlined a numerical scheme to simulate the full nonlinear system. The main contributions were with the decoupling and time-lagging of nonlinear terms; these adjustments to the system allow for linear, decoupled, sequential solvers to be used in approximating solutions.
With our simulations, we verified results from our coarse-grained model and previous stability analyses, namely with deep and shallow convection. However, we additionally observed a type of convection not able to be predicted by marginal stability, shallow-to-deep convection. This kind of convection presents an exciting new direction of work, as there is potential to investigate the dynamic reorganization of the flow. Additionally, the numerical simulations provide opportunities to explore less structured domains, perhaps with more complex interfaces (Allen Reference Allen1984; Han & Wang Reference Han and Wang2016), boundaries that evolve in time (Zhang & Libchaber Reference Zhang and Libchaber2000; Moore Reference Moore2017; Quaife & Moore Reference Quaife and Moore2018; Chiu, Moore & Quaife Reference Chiu, Moore and Quaife2020; Mac Huang & Moore Reference Mac Huang and Moore2022) or with some permeable membrane developing between the regions (Sorribes et al. Reference Sorribes, Moore, Byrne and Jain2019; Eastham et al. Reference Eastham, Moore, Cogan, Wang and Steinbock2020).
Acknowledgements
We would like to thank the reviewers for carefully critiquing the manuscript. With their suggestions, we were able to better convey our arguments and present our results in a clearer manner.
Funding
This work was supported by the National Science Foundation (N.J.M., grant DMS-2012560) and the National Natural Science Foundation of China (X.W., grants 12271237 and 11871159).
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are openly available at http://doi.org/10.5281/zenodo.7320220.
Appendix
We present data from related papers that observed a transition from deep convection to shallow convection as the depth ratio is altered. Each of the papers we consider uses a two-domain approach in modelling the system (so that the interface conditions can be explicitly enforced), and they each use a linear equation of state for the Boussinesq approximation. There are a variety of approaches taken in each work, from analytical methods to numerical simulations to experiments. Table 3 references the approach taken in each paper, and we expand on the analytical approaches here: Chen & Chen (Reference Chen and Chen1988), McKay (Reference McKay1998), Straughan (Reference Straughan2002), Hirata et al. (Reference Hirata, Goyeau, Gobin, Carr and Cotta2007) and Yin et al. (Reference Yin, Fu and Tan2013) each use linear stability theory, Hill & Straughan (Reference Hill and Straughan2009) and McCurdy et al. (Reference McCurdy, Moore and Wang2019) use nonlinear stability arguments and Han et al. (Reference Han, Wang and Wang2020) investigate this bifurcation using centre manifold reduction theory.
Table 3. Data from papers that note a critical depth ratio along with the prediction for 
$\hat {d}^*$ from our theory (3.5) compared with the prediction made using the old model. The approach taken by the authors is noted alongside the citation, with a majority of the analytic approaches taking the form of investigating marginal stability curves and Han et al. (Reference Han, Wang and Wang2020) making use of centre manifold reduction theory to determine where the transition occurs. Notation for the governing equations: (a) Navier–Stokes–Darcy–Boussinesq; (b) Navier–Stokes–Darcy–Brinkman–Boussinesq, with nonlinear Forchheimer term in Darcy; (c) Navier–Stokes–Darcy–Brinkman–Boussinesq; (d) Navier–Stokes–Darcy–Boussinesq with an Oldroyd-B fluid. 
${\dagger}$ From the two-domain approach in Hirata et al. (Reference Hirata, Goyeau, Gobin, Carr and Cotta2007).
Table 3 shows the marked improvements in predicting the critical depth ratio 
$\hat {d}^*$ using our theory from (3.5) compared with predictions made using the old model. Despite several works operating under different assumptions about the system – like different governing equations (e.g. Darcy–Brinkman in lieu of Darcy), assuming the porous medium has a high porosity or using a non-Newtonian fluid – we find our new model produces an accurate prediction for the critical depth ratio.
Our formula in (3.5) can be rearranged to solve for the critical 
$\epsilon _T$ or 
$Da$ values, and two works, Yin et al. (Reference Yin, Fu and Tan2013) and Han et al. (Reference Han, Wang and Wang2020), note how the transition from deep to shallow convection is affected by changing these parameters. Table 4 shows our prediction for the critical values of these parameters alongside their true values. We once again find better agreement between the true values and predictions made with the new model compared with our old formula.
Table 4. Data from related papers that observed a transition from deep convection to shallow convection as the ratio of thermal diffusivities 
$\epsilon _T$ or the Darcy number 
$Da$ is altered. We present our prediction for the critical value from the theory presented in (3.5) (rearranged to solve for 
$\epsilon _T^*$ or 
$Da^*$ instead of 
$\hat {d}^*$) compared with the prediction made using the old model. The approach taken by the authors is noted alongside the citation, with Yin et al. (Reference Yin, Fu and Tan2013) investigating marginal stability curves and Han et al. (Reference Han, Wang and Wang2020) making use of centre manifold reduction theory to determine where the transition occurs. Notation for the governing equations: (i) Navier–Stokes–Darcy–Boussinesq with an Oldroyd-B fluid; (ii) Navier–Stokes–Darcy–Boussinesq.
Open access
