A 3D Numerical Study of Interface Effects Influencing Viscous Gravity Currents in a Parabolic Fissure, with Implications for Modeling with 1D Nonlinear Diffusion Equations

Although one-dimensional non-linear diffusion equations are commonly used to model flow dynamics in aquifers and fissures, they disregard multiple effects of real-life flows. Similarity analysis may allow further analytical reduction of these equations, but it is often difficult to provide applicable initial and boundary conditions in practice, or know the magnitude of effects neglected by the 1D model. Furthermore, when multiple simplifying assumptions are made, the sources of discrepancy between modeled and observed data are difficult to identify. We derive one such model of viscous flow in a parabolic fissure from first principals. The parabolic fissure is formed by extruding an upward opening parabola in a horizontal direction. In this setting, permeability is a power law function of height, resulting in a generalized Boussinesq equation. To gauge the effects neglected by this model, 3D Navier-Stokes multiphase flow simulations are conducted for the same geometry. Parameter variations are performed to assess the nature of errors induced by applying the 1D model to a realistic scenario, where the initial and boundary conditions can not be matched exactly. Numerical simulations reveal an undercutting effect observed in laboratory experiments, but not modeled when the Dupuit-Forchheimer assumption is applied. By selectively controlling the effects placed on the free surface in 3D simulations, we are able to demonstrate that free surface slope is the primary driver of the undercutting effect. A consistent lag and overshoot flow regime is observed in the 3D simulations as compared to the 1D model, based on the choice of initial condition. This implies that the undercutting effect is partially induced by the initial condition. Additionally, the presented numerical evidence shows that some of the flow behavior unaccounted for in the 1D model scales with the 1D model parameters.


Introduction
Multiple mechanisms can be responsible for the movement of fluids. Gravity currents occur when a difference in the density between two or more fluids induces movement. Applications of this type of movement include both natural [1][2][3][4] and industrial [5,6] settings. Examples of such flows are the model outlined in Section 2.2. This approach allows us to isolate the effects of flow assumptions and parameters in a way that would be intractable using physical experiments. Results of the numerical simulations are outlined in Section 3. From these, mechanisms for the discrepancies between the two models is discussed in Section 4.

1D Flow Model
We derive a flow model for a viscous fluid in a parabolic fissure. Consider a fluid element formed by vertically cutting through a parabolic prism, as seen in Figure 1. Conservation of mass dictates that the rate of change of fluid in the elementary volume is equal to the difference in volumetric fluxes through the faces at y and y + dy: where the width of the wedge for height z is given by b(z), u v is the fluid speed in the y direction, and h v (y, t) is the height of the free surface. We consider only the flux of an incompressible high viscosity fluid (e.g., glycerol), disregarding the flux of gaseous atmosphere being displaced. Neglecting the higher order terms of the Taylor expansion in the second term on the right hand side of Equation (1) yields, As seen in Figure 1, the parabolic cross section is oriented against gravity in the x-z plane. The equation of the fissure walls in the x-z plane is, for some constant ω > 0. Then the distance between the walls as a function of height is, Under the assumption of Poiseuille flow between parallel plates of width b(z) for an arbitrary height z, permeability as a power law function of height is obtained: This permeability results in the Darcy velocity, Integrating Equation (2), with substitutions for width (4) and Darcy velocity (6), yields ∂ ∂t which further simplifies to, The change of variables, Υ = h 3 2 , is implemented to obtain a special form of the PME [15]: We are interested in case of the free spreading of a finite pulse of liquid under the force of gravity, which results in the boundary condition, where σ > 0 for an initially empty fissure [11]. The problem defined by Equations (9) and (10) can be solved by introducing dimensionless similarity variables, and where ξ(y, t) is the similarity variable, H(ξ) is a scaling function, m = 7/3 is the nonlinear diffusion exponent in the PME, and The parameter β is the time exponent of the boundary condition, in this case − 1 7/3+1 . Substitution of the similarity variables results in the ordinary differential equation, where λ = β 1+β(m−1) , and H(0) = 1, H(ξ 0 ) = 0. The front position in terms of similarity variables is ξ 0 , and is found in the process of solution. The "mathematician's pressure" change of variables [15], performed on Equation (14) produces an ODE with integer powers: Boundary conditions are given by, and Equations (16) and (18) can be solved using a Runge-Kutta solver [30], or a power series of the form found in [29]:

3D Flow Model
Flow through the parabolic fissure can be modeled with the incompressible 3D Navier-Stokes equations, and where ρ is density, u is velocity, p is pressure, µ is viscosity, and g is the gravity vector. Phases are tracked using a VOF method [31,32]. The indicator function α takes on values in the range of [0, 1] and can be interpreted as the fraction of two phase types within a computational cell. For our purposes, 0 denotes air, and 1 denotes a viscous liquid. For any computational cell in the domain, the average fluid properties for density and viscosity are calculated by, and where α i denotes the volume fraction of phase i. The indicator function evolves under the equation, where the counter gradient transport term ∇ · u c α (1 − α) acts as a limiter based on the maximum velocity u c on interface. Surface tension effects are calculated by, where γ is the surface tension, and K is the curvature −∇ · n. S is added to the momentum Equation (20) to obtain, This method has been successfully used in applications at a variety of scales [33][34][35][36]. Our use of this method is very conservative, in that the flow is laminar and surface effects are not the main driver of flow. We use the interFoam multiphase flow solver included in the finite volume method toolbox openFOAM [37] for all 3D flow simulations. A hexahedral base mesh was constructed with non-uniform grading decreasing cell size towards the y-z plane and base of the parabola. This mesh was snapped to a 3D parabolic geometry using the snappyHexMesh utility. The parabolic geometry was constructed to the specifications given in Section 2.3. The same mesh was utilized for all simulations, with time-step limiting based on a maximum Courant number of 0.25. With the available computing resources, each simulation required under two days to complete 3 s of simulated time.
Boundary conditions are chosen to mimic a lab experiment, with the top of the parabola open to the atmosphere. For velocity, a no-slip condition, is applied to all solid walls. To allow fluxes of air at the top of the experiment, a zero gradient condition is applied with respect to the patch normal vector: A no-flow pressure condition is applied at the solid walls, while at the top of the experiment total pressure P tot is calculated using the gauge pressure p atm , Unless otherwise specified, the indicator function α is fixed at the top of the experiment and zero gradient at solid walls,

Flow Assumptions, Initial Conditions and Parameter Values
The 1D flow model uses the DF assumption, and does not include the effects of air flow, surface tension or contact angle. The 3D model includes a complete velocity field and optional effects at the liquid-gas-wall interface, but requires large amounts of computational power to solve. This is largely due to the aspect ratio of the domain, which requires very small finite volumes between the parabola walls. The parabolic geometry defined by Equation (4) is extruded along the y axis for one meter with a value ω = 2000. This corresponds to a gap width of x = 2 cm at a maximum height of z = 0.2 m.
An initial condition for the 3D simulation was chosen from the 1D analytical solution using the parameter values in Table 1. The 1D viscosity is held constant at µ =1.4 Pa · s for generating the initial condition, though this is effectively varied as part of a (Equation (13)) in the generation of results. Choosing a solution curve from the 1D model at a nonzero starting time is required, as the solution collapses to a singularity at t = 0 s and the height of the solution is above the z ∈ [0, 0.2] m of computational domain at early time. The resulting curve at time t = 0.3 s was chosen as a suitable initial condition containing enough potential energy to produce interesting results. Additional 3D numerical experiments were performed under different values of the parameter a using the same initial condition by exploiting the scaling nature of the similarity solution. After the ODE system Equations (16) and (18) is solved and ξ 0 obtained, the propagation distance y of the 1D solution can be computed as a function of time: To compare between 1D and 3D results, the various moments in time from the 1D model to which the initial condition used in 3D simulations corresponds must be calculated. Let the propagation distance of the aforementioned numerical initial condition be y 0 , and a' be an alternate parameter value calculated from a 3D simulation parameter set. Then the start time of the numerical initial condition under the parameter set that produced a' is, .

Results
In Section 3.1, the models are validated against one another, and the mesh quality of the numerical simulations is shown to be adequate. The mechanisms of discrepancy between the 1D and 3D models is established. Section 3.2 explores the less prominent effects of contact angle and surface tension at the free surface. Additionally, we examine the transition to the similarity flow regime from an initial condition that is not in the similarity solution set. Finally, Section 3.3 examines the behavior of the 3D solution under variations in the parameter a, which is compared against the 1D solution for the analytical initial condition.

Profile Comparison and Validation
A comparison between a baseline 3D simulation with γ = 0 and the 1D similarity solution is shown in Figure 2. Excellent agreement is observed between the two models, both in propagation distance and shape. The poorest fit is seen at the earliest time, Figure 2a. The simulated front lags behind the analytical solution at the base of the parabola, and undercutting is observed at the bottom of the advancing front. At later times, the free surface of the 1D analytical solution lags behind the NSE simulation at the "nose" of the advancing front.
The flow profiles observed in this simulation indicate that a steep free surface is responsible for the undercutting observed at early time. The effects of surface tension and contact angle are not included in the simulation used to produce Figure 2, a situation that would be difficult to replicate with experimental techniques. The 1D flow model presented in Section 2.1 stems from a larger family of nonlinear diffusion equations which have been validated for flow through various fissures. These equations can be derived by substituting an alternative expression for Equation (5). In the case of parallel walls, the resulting Boussinesq equation was experimentally validated by Huppert and Woods [1] using a Hele-Shaw cell and glycerol. Experiments with angled planar walls were validated by Furtak-Cole et al. [29] as well as Ciriello et al. [22] using similar methods. The permeability profile from Equation (5) is bounded on either side by these two classes of experiments, making experimental validation unnecessary.
The interFoam solver presented in Section 2.2 is a type of VOF method [31] which has been used and validated for a variety of applications [33,36]. The cases presented in this paper represent a conservative use of this method, as the maximum Reynolds number is less than one. The role of mesh quality was tested at two levels. Base meshes were constructed with 60 and 80 cells partitioning the width between the parabolic walls. This width is consistent with the x direction shown in the fluid element for the 1D model, see Figure 1. Cell counts in the remaining directions were calculated to maintain a 1-1-1 average aspect ratio. As seen in Figure 3, a consistent rate of propagation was reached at both levels. Propagation distance for the fine mesh lags slightly behind as expected, due to its ability to better resolve the free surface. The finer mesh was used for all simulations presented here. Validation of both 1D and 3D flow models in independent settings, as well as their agreement displayed in Figure 2, implies an acceptable level of accuracy for both solution types.

Effects of Surface Tension, Contact Angle, and Initial Conditions
The similarity solution to the 1D flow problem is a weak, continuous solution. By contrast, the interface resolved by the numerical simulations is double valued at the propagating edge. To quantify this early time behavior observed in experiments and simulations, we measure the height of the undercutting observed in Figure 2a. This is measured vertically from the base of the parabola to the "nose" of the advancing front, where the air-liquid interface is vertical. To further test the assertion that the vertical velocity is primarily responsible for the undercutting behavior, the height of undercutting for simulations with surface tension, dynamic contact angle, and an alternative initial condition are measured. Surface tension for the 3D air-liquid interface is listed in Table 1 and based on the value for glycerol, a viscous fluid commonly used in experiments. A dynamic contact angle θ was modeled using the function, where θ 0 is the static contact angle, u wall the speed of the interface tangent the static wall, u θ is a scaling factor, and θ a and θ r are the advancing and receding angles respectively. Given the hypothetical nature of this study and lack of physical data, we arbitrarily choose a neutral parameter set, θ 0 = 90 • , θ r = 60 • , θ a = 120 • , and scaling parameter u θ = 1.
To perform physical experiments, initial conditions that do not coincide with the similarity solution are often used out of necessity. Such initial conditions have been observed to converge quickly on the similarity solution [5,38]. Numerical experiments offer a previously unobtainable opportunity to compare the convergence of arbitrary initial conditions. A hypothetical mass-conserving initial condition of rectangular cross section was produced by solving, where the upper limit z = 0.18 was chosen arbitrarily. The value for the right hand side was produced by numerically computing the volume integral defined by the analytically-derived initial condition from the 1D flow model in the parabolic geometry. For practical purposes, such an initial condition could be used in an experiment by setting a lock gate to hold a volume of liquid at the initial time.
Undercut heights for simulations with no interface effects (baseline), contact angle, surface tension, and the alternative rectangular initial condition are shown in Figure 4. Propagation distance and speed for the baseline and rectangular initial condition 3D simulations are compared to the 1D model in Figure 5.

Scaling Effects
The previously presented results raise the question of whether or not early time inaccuracies of the 1D model can be predicted knowing only the parameters in Equation (9). Given that the coefficients of the right hand side of Equation (9) can be consolidated in a single constant a (Equation (13)), a single parameter can be varied without changing the ode problem given by Equation (14). A natural choice is to vary µ. This allows for all simulations to be performed with the same computational settings, preventing any false comparison of finite volume mesh-induced errors. For the same initial condition, simulations were conducted at the viscosities µ = 0.8 Pa · s, 1.2 Pa · s, 1.6 Pa · s, 2.0 Pa · s, and 2.4 Pa · s. This approach assumes that the parameter values result in a laminar flow, and that at the boundary y = 0 the fixed parameters of Equation (10) guarantee the liquid phase free surfaces change strictly due to gravity. For each simulation, Equation (33) was used to calculate the simulation start time for the single numerical initial condition chosen in Section 2.3, since the same initial profile for different values of µ corresponds to different initial moments in time. Propagation distances are shown in Figure 6, while undercut heights are presented in Figure 7.

Discussion
The DF assumption is known to perform poorly for steep free surfaces and the data presented from the numerical simulations make it possible to quantify this deficiency. Use of the DF assumption is often preceded with stipulations such as the ratio of horizontal displacement over height being much larger than one [1]. A smaller ratio results in an increasingly steep free surface with more vertical velocity. Our numerical simulations indicate that vertical velocity components are not only important at the advancing front. Figure 8 shows velocity vectors and their magnitude at a cross section of a baseline simulation as given in Figure 2, after 0.4 s of simulation time. Within the fluid mound, velocities are large only near the interface, while fluid at the lower-left corner is nearly at rest. Primarily vertical velocities are seen at both steep and nearly horizontal segments of the free surface. We calculate continuity error numerically using two control volumes given by a modified version of Equation (1). Here, u ν is given by a surface integral of the spatially variable y component of the velocity in the liquid phase of the 3D simulation seen in Figure 8. Control volumes of 0.01 m width are placed 0.02 m from the left wall and 0.01 m to the left of the advancing front. The time derivative for the volume integral on the left hand side of Equation (1) is approximated with a forward derivative calculated over 0.01 s. The left control volume, positioned at a nearly horizontal free surface, results in a relative continuity error of 0.54 %. The right control volume, located on the steep free surface of the advancing front, gives a relative error of 5.06 %. In this sense, the 1D model holds well for vertical motions that do not appear near a vertical interface. This indicates that nose formation at the advancing front is at least partially due to a flow regime where fluid "skims" in a direction parallel to the free surface, perhaps better called "overtopping" than "undercutting". What we have measured as "undercut" is in fact a combination of mechanisms, including viscous drag unaccounted for in the 1D model. In some experiment configurations [1,5] viscous drag is not modeled at the base of the experiment and this effect can not be distinguished from vertical velocity due to the steep free surface. In these cases, it is likely that overtopping is actually compensating viscous drag at the base of experiments, resulting in good agreement between the 1D and experimental results at later times. Indeed, other experiments using vertical wedge configured plates [22,29] experience undercut profiles. No unaccounted drag exists at the base of experiments in these cases, though effects such as surface tension become immense where the plates meet and the cited experiments were performed at an empty initial state with flux boundaries. For the 1D model of a parabolic geometry there is no unaccounted drag at the base of the parabola, though an analogous problem exists: as height decreases to zero, Equation (6) is an increasingly poor approximation of the Darcy velocity because it is based on horizontal width only.  Figure 7. (a) Undercut heights plotted from peak value as a function of time. Shifting of the peaks is largely due to the start times calculated by Equation (33). (b) Undercut heights plotted from peak value as a function of space. In all cases, only µ is varied and no interfacial effects are present.
As shown in Figure 4, the magnitude of measured undercutting can be modified by small perturbations in the initial condition or effects acting at the free surface. The rectangular initial condition's infinite slope dramatically increased the magnitude of undercutting, but the solution rapidly converged on the same trajectory as the analytical initial condition within two seconds. Surface tension caused more extreme long term behavior. The extra force acting parallel to the free surface resists bending, allowing the fluid to override the lowest region of the parabolic wedge for an extended period of time. The Bond number is calculated Bo = ∆ρgL 2 /γ, and represents the interplay between gravitational and surface tension forces. Here, the characteristic length L 2 is determined by the width of the wedge at the height of undercutting in Figure 9a. In the surface tension case, undercut is consistently higher than the baseline case, with the exception of the curves meeting between 1 s and 1.5 s, as seen in the difference plot Figure 9b. This corresponds to large vertical velocities along the leading edge driving the liquid phase downward. Similar residuals are observed at 0.75 s and 2.5 s for an order of magnitude difference in Bo, emphasizing the complex relationship between dynamic velocity and surface tension at the interface. The additional simulation using fairly neutral dynamic contact angle parameters shows little change in undercutting as compared to the baseline case.
This could undoubtedly be altered by selecting more extreme parameters, however the complicated physics and numerical issues of this topic warrant a separate study. As with the undercutting results, the propagation plots of Figure 5 show the numerical simulations converging on one another rapidly from different initial conditions. While there appears to be discrepancy between the simulations and 1D solution at later times, it is important to note that propagation is measured as the furthest point at the nose of the advancing front. Thus the shape of the profile plays an out-sized role in the measurement of propagation distance, which is not an all-inclusive means of comparing a 3D interface to a 1D solution. More importantly, Figure 5b shows convergence between the propagation velocities for the 3D numerical simulations and 1D model. Simulations produced for different viscosities show that a period of velocity deficit and overcompensation occurs whenever a simulation is started, as compared to the 1D model solution. This presents a particular problem for experiments where it would be difficult to provide the initial condition with a velocity field from the 1D solution. The starting point of the initial condition relative to the similarity solution and lack of initial velocity are dominant factors in the lack of agreement between the two models. The simulations show increasingly poor convergence to the 1D model the later they are started in time, see Figure 6. In the case of our specific similarity solution, this corresponds to using a higher value of viscosity; a technique often employed to maintain laminar flow in experiments. Initial profiles close to the singular initial condition for the 1D model, and correct initial velocity, are far from the aspect ratio condition that makes the similarity solution applicable. Our results do not establish a general guideline for how severe this effect is; we have not generated a set of different initial profiles corresponding to the same initial moment in time for the different values of viscosity considered here. Still, Figure 6 does show that it would be easy to unknowingly choose a poor combination of initial condition and viscosity (or other parameters that make up a) prior to performing an experiment. Special attention should be paid to this issue.  Figure 9. (a) Bond number calculated at the height of undercut for the surface tension case. (b) Difference between the surface tension and baseline cases. Bond number and differences from simulations presented in Figure 4. The lack of correlation between Bo and the difference occurs as velocity builds from the initial condition.
The dominant role of the initial condition in generating discrepancy is further underscored by plotting the undercut as a function of space in Figure 7b. For each viscosity the time corresponding to the initial condition is different, though each simulation quickly converges to the same trajectory through space. However, the peak value of undercut and its magnitude is shown to shift in space for variations in the value of a. Further analysis reveals a linear relationship between the peak y displacement and the parameter a, shown in Figure 10a. This relationship suggests realistic effects that are not accounted for in the 1D model are still intimately related to the similarity variables governing the 1D solution, and can be predicted as such. Similarly, the magnitude of undercut is plotted against a in Figure 10b, and shows equally linear behavior. The relationships established are purely empirical and further work is needed to understand their mechanism.

Conclusions
This paper investigated some of the effects that are neglected by 1D models of viscous gravity spreading under the DF assumption. To have the ability to gauge each effect separately, we conducted highly accurate 3D NSE VOF simulations. Effects neglected by the 1D model were included separately to assess the magnitude of their influence on fluid flow. These effects include vertical components of the flow, surface tension, and contact angle. Each simulation involved the spreading of a finite volume of liquid under the action of gravity. Undercutting was observed in both the 3D numerical simulations presented here, as well as in previous work by the authors and others. This effect can not be modeled if the DF assumption is used. Undercutting appears when the free surface of the fluid overruns the position of the wetting front at the bottom of the fissure. In this paper, numerical investigations have shown that the primary cause of this effect is vertical velocity induced by a steep free surface. Effects on the free surface may intensify this behavior. Undercutting is strongly influenced by the choice of initial condition. However, late time overall agreement between the two solutions was observed to be very good. This observation is corroborated by previous experimental work. At late times the similarity solution is known to perform well, and the observed height of undercutting takes a small value that changes very slowly for all initial conditions.
To summarize, the placement of the initial condition for numerical experiments relative to the similarity solution parameters produces a wide dilatation in both undercutting and in the position of the advancing front. Initial conditions for numerical simulations that are close quickly converge on each other. Generally, the 1D model predicted the correct speed of propagation, with shortcomings in modeling the shape of the free surface. The DF assumption and associated vertical velocity near steep segments of the free surface are largely responsible for discrepancies between the two models. Effects included at the interface have a secondary role to the above points for the cases considered here. Acknowledgments: Eden Furtak-Cole thanks the Utah Center for High Performance Computing for the management of computing resources and access to their technical expertise.

Conflicts of Interest:
The authors declare no conflict of interest.