1. Introduction
This paper is devoted to the formation of flows in a stratified fluid at a violation of stable stratification. It is known that the motionless state of a fluid heated strictly from above is stable with respect to any disturbances (Sorokin (Reference Sorokin1954), see also Gershuni & Zhukhovitsky (Reference Gershuni and Zhukhovitsky1976)). A similar situation is observed in solutal convection (in the absence of thermal diffusion), when the concentration depends linearly on the vertical coordinate and the concentration of the denser component is greater at the bottom (Gershuni & Zhukhovitsky Reference Gershuni and Zhukhovitsky1976). Since stable stratification prevents vertical flows associated with the work against gravity, but does not prevent horizontal flows, situations may arise when a weak violation of stable stratification (a weak violation of the homogeneity of the density gradient) can lead to the emergence of flows with a horizontal component of velocity significantly exceeding the vertical one.
Stable stratification can be violated near inclined walls. This problem was studied for the first time by Prandtl (Reference Prandtl1952) who investigated it in application to the downslope wind in the mountains. He considered an inclined sidewall having a constant temperature deviation from the temperature profile corresponding to stable stratification occurring far from the wall and found a flow in the boundary layer near the wall. His analytical solution shows that, if the wall is hotter than the surrounding air, an upward flow along the wall occurs, while in the opposite case a downward flow arises. His work also provides qualitative considerations of the flows at the ends of the wall. For the case of upward flow, a sink is located at the bottom of the mountain into which the valley wind rushes, which can move mainly horizontally due to the stable stratification of the air. The flow that occurs at the top is less noticeable. Wunsch (Reference Wunsch1970) and Phillips (Reference Phillips1970) independently and almost simultaneously investigated a problem similar to that studied by Prandtl, but for ocean flows. They considered two similar tasks: the effect of inclined thermally insulated walls on thermal convection in the case when the fluid is stably stratified far from the wall, and the effect of an impermeable wall on solutal convection in a similar situation. Flows occurring in the boundary layer were investigated. An analytical solution was found. In the second case, mixing occurs near the walls. Obviously, by analogy with Prandtl's problem, there should be sinks and sources at the ends of the wall. Shapiro & Burkholder (Reference Shapiro and Burkholder2012) also performed a study that can be considered as a continuation of Prandtl work. They analytically and numerically investigated the case where a finite-length section on an inclined surface is cooled compared with the equilibrium temperature profile in the environment. Two variants of boundary conditions were considered: the first variant, when the boundary condition is set for the temperature (analogous to the Prandtl condition); and the second, when the boundary condition is set for the normal component of the temperature gradient. In the main part near the wall, they obtained a downward flow, and at the ends of the cooled section, two nearly horizontal jets. This last result confirms Prandtl qualitative considerations.
In Baidulov, Matyushin & Chashechkin (Reference Baidulov, Matyushin and Chashechkin2007), the formation of the diffusion-induced flow around a sphere was studied experimentally and theoretically for the case when the stable stratification was created due to the gradient of solute concentration. Initially, both the sphere and the stratified fluid were at rest, further flows were produced by the small buoyancy force arising due to the inhomogeneity of the horizontal density distribution around a sphere. Numerical simulations have shown that the cellular structures stretched in the horizontal direction are formed around a sphere. In this work, the emphasis was on the investigation of the initial and transient stages of the formation of convective structures, and no special attention was paid to steady flows and their characteristics.
There are also works in the literature devoted to the occurrence of flows in stably stratified binary mixtures, when the stable stratification created by a vertical concentration gradient is violated by horizontal heating. Chen, Briggs & Wirtz (Reference Chen, Briggs and Wirtz1971) experimentally investigated the effect of lateral heating on convective phenomena in a vertical layer containing a salt solution with a vertical concentration gradient providing stable stratification. They found the formation of layered convection when the critical Rayleigh number was exceeded. In this case, a system of horizontal narrow convective cells was observed in the layer. They called the situation unstable when these cells were extended from one wall to the other. Lee & Hyun (Reference Lee and Hyun1991) investigated this problem numerically. They observed the growth of convective cells in the horizontal direction from the heated wall. At a small Rayleigh number, these cells had a small horizontal length and did not reach the other wall. Kamakura (Reference Kamakura1993) investigated a similar problem experimentally and numerically. Common for these studies is the discovery of layered convection in the form of narrow horizontal convective cells, which is associated with the suppression of vertical flows under conditions of stable stratification.
Stable stratification can also be violated by a solid body immersed in the fluid whose thermal conductivity differs from that of the fluid. In this case, different from the case considered by Prandtl (Reference Prandtl1952) and Wunsch (Reference Wunsch1970), the inclusion surface is curved and the inclusion dimensions are limited. This problem was studied in Singh (Reference Singh1977) where a steady convective flow around a sphere in an infinite fluid heated from above was considered. The problem was solved analytically by expanding in a series, in a small parameter, the Grashof number, within the Stokes and Oseen approximations. Using the Stokes approximation cannot give a solution satisfying the boundary conditions at infinity. For this reason, the flow far from the body was considered within the Oseen approximation, and then the solutions found within the two approximations were matched. In fact, the method of matched asymptotic expansion was used. Two cases were considered: (i) when the thermal conductivity of the sphere is low; (ii) when the surface temperature of the sphere is maintained constant and equal to the temperature of the fluid in the central plane. It was found that in the case of low thermal conductivity of the sphere, the flow is similar to that obtained earlier for a stationary sphere in a rotating fluid. In the case of an isothermal sphere, the flow pattern was similar to the flow around a sphere uniformly rotating in a fluid which is at rest at infinity. The results obtained in Singh (Reference Singh1977) are in good agreement with experimental observations of natural convective flow around an isothermal sphere or cylinder immersed in a thermally stratified fluid (Eichhorn, Lienhard & Chen Reference Eichhorn, Lienhard and Chen1974).
Singh's work is limited to the consideration of small Grashof numbers at which an analytical solution to the problem can be obtained, which is why the horizontal structures that could be expected to occur in a stratified fluid at a violation of stable stratification were not found. We study a problem similar to that investigated by Singh, but for a porous medium saturated with fluid. The simple form of the equations of thermal buoyancy convection in a saturated porous medium in the Darcy–Boussinesq approximation should allow for a fairly complete analytical study of the problem (over a wide range of parameters) and a description of the expected horizontal structures.
The importance of the problem under study is related to the following issues. The first issue is related to thermal convection in soils containing water. In permafrost regions, at a certain depth, the soil has been frozen for many centuries. In winter, the upper part of the soil (active layer) above the permafrost freezes completely or partially (in the latter case, a layer of thawed soil remains). In the summer, during the hottest period, the active layer of soil melts. At the bottom, there remains unmelted soil (permafrost). At the same time, in the part of the layer where the temperature is higher than the density inversion temperature, conditions arise for the formation of stable thermal stratification. If there are solid inclusions in the active layer, convection may occur there. Such inclusions can be various utility lines, for example, pipes. Convection in permafrost conditions in the active layer of soil, in particular, was considered in Magnani et al. (Reference Magnani, Musacchio, Provenzale and Boffetta2024). The other issue is associated with the solutal convection, which is very similar to thermal convection. Here one can distinguish two classes of problems. First, soils under certain conditions may contain salt water with density stratification. This may occur when coastal regions containing fresh water aquifers are flooded with sea water. If this water is used for water supply purposes, it is important from a practical point of view to know how sea water interacts with fresh water. Also, sea water intrusion into coastal regions is important for agriculture, if it is carried out in these regions. Density stratification may occur not only in salt water, but also in oil (in oil fields). Density stratification of water may also occur when hot water or water containing dissolved salts is artificially injected through recharge wells or infiltration ponds. Second, there is the problem of waste spreading in groundwater (see, for example, the works Kodešová et al. (Reference Kodešová2023, Reference Kodešová2024)).
The study of convective flows in porous media with solid inclusions is also important for hydrology due to the possibility of such flows in fractured rocks in the presence of a geothermal temperature gradient and a salt concentration gradient. The effect of solid square blocks placed inside a porous medium filled with fluid on the onset of convection was studied by Rees & Nield (Reference Rees and Nield2016). Although their work is devoted to the effect of solid blocks on the onset of convection, in reality they discuss the effect of the blocks on the stability of the ground state, which is not a stationary state. In the ground state, there are vortices at the corners of the square blocks. However, this remark is purely terminological.
The structure of the paper is as follows. In § 2, the governing equations are formulated. Analytical results obtained in the framework of Stokes and Oseen type approximations are presented in §§ 3 and 4. Finally, § 5 presents the results of numerical simulations based on the complete nonlinear equations.
2. Governing equations
Let us consider convective flows near a solid inclusion in a porous medium saturated with fluid. The inclusion has the shape of a horizontal cylinder of circular cross-section (see figure 1). We write down the equations of convective filtration in the framework of the Darcy–Boussinesq approximation,
$$\begin{gather} - \boldsymbol{\nabla} p - {\boldsymbol{v}} + Ra\,T \boldsymbol{\gamma} = 0, \end{gather}$$
$$\begin{gather}\frac{{\partial T}}{{\partial t}} + ({\boldsymbol{v}} \boldsymbol{\cdot} \boldsymbol{\nabla}) T = \Delta T, \end{gather}$$
$$\begin{gather}{\rm{div}}\,{\boldsymbol{v}} = 0. \end{gather}$$
Here, 
${\boldsymbol {v}}$ is the velocity, 
$p$ is the convective addition to the hydrostatic pressure, 
$T$ is the temperature counted from a reference value, 
$\boldsymbol {\gamma }$ is the unit vector directed vertically upwards, 
$Ra = g\beta AKR^2 /(\nu \chi )$ is the analogue of the Rayleigh number for a porous medium (
$g$ is the gravity acceleration, 
$\beta$ is the thermal expansion coefficient, 
$A$ is the value of imposed temperature gradient, 
$K$ is the permeability, 
$R$ is the inclusion radius, 
$\nu$ is the kinematic viscosity of fluid and 
$\chi = \kappa _m /(\rho c_p )_f$ is the effective thermal diffusivity of medium (
$\kappa$ is the thermal conductivity, 
$\rho$ is the density, 
$c_p$ is the specific heat; the quantities with the subscript 
$m$ are related to the porous medium and those with the subscript 
$f$ – to the fluid)). We assume that conditions are uniform in the axial direction. In this case, the problem allows two-dimensional solutions and we restrict our analysis to such solutions. Let us discuss the justification of this approach. Due to the type of boundary conditions for temperature, at small Rayleigh numbers the flow should not depend on the 
$z$ coordinate (the horizontal coordinate along the cylinder axis), i.e. it is two-dimensional. Three-dimensional effects are possible only at large Rayleigh numbers. In this case, the length of the vortices will be a function of the 
$z$ coordinate. But the existence of these vortices is not obvious. The reasons for that are the following. If we were dealing with usual convection between coaxial cylinders with an inner cylinder having a thermal conductivity equal to that of the fluid, heated from below, then, firstly, the motionless state would be possible, and, secondly, at Rayleigh numbers greater than the critical value, convective flows with two opposite circulation directions would be possible, in which case a flow dependent on the 
$z$ coordinate would naturally be included in the problem. In the problem under consideration, it is not so. We have a basic state determined by the boundary conditions that limit the circulation direction. In addition, since we have stable stratification, disturbances in the form of convective flows, similar to the case of usual convection, are suppressed.
Geometry of the problem.
Restricting ourselves to such solutions and introducing the stream function associated with the velocity components by the relations
\begin{equation} {v}_x ={-} \frac{{\partial \psi }}{{\partial y}},\quad {v}_y = \frac{{\partial \psi }}{{\partial x}}, \end{equation}
we obtain the following equations in terms of 
$\psi$ and 
$T$:
$$\begin{gather} \Delta \psi = Ra\frac{{\partial T}}{{\partial x}}, \end{gather}$$
$$\begin{gather}\frac{{\partial T}}{{\partial t}} + J(\psi ,T) = \Delta T. \end{gather}$$
Here, 
$J(\psi,T)$ is the Jacobian with respect to 
$x,y$ (
$x$ is the horizontal coordinate).
To complete the set of governing equations we need to add the energy equation in the inclusion
\begin{equation} \alpha \frac{{\partial \tilde T}}{{\partial t}} = \kappa \Delta \tilde T. \end{equation}
Here, 
$\tilde T$ is the temperature in the inclusion, 
$\kappa = \kappa _s / \kappa _m$, 
$\alpha = (\rho c_p)_s/(\rho c_p)_m$, indices 
$s$ and 
$m$ are related to the inclusion and the porous medium, respectively.
3. Analytical study in the framework of the Stokes-like approximation
Let us at first consider the case where the inclusion is embedded in a saturated porous matrix of infinite extent.
Far from the inclusion, at 
$r \to \infty$, the flow is absent and a constant vertical temperature gradient corresponding to stable stratification (heating from above) is maintained,
\begin{equation} v_{r}={-} \frac{1}{r} \frac{\partial \psi}{\partial \varphi} \to 0,\quad T \to y. \end{equation}
Here, 
$y$ is the vertical coordinate.
On the surface of the inclusion, at 
$r = 1$, we impose the impermeability condition, continuity of temperature and thermal flux:
$$\begin{gather} \psi = 0, \end{gather}$$
$$\begin{gather}T = \tilde T, \end{gather}$$
$$\begin{gather}\frac{\partial T}{\partial r} = \kappa \frac{\partial \tilde T}{\partial r} . \end{gather}$$In the centre of the inclusion we suppose the absence of any singularities.
Let us look for stationary solutions.
Firstly consider the case of small Rayleigh numbers, using the Stokes-like expansion:
$$\begin{gather} \psi= Ra \,\psi_1 + Ra^2 \psi_2 + \cdots , \end{gather}$$
$$\begin{gather}T= T_0 + Ra\, T_1 + \cdots , \end{gather}$$
$$\begin{gather}\tilde{T}= \tilde{T}_0 + Ra \,\tilde{T}_1 + \cdots. \end{gather}$$The equation for the stream function in the first order and the equations for temperatures in the zeroth order in cylindrical system of coordinates have the form
$$\begin{gather} \left[ \frac{1}{r}\frac{\partial}{\partial r}\left(r \frac{\partial}{\partial r}\right) + \frac{1}{r^2} \frac{\partial^2}{\partial \varphi^2} \right] \psi_1= \left( \cos{\varphi} \frac{\partial}{\partial r} - {\frac{\sin{\varphi}}{r} \frac{\partial}{\partial \varphi}} \right) T_0 , \end{gather}$$
$$\begin{gather}\left[ \frac{1}{r}\frac{\partial}{\partial r}\left(r \frac{\partial}{\partial r}\right) + \frac{1}{r^2} \frac{\partial^2}{\partial \varphi^2} \right] T_0= 0 , \end{gather}$$
$$\begin{gather}\left[ \frac{1}{r}\frac{\partial}{\partial r}\left(r \frac{\partial}{\partial r}\right) + \frac{1}{r^2} \frac{\partial^2}{\partial \varphi^2} \right] \tilde{T}_0= 0 . \end{gather}$$The boundary conditions are
$$\begin{gather} r=1:\quad \psi_1=0, \quad T_0=\tilde{T}_0, \quad \frac{\partial T_0}{\partial r} - \kappa \frac{\partial \tilde{T}_0}{\partial r}= 0 , \end{gather}$$
$$\begin{gather}r \to\infty :\quad v_{r1}={-} \frac{1}{r} \frac{\partial \psi_1}{\partial \varphi} \to 0, \quad T_0 \to r \sin{\varphi}. \end{gather}$$The problem (3.8)–(3.12) has the following solution:
$$\begin{gather} \psi_1 ={-} \frac{1}{4} \frac{( \kappa - 1 )}{(\kappa + 1 )} \left(1 - \frac{1}{r^2} \right) \sin{ 2 \varphi } , \end{gather}$$
$$\begin{gather}T_0 = r \sin{\varphi} - {\frac{(\kappa - 1)}{(\kappa + 1)} } \frac{1}{r} \sin{\varphi}, \end{gather}$$
$$\begin{gather}\tilde{T}_0 = r \sin{\varphi} - {\frac{(\kappa - 1)}{(\kappa + 1)} } r \sin{\varphi} . \end{gather}$$The radial component of the velocity is
\begin{equation} v_{r1} ={-} \frac{1}{r} \frac{\partial \psi_1}{\partial \varphi}= \frac{1}{2} {\frac{( \kappa - 1 )}{(\kappa + 1 )}} \left(\frac{1}{r} - \frac{1}{r^3} \right) \cos{ 2 \varphi}. \end{equation}
In the case of high thermal conductivity of the inclusion (
$\kappa \to \infty$) instead of solutions (3.13)–(3.16) we have
$$\begin{gather} \psi_1 ={-} \frac{1}{4} \left(1 - \frac{1}{r^2} \right) \sin{ 2 \varphi } , \quad v_{r1} = \frac{1}{2} \left(\frac{1}{r} - \frac{1}{r^3} \right) \cos{ 2 \varphi }, \end{gather}$$
$$\begin{gather}T_0 = 0, \quad \tilde{T}=0. \end{gather}$$It can be seen that the resulting solution in the considered orders of smallness is uniformly valid at any distance from the inclusion.
Consider now the temperatures in the first order of smallness (for finite 
$\kappa$). In this order the equations for temperatures have the form
$$\begin{gather} \left[\frac{1}{r}\frac{\partial}{\partial r}\left(r \frac{\partial}{\partial r}\right) + \frac{1}{r^2} \frac{\partial^2}{\partial \varphi^2} \right] T_1 = \frac{1}{r} \left( \frac{\partial T_0}{\partial \varphi} \frac{\partial }{\partial r} - \frac{\partial T_0}{\partial r} \frac{\partial }{\partial \varphi} \right) \psi_1 , \end{gather}$$
$$\begin{gather}\left[ \frac{1}{r}\frac{\partial}{\partial r}\left(r \frac{\partial}{\partial r}\right)+ \frac{1}{r^2} \frac{\partial^2}{\partial \varphi^2} \right] \tilde{T}_1= 0. \end{gather}$$The boundary conditions are
$$\begin{gather} r=1:\quad T_1=\tilde{T}_1,\quad \frac{\partial T_1}{\partial r} - \kappa \frac{\partial \tilde{T}_1}{\partial r}= 0 , \end{gather}$$
$$\begin{gather}r \to \infty :\quad T_1=0. \end{gather}$$Equations (3.19)–(3.20) have the following solution:
\begin{gather} \begin{aligned}T_1& = \left(A r + \frac{B}{r} - \frac{1}{4} \frac{\kappa-1}{\kappa+1} r \log{r} - \frac{1}{8} \frac{(\kappa-1)^2}{(\kappa+1)^2} \frac{\log{r}}{r} \right) \sin{\varphi} + \left( \frac{D}{r^3} - \frac{1}{32} \frac{\kappa-1}{\kappa+1} r \right. \nonumber\\ &\quad \left. + \frac{1}{32} \frac{(\kappa-1)(3 \kappa+1)}{(\kappa+1)^2} \frac{1}{r} - \frac{1}{12}\frac{(\kappa-1)^2}{(\kappa+1)^2} \frac{\log{r}}{r^3} \right) \sin{3 \varphi},\end{aligned} \end{gather}
\begin{gather} \tilde{T}_1 = \tilde{A} r \sin{\varphi} + \tilde{C} r^3 \sin{3\varphi} . \end{gather} It can be seen that the solution for 
$T_1$ cannot satisfy the boundary conditions at infinity for arbitrary values of constants 
$A$, 
$B$, 
$D$. The reason for this is the same as in the well-known Whitehead paradox for the problem of the flow around the bodies.
Let us estimate the term 
$(\boldsymbol {v} \boldsymbol {\cdot } \boldsymbol {\nabla }) T$ and terms describing thermal dissipation (terms in equation (2.2)). For this, we use the solutions (3.13) and (3.14). Taking into account that 
$T \approx {T_0}$, 
$\psi \approx Ra \psi _1$, we obtain the following order of magnitude for the nonlinear term at large distances, 
$Ra/r$. A typical dissipative term has the form
\begin{equation} \left( {\frac{1}{r}\frac{\partial }{{\partial r}} + \frac{1}{{{r^2}}}\frac{{{\partial ^2}}}{{\partial {\varphi ^2}}}} \right)T = \frac{{\kappa - 1}}{{\kappa + 1}}\frac{1}{{{r^3}}}, \end{equation}
i.e. it has the order 
$1/{r^3}$. The ratio of these two members is
\begin{equation} \frac{\hbox{Nonlinear term}}{\hbox{Nonlinear thermal dissipation term}} = O(Ra\,r^2) \quad \text{at} \ r \to \infty, \end{equation}
can be seen from this relation, at 
$r \to \infty$ this ratio tends to infinity and nonlinear terms cannot be neglected. This means that nonlinear terms can only be neglected for 
$r \ll 1/\sqrt {Ra}$. In dimensional form, this condition has the form 
${r_d} \ll l \equiv \sqrt {{{\nu \chi }}/{{g\beta AK}}}$, where 
${r_d}$ is the dimensional radial coordinate. Only for such small 
${r_d}$ can the Stokes approximation be used. For large 
$r_d$, nonlinear terms must be taken into account. In this connection, the Oseen approximation is used.
To resolve this problem the Oseen-like approximation can be used.
4. Analytical study in the framework of the Oseen-like approximation
4.1. Case of the inclusion of high thermal conductivity
Let us implement the Oseen-like approximation but, different from the conventionally used one, apply the quasilinearization to nonlinear terms in the energy equation and not in the momentum equation.
Neglecting nonlinear terms in the momentum equation for porous media is a common practice, since the presence of a porous skeleton greatly slows down the motion of the fluid. In this case, the nonlinear terms are much smaller in order of magnitude than the friction force. Darcy law of friction was originally suggested neglecting nonlinear terms. If we consider taking into account nonlinear terms, then this is primarily the Forchheimer term, 
$- {c_F}{K^{ - 1/2}}{\rho _F}| \boldsymbol {v} |\boldsymbol {v}$, and not the term 
$(\boldsymbol {v} \boldsymbol {\cdot } \boldsymbol {\nabla })\boldsymbol {v}$. One of the reasons for this is that in most situations, the term 
$- {c_F}{K^{ - 1/2}}{\rho _F}| \boldsymbol {v} |\boldsymbol {v}$ is larger in order of magnitude than the term 
$(\boldsymbol {v} \boldsymbol {\cdot } \boldsymbol {\nabla })\boldsymbol {v}$. This issue is discussed in more detail in the book by Nield 2006.
The classical Oseen approximation is valid for sufficiently small Reynolds numbers. In this case, the nonlinear terms at small distances from the inclusion are much smaller than the friction force. Quasilinearization is required only because at large distances the nonlinear term exceeds the friction force by an order of magnitude. In the case of large Reynolds numbers, the use of boundary layer concepts is required. This is the opposite limiting case (see Dyke Reference Dyke1964).
Note that if we start from the nonlinear version of the momentum equation, then within the framework of the Oseen approximation, we will still have to omit the nonlinear term (at least in the leading order of smallness). This is due to the fact that when using the Oseen approximation, we replace the nonlinear term 
$(\boldsymbol {v} \boldsymbol {\cdot } \boldsymbol {\nabla })\boldsymbol {v}$ by 
$(\boldsymbol {U} \boldsymbol {\cdot } \boldsymbol {\nabla } )\boldsymbol {v}$, where 
$\boldsymbol {U}$ is the velocity at infinity. In our case, the velocity vanishes at infinity, so we do not take into account the nonlinear terms in the momentum equation (this is also true with respect to the Forchheimer term). In our case, the temperature gradient at infinity is non-zero, so the quasinonlinear term is present in the energy equation.
A similar approach was used in Balasubramaniam & Dill (Reference Balasubramaniam and Dill1991) for the investigation of the thermocapillary bubble migration.
We start with the case of the inclusion of high thermal conductivity. Then, we have the equations
$$\begin{gather} \Delta \psi = Ra\frac{{\partial T}}{{\partial x}}, \end{gather}$$
$$\begin{gather}\Delta T = \frac{{\partial \psi }}{{\partial x}}, \end{gather}$$with the previous boundary conditions.
It is convenient to introduce the temperature deviation from the distribution at infinity (
$T=y+\vartheta$). Equations (4.1)–(4.2) retain their form (though 
$T$ is replaced by 
$\vartheta$), and the boundary conditions for 
$\vartheta$ are
$$\begin{gather} r = 1: \quad \vartheta ={-}y, \end{gather}$$
$$\begin{gather}r \to \infty : \quad \vartheta = 0. \end{gather}$$We arrive at the linear inhomogeneous problem, whose solutions have certain symmetry
\begin{equation} \left.\begin{array}{c@{}}
\psi ( {x, - y} ) ={-} \psi ( {x,y} ),\quad \vartheta ( {x,
- y} ) ={-} \vartheta ( {x,y} ) \quad \theta ( {x, - y} )
={-}\theta ( {x,y} ), \\ \psi ( { - x,y} ) ={-} \psi (
{x,y} ),\quad \vartheta ( { - x,y} ) = \vartheta ( {x,y} ).
\end{array}\right\} \end{equation}
We seek a solution to the problem (4.1)–(4.4) by renormalizing the stream function
\begin{equation} \psi \to \rho \psi ,\quad \rho = \sqrt{Ra} \end{equation}(Rayleigh number is positive for heating from above). Then the equations take the form
$$\begin{gather} \Delta \psi = \rho \frac{{\partial \vartheta}}{{\partial x}}, \end{gather}$$
$$\begin{gather}\Delta \vartheta = \rho \frac{{\partial \psi }}{{\partial x}}. \end{gather}$$We can simplify equations (4.7)–(4.8). For that, we introduce two auxiliary functions,
$$\begin{gather} \varPhi = (\psi + \vartheta )/2, \end{gather}$$
$$\begin{gather}\varPsi = (\psi - \vartheta )/2. \end{gather}$$
For functions 
$\varPhi$ and 
$\varPsi$, the decoupled equations are obtained:
$$\begin{gather} \Delta \varPhi = \rho \frac{{\partial \varPhi }}{{\partial x}}, \end{gather}$$
$$\begin{gather}\Delta \varPsi ={-} \rho \frac{{\partial \varPsi }}{{\partial x}}. \end{gather}$$
Eliminate the first derivative in (4.11) by representing 
$\varPhi$ in the form 
$\varPhi = \textrm {e}^{\alpha x} F$. For the function 
$F$ we obtain the equation
\begin{equation} \Delta F + 2\alpha \frac{{\partial F}}{{\partial x}} + \alpha ^2 F = \rho \frac{{\partial F}}{{\partial x}} + \alpha \rho F. \end{equation}
It is seen that, at 
$\alpha = \rho /2$, the terms with the first derivative cancel each other, and the equation for 
$F$ takes the form
\begin{equation} \Delta F = \tfrac{1}{4}\rho ^2 F. \end{equation}
Similarly, representing 
$\varPsi$ in the form 
$\varPsi = \textrm {e}^{ - \alpha x} G$ with the same 
$\alpha$, we arrive at the same equation for 
$G$ as for 
$F$:
\begin{equation} \Delta G = \tfrac{1}{4}\rho ^2 G. \end{equation}From relations (4.9)–(4.10), we get
$$\begin{gather} \psi = \varPhi + \varPsi = {\rm e}^{\alpha x} F + {\rm e}^{ - \alpha x} G, \end{gather}$$
$$\begin{gather}\vartheta = \varPhi - \varPsi = {\rm e}^{\alpha x} F - {\rm e}^{ - \alpha x} G. \end{gather}$$
We can rewrite the boundary conditions (3.2)–(3.4) for 
$r=1$ in the form
\begin{equation} F ={-} \tfrac{1}{2} {\rm e}^{-\alpha \cos{\varphi}} \sin{\varphi}, \quad G ={-} \tfrac{1}{2} {\rm e}^{\alpha \cos{\varphi}} \sin{\varphi}. \end{equation}
At the infinity (
$r \to \infty$) we have the boundary conditions
\begin{equation} F \,{\rm e}^{\alpha \cos{\varphi}} \to 0, \quad G\,{\rm e}^{-\alpha \cos{\varphi}} \to 0. \end{equation} The functions 
$F$ and 
$G$ are odd in 
$y$. For symmetry with respect to 
$x$, we have the relations
$$\begin{gather} {\rm e}^{\alpha x} F(x) + {\rm e}^{ - \alpha x} G(x) ={-} {\rm e}^{ - \alpha x} F( - x) - {\rm e}^{\alpha x} G( - x), \end{gather}$$
$$\begin{gather}{\rm e}^{\alpha x} F(x) - {\rm e}^{ - \alpha x} G(x) = {\rm e}^{ - \alpha x} F( - x) - {\rm e}^{\alpha x} G( - x). \end{gather}$$Subtracting (4.21) from (4.20) yields
\begin{equation} G(x) ={-} F( - x). \end{equation}Thus,
$$\begin{gather} \psi = {\rm e}^{\alpha x} F(x) - {\rm e}^{ - \alpha x} F( - x), \end{gather}$$
$$\begin{gather}\vartheta = {\rm e}^{\alpha x} F(x) + {\rm e}^{ - \alpha x} F( - x). \end{gather}$$
Now we need to find a function 
$F$ satisfying the equation of membrane type. The eigenfunctions of this equation are written in the form
\begin{equation} F_n = f_n (r)\sin n\varphi, \end{equation}
where the angle 
$\varphi$ is counted from the axis 
$x$. For 
$f_n$, the following equation is obtained:
\begin{equation} f''_n + \frac{1}{r}f'_n - \frac{{n^2 }}{{r^2 }}f_n - \frac{{\rho ^2 }}{4}f_n = 0. \end{equation}
By introducing a new variable 
$\xi = \rho r/2$, we rewrite the (4.26) in the form
\begin{equation} f''_n + \frac{1}{\xi }f'_n - \left( {\frac{{n^2 }}{{\xi ^2 }} + 1} \right)f_n = 0. \end{equation}
The solution of (4.27) involves the modified Bessel functions of the second kind 
$K_n(\xi )$. Thus, we have
\begin{equation} F = \sum_{n = 1}^\infty {C_n } K_n (\xi )\sin n\varphi. \end{equation}Computations yield
\begin{equation} C_n = n\frac{{( - 1)^{n } }}{{\alpha K_n (\alpha )}}I_n (\alpha ), \end{equation}
where 
$I_n (\alpha )$ are the modified Bessel functions of the first kind.
For large values of the argument, the Bessel functions of the second kind have the asymptotics
\begin{equation} K_n (\xi ) = \sqrt {\frac{{\rm \pi} }{{2\xi }}} \,{\rm e}^{ - \xi } \left( {1 + \frac{{4n^2 - 1}}{{8\xi }} + \cdots} \right). \end{equation}
This means that for large 
$r$
\begin{equation} F \sim {\rm e}^{-({\rho}/{2})r}. \end{equation}
Therefore, the functions 
$\psi$ and 
$\vartheta$ have the asymptotics (for 
$x > 0$)
\begin{equation} \psi ,\vartheta \sim {\rm e}^{({\rho}/{2})r(\cos \varphi - 1)}, \end{equation}
i.e. for any directions, except the ones close to horizontal, 
$\psi$ and 
$\vartheta$ decay exponentially and, for a nearly horizontal direction, the power-law asymptotic behaviour takes place.
Let us find out a more precise asymptotic formula for the horizontal velocity 
$v_x$ at 
$y = 0$, i.e. at 
$\varphi = 0$. Since 
$v_x = - {{\partial \psi }}/{{\partial y}}$, for 
$y = 0$, we have
\begin{equation} v_x |_{y = 0} ={-} {\rm e}^{\alpha x} \left. {\frac{{\partial F(x)}}{{\partial y}}} \right|_{y = 0} + {\rm e}^{ - \alpha x} \left. {\frac{{\partial F( - x)}}{{\partial y}}} \right|_{y = 0}. \end{equation}Since
\begin{equation} \frac{\partial }{{\partial y}} = \sin \varphi \frac{\partial }{{\partial r}} + \frac{1}{r}\cos \varphi \frac{\partial }{{\partial \varphi }}, \end{equation}then
\begin{equation} \left. {\frac{{\partial F}}{{\partial y}}} \right|_{y = 0} = \sum_{n = 1}^\infty {\frac{n}{x}C_n } K_n (\xi ), \end{equation}
where now 
$\xi = \rho x/2$.
Substituting (4.35) into (4.33) leads to
\begin{equation} v ={-} \frac{1}{{x^{3/2} }}\sqrt {\frac{{\rm \pi} }{\rho }} \sum_{n = 1}^\infty {nC_n } + \cdots. \end{equation}Thus, for the velocity, we got the power-law asymptotics with a scaling exponent 3/2.
The flow has the form of horizontal vortices directed along the 
$x$ axis, and the jet velocity decreases according to the law 
$x^{ - 3/2}$.
Thus, we have obtained the asymptotic representation of the solutions of (4.14), satisfying the boundary conditions and for small Rayleigh number, uniformly valid at any distances from the inclusion. The analysis of the structure of the solution allows us to conclude that, even at low values of the Rayleigh number, the convective flow at the distances from the inclusion larger than the inclusion size has the form of horizontal vortices directed away from the inclusion. As an example, the streamlines of stationary flows at 
$Ra = 10$ and 
$Ra=100$ are presented in figures 2 and 3.
Streamlines of the stationary flow (a) and isolines of temperature deviation from the linear profile (b) found in the framework of the Oseen approximation for 
$Ra = 10$, 
$\kappa \to \infty$.
Streamlines of the stationary flow (a) and isolines of temperature deviation from the linear profile (b) found in the framework of the Oseen approximation for 
$Ra = 100$, 
$\kappa \to \infty$.
The distributions of the radial component of velocity in the radial direction at fixed 
$\varphi$ and in the azimuthal direction at fixed 
$r$ are shown in figure 4. It can be seen that with an increase in the Rayleigh number, the maximum value of the radial component of the fluid velocity (and the stream function) increases, and the thickness of the vortices decreases.
Distributions of the radial component of velocity in the radial direction at 
$\varphi =0$ (a) and in the azimuthal direction at 
$r=5$ (b) (
$\kappa \to \infty$, 
$Ra=1$, 10, 100).
In figure 5 we compare the distributions of the radial velocity component in the radial direction at 
$\varphi =0$, 
$\kappa \to \infty$ for 
$Ra=0.01$ and 
$Ra=1$ obtained in the Oseen- and Stokes-type approximations. There is good agreement for 
$Ra=0.01$ and some worse for 
$Ra=1$. This means that the Oseen and Stokes approximations give close results only for sufficiently small values of the Rayleigh number.
Distributions of the radial component of velocity in the radial direction at 
$\varphi =0$, 
$\kappa \to \infty$ for 
$Ra=0.01$ (a) and 
$Ra=1$ (b) obtained in the Stokes-like (i) and Oseen-like (ii) approximations.
4.2. Case of the inclusion of finite thermal conductivity
Now consider the case for finite 
$\kappa$. In this case the equation for 
$\tilde T$ is additionally required. The temperature deviation from the distribution at infinity (
$\tilde T=y+\theta$) is also used. Equation (2.7) retains its form (though 
$\tilde T$ is replaced by 
$\theta$). The thermal boundary conditions in this case are more complex and can be written as follows:
$$\begin{gather} r = 1: \quad \vartheta = \theta, \quad \frac{\partial \vartheta}{\partial r} - \kappa \frac{\partial \theta}{\partial r}=(\kappa-1) \sin{\varphi}, \end{gather}$$
$$\begin{gather}r \to \infty: \quad \vartheta = 0. \end{gather}$$ In terms of 
$F$ and 
$G$, the boundary conditions have the form
$$\begin{gather} F = \frac{1}{2} {\rm e}^{-\alpha \cos{\varphi}} \theta, \quad G ={-} \frac{1}{2} {\rm e}^{\alpha \cos{\varphi}} \theta, \end{gather}$$
$$\begin{gather}\frac{\partial F}{\partial r} - {\rm e}^{- 2 \alpha \cos{\varphi}} \frac{\partial G}{\partial r} -\kappa \,{\rm e}^{- \alpha \cos{\varphi}} \frac{\partial \theta}{\partial r}= (\kappa-1) \sin{\varphi} . \end{gather}$$
At infinity (
$r \to \infty$) we have the boundary conditions
\begin{equation} F \,{\rm e}^{\alpha \cos{\varphi}} \to 0, \quad G \,{\rm e}^{-\alpha \cos{\varphi}} \to 0. \end{equation}
Now we write the solution for 
$\theta$ in the form
\begin{equation} \theta = \sum_{k = 1}^\infty {E_k } r^{2k-1} \sin{(2k-1)\varphi}. \end{equation}
The coefficients 
$C_n$ and 
$E_k$ should be found from the boundary conditions. Using the boundary condition (4.39), we find the following expressions for the coefficients 
$C_n$:
\begin{equation} C_n =\sum_{k=1}^{\infty} c_{nk} E_k, \quad c_{nk}=\frac{({-}1)^{n+1}}{2 K_{n}(\alpha)} [I_{|n-2k+1|}(\alpha) - I_{n+2k-1}(\alpha)]. \end{equation}Instead of the boundary condition (4.40) we can write
\begin{equation} \sum_{k=1}^{\infty} a_{pk} E_k=b_p, \end{equation}where
\begin{gather} \begin{aligned}a_{pk}&={-} \frac{\alpha}{2} \sum_{l=1}^{\infty} \frac{({-}1)^l}{K_{l}(\alpha)}[K_{l-1}(\alpha)+K_{l+1}(\alpha)] [I_{|l-2k+1|}(\alpha)-I_{l+2k-1}(\alpha)] \nonumber\\ &\quad \times [I_{|l-2p+1|}(\alpha)-I_{l+2p-1}(\alpha)] - (2p-1)\kappa \delta_{pk},\end{aligned} \end{gather}
\begin{gather} b_p = (\kappa-1) \delta_{1p}. \end{gather}
This boundary condition can be used for finding of 
$E_p$, but previously we should cut all sums. Then, the boundary condition (4.40) will be
\begin{equation} \sum_{k=1}^{N} a_{pk} E_k=b_p, \end{equation}where
\begin{align} a_{pk}&={-} \frac{\alpha}{2} \sum_{l=1}^{2N} \frac{({-}1)^l}{K_{l}(\alpha)}[K_{l-1}(\alpha)+K_{l+1}(\alpha)] [I_{|l-2k+1|}(\alpha)-I_{l+2k-1}(\alpha)] \nonumber\\ &\quad \times [I_{|l-2p+1|}(\alpha)-I_{l+2p-1}(\alpha)] - (2p-1)\kappa \delta_{p k}, \end{align}
\begin{equation} b_p = (\kappa-1) \delta_{1p}. \end{equation}
We also cut sums in the expressions for 
$F$ and 
$\theta$:
\begin{equation} F = \sum_{n = 1}^{2N} {C_n } K_n (\xi )\sin n\varphi, \quad \theta = \sum_{k = 1}^N {E_k } r^{2k-1} \sin{(2k-1)\varphi}. \end{equation}Now we can write the solution for (4.47) in the form
\begin{equation} \boldsymbol{E}= \boldsymbol{\mathsf{A}}^{{-}1} \boldsymbol{B}, \end{equation}
where we used the matrix notations
\begin{equation}
\boldsymbol{E}=\left(\begin{array}{@{}c@{}} E_1 \\ E_2 \\ \vdots \\ E_N
\end{array}\right),\quad \boldsymbol{B}=\left(\begin{array}{@{}c@{}} b_1 \\
b_2 \\ \vdots \\ b_N \end{array}\right),\quad
\boldsymbol{\mathsf{A}}=\left(\begin{array}{@{}cccc@{}} a_{11} & a_{12} & \ldots a_{1N}
\\ a_{21} & a_{22} & \ldots a_{2n} \\ \vdots & \vdots &
\ddots \vdots \\ a_{N1} & a_{N2} & \ldots a_{NN}
\end{array}\right).
\end{equation}
The wxMaxima package was used for finding the solution.
The results for the specific values of parameters are presented in figures 6 and 7 (
$N=20$). It is seen that, similar to the case of high thermal conductivity of the inclusion, with the increase of the Rayleigh number, the flow gradually takes the form of horizontal vortices, in which the fluid moves in the horizontal plane away from the inclusion. With a further increase in the Rayleigh number, the maximum value of the stream function grows and the thickness of the vortices decreases. Isolines of temperature deviation from the linear profile with increasing Rayleigh number form less pronounced horizontal structures (thermal plumes) than fluid flows. The temperature deviation decreases with increasing distance from the inclusion faster than the stream function.
Streamlines of the stationary flow (a) and isolines of temperature deviation from the linear profile (both in fluid and solid inclusion) (b) obtained in the framework of the Oseen approximation at 
$Ra = 1$, 
$\kappa =1.5$.
Streamlines of the stationary flow (a) and isolines of temperature deviation from the linear profile (both in fluid and solid inclusion) (b) obtained in the framework of the Oseen approximation at 
$Ra = 10$, 
$\kappa =1.5$.
In figure 8 we compare the distributions of the radial component of velocity in the radial direction at fixed 
$\varphi$ and in the azimuthal direction at fixed 
$r$ for 
$\kappa \to \infty$, 
$\kappa =5$ and 
$\kappa =1.5$. It can be seen that as 
$\kappa$ decreases (but 
$\kappa > 1$), the maximum value of the radial component of the fluid velocity decreases.
Distributions of the radial component of velocity in the radial direction at 
$\varphi = 0$ (a) and in the azimuthal direction at 
$r = 5$ (b) for 
$Ra=1$, 
$\kappa \to \infty$ (i), 
$\kappa = 5$ (ii), 
$\kappa = 1.5$ (iii).
At 
$\kappa =1$ there is no flow, since in this case the conductivity of the inclusion is equal to the conductivity of the fluid and the basic temperature profile is not disturbed.
The case of low thermal conductivity (
$\kappa \to 0$) can be formally obtained from the solution for finite value of 
$\kappa$ by setting 
$\kappa = 0$. The results for this case can be seen in figures 9 and 10. For comparison the results for 
$\kappa =0.5$ are shown in figure 11. It can be seen that for 
$\kappa <1$ (then the thermal conductivity of the solid inclusion is less than the thermal conductivity of the fluid) the direction of circulation changes to the opposite (in comparison with the case of high thermal conductivity). The remaining results are similar to those discussed above.
Streamlines of the stationary flow (a) and isolines of temperature deviation from the linear profile (b) obtained in the framework of the Oseen approximation at 
$\kappa =0$, 
$Ra=10$.
Distributions of the radial component of velocity in the radial direction at 
$\varphi =0$ (a) and in the azimuthal direction at 
$r=5$ (b) for 
$\kappa =0$, 
$Ra=1$, 10, 100.
Distributions of the radial component of velocity in the radial direction at 
$\varphi =0$ (a) and in the azimuthal direction at 
$r=5$ (b) for 
$\kappa =0.5$, 
$Ra=1$, 10, 100.
5. Numerical study for finite values of parameters
For finite values of the Rayleigh number and finite size of the container, the problem (2.5)–(2.6) was solved numerically using the finite difference method. Computations were carried out for a cylindrical inclusion, surrounded by a porous medium saturated with the fluid filling the cylindrical container with a radius larger than the inclusion radius (the ratio of the container radius to the radius of the inclusion was varied in the range 1–30). The calculations were carried out in the quarter of cylindrical gap using the symmetry conditions at 
$\varphi = 0$ and 
$\varphi = {\rm \pi}/2$.
Two home-made numerical codes were used in the calculations. The first code was based on an unsteady approach with the explicit numerical scheme for time integration. In this case, the calculations were carried out until the stationary values of the integral characteristics of the stream function and temperature fields were attained with relative accuracy 
$10^{-5}$. The spatial discretization of the derivatives was made with the second order. The Poisson equation for the stream function was solved by the iteration method. Most of the calculations were performed using the grid with 150 nodes in the radial direction and 150 nodes in the azimuthal direction.
At large enough values of the Rayleigh number (
$Ra>200$), thin boundary layers are formed near the rigid boundaries and near the symmetry axes. In this case, numerical calculations were performed using a home-made code (Khlybov Reference Khlybov2022) which is aimed at solving steady state and transient problems on two- and three-dimensional rectilinear grids employing semiautomatic (guided) and automatic discretization of partial differential equations using finite difference and finite volume methods, which transform original differential–algebraic equations into corresponding discrete forms. Implemented samplers perform sampling on Cartesian grid with unit spatial step; the mapping of a physical computational domain to a computational one is carried out with a user-specified coordinate transformation, which can either be specified analytically or be a result of the operation of an external mesh generator. In the present paper the following coordinate transformations were used:
$$\begin{gather} r = 1 + ({R_{ext}} - 1)\bar \xi ,\quad \bar \xi = \frac{1}{2}\left\{ {1 + \frac{{\tanh\left( {{\delta _r}\left[ {\dfrac{{\xi - 1}}{{{R_{ext}} - 1}} - \dfrac{1}{2}} \right]} \right)}}{{\tanh({\delta _r}/2)}}} \right\}, \end{gather}$$
$$\begin{gather}\varphi = \frac{{\rm \pi} }{2}\bar \eta ,\quad \bar \eta = \frac{1}{2}\left\{ {1 + \frac{{\tanh\left( {{\delta _\eta }\left[ {\dfrac{{2\eta }}{{\rm \pi} } - \dfrac{1}{2}} \right]} \right)}}{{\tanh({\delta _\eta }/2)}}} \right\}. \end{gather}$$ An arbitrary partition of a discretized system of nonlinear algebraic equations into subsystems is allowed to obtain various solution schemes ranging from the fully implicit strategy (where all unknown fields are solved simultaneously) to the segregated one (for each field a specific solution procedure is generated, there the remaining fields included in equations are not considered a part of solution). Discretized systems of (non)linear algebraic equations are solved by Newton's method using third-party SLAE solvers; one can choose between sequential and parallel solution methods. In the latter case, it is possible to use multiprocessor systems with shared memory or parallel systems (clusters), or even in heterogeneous environments. The partial differential equation system under discussion was discretized with a second-order finite-difference method with analytical coordinate transformation mapping of the axisymmetric physical domain with non-uniform grid into the uniform Cartesian computational grid. The employed coordinate transformation generates a finer grid near the domain borders at the expense of the bulk of the domain in order to attain boundary layer resolution. The obtained system of nonlinear algebraic equations is solved with the fully implicit solution strategy with nonlinearity and is treated by the Newton's method, and the underlying linear system solved by the MUMPS direct sparse SLAE solver. Most of the calculations were performed using the grid with 300 nodes in the radial direction and 300 nodes in the azimuthal direction. Parameters 
$\delta _r$ and 
$\delta _\eta$ were taken equal to 5. For this mesh at 
$Ra=1000$, 
$\kappa =0$ there are 27 nodes in the boundary layer in the radial direction at 
$\phi =0$ near the inclusion surface, and 111 nodes in the boundary layer in the azimuthal direction at 
$r=5$ near symmetry line 
$\phi =0$. Test calculations using the grids 
$100\times 100$, 
$200\times 200$, 
$400\times 400$ and 
$500\times 500$, confirmed very good accuracy of calculations with grid 
$300\times 300$ (see table 1).
Dependence of the maximum value of the stream function on the number of nodes.
Figure 12 shows how the radial size of the numerical domain affects the accuracy of the results. It can be seen that increasing the size of the numerical domain has virtually no effect on the accuracy of the results.
Distributions of radial component of velocity in the radial direction at 
$\varphi =0$ (a) and in the azimuthal direction at 
$r=5$ (b) for 
$\kappa \to \infty$, 
$Ra=1000$, 
$r_{ext} = 15$ for dashed line, 
$r_{ext} = 30$ for solid line (grid 
$500\times 500$).
5.1. Case of the inclusion of high thermal conductivity
Let us first consider the case of the inclusion of high thermal conductivity. In this case, at the container wall, we imposed the temperature distribution corresponding to the uniform vertical temperature gradient (heating from above).
The calculations show that at low Rayleigh numbers there is a single vortex occupying the whole computational domain (remember that the computational domain is a quarter of the whole cylindrical gap) and symmetric about the diagonal (figure 13a).
Streamlines of the stationary flow calculated numerically using the full nonlinear approach at 
$Ra= 1$ (a), 
$Ra= 100$ (b) and 
$Ra= 1000$ (c).
With the increase in the Rayleigh number, the flow gradually takes the form of horizontal vortices, in which the fluid moves in the horizontal plane away from the inclusion (figure 13b). With a further increase in the Rayleigh number, the maximum value of the stream function grows and the thickness of the horizontal vortices decreases (figure 13c).
This behaviour is also illustrated in figure 14, where the distributions of the radial component of velocity in the azimuthal direction at fixed 
$r$ and in the radial direction at fixed 
$\varphi$ are presented.
Distributions of the radial component of velocity in the radial direction at 
$\varphi =0$ (a) and in the azimuthal direction at 
$r=5$ (b) obtained numerically for the case of inclusion of high thermal conductivity.
5.2. Case of the inclusion of low thermal conductivity
Numerical solution was carried out not only for the high thermal conductivity of the inclusion, but also in the case of the inclusion whose thermal conductivity is much less than the thermal conductivity of the surrounding liquid. The following thermal boundary condition at the inclusion surface was used:
\begin{equation} \frac{\partial T}{\partial r}=0. \end{equation}It follows from the calculations that in this case the sign of the convective circulation changes, so that in the horizontal plane passing through the axis of the cylinder the fluid is leaking to the inclusion (figure 15a–c).
Streamlines of the stationary flow obtained numerically using the full nonlinear approach in the case of inclusion of low thermal conductivity at 
$Ra = 1$ (a), 
$Ra = 100$ (b) and 
$Ra = 1000$ (c).
In figure 16, the distributions of the radial component of velocity in the radial direction at fixed 
$\varphi$ and in the azimuthal direction at fixed 
$r$ in the case of low thermal conductivity of the inclusion are presented.
Distributions of the radial component of velocity in the radial direction at 
$\varphi = 0$ (a) and in the azimuthal direction at 
$r = 5$ (b) obtained numerically for the case of inclusion of low thermal conductivity.
The evolution of the maximal value (in absolute value) of the stream function on the Rayleigh number for the limit cases of high and low thermal conductivities is presented in figure 17(a,b). It can be seen that for the case of low conductivity, the sign of the stream function differs from the sign for the case of high thermal conductivity.
Dependencies of the maximal value (in absolute value) of the stream function on the Rayleigh number obtained numerically for high (a) and low (b) thermal conductivities.
6. Comparison of analytical and numerical results
Figures 18, 19, and figures 20, 21, show the comparison of the velocity profiles in radial and azimuthal directions obtained at different Rayleigh number values for high and low thermal conductivities of the inclusion in the framework of the Oseen approximation and numerically in the framework of the full nonlinear approach. It is seen that there exists not only qualitative but also good quantitative agreement for 
$Ra=100$.
Distributions of the radial component of velocity in the radial direction at 
$\varphi =0$, 
$\kappa \to \infty$ obtained analytically in the framework of the Oseen approximation (i) and numerically using full nonlinear approach (ii) for 
$Ra=100$ (a), 
$Ra=1000$ (b).
Distributions of the radial component of velocity in the azimuthal direction at 
$r=5$, 
$\kappa \to \infty$ obtained analytically in the framework of the Oseen approximation (i) and numerically using full nonlinear approach (ii) for 
$Ra=100$ (a), 
$Ra=1000$ (b).
Distributions of radial component of velocity in the radial direction at 
$\varphi =0$, 
$\kappa =0$ obtained analytically in the framework of the Oseen approximation (i) and numerically using full nonlinear approach (ii) for 
$Ra=100$ (a), 
$Ra=1000$ (b).
Distributions of radial component of velocity in the azimuthal direction at 
$r=5$, 
$\kappa =0$ obtained analytically in the framework of the Oseen approximation (i) and numerically using full nonlinear approach (ii) for 
$Ra=100$ (a), 
$Ra=1000$ (b).
Thus, the analytical data computed with Oseen quasilinearization describe well the arising flow within a wide enough range of the Rayleigh number values. At higher values of the Rayleigh number the difference in analytical and numerical results becomes more substantial (see figures 18b and 21b). This shows the validity range of the Oseen approximation.
7. Conclusions
Formation of flows in stably stratified fluid at weak violation of stable stratification has been studied using the problem on thermal buoyancy convection near a solid inclusion embedded in a fluid-saturated porous medium heated from above as an example. The simple form of the equations of thermal buoyancy convection in a fluid-saturated porous medium in the Darcy–Boussinesq approximation allowed for a fairly complete analytical study of the problem. The equations are written in terms of the stream function and temperature deviations from a linear vertical distribution and possessed the internal symmetry, which allowed us to reduce the system of equations to one equation without increasing the order.
An analytical solution of the problem has been obtained for the limit case of large thermal conductivity of the inclusion, when the ratio 
$\kappa = \kappa _s / \kappa _f$ tends to infinity, in the framework of the Oseen-like approximation where, different from the conventional Oseen approximation, the quasilinearization is applied to the nonlinear terms in the energy equation and not in the momentum equation. This solution satisfies the given boundary conditions and is uniformly valid at any distance from the inclusion.
Analysis of this solution shows that, even at low values of the Rayleigh number, the convective flow at a distance from the inclusion greater than the inclusion size has the form of four horizontal vortices; besides, in a horizontal plane passing through the inclusion axis the fluid moves away from the inclusion. Thus, the formation of the horizontal structures in the laminar flow of a stably stratified fluid has been discovered and described for the first time. The flow velocity in a horizontal plane passing through the inclusion axis decreases with the distance from the inclusion as 
$x^{-3/2}$.
For finite values of the thermal conductivity ratio 
$\kappa$, a semianalytical solution has been obtained within the framework of Oseen-like approximation with the truncation of a system of linear equations. Analysis of this solution shows that the convective flow in the form of horizontal vortices arises at any values of 
$\kappa$ except for the case 
$\kappa = 1$ when the inclusion thermal conductivity is the same as that of the fluid and the conditions of stable stratification are not violated. However, the directions of the flow at 
$\kappa > 1$ and 
$\kappa < 1$ are different: in the first case the fluid moves in a horizontal plane passing through the inclusion axis away from the inclusion and in the second case – towards the inclusion. The flow velocity magnitude increases with the increase of deviation of 
$\kappa$ from unity. The results for the limit case of low conductivity of the inclusion are directly obtained from semianalytical solution obtained for finite 
$\kappa$.
Analytical results obtained in the framework of the Oseen-like approximation have been confirmed by the comparison with the numerical results on the velocity profiles and flow structure obtained in the framework of full nonlinear approach by finite difference method.
Good qualitative and quantitative agreement for the Rayleigh number values 
$Ra < 100$ identifies the range of the validity of the Oseen-like approximation.
The phenomena considered in the present paper have general character. Particularly, similar phenomena should be observed in the solutal convection. A general feature of these phenomena is that weak violation of stable stratification in stably stratified fluid leads to the emergence of laminar convective flows in the form of horizontal structures with the horizontal component of the velocity significantly exceeding the vertical one. The reason for the arising of these flows is that stable stratification prevents vertical flows associated with the work against gravity, but does not prevent horizontal flows.
Acknowledgements
Authors thank O. Khlybov for the numerical code implementation.
Funding
The work was carried out as part of a major scientific project funded by the Ministry of Science and Higher Education of the Russian Federation (agreement no. 075-15-2024-535 dated 23 April 2024).
Declaration of interests
The authors report no conflict of interest.
