Assessment of a symmetry preserving JFNK method for atmospheric convection
AAssessment of a symmetry-preserving JFNK method foratmospheric convection
M. Alamgir Hossain b , Jahrul M Alam a, ∗ a Department of Mathematics and Statistics, Memorial University, Elizabeth Ave, NL A1C 5S7, Canada b Department of Mathematics and Statistics, Simon Fraser University, 8888 University Dr, Burnaby, BC V5A 1S6,Canada
Abstract
Numerical simulations of nonhydrostatic atmospheric flow, based on linearly decoupled semi-implicit or fully-implicit techniques, usually solve linear systems by a pre-conditioned Krylovmethod without preserving the skew-symmetry of convective operators. We propose to performatmospheric simulations in such a fully-implicit manner that the difference operators preserve boththe skew-symmetry and the tightly nonlinear coupling of the differential operators. We demon-strate that a symmetry-preserving Jacobian-free Newton-Krylov (JFNK) method mimics a balancebetween convective transport and turbulence dissipation. We present a wavelet method as an effec-tive symmetry preserving discretization technique. The symmetry-preserving JFNK method forsolving equations of nonhydrostatic atmospheric flows has been examined using two benchmarksimulations of penetrative convection – a) dry thermals rising in a neutrally stratified and stablystratified environment, and b) urban heat island circulations for effects of the surface heat flux H varying in the range of 25 ≤ H ≤
930 W m − . The results show that an eddy viscosity model pro-vides the necessary dissipation of the subgrid-scale modes, while the symmetry-preserving JFNKmethod provides the conservation of mass and energy at a satisfactory level. Comparisons of theresults from a laboratory experiment of heat island circulation and a field measurement of poten-tial temperature also suggest the modelling accuracy of the present symmetry-preserving JFNKframework. Keywords:
JFNK; skew-symmetry; atmospheric convection; physics-based preconditioning;
1. Introduction
Even with increasing power of computers and advances in numerical methods, it is a chal-lenging endeavour to resolve the important physics of convective motion and cascade of turbu-lence kinetic energy (TKE) in atmospheric simulations. The large-scale physics cannot reacha near equilibrium of the interplay between convective transport and diffusive dissipation. Wehave to cope with the formidable problem of the subgrid-scale parameterization of convective ∗ Corresponding author
Email address: [email protected] (Jahrul M Alam)
Preprint submitted to Journal July 22, 2020 a r X i v : . [ phy s i c s . c o m p - ph ] J u l rocesses (Ferziger and Peric, 1997; Pielke, 2002; Pletcher et al., 2013). Mathematically, theconvection ( u · ∇ u ) and the diffusion ( ∇ · τ ) are governed by skew-symmetric and symmetricpositive-definite operators, respectively, which are not fully preserved in many operational at-mospheric modelling codes (see Wicker and Skamarock, 1998; Skamarock et al., 1997; Pielke,2002). A parameterization of the subgrid-scale stress τ would provide a subtle balance of the twooperators, which is broken by the non-symmetric discretization of the skew-symmetric convec-tive operator. Most important reasons for preserving the symmetry of convective operator are: i)improved forecasting skill for meso-scale phenomena; and ii) reduced cost for highly complexforecasting systems when unified for both the boundary layer and the meso-scale phenomena. Asthe discretization breaks the skew-symmetry of the convective operator, the underlying conserva-tion law is not globally satisfied at discrete times (Verstappen and Veldman, 2003). The stabil-ity and conservative properties of the existing non-symmetric schemes are thoroughly reviewedby Steppeler et al. (2003) and Klemp et al. (2007). Studies observed that higher-order linearlyconsistent schemes are sometimes unable to deal with the contamination of poorly resolved short-wavelength perturbations, usually triggered by atmospheric convection (Skamarock et al., 1997;Bryan and Fritsch, 2002; Pielke, 2002; Steppeler et al., 2003).In this article, we present the Jacobian-free Newton-Krylov (JFNK) method (Knoll and Keyes,2004; Chisholm and Zingg, 2009) for studying mesoscale penetrative convection and thermal dy-namics of the atmospheric boundary layer flow (Carpenter et al., 1990; Skamarock et al., 1997;Bryan and Fritsch, 2002; Lane, 2008). The JFNK method is increasingly considered in manybranches of computational fluid dynamics (CFD). However, it has not been a choice in major at-mospheric flow solvers except in a few academic studies (e.g. Mousseau et al., 2002; Alam andIslam, 2015). The lack of a broad acceptance of the JFNK method by the atmospheric sciencecommunity is somewhat related to the known challenge of constructing appropriate precondition-ers. If the nonlinear convection and other physical effects such as turbulence, radiation, or latentheat release are included within the matrix to be inverted at each time step, the construction ofa preconditioned JFNK method for atmospheric modelling is not fully clear from the existingliterature. The present study fills in the research gap in atmospheric modelling, where we demon-strate that preserving the symmetry of operators in their discretization can partially circumvent thepreconditioning challenge through a physics-based nonlinear preconditioning approach.To develop a JFNK solver, we discretize the tightly coupled residual of mass, momentum, andenergy of the nonhydrostatic atmospheric model equations (see Bryan and Fritsch, 2002) in a non-linearly consistent manner. To account for the lack of implicit dissipation by the skew-symmetricdiscretization of the convective operator, we show that a subgrid model is capable of dissipatingthe short-wavelength perturbations triggered by atmospheric convection. For the proposed JFNKmethod, we consider a wavelet-based approximation of differential operators, which filters theshort-wavelength perturbation. We study penetrative convection and convective boundary layer(CBL) flow over a heterogeneously heated surface. Although convection may not illustrate all ofthe computational issues of atmospheric modelling, the study of thermal dynamics indicates thatthe JFNK method offers much insight into the more complicated dynamics of atmospheric con-vection. We demonstrate that a rising thermal penetrating into a stably stratified atmosphere willeventually overshoot its level of neutral buoyancy, a crucial component of which is the generationof internal gravity waves. This overshooting involves entertainment and detrainment, which plays2 key role in atmospheric mixing and convective redistribution of heat and other scalars. Capturingsuch phenomena of atmospheric convection illustrates our understanding of the JFNK methodol-ogy in dealing with the coupled nature of atmospheric multiphysics and the fascinating nonlinearcascade of scales of atmospheric dynamics.Section 2 presents the JFNK methodology for solving the governing equations for compress-ible nonhydrostatic atmospheric flows, where a technical details of the wavelet-based discretiza-tion is outlined briefly. Section 3 summarizes the numerical results of penetrative turbulent convec-tion in the atmospheric boundary layer. We have discussed the results with respect to neutrally-and stably-stratified configurations. The test cases considered in this article are representativecases for the verification of atmospheric modelling. Results of other numerical models and fieldmeasurements have been utilized to validate the symmetry preserving JFNK methodology. Finally,Section 4 discusses the present findings and outlines how the presented methodology may furtherbe extended to advance the field of atmospheric modelling.
2. Methodology
Let us consider the dynamics of idealized dry thermals without condensation, evaporation,or any background wind shear, where the continuity, momentum, and energy equations take thefollowing form in Cartesian coordinates (see Skamarock et al., 1997; Pielke, 2002; Bryan andFritsch, 2002), ∂ p ∂ t + u j ∂ p ∂ x j = − c p c v (cid:32) ∂ u i ∂ x i (cid:33) p , (1) ∂ u i ∂ t + u j ∂ u i ∂ x j = − θ ∂ p ∂ x i − θ (cid:32) ∂ p ∂ x δ i + ∂ p ∂ y δ i (cid:33) + ∂τ i j ∂ x j + θθ g δ i , (2) ∂θ∂ t + u j ∂θ∂ x j + u j δ j β = ∂τ θ j ∂ x j . (3)The equations (1-3) are nonlinearly coupled by the convective operator. The conservative proper-ties and stability of this model are directly related to the energy contribution and the propagationspeed of atmospheric waves (Steppeler et al., 2003). The velocity u i is coupled with the non-dimensional pressure ( p ) that is referred to as the Exner function and related to the dimensionalpressure P : p = c p (cid:32) Pp (cid:33) R / c p . In Eq (3), a splitting the total potential temperature is considered, such as θ + ¯ θ + θ , where θ is a constant background temperature, β = ∂ ¯ θ/∂ z , ¯ θ ( z ) denotes a mean vertical distributionof the temperature, and θ is a temperature perturbation. This decomposition is often useful forimplementing the heat flux boundary condition at the ground (Dubois and Touzani, 2009). Note, p = R =
287 J kg − K − is the gas constant, c p = − K − is the specific heat at constant pressure, and c v =
717 J kg − K − specific heat at constantvolume. In Eqs (1-3), x i denotes the Cartesian coordinate x = ( x , y , z ), δ i j is the Kronecker delta,3 i j is the turbulent momentum flux, and τ θ j is the turbulent heat flux. Note that u i and u j denotesthe velocity components u , u , u ; however, as mentioned below, the bold-face u = [ u i , p , θ ]represents the numerical solution vector of the system (1-3). To illustrate how the symmetry of underlying physics is preserved numerically by the JFNKmethod, let u = [ u k ] be a column vector [ u i , p , θ ] k , i.e. the numerical solution of Eqs (1-3) ateach spatial grid point x k , for k = . . . N −
1, and F represent the discretization of all spatialdifferential operators involved in the system (1-3). Then, the following dynamical system ∂ u ∂ t = F u (4)represents the spatially discretized form of Eqs (1-3). Considering a time centred implicit (trape-zoidal) scheme for the dynamical system (4), we get2 u n + ∆ t − F ( u n + ) = u n ∆ t + F ( u n ) , which is a nonlinear system of algebraic equations of the following compact form L ( u n + ) = f ( u n ) . (5)For a nonlinear problem, the time centred scheme (5) with a fixed positive step size ∆ t leads to abounded error for n → ∞ , which is equivalent to A-stability of the scheme when it is applied to alinear problem. Since the order of an A-stable linear multistep method cannot exceed 2, the timecentered method is the best choice to deal with waves not contributing to energy conservation inthe solutions of atmospheric model equations (1-3) (see, e.g. LeVeque, 1990).The nonlinear convective operator u j ∂ u i /∂ ( · ) is skew-symmetric because of the property of thetrilinear form that (cid:104) u j ∂ v i , w i (cid:105) + (cid:104) w i , u j ∂ v i (cid:105) . In Eq (4), the operator F is said to be skew-symmetricwith respect to the inner product (cid:104)·(cid:105) if we have (cid:104)F u , u (cid:105) + (cid:104) u , F u (cid:105) = u . In otherwords, an anti self-adjoint operator is skew-symmetric. If the operator F is a matrix, the skew-symmetry is equivalent to F = −F T . Now, taking the inner product of u with both the sides ofEq (4) and ignoring the effects of boundary conditions, we find that ∂∂ t (cid:104) u , u (cid:105) = (cid:104)F u , u (cid:105) + (cid:104) u , F u (cid:105) . Clearly, the dynamical system (4) conserves the inner product (cid:104) u , u (cid:105) if the corresponding operator F satisfies the above skew-symmetric property. In order to satisfy the conservation of the innerproduct (cid:104) u , u (cid:105) at discrete level in the context of the dynamical system (4), we must have the innerproduct satisfying (cid:104) u n + , u n + (cid:105) = (cid:104) u n , u n (cid:105) for two consecutive time steps. It can be shown thatsuch a requirement at each time step is satisfied, subject to a truncation error of O ( ∆ t ), by thetrapezoidal time integration scheme (5) considered above if the operator F is skew-symmetric.Consider a higher-order upwind-biased discretization of convection in Eq (4), which minimizesthe local truncation error. An upwind discretization does not retain the skew-symmetry of the4onvective operator. With classical upwind methods, the eigenvalues of the operator F will havenegative real parts (see Klemp et al., 2007). The negative real part of the eigenvalues of theoperator F help ensure the stability of the system Eq (4). However, they artificially dampen theenergy (cid:104) u , u (cid:105) , and thus, the conservation of energy cannot be satisfied globally (see Klemp et al.,2007). To preserve the skew-symmetry of the convective differential operator u j ∂ u i /∂ x j with asecond order finite difference method, momentum transport equation can be discretized in thefollowing skew-symmetric form F u = ∂ u i u j ∂ x j + u j ∂ u i ∂ x j , where one considers the arithmetic mean of flux- and convective-forms of the convective operator.Notice that the skew-symmetry of the convective operator is directly related to the conservationof the convective variable. It is worth mentioning that the nonlinear stability of numerical schemesis often easier to establish if the nonlinear convective term is expressed in the skew-symmetricform. First, preserving skew-symmetry results in reduced levels of artificial dissipation, which isdesired in atmospheric simulations. Second, it eliminates the convective instability associated withspurious transfer of kinetic energy on grids that are not fine enough to resolve short-wavelengthperturbations caused by convective transport. For example, it was reported by previous researchersthat due to the artificial numerical damping of the shorter-wavelength, the Weather Research andForecasting (WRF-ARW) model is unable to adequately resolve the capping inversions. Third, itensures that the numerical dissipation of resolved kinetic energy does not overwhelm the dissipa-tion provided by a subgrid-scale model.Preserving the skew-symmetry of the convective transport by the wavelet method, in additionto the tightly nonlinear coupling of the JFNK method, brings multifold benefits discussed above. Deslaurier-Dubuc interpolating wavelets (see Deslauriers and Dubuc, 1989; Mallat, 2009) aredefined on a sequence of nested grids G j = { x k } which are embedded over nested approximationspaces V j ⊂ V j + . An element of the basis { ψ k ( x ) } of V j is presented in Fig 1 a , b . The waveletcollocation method finds an approximation u N ( x ) ∈ V j of u ( x ) ∈ L ( Ω ) such that (cid:104)L ( u N ) − f , ˜ ψ k (cid:105) = , where N is the number of grid points, L is a differential operator including the boundaryconditions, and { ˜ ψ k } is a dual basis corresponding to the approximation space V j . For the givenbasis { ψ k ( x ) } of V j , there exist a dual approximation space equipped with a basis { ˜ ψ k ( x ) } such that˜ ψ j ( ψ i ) = δ i j .The wavelet-based approximation u N ( x ) = N− (cid:80) k = u ( x k ) ψ k ( x ) (6)projects the coefficients { u ( x k ) } into V j in which the projection u N ( x ) does not oscillate at wave-lengths smaller than the grid-spacing. The discretization of differential operators are performedthrough projection of derivatives into V j . Without the details of wavelet theory, the wavelet-based5 a ) ( b ) Figure 1: ( a ) A wavelet function ψ k ( x ) satisfying ψ k ( x l ) = δ kl ; it takes a value of 1 on a given grid point x k ∈ G j and 0 on all other grid points x l ∈ G j . Then, ψ k ( x ) is extended to all grid points x k ∈ G j + by Deslauriers and Dubuc (1989) interpolation, and subsequent iterations forms a continuousfunction in G j as j → ∞ . ( b ) A restriction of ψ k ( x ) on y = ψ k ( x ).projection (see Deslauriers and Dubuc, 1989, for a technical details) of the first derivative withrespect to x is given by ∂∂ x u N ( x ) = N (cid:88) k = u (cid:48) ( x k ) ψ k ( x ) = j + p − (cid:88) k = j − p + u ( x k ) ∂∂ x ψ k ( x ) . The symmetry of differential operators are preserved due to the symmetry of ψ k ( x ). On a uniformlyrefined grid having a grid-spacing of ∆ x in all directions, the local truncation error is O ( ∆ x p ) forthe above approximation of derivatives (Alam et al., 2014). The derivative of u ( x ) is exactlyrepresented by this wavelet method if u ( x ) is a polynomial of degree 2 p −
1. The subgrid-scalemodes of half the wavelength of the resolved scale modes, which are contributed by convection u ∂ u /∂ x contributes, are explicitly filtered and parameterized by the subgrid model. In atmospheric modelling (see Pielke, 2002), the sugrid-scale schemes assume that turbulenceproduces vertical mixing in the real atmosphere, and that the role of the horizontal componentsof the subgrid scale stress is to control nonlinear aliasing errors (Pielke, 2002). Such schemes arebased on the momentum exchange coefficient (e.g. Deardorff, 1970, 1980), ν τ = ( ∆ LES C s ) (cid:113) S i j S i j , where the stresses τ i j = ν τ S i j − τ kk δ i j S i j = (cid:16) ∂ u i ∂ x j + ∂ u j ∂ x i (cid:17) of the resolved flow. Computing resources limit atmo-spheric simulations on coarse grids, where the subgrid model acts on turbulent motions that areanisotropic and intermittent (Moeng and Sullivan, 2015). Moreover, the vertical dissipation re-mains stronger than the horizontal dissipation, for example in penetrative convection (Bartello andTobias, 2013). To examine the symmetry preserving JFNK solver with respect to a basic subgridmodel, we follow the dimensional reasoning outlined by Deardorff (1970) (see Eq 3.4 therein) toestimate the horizontal and the vertical exchange coefficients separately, K m = C / s ε / ( ∆ x ∆ y ) / and K M = C / s ε / ∆ / z . (7)Here, C s is a dimensionless constant and the rate of dissipation of turbulent kinetic energy is ε = ν τ S i j S i j . According to Deardorff (1980), the horizontal eddy diffusivity is K h = (1 + l / (cid:112) ∆ x ∆ y ) K m and the vertical eddy diffusivity is K H = (1 + l / ∆ z ) K M , where l = . √ τ kk / N is a subgrid scalemixing length and N is the Brunt-V¨ais¨al¨a frequency. To ensure a minimal technical details of present contribution, we closely follow the precon-ditioned Krylov method considered by Skamarock et al. (1997) in their semi-implicit scheme forsolving the linearized nonhydrostatic atmospheric model Eqs (1-3). The JFNK method is a classof practical iterative methods for finding the solution u ∗ to the nonlinear system L ( u ) = f , Eq (5),when an initial approximation, u , is known. The nonlinear function L : R N → R N is assumeddifferentiable and J ( u ) denotes the Jacobian of L at the point u . For the nonlinear system (5), let us find variations δ u k of the solution vector u n + iterativelysuch that u n + , k + = u n + , k + δ u k for k ≥
0. This outer loop of iterations forms the Newton’smethod. The solution from the previous time step provides the first iteration u n + , = u n . At k -th Newton iteration, we minimize the residual vector R ( u n + , k ) ˙ = L ( u n + , k ) − f ( u n ) using thegeneralized minimal residual (GMRES) method of Saad and Schultz (1986). Thus, we look forthe variation δ u k satisfying the linear system J ( u ) δ u k = −R ( u n + , k ) . (8)Only a few Krlov iterations of the inner loop solve Eq (8). This ‘inexact’ Newton-Krylov methodis equivalent to solving the ordinary differential equation d u n + , k dk = −J − ( u ) R ( u n + , k )by the Euler-explicit method with a step size of one. Therefore, R ( u n + , k ) = e − k R ( u n ) if Eq (8) issolved exactly. This shows the fast rate of convergence of the inexact Newton’s method.7 .5.2. Physics-based nonlinear preconditioning Two families of preconditioning method are usually considered for the JFNK method. Thelinear preconditioning is quite similar to the Krylov method presented by Skamarock et al. (1997).For solving the linear system (8) by a preconditioned Krylov method, linear preconditioning isclassified as right- and left-preconditioning. In contrast, a physics-based nonlinear preconditioningis cost-effective thanks the wavelet-based discetization. Each diagonal block of the operator L inEq (5) means a physical field that couples with itself and each off-diagonal block means a physicalfield that couples with another field. Consider a physics-based preconditioner in which the weakcoupling of off-diagonal terms in the operator L is ignored. Physics-based preconditioner gathersthe eigenvalues of the preconditioned system in small areas, thereby increasing the convergencerate.The physics-based nonlinear preconditioning approach constructs an equivalent system of non-linear equations which provides faster rate of convergence with respect to the original system. Itcan be shown that if a fixed point iteration converges for the system (5), the eigenvalues of theJacobian ˜ J ( u ) = ∂ ˜ R ( u ) /∂ u for the preconditioned nonlinear system ˜ R ( u ) = u − P − ( u )[ f ( u n ) −L ( u ) + P ( u )] are gathered in a small area, where P ( u ) is the nonlinear preconditioning matrix.The most attractive feature of nonlinear preconditioning is that faster convergence rate of Kryloviteration is achieved with minimal mathematical and coding effort. For example, in case of imple-menting the JFNK method within an existing atmospheric modelling code, a fixed point iterationcan be performed with only a few code modification.For the matrix-vector product on the left side of Eq (8), the JFNK method needs to computethe action of the linear map J : V l ⊂ R N → R N on the variation δ u k of the solution vector u n + .To compute this action with a complexity of O ( N ), consider the Fr´echet derivative of the operator R ( u ) defined by J ( u ) δ u k = lim η → ∂∂η R ( u + ηδ u k ) . (9)The limit in Eq (9) exists, and thus, R ( u ) is Fr´echet differentiable, wherelim η → ||L ( u n + , k + ηδ u k ) − L ( u n + , k ) − J δ u k |||| ηδ u k || = . In other words, the same algorithm that provides the differentiation matrix L is applied to calcu-late the action of J on δ u k without requiring additional technical development a preconditioner.This observation suggest that the implementation of JFNK method within an existing atmosphericmodelling code is straight forward. Moreover, the complexity of the JFNK method scales like thecomplexity of the algorithm used for the discretization of Eqs (1-3), which is O ( N ) for the waveletmethod.
3. Numerical results and discussion
We report primary results on the accuracy, efficiency, and efficacy of the symmetry-preservingJFNK method as a potential candidate for problems of meteorological interest. We have stud-ied two categories of convective phenomena to test the tightly nonlinear strategy of the JFNK8 M (m s − ) K H (m s − ) Re Ra PrCase A 10 14.1 2 . × . × . × . × . × . × . × . × . × . × . × . × . × . × Table 1: List of representative parameters for the simulation of two cases in a neutral environment. In case A, theRayleigh number R a is varied with fixed Prandtl number P r = K M / K H . In case B, both the Rayleigh number R a andthe Prandtl number P r are varied. In both cases, the bulk Richardson number is fixed at R i b = .
1. The correspondingReynolds number R e is listed in this table. method. Comparisons of the present results among experimental and numerical data collectedfrom the literature have been considered. These numerical exercises indicate that the tightly non-linear physics-based coupling of all physical processes considered within the JFNK method hasthe potential to be a scale-adaptive frame-work for modelling the transition from the near-surfacesmall-scale 3D physics to the outer-layer meso-scale meteorology. We have compared the JFNK simulation of penetrative convection with the results providedby Bryan and Fritsch (2002), Wicker and Skamarock (1998), and Carpenter et al. (1990). Bryanand Fritsch (2002) examined a time-split segregated algorithm in which the convective operatorwas discretized in its flux-form without preserving its skew-symmetry. They adopted a divergencedamping term to help maintain the quality of the scheme. Carpenter et al. (1990) noted that thechoice of not preserving skew-symmetry of convective operators by the positive definite upwindschemes is to help control the Gibbs phenomenon. They also reported upwind schemes tend tosmear sharp gradient of penetrative thermals.We consider that a warm perturbation of θ is prescribed at the horizontal midpoint and at aheight of 2 km in the domain of [ − , × [0 ,
10] km , where the surrounding environment isneutrally stratified with lapse rate of 10 o K km − and θ =
300 K. The initial thermal has a radiusof 2 km (e.g. Bryan and Fritsch, 2002). As mentioned in Table 1, momentum- and heat-exchangecoefficients are varied for Case A at a fixed Prandtl number of P r = .
71. In Case B, the Prandtlnumber is varied between 0 . .
0. This test confirms hypothesis that short-wavelength modes,triggered by the nonlinear convection process in a period of time evolution of the thermal, areaccurately filtered by the subgrid model.
In Case A-B, the effect of the horizontal momentum exchange coefficient, e.g. K M = , , . , and 1 . s − , with respect to the skew-symmetry of convective operator is studied for atmo-spheric convection in a neutral environment ( i.e. the Buoyancy frequency N = M
10 m s − s − . s − . s − B & F θ min (K) -0.000632 -0.003814 -0.009359 -0.133971 -0.144409 θ max (K) 1.408749 1.629635 1.843659 2.138108 2.02178 u min (m s − ) -9.511412 -10.058257 -10.636190 -11.667357 u max (m s − ) 9.512040 10.059020 10.637147 11.668235 w min (m s − ) -6.360285 -6.527770 -6.596165 -6.627753 -8.58069 w max (m s − ) 15.352078 15.599483 15.833058 16.018170 14.5396 ω min (s − ) -0.065906 -0.095523 -0.137061 -0.188188 ω max (s − ) 0.065906 0.095523 0.137049 0.188187 Table 2: Max/min values of θ , u , w , and ω for the dry thermal simulation (Case A) at t = K M
10 m s − s − . s − θ min (K) -0.0061 -0.0062 -0.0063 θ max (K) 1.7138 1.7268 1.7371 u min (m s − ) -9.7735 -10.1486 -10.4817 u max (m s − ) 9.7744 10.1494 10.4824 w min (m s − ) -6.5080 -6.5204 -6.5481 w max (m s − ) 15.4442 15.6334 15.8046 ω min (s − ) -0.0849 -0.1053 -0.1271 ω max (s − ) 0.0849 0.1053 0.1271 Table 3: Max/min values of θ , u , w , and ω for the dry thermal simulation (Case B, see Table 1), where θ =
300 K and t = plots of the potential temperature θ in Fig 2 shows the development of two ‘rotors’ around therising thermal, which replicates the corresponding dynamics predicted by the mesoscale modelsof Wicker and Skamarock (1998) and Bryan and Fritsch (2002). In comparison to Fig 1 of Bryanand Fritsch (2002), one finds that the dynamics of penetrative thermals in a neutrally stratified dryatmosphere has been accurately simulated by the tightly nonlinear coupling strategy of the JFNKmethod. In particular, the sensitivity of the simulated dynamics on the values of K M (Table 1) isconsistent with the similar results that appeared in the literature (e.g. Carpenter et al., 1990). For aquantitative comparison, we note that the minimum and maximum potential temperature reportedby Bryan and Fritsch (2002) are θ min = − . θ max = . i b F r N (s − ) ω/ N α . × − .
999 2 . o . × − .
997 4 . o . × − .
985 9 . o . × − .
977 12 . o . × − .
962 15 . o . × − .
929 21 . o Table 4: Relation between the wave frequency and buoyancy frequency for penetrative convection in a stably stratifiedenvironment. sidered values of K M = . s − . However, the vertical velocity is predicted more accurately witha higher value of K M = . s − . It is also evident from Table 2 that the numerical predictionsare not noticeably sensitive to changes in momentum exchange coefficients. One observes that theupper surfaces of thermals at t = . , . , and 8 .
15 kmfor K M = , , and 2 . s − , respectively. From the predicted maximum vertical velocity inTable 2, we see that the rate of vertical momentum transfer in a turbulent penetration of dry ther-mals may be weakly sensitive to the subgrid-scale mixing length l provided by the eddy diffusivitymodel of Deardorff (1980).Color-filled contour plots of the horizontal and the vertical velocities are presented in Figure 3.Notice that the velocity field is symmetric about x = P r was also varied, are similar to that for Case A, as reported in Table 3. Corresponding to the neutral simulations (Case-A and Case-B), penetrative convection in astably stratified environment is simulated, where the model was initialized for the potential tem-perature θ + β z + θ ( x , z ). In this situation, the frequency of internal wave ( ω ) correlates with thebuoyancy frequency ( N ) through the dispersion relation ω = N cos α . Values of the buoyancyfrequency N and α are listed in Table 4.The linear theory suggests the existence of evanescent and vertically propagating waves (seeLin, 2007, p. 187) if 2 π N > LU and 2 π N < LU , respectively. Moreover, if π N >> LU , the buoyancy force becomes extremely weak so that thevertical velocity can be estimated by (see Lin, 2007) w ( x , z ) = W ( x ) e − k | z | . Here, U and W are horizontal and vertical velocity scales, respectively, and L denotes a horizontallength scale. Simulated vertical velocity w (0 , z , t ) and potential temperature θ (0 , z , t ) at t = a ) ( b )( c ) Figure 2: Contour plots of the potential temperature perturbation ( θ ) for Case A with a neutral environment at P r = .
71 and t = K M =
10 m s − ; (b) K M = s − , and (c) K M = . s − . a ) ( b )( c ) ( d )( e ) ( f ) Figure 3: Color filled contour plots of the horizontal velocity ( u , left column) and the vertical velocity ( w , rightcolumn) for Case A with a neutral environment at t = K M =
10 m s − , (c,d) K M = s − and (e,f) K M = . s − . Red, blue and yellow colors represent positive, negative and zero values, respectively. U ∼
10 m s − (based on Table 2) and the buoyancy frequency, N ∼ . × − s − , one finds a vertical length scale of L z = π UN ∼ .
25 km. Thus, we recordedthe vertical velocity at a location on the center of the horizontal domain at z = .
25 km for everytime step with respect to four values of the buoyancy frequency, N = . × − , . × − ,1 . × − , . × − , and the result is shown in Fig 5. For each N , the corresponding bulkRichardson numbers R i b , as well as the Froude number F r are presented in Table 4. For the resultin Fig 5, the ratio of the wave frequency to the buoyancy frequency ( ω/ N ) at z = .
25 km arereported in Table 4. The angle ( α ) between the phase velocity vector and the horizontal directionare also presented in Table 4, suggesting the dispersion relation: ω = N cos α. Table 4 also suggests that the wave frequency satisfies ω < N , indicating the maximum possiblefrequency of internal waves in a stratified fluid is N . The angle ( α ) increases as buoyancy fre-quency decreases. To illustrate internal wave propagation, Fig 6 presents the time series of thevertical velocity corresponding to seven locations (0 , . , . , . , . , . , .
75) for three values of N . These results show a good agreement of the phenomena sim-ulated by the JFNK method with the corresponding findings reported in the literature (e.g. seeMorton et al., 1956), which means that the wavelet-based JFNK model accurately predicts thepenetrative convection of thermals in a stably stratified environment. The governing equations Eq. (1-3) leads to the following energy balance laws (e.g. Wintersand Young, 2009), where the kinetic and potential energies E k = (cid:90) Ω ( u + w )d V , and E p = (cid:90) Ω ( z max − z ) θ d V , respectively, satisfy the following energy equations (see also Carpenter et al., 1990, e.g. Fig 19therein): d E k d t = (cid:90) Ω w θ d V − K M (cid:15) , (cid:15) = (cid:90) Ω | ∆ u | + | ∆ w | d V and d E p d t = − (cid:90) Ω w θ d V + K H θ max − θ min z max − z min . These energy equations quantify the rate of production of E p , the conversion from E p to E k , andthe rate of kinetic energy dissipation, (cid:15) , thereby making a steady energy balance for the isolatedthermal in the neutral environment.The time evolution of E p , E k , and E = E p + E k for a rising thermal in the neutral environ-ment have been reported in Fig 7a. Clearly, the potential energy E p , decreases with time as aresult of the potential energy conversion into kinetic energy E k . Also, the total energy, E , re-mains approximately constant. In order to conserve energy, Carpenter et al. (1990) considered14 a ) vertical velocity( b ) perturbation potential temperature Figure 4: Plots of (a) vertical velocity (b) potential temperature probed along the vertical line x = t = K M =
10 m s − and P r = . a) N = . × − s − (b) N = . × − s − (c) N = . × − s − (d) N = . × − s − Figure 5: Evolution of vertical velocity at x = z = .
25 km for the stable case (a) N = . × − s − , (b) N = . × − s − , (c) N = . × − s − and (d) N = . × − s − . a) Ri b = . , N = . × − s − (b) Ri b = . , N = . × − s − (c) Ri b = . , N = . × − s − Figure 6: Time series of vertical velocity at x = z = .
25, 2 .
5, 3 . , . . , . , and 8 .
75 km for three values of N a) (b) Figure 7: (a) Time evolution of energy for Case A, showing that the total energy is approximately conserved, wherepotential energy is converted to kinetic energy. Note, θ =
300 K, K M = s − , t ∈ [0 , s ]. (b) The effect ofstratification on the time evolution of energy for Case A, but in the stable environment with N = . × − s − for t ∈ [0 , s ]. the piece-parabolic method for the discretization of convective operators. The present result ofthe energy conservation depicted in Figure 7a has an excellent agreement with the correspondingresult reported by Carpenter et al. (1990).However, the time evolution of the potential and kinetic energies feature more complex andoscillating behaviour when the environment is stability stratified. It is because a rising thermalfinds itself in an environment with a higher potential temperature, where the buoyancy force pushesit downwards. In Fig 7b, the energy curves for N = . × − s − have been displayed. Theresults clearly indicates the overshooting of thermals beyond their level of buoyancy in the stablystratified environment (see also Lane, 2008). There is a growing trend of land-surface modification through urbanization because over halfthe world’s population lives in urban areas. Urban Heat Island (UHI) is the source of the mesoscaleresponse of the atmosphere to horizontal variations in temperature associated with dry convec-tion (Grimmond and Oke, 2002). Urban Heat Island (UHI) is a potential atmospheric model toinvestigate the influence of land-surface modification on the health and welfare of urban residents.UHI simulations help quantify mesoscale circulation triggered by surface heterogeneity of urban-ization (see Zhang et al., 2014).
Using the LES mode of the Weather Research and Forecasting (WRF-LES) model, Zhanget al. (2014) investigated the influence of the UHI circulation over an isolated urban area that is18 able 5: The effect of surface heat flux variation on the potential temperature θ min , horizontal velocity u max , verticalvelocity w max , and vorticity ω max has been predicted accurately by the present model in comparison to the results ofDubois and Touzani (2009), i.e. D&T H .
93 W m − .
87 W m − .
74 W m − Present D & T Present D & T Present D & T θ min -0.023537 -0.024823 -0.064457 -0.071289 -0.167264 -0.166316 u max w min -0.030134 -0.030470 -0.037337 -0.039291 -0.085519 -0.079265 ω max y -direction, where the atmosphere is dry and the terrain is flat. The WRF-LESmodel of Zhang et al. (2014) is similar to the simulation of Dubois and Touzani (2009) (hereinafterD & T). Similarly, Kimura (1975) provides a laboratory model of an equivalent UHI circulation. Inthis section, these three reference models are considered to understand the accuracy of the wavelet-based simulation, where the model domain extends 100 km horizontally and 2 km vertically.Following Zhang et al. (2014) a constant heat flux ( H surface [W m − ]) and u = w = z =
0, such that H surface − H rural = (cid:40) x < − l or x > lH tanh (cid:16) x + ξ (cid:17) − H tanh (cid:16) x − ξ (cid:17) if − l ≤ x ≤ l . , where H rural is the surface heat flux over the rural area and H is the surface heat flux over theurban area (see Fig 8 a ).Values of H = .
93, 57 . , . , . , .
96, and 926 .
92 [W m − ] were tested. Thesesix values of H correspond to six values of Rayleigh numbers R a = , 10 , , , , and 10 . Fig 8 demonstrates the vertical profiles of the vertical velocity and the potential temperaturecomputed at the centre of the heat island for H = . , .
87 and 115 .
74 Wm − . In dimen-sionless variables, these plots are found in a very good agreement with the corresponding plotsof Dubois and Touzani (2009). The temperature, vertical velocity, and horizontal velocity de-cay rapidly with respect to the elevation z . In Figure 8b and 8c, it is clear that the mixed layerheight appears between z = . . H increases. Thevalues of the horizontal velocity u , the vertical velocity w , the potential temperature θ , and thevorticity ω y = ∂ u /∂ z − ∂ w /∂ x were compared between the results of the present model and thatof Dubois and Touzani (2009). The results presented in Table 5 indicates that the accuracy of thewavelet-based JFNK model is within 5% to 10% of the results of D & T for the test case of UHI. It was observed in the experimental study of Kimura (1975) that the centre of the heat islandcirculation is located near the edges of the heat island when the surface heat flux is relatively weak,19 a ) θ Stable atmosphere θ ¯ θ = θ + β z Heated region − L x / km L x / km − l km l km z = L z km z = h ~ 1 km z = km m i x ed l a y e r i n i t i a l t e m pe r a t u r e ( b ) ( c )( d ) Figure 8: (a) A sketch of heat island circulation and the initial profile of the total potential temperature. Profiles of (b)temperature variation θ at x = x = x = . H = . , . .
74 W m − , respectively. )c) d)b) Figure 9: (Left) Two types of flow regimes observed in a laboratory experiment (Kimura, 1975). (a) Center of thecirculation is located near the edge of the heating element; (c) the center of the circulation is near the center of theheating element. (Right) Numerically experimented results: (b) H = .
93 W m − , β = − (d) H = .
74W m − , β =
10 K km − . Thick black horizontal line represents the heating element. Note, the center of the circulationapproximately corresponds to each other in ( a − b ) and ( c − d ). The thick arrow on the left panel is about z = and the up-draft prevails all over the urban area (e.g Fig. 9a). On the other hand, a strong narrowup-draft is concentrated above the centre of the island, when the heat flux is relatively strong (e.g.Fig. 9c).The simulation results in Fig 9b and 9d are in a very good agreement with the correspondingexperimental results. It was found that if H increases, the center of the circulation moves towardthe center of the urban region. In Fig. 9b, the surface heat flux is 28 .
93 W m − , where the centerof the circulation is away from the center of the heat island, and in Fig. 9d, the surface heat flux is115 .
74 W m − , where the center of the circulation is at the center of the heat island. To provide a primary assessment of the model for predicting the structure of a convectiveboundary layer (CBL), an idealized dry case is studied, which is driven solely by a surface heatflux. The result is analyzed with respect to the field measurements of the day-33 Wangara experi-ment (e.g. Moeng, 1984), where the temporal evolution of the mixed layer is studied. In a similarstudy, Moeng (1984) considered LES to reproduce Wangara data when the numerical model in-corporates a prameterization of the moisture field, radiation effect, Coriolis effect, and surfaceroughness in addition to having the wind shear consistent with the field measurements. Mukherjeeet al. (2016) provides an idealized numerical study of the CBL using a turbulence-resolving LESin the domain 3 . × . × . .In comparison to the reference work mentioned above, the present study considers a simplifiedcase in a two-dimensional domain. Here, the horizontal length of 100 km is divided into 512segments ( ∆ x ≈
195 m), and the vertical length of 2 km is divided into 256 segments ( ∆ z ≈ ) Z ( k m ) potential temperature, ⟨ θ ⟩⟨ θ ⟩ Observed, ( C ) Figure 10: The vertical profile of the average potential temperature (cid:104) ¯ θ (cid:105) . The temporal average is obtained in the timeinterval of t ∈ [0 h, 6 h]. The prediction of the JFNK model at H = .
92 W m − is compared with the fieldmeasurement (by digitally extracting the data from Fig 3b of Moeng (1984)), which represents the vertical profile of (cid:104) ¯ θ (cid:105) at 12-th hour of day 33. The hourly evolution of the vertical temperature profile is also shown for clarity. Note forthe unit conversion to o C, which aims for being consistent with Moeng (1984).
Horizontally (cid:104)·(cid:105) and temporally ( · ) averaged vertical profile of the total potential temperature, i.e. (cid:104) ¯ θ (cid:105) = (cid:104) θ + β z + θ (cid:105) , for H = .
96 W m − and H = .
92 W m − are analyzed. We observedthat the estimated mixed layer height, for H = .
96 W m − and H = .
92 W m − , are about0 . . . − . . − . (cid:104) ¯ θ (cid:105) ( z ) of the present simulation with the sim-ilar profile observed in the Wangara day-33 experiment ( e.g. , Moeng, 1984). The displayed datacorresponds to the surface heat flux of H = .
92 W m − . We can see the development of awell-mixed layer within 2 h. The mixed layer in the turbulent region is capped with the inversionlayer approximately at z ∼
800 m. At the bottom of the boundary layer, the unstable surface layerappears at ∼
50 m. The mixed layer height did not fully agree between the simulation and themeasurement because the idealized simulation ignored the realistic meteorological features, suchas wind-shear, Coriolis effect etc. However, the vertical structure of the daytime boundary layerhas been predicted with a good quality. 22 . Summary and future developments
A symmetry preserving JFNK method for accurate simulations of nonlinearly coupled atmo-spheric physics has been illustrated in this article. The method has been applied to simulate at-mospheric convection using the nonhydrostatic atmospheric model equations. Unlike linearizedmethods commonly adopted in atmospheric flow solvers, the JFNK method emphasizes that non-linear coupling of physics to be modelled without any compromise. In atmospheric flow simula-tions, capturing the tight nonlinear coupling may lead to a next generation atmospheric model tocorrectly forecast weather events. In this work we have shown that preserving the skew-symmetryof the convective operator by the wavelet method brings two modelling benefits. First, short-wavelength modes are properly cascaded toward the subgrid scale dissipation mechanism. Thisis done by ensuring the role of the convective operator it would play on the physics of the flow.Second, the subgrid-scale modes are properly transferred to the subgrid model for being dissipatedat a rate offered by the subgrid model.In mathematical terms, an interplay between the skew-symmetric convection and symmetric,negative-definite diffusion leads to small-scale motion in a turbulent flow. With this hypothesisof energy cascade in mind, we have combined the JFNK method with a symmetry-preservingdescretization that is based on the wavelet method. This article presents the efficiency and reliabil-ity of the JFNK methodology for numerical simulation of nonhydrostatic atmospheric flows in thecontext of penetrative convection. We have chosen to simulate idealized dry convection for the pre-sentation of the JFNK method because the evolution of plumes and thermals offers much insightinto the more complicated dynamics of atmospheric motion. Because of the tightly nonlinear cou-pling of atmospheric motions, present authors envision to perform simulations in such a mannerthat the discretized operators preserve the same symmetry properties and the same nonlinear cou-pling as the underlying differential operators. The main question becomes whether the symmetry-preserving nonlinearly coupled discretization is appropriate for atmospheric simulations since theatmospheric modelling community has accepted the discretization of the skew-symmetric con-vective opterator to a positive-definite convective operator in their publicly available codes (seeKlemp et al., 2007; Smolarkiewicz et al., 2014, 2017, and the refs therein). This question has beenaddressed by simulations for which conservation of energy is highly desired to correctly forecastintensification of particular weather events (see Carpenter et al., 1990; Klemp et al., 2007).The tight nonlinear coupling of physics offered by the JFNK solver will bring full benefit toatmospheric simulations if a scale-adaptive subgrid model is considered (see Alam and Fitzpatrick,2018). Transition of scales – often labelled as the ‘gray-zone’ – is currently one of the mostchallenging problems in the field of meteorology (Wyngaard, 2004; Kurowski and Teixeira, 2018).The findings of this article encourages to further test the JFNK method using a more appropriatescale-adaptive subgrid model that adapts the cut-off scale dynamically as the characteristic scaleexhibits transitions. For example, a balance between the local production and the dissipation ofturbulence occurs in the surface-layer at much smaller scales than the characteristic scale of eddiesin the outer layer. In the future research, we plan on developing wavelet-based preconditionersbased on lifting schemes, which would discretize the differential operators on to a hierarchy of‘details’ wavelet space V j + \V j instead of the approximation space V j considered in the presentwork. 23 cknowledgements References
J. Ferziger, M. Peric, Computational Methods for Fluid Dynamics, Springer, New York, 1997.R. A. Pielke, Mesoscale Meteorological Modeling, Academic Press, 2nd edition, 2002.R. H. Pletcher, J. C. Tannehill, D. A. Anderson, Computational Fluid Mechanics and Heat transfer, Taylor & Francisgroup, third edition, 2013.L. J. Wicker, W. C. Skamarock, A time-splitting scheme for the elastic equations incorporating second-order runge-kutta time differencing., Mon. Wea. Rev. 126 (1998) 1992–1999.W. C. Skamarock, P. K. Smolarkiewicz, J. B. Klemp, Preconditioned conjugate-residual solvers for helmholtz equa-tions in nonhydrostatic models, Monthly Weather Review 125 (1997) 587–599.R. Verstappen, A. Veldman, Symmetry-preserving discretization of turbulent flow, Journal of Computational Physics187 (2003) 343 – 368.J. Steppeler, R. Hess, U. Sch¨attler, L. Bonaventura, Review of numerical methods for nonhydrostatic weather predic-tion models, Meteorology and Atmospheric Physics 82 (2003) 287–301.J. B. Klemp, W. C. Skamarock, J. Dudhia, Conservative split-explicit time integration methods for the compressiblenonhydrostatic equations, Monthly Weather Review 135 (2007) 2897–2913.G. H. Bryan, J. M. Fritsch, A benchmark simulation for moist nonhydrostatic numerical model, Mon. Wea. Rev. 130(2002).D. A. Knoll, D. E. Keyes, Jacobian-free newton-krylov methods: a survey of approaches and applications, J. Comput.Phys. 193 (2004) 357–397.T. T. Chisholm, D. W. Zingg, A jacobian-free newtonkrylov algorithm for compressible turbulent fluid flows, Journalof Computational Physics 228 (2009) 3490 – 3507.R. L. J. Carpenter, K. K. Droegemeier, P. R. Woodward, C. E. Hane, Application of the piecewise parabolicmethod(ppm) to meteorological modelling., Mon. Wea. Rev. 118 (1990) 586–612.T. P. Lane, The vortical response to penetrative convection and the associated gravity-wave generation, Atmos. Sci.Lett. 9 (2008) 103–110.V. A. Mousseau, D. A. Knoll, J. M. Reisner, An implicit nonlinearly consistent method for the two-dimensionalshallow-water equations with coriolis force, Monthly Weather Review 130 (2002) 2611–2625.J. Alam, M. R. Islam, A multiscale eddy simulation methodology for the atmospheric ekman boundary layer, Geo-physical & Astrophysical Fluid Dynamics 109 (2015) 1–20.T. Dubois, R. Touzani, A numerical study of heat island flows: Stationary solutions, International Journal forNumerical Methods in Fluids 59 (2009) 631–655.R. J. LeVeque, Numerical Methods for Conservation Laws, Birkhauser Verlag, Boston, 1990.G. Deslauriers, S. Dubuc, Symmetric iterative interpolation process., Constructive Approximation 5 (1989) 49–68.S. Mallat, A wavelet tour of signal processing, Academic Press, 2009.J. M. Alam, R. P. Walsh, M. Alamgir Hossain, A. M. Rose, A computational methodology for two-dimensional fluidflows, International Journal for Numerical Methods in Fluids 75 (2014) 835–859.J. W. Deardorff, A three-dimensional numerical investigation of idealized planetary boundary layer, Geophys. FluidDyn. 1 (1970) 377–410.J. W. Deardorff, Stratocumulus-capped mixed layers derived from a three-dimensional model, Boundary-LayerMeteorology 18 (1980) 495–527.C.-H. Moeng, P. Sullivan, Encyclopedia of Atmospheric Sciences, 2nd Edition, volume 4, Elsevier Ltd, AcademicPress, pp. 232–240. . Bartello, S. M. Tobias, Sensitivity of stratified turbulence to the buoyancy reynolds number, Journal of FluidMechanics 725 (2013) 122.Y. Saad, M. H. Schultz, Gmres: a generalized minimal residual algorithm for solving nonsymmetric linear systems,SIAM J. Sci. Stat. Comput. 7 (1986) 856–869.Y.-L. Lin, Mesocale Dynamics, Cambridge University Press, 2007.B. R. Morton, G. Taylor, J. S. Turner, Turbulent and gravitational convection from maintained and instaneous sources,Proceedings of the Royal Society of London A234 (1956) 1–23.K. B. Winters, W. R. Young, Available potential energy and buoyancy variance in horizontal convection, Journal ofFluid Mechanics 629 (2009) 221–230.C. S. B. Grimmond, T. R. Oke, Turbulent heat fluxes in urban areas: observations and a local-scale urban meteoro-logical parameterization scheme (LUMPS), J. Appl. Meteor 41 (2002) 792–810.N. Zhang, X. Wang, Z. Peng, Large-eddy simulation of mesoscale circulations forced by inhomogeneous urban heatisland, Boundary-Layer Meteorology 151 (2014) 179–194.R. Kimura, Dynamics of steady convections over heat and cool islands, Journal of the Meteorological Society ofJapan 53 (1975) 440–457.C.-H. Moeng, A large-eddy-simulation model for the study of planetary boundary-layer turbulence, Journal of theAtmospheric Sciences 41 (1984) 2052–2062.S. Mukherjee, J. Schalkwijk, H. J. J. Jonker, Predictability of dry convective boundary layers: An les study, Journalof the Atmospheric Sciences 73 (2016) 2715–2727.P. K. Smolarkiewicz, C. K¨uhnlein, N. P. Wedi, A consistent framework for discrete integrations of soundproof andcompressible pdes of atmospheric dynamics, J. Comput. Phys. 263 (2014) 185–205.P. K. Smolarkiewicz, C. K¨uhnlein, W. W. Grabowski, A finite-volume module for cloud-resolving simulations ofglobal atmospheric flows, Journal of Computational Physics 341 (2017) 208 – 229.J. M. Alam, L. P. J. Fitzpatrick, Large eddy simulation of urban boundary layer flows using a canopy stress method,Computers & Fluids 171 (2018) 65–78.J. C. Wyngaard, Toward Numerical Modeling in the ”Terra Incognita”, Journal of the atmospheric sciences 3 (2004)1816–1826.M. J. Kurowski, J. Teixeira, A scale-adaptive turbulent kinetic energy closure for the dry convective boundary layer,Journal of the Atmospheric Sciences 75 (2018) 675–690.. Bartello, S. M. Tobias, Sensitivity of stratified turbulence to the buoyancy reynolds number, Journal of FluidMechanics 725 (2013) 122.Y. Saad, M. H. Schultz, Gmres: a generalized minimal residual algorithm for solving nonsymmetric linear systems,SIAM J. Sci. Stat. Comput. 7 (1986) 856–869.Y.-L. Lin, Mesocale Dynamics, Cambridge University Press, 2007.B. R. Morton, G. Taylor, J. S. Turner, Turbulent and gravitational convection from maintained and instaneous sources,Proceedings of the Royal Society of London A234 (1956) 1–23.K. B. Winters, W. R. Young, Available potential energy and buoyancy variance in horizontal convection, Journal ofFluid Mechanics 629 (2009) 221–230.C. S. B. Grimmond, T. R. Oke, Turbulent heat fluxes in urban areas: observations and a local-scale urban meteoro-logical parameterization scheme (LUMPS), J. Appl. Meteor 41 (2002) 792–810.N. Zhang, X. Wang, Z. Peng, Large-eddy simulation of mesoscale circulations forced by inhomogeneous urban heatisland, Boundary-Layer Meteorology 151 (2014) 179–194.R. Kimura, Dynamics of steady convections over heat and cool islands, Journal of the Meteorological Society ofJapan 53 (1975) 440–457.C.-H. Moeng, A large-eddy-simulation model for the study of planetary boundary-layer turbulence, Journal of theAtmospheric Sciences 41 (1984) 2052–2062.S. Mukherjee, J. Schalkwijk, H. J. J. Jonker, Predictability of dry convective boundary layers: An les study, Journalof the Atmospheric Sciences 73 (2016) 2715–2727.P. K. Smolarkiewicz, C. K¨uhnlein, N. P. Wedi, A consistent framework for discrete integrations of soundproof andcompressible pdes of atmospheric dynamics, J. Comput. Phys. 263 (2014) 185–205.P. K. Smolarkiewicz, C. K¨uhnlein, W. W. Grabowski, A finite-volume module for cloud-resolving simulations ofglobal atmospheric flows, Journal of Computational Physics 341 (2017) 208 – 229.J. M. Alam, L. P. J. Fitzpatrick, Large eddy simulation of urban boundary layer flows using a canopy stress method,Computers & Fluids 171 (2018) 65–78.J. C. Wyngaard, Toward Numerical Modeling in the ”Terra Incognita”, Journal of the atmospheric sciences 3 (2004)1816–1826.M. J. Kurowski, J. Teixeira, A scale-adaptive turbulent kinetic energy closure for the dry convective boundary layer,Journal of the Atmospheric Sciences 75 (2018) 675–690.