Direct numerical simulation of “fountain filling box” flow with a confined weak laminar plane fountain

A “fountain filling box” flow produced by discharging a weak laminar plane fountain in a confined open channel is studied numerically. Two‐dimensional direct numerical simulations were performed for weak plane fountains. The development of the fountain flow experiences five stages; the initial upflow and the subsequent downflow after the fountain penetrates to the maximum height, followed by the outward movement of the intrusion of the fallen fountain fluid on the channel bottom, and then the wall fountain formed by the impingement of the intrusion on the vertical sidewall, which results in the reversed flow, and finally the gradual stratification of the fluid. The behavior of the intrusion can be approximately described with the plane gravity current theory. The period for the intrusion to reach the bounded side wall increases with increasing Re or decreasing Fr. Three regimes are found for the wall fountain behavior; “no‐falling,” “slumping down,” and “rolling down” behavior. Convection, mixing, conduction, and filling all contribute to the formation and development of stratification, but their effects vary at different stages. For the initial stages, convection and mixing play a key role, resulting in an increasing bulk entrainment rate, while conduction and filling are dominant after quasi‐steady stratification is created, presenting a decreasing bulk entrainment rate.

For a typical confined fountain flow, it is seen from Figure 1 that the fountain fluid falls back after the fountain attains its maximum penetration height, subsequently impinges on the base, and then moves outwards as an intrusion. Due to the relatively narrow confinement of the ambient, the intrusion will impinge on the vertical sidewall and move subsequently upwards along the wall to form a plane wall fountain. The wall fountain, which is usually much weaker than the source fountain, will fall back after that due to its larger density to form a reversed flow that moves towards the source, interacting with the intrusion and the ambient fluid on the way. The continuous discharge of the fountain fluid from the source will gradually generate a density stratification of fluid in the confined space. Baines et al. 14,15 denoted such a confined fountain flow as the "fountain filling box model," which is also adopted in the present study.
The characteristics of the "fountain filling box" model share some similarities with those of the "plume filling box" model 15 with overturning, as studied by Kaye and Hunt,15 particularly the intrusion and the wall fountain.
The "plume filling box model," produced by a plume in a confined region, has received considerable attention because of its fundamental and application significance (e.g., a fire plume in buildings). [16][17][18][19] The research has been focused on the behavior of turbulent plumes and the time evolution of density stratification. It has been found that a turbulent plume in a confined space can cause strong flows such as shear flow (i.e., intrusion) resulting from the plume outflow and the overturning structure (i.e., wall fountain) created by the impingement of the intrusion on side boundaries, in the early development stage. 16 Kaye and Hunt 16 developed a theoretical model for an axisymmetric constant-buoyancy-flux turbulent plume in confined cylinders, in which the outflow from the plume was treated as a forced gravity current with constant buoyancy flux, whereas the flow along the sidewall after the impingement of the gravity current is modeled as a wall fountain. Two regimes are identified in terms of the aspect ratio H/W, where W is the cylinder radius and H is its height. For a room with a larger aspect ratio (H/W > 1.5), the intrusion is a pure gravity current when impinging on the sidewall, and the rise height of the wall fountain only depends on H. However, for H/W ≤ 1.5, the intrusion is not fully developed into pure gravity before the impingement, resulting in the wall fountain height depending on both H and W. For H/W > 4.0, Barnett 20 found that the upward plume was prevented from impinging on the ceiling due to the downflow in the ambient, and named these "blocked" regimes. Apparently, the secondary flows can in turn affect the behavior of the plumes. In an experimental study on the wall fountain created by a ceiling jet turning downward at the corner of a compartment, Jaluria and Kapoor 21 found the penetration depth of the wall fountain only depends on Fr of the outflowing buoyant jet, with a fixed distance between the jet source and the corner. Therefore, the secondary flows and the development of stratification are significantly influenced by the source conditions and the geometry of the confinement.
Due to the similarity to the "plume filling box" model, it is reasonable to expect that the intrusion, wall fountain, and time-dependent density stratification in the "fountain filling box" flow are also influenced by the source condition (via Fr and Re) and the confinement conditions. In their analysis of "fountain filling box" flows with turbulent confined fountains, Baines et al. 14 considered the entrainment but excluded the confinement size. There is a similar knowledge gap for plumes impinging on a density interface. To our best knowledge, only Shrinivas and Hunt 22 took into account the confinement influence on the entrainment of the fountain formed by the impinging plume. In their study, the confinement parameter λ i is characterized by the ratio of the interfacial turbulence length scale to the depth of the upper layer of the two-layer ambient. When λ i is small, a weak secondary flow does not influence the entrainment significantly, with the scaling relation E i ∝ Fr i 2 , where E i and Fr i represent the dimensionless entrainment flux over the interface and the Froude number on the interface respectively, while for large λ i , E i ∝ Fr i 3 .
Recently, there have been some studies focusing on the behavior of fountains under confinement. Debugne and Hunt 23 investigated the impact of the spanwise confinement on turbulent round fountains and identified four flow regimes over 0.5 ≤ Fr ≤ 96 and 2 ≤ λ ≤ 24, where λ = W/R 0 is the dimensionless confinement length, with W the dimensional confinement width and R 0 is the source radius. To account for the effects of confinement, a "confined" Froude number Fr c ≡ Frλ -5/4 was introduced as the governing parameter for confined fountains. Notably, the fountain is only spanwise confined, with no confinement in the lateral direction. Xue et al. 24 characterized the development of round fountains in a bounded container at the initial stage, in terms of Re and Fr over Re < 500 and 5 ≤ Fr ≤ 35. Their results show that for a given Fr, the fountain volume entrainment flux ratio attains a local peak at Re ≈ 200, which is also supported by the measurement of the fountain penetration heights. Their experiments showed that the confinement strengthens the horizontal mixing. However, only two square tanks were tested in their experiments, which makes a systematic study on the impact of the weak confinement impossible. Lippert and Woods 25 examined theoretically and experimentally a particle fountain in a confined environment, and identified four regimes for the flow with different source conditions. Very recently, Dong et al. 26 studied experimentally the behavior of confined laminar and turbulent round "fountain filling box" flow over 1.0 ≤ Fr ≤ 20.0, 102 ≤ Re ≤ 1502, and 27.9 ≤ λ ≤ 48.75. We showed that the confinement significantly changes the transient behavior of round fountains in the confined containers, particularly the intrusion, reversed flow, and stratification. We also used three-dimensional direct numerical simulation (DNS) over 0.25 ≤ Fr ≤ 3.0, 5 ≤ Re ≤ 800, and 10 ≤ λ ≤ 35 to analyze the intrusion, reversed flow, and stratification of confined weak round fountains 27 and identified three development stages of the bulk entrainment rate. In the present study, our study 27 is extended to the confined weak laminar plane fountains.
The rest of this paper is structured as follows. In Section 2, the details of the DNS runs are provided, including a brief introduction of the physical system. The governing equations, initial and boundary conditions, numerical solution techniques, construction of computational mesh, and the mesh and time-step independence testing are detailed in the accompanying Supporting Information Materials. The snapshots of the contours of temperature obtained from the numerical results are presented in Section 3 to show the development of the transient flow of typical laminar planar fountains in confined open channels when Fr, Re, and λ vary. Additionally, the influence of these governing parameters on the confined laminar planar fountains is then analyzed and discussed quantitatively in Section 3, including the movement of the intrusion and wall fountain fronts, the development of stratification, as well as the associated bulk entrainment. The conclusions are finally summarized in Section 4.

| NUMERICAL METHODOLOGY
The details of this numerical methodology section are presented in the accompanying Supporting Information Materials. Here only key information is presented. The physical system considered is a two-dimensional (2D) rectangular box of height H and half-width L. The vertical sidewalls of the box are no-slip and adiabatic, and the box top is open. On the center of the bottom floor, a slot of half-width X 0 serves as the source for the 2D plane fountain. The remaining floor area is rigid, no-slip, and adiabatic. Initially, a quiescent homogeneous Newtonian fluid at uniform temperature T a is filled in the box. At time t = 0, a dense jet at temperature T 0 (T 0 < T a ) is ejected upward into the box at a velocity W 0 and this discharge is maintained thereafter.
The flow of the 2D laminar plane fountain in a confined homogeneous environment is governed by the 2D incompressible Navier-Stokes equations together with the temperature equation. These governing equations together with the initial and boundary conditions are presented in the Supporting Information Materials.
The main governing parameters are Re, Fr, and dimensionless length of the rectangular box (λ), which are defined as follows: where ρ 0 and ρ a are the densities of the jet fluid at the source and the ambient fluid, respectively, g is the gravitational acceleration, ρ, β, and ν are the density, volumetric expansion coefficient, and kinematic viscosity of the fluid, respectively. For fountains resulting from the temperature difference between the jet and ambient fluid, Fr can also be calculated with the temperature difference using the Oberbeck-Boussinesq approximation for buoyancy, as shown in the second expression in Equation (2), which requires (ρ 0 − ρ a )/ρ a to be significantly less than unity.
The discretization and integration of the governing equations are described in the Supporting Information Materials. All direct numerical simulation runs were carried out using ANSYS Fluent 17.0.
A total of 52 DNS runs have been carried out over 0.1 ≤ Fr ≤ 3.0, 5 ≤ Re ≤ 800, and 10 ≤ λ ≤ 35, all at Pr = 7, with their main data listed in Table 1.
The meshes are nonuniform, with fine, uniform grids in the bottom region and a relatively coarse vertically stretched mesh in the top region. Extensive grid and time-step independence tests were conducted to ensure accurate simulation results are produced. Additionally, to ensure the DNS results are accurate, the simulation data for the plane fountain of 0.1 ≤ Fr ≤ 2.5 and Re = 200 is benchmarked against the existing numerical and experimental data summarized by Hunt and Coffey. 28 These are presented in the Supporting Inforation Materials. The meshes used for each DNS run are also presented in Table 1 where h is the dimensionless height of the rectangular box, which is nondimensionalized by h = H/X 0 .  where τ is the dimensionless time as presented in Supporting Inforation Materials. After the inception of the fountain flow, the fountain will rise until it attains the maximum height as illustrated in Figure 2A. Since the fountain fluid remains heavier than the ambient fluid, it will then descend to and spread outwards along the bottom, resulting in a layer of denser fluid (intrusion), as illustrated in Figure 2B,C. The intrusion can be treated as a gravity current, with behavior that may be characterized by different regimes in terms of the governing forces, further discussed in Section 3.2.1. Since the evolution of the initial fountain formation and the intrusion was described in Lin and Armfield, 29 the account of these stages is omitted here.
The influence of the bounded sidewall on gravity intrusion becomes significant as the intrusion flow approaches the sidewall. Figure 2D,E show that the circulation above the intrusion head is stretched and spread upwards along the sidewall with the current impinging on the sidewall. Notably, the density of the leaf-like head of the upward flow is significantly smaller than the intrusion current due to the ambient fluid entrained as the intrusion spreads and then impinges on the sidewall. Therefore, the flow could be divided into two regions, that is, the lighter leaf-like region evaluated from the circulation on the top and the heavier flow from the continuous intrusion current at the bottom, as shown in Figure 2E. Since the flow is still heavier than the ambient fluid, the upward flow along the sidewall can be treated as a wall fountain, which keeps moving upwards until a certain height and then slumps back due to the negative buoyancy ( Figure 2F). The wall fountain flow behavior will be further discussed in Section 3.2.2. The fallen flow of the wall fountain interacts with the intrusion, resulting in an increase in the thickness of the dense current close to the sidewall region. Additionally, fingerlike structures are observed to appear and disappear, as shown in Figure 2G,H. Due to negative buoyancy, a reversed flow is formed, moving from the sidewall to the fountain core as shown in Figure 2G-I. The reversed flow interacts with the sidewall, the intrusion, and the fountain core. As a result, the height of the fountain core experiences a significant increase, until it reaches a DONG ET AL.
| certain height, then falls down again due to gravity and creates a stronger intrusion. Figure 2J shows the reversed flow fronts from the sidewall colliding at the center and being projected to a higher position. Due to the negative buoyancy, the flow drops down and separates into the two reversed flow fronts as shown in Figure 2K. This swaying back-and-forth process repeats several times, behaving like a seiche with a decreasing magnitude. After a certain time, a thermal stratification is formed in the domain with the fountain flow submerged, as can be seen in Figure 2L, which grows gradually with time resulting from convection, mixing, thermal conduction, and filling. Particularly, after the fountain is immersed in the stratified fluid, the increase in the stratification height is mainly due to the continual filling of the denser fluid from the source and thermal conduction. (the distance between the blue part and the red part) increases with decreasing Re, which indicates thermal conduction plays a more significant role.

| Intrusion
The left column in Figure 6 presents the time series of the movement of the intrusion front for Fr, Re, and λ over 0.1 ≤ Fr ≤ 3.0, 5 ≤ Re ≤ 800, and 10 ≤ λ ≤ 35, respectively. The intrusion front is defined as the x-location at which T(x) = T a − 1%(T a − T 0 ) within the right half of the computation domain from y = 0 to H. The kink and the endpoints of the time series profile are the instant when the intrusion is created and the instant when the intrusion front impinges on the sidewall, respectively. The period for the intrusion from creation to impinge on the sidewall is denoted as Δτ w , as shown in Figure 6A for the Fr = 1.0 result. The results in Figure 6A,B show that it takes a longer time for the intrusion front with larger Fr or smaller Re to impinge on the sidewall. However, the Re dependency reduces with increasing Re, with the Re = 200 to Re = 800 overlappings, indicating a minimal dependency for Re ≥ 200. Figure 6C presents the intrusion profiles for the cases with various confinement sizes λ. It is seen that the effect of the confinement size becomes noticeable only when the intrusion approaches the sidewall.
The impacts of Fr, Re, and λ on Δτ w are shown in the right column in Figure 6. From Figure 6D, three ranges can be distinguished with two critical numbers at Fr = 1.0 and Fr = 2.25. Three corresponding correlations are obtained as follows: The regression constants for the three correlations are R 2 = 1, 0.998 and 0.974, respectively. Similarly, the influence of Re on Δτ w is presented in Figure 6E. Re = 150 is found to distinguish the range into two parts, that is, for 5 ≤ Re ≤ 150, Δτ w is well approximated by a power law relation, while for Re ≥ 200, the influence of Re on Δτ w is negligible. The corresponding correlations are determined with the numerical results as follows: with a regression constant of R 2 = 0.991 for the power relation. Figure 6F shows that Δτ w has a linear relation with the confinement size of the channel as shown below where the regression constant is R 2 = 0.999. During the development of a pure gravity current, the flow may experience three regimes, that is, the wall jet regime (W-J), the buoyancy-inertial regime (B-I), and the buoyancy-viscosity regime (B-V), which are determined by the dominating forces. 30 For the initial stage, the flow is dominated by momentum as a plane wall jet, so a scaling relation X(t)~M 1/3 t 2/3 is expected, where M is the momentum flux. After that, the driving force of the current becomes buoyancy (gravity) which is balanced by the inertial force, thus the current is in a B-I regime. The balance between the gravity and the inertial force is maintained until the inertial force is small compared to the total viscous drag force resulting from the interfacial shear stress between the current and the ambient fluid and the bottom shear stress, leading to the third regime where the buoyancy force is balanced by the viscous drag force. 30 Two corresponding scaling relations, that is, X(t)~B 1/3 t and X(t)~(BQ 2 /ν) 1/5 t 4/5 , were obtained for the gravity current in the B-I and B-V regimes, respectively, 30 where Q = X 0 W 0 , B = g(ρ 0 − ρ a )/ρ a Q are volume and buoyancy fluxes, respectively. The correlations above may be written in dimensionless form as follows: x τ , for the W-J regime, x Fr τ ,
For a free gravity current, Equation (7) indicates that the intrusion front location of the W-J regime is dependent on time, but independent of both Fr and Re. While for the B-I regime described by Equation (8), the intrusion front presents a linear correlation with time, and only depends on the parameter of Fr. In the B-V regime, the correlation follows a power law with time again and is dependent on both Fr and Re, as shown in Equation (9).
For the fountains with 10 ≤ λ ≤ 35, Δτ w follows a correlation of Δτ w~λ as shown in Equation (6). Equation (8) can be transformed into τ~Fr 2/3 x. For the fountain with the same Fr, the formula can be further written as τ~x, which is consistent with Δτ w~λ . This indicates that for the cases of 10 ≤ λ ≤ 35 with Fr = 1.0 and Re = 200, the intrusion flows all fall into the B-I regimes.

| Wall fountain
As noted above, the bulk behavior of the wall fountain can be characterized as "No-falling" (i.e., No overturning), "Slumping down," and "Rolling down," as summarized in Figure 7. From this figure, it is found that the overturning behavior of the wall fountain is mainly under the influence of Re; for small Re values (5 ≤ Re ≤ 20), no overturning is observed; for 35 ≤ Re ≤ 100, the wall fountain slumps down after reaching the maximum height; and the fountain rolls down for higher Re values (Re ≥ 150). Notably, it is found that the wall fountain for Fr = 1.0, Re = 200, and Pr = 7.0 slumps down when the confinement size reduces to λ = 10, in contrast to the rolling down behavior for its counterpart of λ ≥ 15. The influence of Fr on the overturning behavior is minimal.
To further investigate the influence of these parameters on the wall fountain, the time series of the wall fountain front along the sidewall are presented in Figure 8 (top row). Again, a part of the time series is omitted for clarity. The wall fountain front is defined as the y-location F I G U R E 7 Behaviors of the wall fountain for 0.1≤ Fr ≤ 3.0, 5 ≤ Re ≤ 800 at Pr = 7 and λ = 20: "No falling," "Slumping down," and "Rolling down." where T(y) = T a − 1%(T a − T 0 ) on the sidewall. As shown in Figure 8A, the maximum height of the wall fountain on the sidewall, y m , is the maximum value of y in the time series and the corresponding, Δτ m , is the period taken for the wall fountain to attain y m from the creation of the intrusion flow. It is found that y m increases with Fr and λ but decreases with Re. Similarly, Δτ m is larger when Fr or λ increases, but becomes smaller when Re increases.
The impact of Fr, Re, and λ on Δτ m is shown in Figure 8 (middle row). Similar to Δτ w , the impact of Fr on Δτ m can be divided into three ranges as identified in Figure 8D, with Fr = 1.0 and 2.25 as the critical values, and the following correlations are found  As there is no falling down process for Re ≤ 20, no Δτ m exists for these cases. As shown in Figure 8E, it is seen that the impact of Re on Δτ m can be divided into two different ranges, with Re = 150 as the critical value, although the effects of Re on Δτ m are not significant.
A linear relation is found for the influence of λ on Δτ m , as shown in Figure 8F, which can be approximated by, with the regression constant of R 2 = 0.996. Figure 8G demonstrates that similar to the impact of Fr on Δτ w and Δτ m , Fr = 1.0 and Fr = 2.25 distinguish the impact of Fr on y m into three different ranges. The corresponding correlations for the three individual ranges are obtained as follows: with R 2 = 0.999, 0.995 and 0.978, respectively. From Figure 8H, it is shown that before and after Re = 150 Re has different influences on y m , which is most likely caused by the different overturning behaviors, that is, slumping down for 50 ≤ Re ≤ 100 and rolling down for 150 ≤ Re ≤ 800.
The influence of λ on y m can be approximated by, with R 2 = 0.984, which is similar to the impact of λ on Δτ m .

| Stratification
The height of the thermally stratified fluid within the domain is defined as the vertical location where T(y) = T a − 1%(T a − T 0 ), which is the height of the interface between the stratified fluid produced by the filling of cold fluid through the fountain flow and the ambient fluid. The time series of the maximum, minimum, and average heights of the thermally stratified fluid after the intrusion reaching the sidewall and the establishment of reversed flow are depicted in Figure 9 (left column) for various Fr, Re, and λ. The magnitudes of the differences among these heights are initially significant, due to the key roles played by convection and mixing. After that, over a relatively long time, the differences decrease and the time series of these heights follow essentially the same trend when the filling becomes dominant in the subsequent stratification formation. From this figure, it is seen that the development rate of the quasi-steady stratification can be approximately described by the rate of the averaged stratification height profile. As shown in Figure 9A, when Fr increases, the magnitudes of the differences among the maximum, average and minimum height profiles increase and a longer time is required for the confined fountain to reach a quasi-steady stratification. Similarly, the differences among these heights increase when Re becomes larger, however, the time to reach a quasi-steady stratification is decreased. Additionally, the development rate of the stratification decreases with the increase of Re.
For a pure filling box without shear flow, convection, mixing, and heat conduction between the filled fluid and the ambient fluid, the velocity of the stratification height, v s , is the reciprocal of λ (v s = 1/λ), based on the conservation of mass. The velocity of the averaged stratification height profile and the corresponding pure filling rate is plotted in Figure 9 (right column) to illustrate the impact of Fr, Re, and λ on v s for the fountains filling box. Figure 9D indicates that the impact of Fr on v s , again, can be divided into three ranges by Fr = 1.0 and 2.25 as the critical values, although no suitable correlations are obtained. The influence of Re on v s , as shown in Figure 9E, can be divided into two different regimes, with Re = 200 as the critical value. The results show that, As it can be seen, v s becomes almost independent of Re for Re larger than 150. The relation between v s and λ is shown in Figure 9F, and can be approximated by the following correlation: with the regression constant of R 2 = 1.0.

| Bulk entrainment or dilution
Since no assumptions are required for the interactions between the upflow and downflow, or the entrainment process between fountains and surroundings, the bulk entrainment is selected to estimate the mean dilution of the buoyancy scalar over the fountain as a whole instead of resolving the local entrainment rate. The bulk entrainment rate is defined as η = Q E /Q 0 , where Q E is the bulk entrainment and Q 0 is the source volume flux. In this study, Q E is calculated by Q E = Q s − Q 0 , where Q s represents the volume flux of the stratified fluid. Q s could be obtained by integrating the area under the thermal stratified surface where the temperature is at T (y) = T a − 1%(T a − T 0 ) as defined in Section 3.2.3. However, processing the integration for hundreds of thousands of time steps is an excessive workload. Hence, Q s here is approximated by the product of the x-location and the average thermal stratification height, although a certain error exists for the initial formation of the fountain. The values of η are calculated for fountains over 0.1 ≤ Fr ≤ 3.0, 5 ≤ Re ≤ 800, and 10 ≤ λ ≤ 35, from the intrusion stage to the filling stage. The time series of η, as shown in Figure 10 (top row), can be approximately divided into three stages. First, η increases monotonically until it attains the first peak point which is denoted as η 1 . Similarly, the period from the formation of intrusion flow to the instant to attain η 1 is denoted as Δτ e . This stage corresponds to the intrusion development and the evolution of the wall fountain. During the intrusion moving towards the sidewall, the ambient fluid is entrained mainly by the eddy over the intrusion head, resulting in the monotonic increase of η.
After that, η shows a fluctuating decrease with a series of subpeak points. The increases of Fr, Re, and λ strengthen the fluctuations as well as the subpeak points. For most of the cases in this study, η 1 is also the maximum value of η. However, the first peak point decreases with increasing Fr, leading to the role reversal between the first peak and the subpeaks, that is, the "subpeak" point exceeding the first peak point to become the maximum value of η for the case of Fr = 3.0, as shown in Figure 10A. Additionally, no obvious fluctuation is observed after the first peak point for the cases of 5 ≤ Re ≤ 100 with Fr = 1.0 and λ = 20, as shown in Figure 10B. However, for the counterpart stage of the cases of 150 ≤ Re ≤ 800, the interactions among the intrusion, the reversed flow, and the ambient fluid induce a fluctuation into the time series until a quasi-steady stratification is established.
With the formation of the quasi-steady stratification, the time series of η finally enters a smooth decline stage. In this stage, entrainment is mainly by thermal conduction, resulting in a decrease of η with the increase of Re and the decreases of Fr and λ. Thermal conduction keeps influencing the dilution process at all stages. The impact of Fr, Re, and λ on η is found to be consistent with the results presented in Section 3.2.3 for the stratification rate.
The influence of Fr, Re, and λ on Δτ e and η 1 is presented in Figure 10D-F, respectively. Notably, for all the cases except for the one of Fr = 3.0 in this study, η 1 could be treated as the maximum value of η, therefore the corresponding Δτ e can be also determined as the time-scale for the filling flow to reach its maximum entrainment rate. Similarly, Fr = 1.0 and Fr = 2.25 are determined as the approximate critical values to separate the impact of Fr into three different ranges, with the following correlations obtained,  with R 2 = 1.0, 0.985 and 0.985, respectively. From Figure 10E, Re = 150 is determined as the approximately critical value to divide the impact of Re into two ranges, with a power-law correlation obtained for 5 ≤ Re ≤ 100 as follows:  Figure 10F, the following linear correlation is obtained, with R 2 = 0.987. η 1 is plotted against Fr, Re, and λ in Figure 10G-I, respectively to depict the impact of these parameters on η 1 . The value of η 1 presents a fluctuant decreasing trend with the increase of Fr as shown in Figure 10G. Again, Re = 150 is determined as the approximate critical value to divide the impact of Re into two ranges, as shown in Figure 10H. For 5 ≤ Re ≤ 100, η 1 decreases with the increase of Re. However, a fluctuant increasing η 1 is observed for 150 ≤ Re ≤ 800. The influence of λ on η 1 presented in Figure 10I can be quantified by, with R 2 = 0.987.

| CONCLUSIONS
The "fountain filling box" flow with a confined weak laminar plane fountain in a confined open channel with a homogeneous ambient fluid is studied using 2D DNS over the ranges 0.1 ≤ Fr ≤ 3.0, 5 ≤ Re ≤ 800, and 10 ≤ λ ≤ 35 at Pr = 7. The conclusions from the detailed qualitative and quantitative analysis of the behavior of these confined weak laminar plane fountains can be summarized as follows: • An intrusion current results from the downflow of the fountain impinging on the channel bottom. The behavior of the intrusion flow can be approximately described as a gravity current. 30 The decrease of Fr or the increase of Re can decrease τ w , that is, the time scales for the intrusion front to impinge on the sidewall. • Three mechanisms are observed for the behavior of the wall fountain, that is, no-falling (no overturning), slumping down, and rolling down. The maximum penetration height y m of the wall fountain increases with the increase of Fr and λ, due to the reduction of buoyancy flux. • For the fountains with 0.1 ≤ Fr ≤ 1.0, Δτ w~F r 0.61 and Δτ m~F r 0.69 are obtained, where Δτ m is the period for the wall fountain front to reach its maximum penetration height from the creation of the intrusion flow. For the space with the same dimension, only the buoyancy flux keeps constant, and the time-scale will follow a 2/3 power law with Fr. For Fr ≤ 1.0, the intrusion current spreads fast and remains laminar.
• Convection, filling and conduction all contribute to the formation of thermal stratification. In the initial stage, convection and mixing play a key role. After a quasi-steady stratification is formed, filling and thermal conduction become dominant. The behavior of intrusion and the wall fountain for the fountain at 5 ≤ Re ≤ 100, due to its conduction-dominant nature, is significantly different from that at larger Re values considered (150 ≤ Re ≤ 800), where convection plays a more significant role. Additionally, Fr = 1.0 and Fr = 2.25 are found to distinguish the influence of Fr into three ranges. With a smaller Re, the influence of thermal conduction becomes more significant. • Compared with the previous DNS results of confined weak round fountains, 27 it is found that it takes a shorter time for the plane fountain than the round one to fill the box of the same dimension, that is, the stratification rate of the plane fountain is larger than that of the round fountain with the same Fr and Re. Significant instability is observed for the confined weak round fountains during the filling process, however, there is no bobbing or flapping instability observed for the corresponding confined plane fountain.