An Analytical Verification Test for Numerically Simulated Convective Flow above a Thermally Heterogeneous Surface

An analytical solution of the Boussinesq equations for the motion of a viscous stably stratified fluid driven by a surface thermal forcing with large horizontal gradients (step changes) is obtained. This analytical solution is one of the few available for wall-bounded buoyancy-driven flows. The solution can be used to verify that computer codes for Boussinesq fluid system simulations are free of errors in formulation of wall boundary conditions and to evaluate the relative performances of competing numerical algorithms. Because the solution pertains to flows driven by a surface thermal forcing, one of its main applications may be for testing the no-slip, impermeable wall boundary conditions for the pressure Poisson equation. Examples of such tests are presented .


Introduction
Thermal disturbances associated with variations in underlying surface properties can drive local circulations in the atmospheric boundary layer (Atkinson, 1981;Briggs, 1988;Hadfield et al., 1991;Segal and Arritt, 1992;Simpson, 1994;Mahrt et al., 1994;Pielke, 2001;McPherson, 2007;Kang et al., 2012) and affect the development of the convective boundary layer (Patton et al., 2005;van Heerwaarden et al., 2014).Computational fluid dynamics (CFD) codes for modelling such flows commonly solve the Boussinesq equations of motion and thermal energy for a viscous/diffusive stably stratified fluid.In this paper we present an analytical solution of the Boussinesq equations for flows driven by a surface thermal forcing with large gradients (step changes) in the horizontal.The solution can be used to verify that CFD codes for Boussinesq fluid system simulations are free of errors, and to evaluate the relative performances of competing numerical algorithms.Such verification procedures are important in the development of CFD models designed for research, operational, and classroom applications.
We solve the linearized Navier-Stokes and thermal energy equations analytically for the case where the surface buoyancy varies laterally as a square wave (Fig. 1).Attention is restricted to the steady state.No boundary-layer approximations are made; the solution is non-hydrostatic, and both horizontal and vertical derivatives are included in the viscous stress and thermal diffusion terms.The solution is similar to that of Axelsen et al. (2010) for katabatic flow above a cold strip but is easier to evaluate (no slope present) and applies to the more general scenario where the viscosity and diffusivity coefficients can differ.The flow is also similar to a special case (no slope) considered by Egger (1981), although a final analytical solution was not provided in that study.Strictly speaking, the linearized Navier-Stokes equations apply to a class of very low Reynolds number motions known as creeping flows.Such flows appear in studies of lubrication, locomotion of microorganisms, lava flow, and flow in porous media.Of course, for the task at hand, if our linear solution is to serve as a benchmark for a nonlinear numerical model solution, it is essential that the parameter space be restricted to values for which the model's nonlinear terms are negligible.
Because the solution pertains to flows driven by a surface thermal forcing, one of its main applications may be as a test for surface boundary conditions in the pressure Poisson equation.In models of atmospheric boundary-layer flows, the buoyancy is a major contributor to the forcing term in the Poisson equation and also appears in the associated Published by Copernicus Publications on behalf of the European Geosciences Union.surface boundary condition.The pressure boundary condition on a solid boundary in incompressible (Boussinesq) fluid flows is an important and complex issue that has long been fraught with technical difficulties and controversies (Strikwerda, 1984;Orszag et al., 1986;Gresho and Sani, 1987;Gresho, 1990;Temam, 1991;Henshaw, 1994;Petersson, 2001;Sani et al., 2006;Rempfer, 2006;Guermond et al., 2006;Nordström et al., 2007;Shirokoff and Rosales, 2011;Hosseini and Feng, 2011;Vreman, 2014).Typical fractional-step solution methodologies and associated pressure (or pseudo-pressure) boundary-condition implementations are often verified using various prototypic flows such as Poiseuille flows, lid-driven cavity flows, flows over cylinders or bluff bodies, viscously decaying vortices, and dam-break flows.We are unaware of verification tests in which flows were driven by a heterogeneous surface buoyancy forcing.Our solution is designed to fill this gap.
The analytical solution is derived in Sect. 2. In Sect.3, this solution is compared to numerically simulated fields in a steady state.Two versions of a numerical code are run: a version in which the correct surface pressure boundary condition is applied, and a version in which the pressure condition is mis-specified.A summary follows in Sect. 4.

Analytical solution
We derive the solution for steady flow over an underlying surface along which the buoyancy varies laterally as a single-harmonic function.This single-harmonic solution is then used as a building block in a Fourier representation of the square-wave solution.

Governing equations
Consider the flow of a viscous stably stratified fluid that fills the semi-infinite domain above a solid horizontal surface (placed at z = 0).This surface undergoes a steady thermal forcing that varies periodically in the right-hand Cartesian x direction but is independent of the y direction.The two-dimensional (x, z) flow is periodic in x and satisfies the linearized (assuming the disturbance is of small amplitude) governing equations under the Boussinesq approximation, Apart from notational differences, Eqs. ( 1)-( 4) are the two-dimensional steady state versions of Eqs. ( 55)-(57) of Sect.II of Chandrasekhar (1961).Equations ( 1) and ( 2) are the horizontal (x) and vertical (z) equations of motion, respectively, Eq. ( 3) is the thermal energy equation (differential form of the first law of thermodynamics) expressed in terms of the buoyancy variable (defined below), and Eq. ( 4) is the incompressibility condition.Here u and w are the horizontal and vertical velocity components, ≡ [p − p e (z)]/ρ w is the kinematic pressure perturbation [p is pressure, p e (z) is pressure in a hydrostatic environmental state in which the density profile is ρ e (z), ρ w is a constant reference density, say, ρ e (0)], and b ≡ −g[ρ − ρ e (z)]/ρ w is the buoyancy, where ρ is the actual density, and g is the acceleration due to gravity.The Brunt-Väisälä frequency N ≡ √ −(g/ρ w )dρ e /dz of the ambient fluid (Kundu, 1990), kinematic viscosity ν, and thermal diffusivity α are taken constant.
We obtain our solution using a standard vorticity/streamfunction formulation.
(3) and (5) yields Introducing a streamfunction ψ defined through guarantees that Eq. ( 4) is satisfied and transforms Eq. ( 6) into a single equation for ψ, The dependent variables are assumed to vanish far above the surface (z → ∞).On the surface we apply no-slip (u = 0) and impermeability (w = 0) conditions, and specify a periodic (in x) buoyancy distribution.As we will now see, restricting the dependent variables to steady periodic forms that vanish as z → ∞ also restricts acceptable distributions of the surface buoyancy.The restriction was first noted by Egger (1981, Sect. 3c), though without details.Averaging Eq. ( 3) over one period (using w = −∂ψ/∂x) yields d 2 b/dz 2 = 0, which integrates to b = A + B z ( b is the average of b; A and B are constants).Taking b → 0 as z → ∞, implies that b → 0 as z → ∞, in which case A = B = 0, and b(z) = 0.In particular, at the surface, b(0) = 0.If a surface distribution b(x, 0) violates this condition, the ground acts as a net heat source/sink.In an unsteady model, such a source/sink would force a continually upward-developing disturbance and a steady state could never be attained.

Single-harmonic forcing
For a surface buoyancy of the form b(x, 0) ∝ sin kx, Eq. ( 3) indicates that ψ is of the form Application of Eq. ( 9) in Eq. ( 8) yields which has solutions of the form A ∝ e Mz for M satisfying Taking the one-third power of Eq. ( 11) yields a useful intermediate result: where n is an integer.Rearranging Eq. ( 12) and taking the square root yields Equation ( 13) furnishes six roots, two for each of n = 0, 1, 2.
To ensure that A(z) → 0 as z → ∞, we reject the roots with a positive real part.With the radicand of Eq. ( 13) expressed in polar form, the physically acceptable roots are where the subscript on M denotes the associated value of n, and r and φ are defined by While solving Eq. ( 16) for φ, care must be taken when evaluating arcsin or arccos functions that φ appears in the correct quadrant (φ should be in quadrant I or II so φ/2 should always be in quadrant I).Also note from Eq. ( 14b) and (14c) that M 2 is the complex conjugate of M 1 (M 2 = M * 1 ), a fact that will often be used below.
With the general solution for ψ written as where B, C, and D are constants, the vorticity becomes, and the buoyancy follows from Eq. ( 3) as where ∇ 2 b h = 0.In view of Eq. ( 12), Eq. ( 19) becomes Applying Eqs. ( 18) and (20) in Eq. ( 5) yields an equation for ∂b h /∂x which, upon use of Eq. ( 12) and M 2 = M * 1 , reduces to ∂b h /∂x = 0.So b h is, at most, a function of z.Since ∇ 2 b h = 0, b h is, at most, a linear function of z, and since b should vanish as z → ∞, that linear function must be 0. Thus, b h = 0.

Piecewise constant (square wave) forcing
Next, consider the case where the surface buoyancy varies horizontally as a square wave, with a distribution over one period L given by b Such a distribution can be expressed as the Fourier series: Application of Eq. ( 34) in Eq. ( 36) yields The solutions for b, ψ, u, and w can then be written as summations over the single-harmonic solutions ( 29), ( 30), (32), and ( 33), with k related to n by and with b 0 replaced by b n :

Verification tests
A solution of the linearized equations may be used to verify a nonlinear code if the nonlinear terms are sufficiently small.Unfortunately, a priori estimates of such terms expressed, for example, through a Reynolds number, are not straightforward since the relevant velocity and length scales in our problem are only evident after a solution has been obtained.We thus seek an appropriate set of test parameters through trial and error, guided by a posteriori linear solution estimates of the terms u • ∇b and u • ∇η[u = (u, w)] present in nonlinear versions of Eqs.
(3) and ( 5), respectively.Specifically, for any computed candidate solution, we formed the ratios of the largest values of those nonlinear terms to the largest values of the corresponding linear terms, that is, the terms actually present in Eqs.
(3) and ( 5).We need only consider one such linear term per ratio since Eqs.
(3) and ( 5) are comprised of two terms of equal magnitude.A solution was deemed to be sufficiently linear if where ε ( 1) is a prescribed threshold.The suitability of this approach was confirmed by the very close agreement between the analytical solutions and the numerical solutions obtained with the correct surface pressure condition.
The numerical model employed in our tests is a variant of a direct numerical simulation (DNS) code used in the boundary-layer and slope-flow studies of Fedorovich et al. (2001), Fedorovich andShapiro (2009a, b), and Shapiro andFedorovich (2013, 2014).The model solves the Boussinesq governing equations on a staggered (Arakawa C) grid.Although designed for three-dimensional simulations, the model was run in a two-dimensional (x, z) mode.The overall solution procedure is patterned on a fractional step method proposed by Chorin (1968).In our version, the prognostic equations are integrated using a filtered leapfrog scheme with explicit treatment of the viscous term.The pressure is diagnosed from a Poisson equation (Eq.A3b, discussed in the Appendix), which is solved using a fast Fourier transform technique in horizontal planes, and a tridiagonal matrix inversion in the vertical.The surface condition on pressure is the inhomogeneous Neumann condition (INC) that arises from projecting the vertical equation of motion into the vertical and imposing the impermeability condition (Vreman, 2014; also see the Appendix).We also run a version of the code in which the surface pressure condition is mis-specified as a homogeneous Neumann condition (HNC).We hasten to add, however, that our implementation of the HNC may be quite different from implementations described in the literature.We elaborate on these technical differences and review general aspects of the problem of surface pressure specification in the Appendix.
The analytical solution was evaluated on an un-staggered (x, z) grid extending over one period of the square wave (x = 0 to x = L).The series were truncated at 50 000 terms.The governing parameters were adjusted so that the linearity criteria were satisfied in comparisons with ε = 5 × 10 −3 .
In the first test, we set ν = α = 0.001 m 2 s −1 , N = 0.02 s −1 , L = 5.12 m, and b max = 1×10 −5 m s −2 .For the analytical solution A-1, the (x, z) grid consisted of 513 points in the x direction and 1025 points in the z direction, with grid spacings x = z = 0.01 m.The linearity criteria (Eq.43) were satisfied with R η ∼ = 8.2 × 10 −5 and R b ∼ = 2.8 × 10 −3 .The analytical b and w fields shown in Fig. 2 depict a broad zone of ascent above the warm surface and a compensating zone of descent over the cold surface, roughly for z < 1.8 m.In the upper part of these zones (at roughly 0.9 m < z < 1.8 m), adiabatic expansion/compression has reversed the senses of the buoyancy fields.Surprisingly, the numerical fields in the inhomogeneous INC-1 and homogeneous HNC-1 cases are very similar to each other and to the A-1 fields.The u fields from A-1, INC-1, and HNC-1 shown in Fig. 3 are visually indistinguishable from one another.
To understand why the INC-1 and HNC-1 simulations are so similar, and to identify simulation parameters that might evince more substantial differences, we consider the idealized problem in which a specified buoyancy b = b 0 e −γ z sin kx (γ = h −1 , where h is the e-folding depth scale) is the only forcing term in the Poisson equation ∇ 2 = ∂b/∂z, with Neumann surface condition ∂ /∂z| 0 = b(x, 0).This idealized problem is solved as The corresponding solution obtained with the homogeneous Neumann condition, ∂ /∂z| 0 = 0, is The relative error (RE) in the vertical pressure gradient force associated with Eqs. ( 44) and ( 45), defined as the local absolute error in that force divided by the local buoyancy, is calculated as where a ≡ γ /k.Written in terms of the depth scale h and wavelength λ = 2π/k, a can be interpreted as an aspect ratio characterizing the width to depth scales of the disturbance, a = λ/(2π h) ∝ λγ .From Eq. ( 46) we see that RE decreases exponentially with z for disturbances characterized by small aspect ratios, a < 1 (which we refer to as deep disturbances), and increases exponentially with z for disturbances characterized by large aspect ratios, a > 1 (which we refer to as shallow disturbances).The buoyancy in Fig. 2 is suggestive of a < 1, which indicates that the first test could be classified as a deep (error-forgiving) simulation.
The preceding analysis suggests that simulations with shallow thermal disturbances (a > 1) might yield large differences between cases with inhomogeneous and homogeneous Neumann conditions.There did not appear to be a straightforward way to increase the effective a by systematically varying the parameters (e.g.increasing L tended to increase the effective h), but a set of suitable parameters were identified through trial and error and were used as the basis for the second test case.

HNC-2
x (m) z (m)  and 2: the low-level thermal disturbance in the second test is much shallower than the disturbance in the first test (and is suggestive of a > 1).In this second test case we find dramatic differences between the inhomoge-neous INC-2 and homogeneous HNC-2 cases.Specifically, while the INC-2 and A-2 fields are in excellent agreement, the HNC-2 fields showed no signs of even approaching a steady state.Long after the INC-2 simulation had reached a steady state, the HNC-2 fields continued to amplify and develop asymmetric structures associated with flow nonlinearities.The very close agreement between the A-2 solution and the steady state in the INC-2 simulation is shown for the u field in Fig. 5.The u field in the disastrous HNC-2 simulation, at a time when a steady state had already been attained in the INC-2 simulation, is shown in Fig. 6.

Summary
The linearized Boussinesq equations for the motion of a viscous stably stratified fluid are solved analytically for a surface buoyancy that varies laterally as a square wave.The solution describes two-dimensional laminar convective structures such as thermal convective rolls and updraft/downdraft pairs.The main applications of the solution may be in code verification and the evaluation of different implementations of the surface pressure condition for the pressure Poisson equation.Tests have been conducted for cases where the aspect ratios of the thermal disturbance have been large and small.With attention restricted to disturbances of sufficiently small amplitude, the linear solution and numerically simulated fields with the inhomogeneous Neumann condition for pressure (which is appropriate in the context of the particular fractional step procedure adopted in our DNS code) have been found to be in excellent agreement for both tests.However, in tests with a mis-specified Neumann condition, an excellent agreement with the analytical solution has been found only for the deep (small aspect ratio) disturbance case; errors in the shallow (large aspect ratio) disturbance case have been catastrophic.associated with the conflation of the two physically unjustifiable specifications (homogeneous Neumann condition for pressure, and u| 0 = 0) cancel out.The homogeneous Neumann condition for pressure can be the source of confusion if the context in which the condition is applied is not made clear: it would be a correct condition if u| 0 is set to 0 (per the equivalence described above), but it would be an incorrect condition if the explicit model-computed values of u| 0 are used.In the experiments with the mis-specified condition described in Sect.3, the homogeneous condition is imposed in the latter context.Unfortunately, in many numerical model descriptions, the nature of the surface pressure condition is left vague, for example, by not indicating whether a Neumann condition is homogeneous or inhomogeneous or, if a homogeneous Neumann condition is indicated, not mentioning how u| 0 is treated.
Finally, we note that in fractional step procedures that treat the viscous term implicitly (e.g.Kim and Moin, 1985;Gresho, 1990;Armfield and Street, 2002;Guermond et al., 2006, and many others), the homogeneous Neumann con-dition is often applied as a surface condition for a Poisson equation, but it is again different from our implementation described in Sect.3. In the implicit treatments, the provisional velocity is obtained as the solution of a boundary value problem ( u| 0 should be specified; often it is set to 0) in which the relevant Poisson equation resembles Eq. (A5) but applies to a scalar function (sometimes called a pseudo-pressure) that is not the real pressure.Temam (1991) refers to this scalar as, ". . .a technical quantity, a mathematical auxiliary . . ." and advocates that it should not even be considered as an approximation of the pressure.Interestingly, in the context of implicit treatments, the homogeneous Neumann condition on the pseudo-pressure has sometimes been implicated as corrupting solution accuracy through the development of spurious numerical boundary layers adjacent to solid boundaries (Gresho, 1990;Guermond et al., 2006;Hosseini and Feng, 2011).

Figure 1 .
Figure1.Schematic of two-dimensional (x, z) thermal convection induced by a surface buoyancy that varies horizontally (x) as a square wave.Red denotes positive surface buoyancy, blue denotes negative surface buoyancy.

Figure 2 .
Figure 2. Vertical cross section of the analytical (A-1) buoyancy b and vertical velocity w fields from the first test case.Colour bar units are m s −2 for b, and m s −1 for w.

Figure 3 .Figure 4 .
Figure 3. Vertical cross section of u from the first test case.A-1 is the analytical solution.INC-1 is the numerical simulation with inhomogeneous Neumann condition for pressure.HNC-1 is the numerical simulation with the homogeneous Neumann condition for pressure.Colour bar units are m s −1 .

Figure 5 .
Figure 5. Vertical cross section of u from the second test case.A-2 is the analytical solution.INC-2 is the numerical simulation with inhomogeneous Neumann condition for pressure.Colour bar units are m s −1 .

Figure 6 .
Figure 6.Vertical cross section of u from HNC-2, the numerical simulation with homogeneous Neumann condition for pressure in the second test case.Colour bar units are m s −1 .