Response to Topical editor. Implementation of an Immersed Boundary Method in the Meso-NH model: Applications to an idealized urban-like environment

This study describes the numerical implementation, verification and validation of an Immersed Boundary Method (IBM) in the atmospheric solver Meso-NH for applications to urban flow modelling. The IBM represents the fluid-solid interface by means of a LevelSet Function and models the obstacles as part of the resolved scales. The IBM is implemented by means of a three-steps procedure: first, an explicit-in-time forcing is developed based on a novel 5 Ghost-Cell Technique that uses multiple image points instead of the classical single mirror point. The second step consists in an implicit step projection where the right-hand side of the Poisson equation is modified by means of a Cut-Cell Technique to satisfy the incompressibility constraint. The condition of non-permeability is achieved at the embedded fluid-solid interface by an iterative procedure applied on the modified Poisson equation. In the final step, the turbulent fluxes and the wall model used for Large-Eddy-Simulations (LES) are corrected and a wall model is proposed to ensure consistency of the subgrid scales with 10 the IBM treatment. In the second of part of the paper, the IBM is verified and validated for several analytical and benchmark test cases of flows around single bluff bodies with increasing level of complexity. The analysis showed that MNH-IBM reproduces the expected physical features of the flow, which are also found in the atmosphere at much larger scales. Finally the IBM is validated in the LES mode against the Mock Urban Setting (MUST) field experiment, which is characterized by strong roughness caused by 15 the presence of a set of obstacles placed in the atmospheric boundary layer in nearly-neutral stability conditions. The Meso-NH IBM-LES reproduces with reasonable accuracy both the mean flow and turbulent fluctuations observed in this idealized urban environment.


Introduction
Once more we thank the Topical Editor for his interest giving some suggestions to improve the manuscript. We kindly revise the manuscript accordingly. A point-by-point reply to the comments is done in the next section. The reply is supported by the revised and marked-up manuscript versions. The red color is used in the marked-up manuscript version to highlight the corrections.  15 presented in the core of the paper. The classification of these input files (preparation and execution) respect the one of the paper sections.
In addition, the tar ball https://cerfacs.fr/MNHIBM/Auguste-GMD-2019.tar has been created. This file can be directly download on the webpage. The size of the tar ball (1Mo) is reasonable and therefore can be a supplementary material of this article.
Implementation of an Immersed Boundary Method in the Meso-NH v5.2 model: Applications to an idealized urban environment.

Introduction
Urbanization impacts the physical and dynamical structure of the atmospheric boundary layer, influencing both the local 20 weather and the concentration and residence time of pollutants in the atmosphere, which in turn impact air quality. While the physical mechanisms driving these interactions and their connections to climate change are well understood (the Urban Heat Island effect, anthropical effects), their precise quantification remains a major modelling challenge. Accurate predictions of these interactions require modelling and simulating the underlying fluid mechanics processes to resolve the complex terrain featured in large urban areas, including buildings of different sizes, street canyons, parks, etc. For example, it is well known that pollution originates from traffic and industry in and around cities, but the actual dispersion mechanisms are driven by the local weather. Furthermore, fine-scale flow fluctuations influence nonlinear physicochemical processes. The present study 5 addresses these issues focusing on the numerical aspects of the problem.
With the progress in metrology, it is now possible to obtain reliable measurements of the atmospheric conditions over a city.
For example, during the Joint Urban experiment (JU2003), scalar dispersion was measured experimentally over the streets of Oklahoma City (Clawson et al., 2005;Liu et al., 2006). Similarly, the CAPITOUL experiment (2004-2005) conducted in Toulouse, analyzed the turbulent boundary layer developed over the urban topography and evaluated the energy exchanges 10 between surface and atmosphere Hidalgo et al., 2008). More recently, the multiscale field study by Allwine et al. (2012) provided meteorological observations and tracer concentrations data in Salt Lake City. Other studies analyzed reduced-scale and/or idealized models to understand urban climate features as in the COSMO (Comprehensive Outdoor Scale Model Experiment for Urban Climate) and Kugahara projects (Moriwaki and Kanda, 2004). For example, Kanda et al. (2007) and Wang et al. (2015) used, respectively, an array of cubic bodies and stones fields as small-scale models. 15 In order to use in the future these experimental data for model validation, the numerical models need first to be verified for academic test cases and simplified scenarios representative of atmospheric turbulent boundary layer flows. In particular, the flow interaction with buildings or any generic obstacles plays a crucial role in urban flow modelling. The range of scales of objects acting as obstacles is huge in urban setting, encompassing large buildings and small vegetation scales and so is the range of the corresponding flow-obstacles interactions. Covering all possible cases is obviously impossible but we can rely on 20 the invariance of certain flow characteristics at different scales. For example, the von Kármán streets are observed in the wake of a centimeter-scale cylinder as well as in the cloud layout behind an island. Following this, a wide selection of benchmark flows can be analyzed to verify and validate the numerical treatment of fluid-obstacle interaction with a view to atmospheric applications.
The physical application considered in this work is the atmospheric mesoscale reaction to perturbations induced by urban 25 areas and the more the obstacles are considered as a part of the scales numerically resolved the more the results accuracy is. To access this resolution, this study presents the development, implementation, verification and validation of an Immersed Boundary Method IBM (Mittal and Iaccarino, 2005) in the Meso-NH model MNH (Lafore et al., 1998;Lac et al., 2018) for applications to urban flow modelling 1 . This choice was dictated by the fact that numerical solvers in MNH enforce conservation on structured grids and hence cannot handle body fitted grids with steep topological gradients. The main idea behind IBM is 30 the detection of an interface separating a fluid region (where conservations laws hold) from a solid region (corresponding to the obstacle volume) using different techniques (e.g. markers, LevelSet functions, local volume fraction, etc). Two main classes of IBM exist based on the continuous and discrete forcing approaches, respectively. The continuous forcing approach was developed by Peskin (1972) for biomechanics applications and consists in the addition of a continuous artificial force (acceleration indeed) in the momentum conservation equation that mimic the effect of the obstacles (heart linings) and drive the flow to relax to no-slip conditions at the wall of the obstacles. This approach and its variant developed by Goldstein et al. (1993) for a rigid interface can suffer from the lack of stiffness (fluid-solid interface is generally spread over few cells) and the time step restriction (spring and damper model with large natural frequency). Nevertheless, the continuous forcing approaches are very successful in many applications (penalization method as in Angot et al. (1999), fictitious domain method, etc). In the 5 second IBM class, the discrete approach, the boundary conditions are specified at the immersed interface. To simulate flows around non moving and rigid bodies, two sub-classes of discrete approaches can be defined as in Mittal and Iaccarino (2005): direct or indirect approaches, depending on the forcing location (Pierson, 2015). Many types of discrete forcing exist, e.g. direct forcing in the fluid region near the interface as in Mohd-Yusof (1997), immersed interface method (Leveque and Li, 1994), Cartesian grid method (Clarke et al., 1986). Depending on how to resolve the partial differential equations, Cartesian 10 grid methods (Ye et al., 1999) are written for finite-volume discretizations (Cut-Cell Technique, CCT) and for finite-difference discretizations (Ghost-Cell Technique, GCT) as in Tseng and Ferziger (2003). CCT reshapes the cell cut by the interface to preserve mass, momentum and energy. Using GCT, the local spatial reconstruction is done inside the solid region. Note that the latter technique has been successfully implemented in Weather Research and Forecasting WRF model (Lundquist et al., 2010(Lundquist et al., , 2012. 15 In this study, a discrete forcing approach is adopted where the fluid-solid interface is modelled by means of a LevelSet function (Sussman et al., 1994). The motivation behind this choice is that we are primality interested in modelling explicitly rigid and non-moving bodies in a turbulent flow, and with sufficiently fine resolution to avoid the large dissipation inherent to the presumed spread interface. The GCT does not introduce source terms in the conservation equations modelling the fluid region so that boundary conditions are imposed at the interface and/or in the solid region. The only corrections to the physical 20 model in the fluid region come from subgrid turbulent parameterizations. The idea is that in future mesoscale applications, IBM will be used to resolve large obstacles (in the solid region) such as buildings but also mountains, whereas a subgrid drag model will be used to handle unresolved obstacles such as vegetation (Aumond et al., 2013).
The paper is organized as follows. Section 2 briefly describes the general features of MNH. Section 3 details the numerical implementation of the IBM. Inspired by the works of Bredberg (2000), Piomelli and Balaras (2002), Craft et al. (2002) and to 25 close the turbulence problem, an immersed wall model is proposed in Section 3.3. Sections 4.1 and 4.2 describe the validation of the method for academic flows, respectively potential (Lamb, 1932;Batchelor, 2000) and viscous flows. Finally, Sections 5.1 and 5.2 discuss the results of turbulent flows simulations and comparisons with data from field experiments. Conclusions are drawn in Section 6. Additional tests and validations on potential and inviscid flows are respectively documented in the Appendix A and B. The study of a viscous and thermodynamic case (Straka et al., 1993) is given in supplementary materials. 30 2 The Meso-NH code at a glance MNH is an atmospheric non hydrostatic research model. Its spatio-temporal resolution ranges from the large meso-alpha scale (hundred of kilometers and days) down to the micro-scale (meters and seconds). It is massively parallel on nested and structured grids, adapted on most of international hosting computer platforms. Several parametrizations are available: radiation, turbulence, microphysics, moist convection with phase change, chemical reactions, electric scheme, externalized surface scheme. In the present study, only two subgrid parametrizations are approached: turbulence and surface schemes.

The conservation laws
The spatial discretization x is based on the terrain-following coordinates (Gal-Chen and Somerville, 1975). The staggered 5 mesh is regular ∆x = ∆y = ∆ in the horizontal directions and a transformation of the vertical one is available in order to fit a non-plane surface. The release of the vertical space step is available where a fine resolution is unnecessary. In the current study, only flat problems are considered with a ∆z = ∆ restriction for altitudes in presence of immersed obstacles.
The core of the MNH dynamic in its dry version is based on the resolution of the Euler and thermodynamic equations (energy preserving). The anelastic approximation (Lipps and Hemler, 1982;Durran, 1989) is assumed; the reference state is stratified 10 and the density deviation to the hydrostatic case in the buoyancy term is considered. The system can be simplified into the Boussinesq approximation when considering an uniform reference state. The tendencies of each prognostic variable ψ satisfying the usual conservation laws in MHN are expressed as ∂ψ ∂t csl where subscript csl is used to distinguish these tendencies from those arising from the application of the IBM procedure (Sect. 3.1). The prognostic variables are the resolved momentum, the potential temperature and if necessary an arbitrary passive scalar. The prognostic variable is decomposed into a resolved 15 component, ψ, and an unresolved component, ψ ′ (ψ ′ = 0 in a Direct Numerical Simulation, DNS). An additional prognostic equation on the subgrid turbulent kinetic energy is solved for a Large Eddy Simulation (LES, Sect. 2.3). The potential temperature is defined through the Exner function Π and the absolute temperature is the density of the reference state, P is the local pressure, P r0 the reference ground level pressure, R d the gas constant and C p the specific heat capacity at pressure constant for dry air. The thermodynamic equations and an additional passive scalar 20 equation are: 25 where F Π corresponds to pressure effects. The transport of each prognostic scalar in Equations (1), (2) and (6) is made by a Piecewise-Parabolic Method (PPM) with undershoots and overshoots limitation (Colella and Woodward, 1984;Lin and Rood, 1996). The temporal algorithm of the advection term in these scalar transport is a forward-in-time scheme (noted FT). The momentum equations are: 30 where u is the resolved wind, g the acceleration due to the gravity appearing in the buoyancy term, θ r is the potential temperature of the reference state, µ f the dynamic viscosity, ∇.(ρ r u ′ ⊗ u ′ ) the Reynolds stresses. The spatial discretization of ∇.(ρ r u ⊗ u) in Equation (3) can be done by either second-or fourth-order centered schemes, third-or fifth-order Weighted-Essentially-Non-Oscillatory schemes (Jiang and Shu, 1996). The temporal evolution of the resolved wind is achieved by a fourth-order ERK Explicit Runge-Kutta algorithm (Shu and Osher, 1989;Lunet et al., 2017). In the present study ∆t is fixed 5 to respect the Courant number |u n |∆t ∆ < 1 (n, the time step index) and no additional time splitting is implied. The temporal viscous stability condition O(ν f /∆ 2 ) (ν f the kinematic viscosity) imposes an additional restriction when viscous term is explicitly resolved in time.
The bottom, lateral wall and top surfaces take a free-slip, impermeable and adiabatic behaviours without the call of an externalized surface scheme. The open boundary condition is a Sommerfeld equation defined as a wave-radiation (Carpenter, 10 1982) to enforce the large scales and allow the reflection wave damping.

The incompressibility condition
The wind of the resolved scales has to satisfy the continuity equation ∇.(ρ r u n+1 ) = 0. The method consists in the projection of the predicted velocity field u * (solution of Eq. (3) without the pressure term) into the null-divergence subspace. This projection estimates the irrotational correction to apply to u * through a potential scalar Ψ * : Ψ * is obtained with the resolution of the pseudo-Poisson equation written as: The horizontal part of the operator to invert in the elliptic problem is treated in the Fourier space (Schumann and Sweet, 1988) and its vertical part brings to the classical tridiagonal matrix. The mathematical operator to invert ∇.(∇) is exact for flat problems (Bernadet, 1995). When the mesh is built with a terrain-following coordinates over a flat surface, the solution of the pressure problem becomes inaccurate. In this orography-presence case, an iterative procedure is employed such as a Richardson, a Conjugate-Gradient (Young and Jea, 1980) or a Residual Conjugate-Gradient (Skamarock et al., 1997) algorithms. 25

The turbulent subgrid scales
To execute LES, the Reynolds stresses ∇.(ρ r u ′ ⊗ u ′ ) appearing in Equation (3) are estimated. The LES closure is done by an eddy-diffusivity approach called 1.5TKE with a 1.5 order closure scheme (Cuxart et al., 2000). The isotropic part of the subgrid turbulence is given by the prognosis of the subgrid turbulent kinetic energy e = 1 2 (u ′2 + v ′2 + w ′2 ): 30 where K e and K ϵ are constants prescribed in the turbulence scheme, l m and l ϵ the length scales defining the turbulent viscosity. The dissipation term is directly estimated from e and l ϵ (the left-hand term in Equation 6). The anisotropic part of the subgrid turbulence is diagnosed from the ψ gradient and e. The diagnosis of the anisotropic part of the subgrid turbulence is obtained using ψ and e: where Einstein summation convention is applied and Φ i,m are atmospheric stability functions (Redelsperger and Sommeria, 1981). The ground condition can be modelled by the externalized surface scheme SURFEX . In the dryversion of MNH and with the hypothesis of a zero thermal flux at the ground and buildings, only the turbulent friction is used.
To compute the non-zero values of u ′ i u ′ j at the ground, the SURFEX call employed in this paper consists in a simple activation 15 of a dynamic wall model related to the Prandtl theory (eddy viscosity concept). The form of the surface turbulent fluxes are Defining a friction velocity u * proportional to the turbulent wall shear and a roughness length z 0 , the vertical gradient of u is recovered by specifying a logarithmic profile (Kármán, 1930) as u(z) = u * k ln(1 + z z0 ) (note that the atmospheric stability conditions are neutral or near-neutral in this manuscript, therefore the additional Monin and Obukhov (1954) term is neglected and Businger et al. (1971) functions do not appear in the previous formulation). SURFEX is employed 20 in Sect. 5.2.

The IBM forcing in the Meso-NH code
The numerical domain is divided in two regions: the equations of continuum mechanics hold and a solid region embedding the obstacle where they do not. After comparing the methods ( Fig. 1-a) based on a local volume fraction function and the (LSF) LevelSet Function (Sussman et al., 1994;Kempe and Fröhlich, 2012), it was decided to use the LSF as it was able to capture 25 the interface between the fluid and solid regions more accurately : The | ϕ | distance informs about the minimal distance to the fluid-solid interface and the ϕ-sign about the region nature: sgn(ϕ) > 0 for the solid one; otherwise sgn(ϕ) < 0. The vector n normal to the interface and its local curvature σ are defined such as n = ∇ϕ |∇ϕ| and σ = −∇.n. Fig. 1-a illustrates the continuous variation of LSF for an arbitrary bell shape interface. The LSF is estimated at the seven available points type per cell to limit the discretization errors ( Fig. 1-b): at the mass point P where prognostic scalar variables are localized, at the three velocity nodes U/V /W where are characterized each projection u, at the A/B/C vorticity nodes employed by turbulent variables. The points of the solid region acts as external points of the computational grid (as acts external points in a boundary-fitted method at the grid limit). An intensive study had been done to estimate the well-modelling of the vector normal to the interface and the local curvature using LSF. The forcing based on a GCT Ghost-Cell Technique (resp. CCT Cut-Cell Technique) is applied to the explicit-in-time schemes (resp. the pressure solver) and detailed in Sect. 3.1 (resp. Sect. 3.2).

Ghost-Cell Technique and explicit-in-time schemes
The ψ n value is estimated at the time n∆t (∆t, the time step). The tendencies of the prognostic variables ψ = [u, θ, s, (e)] (Sect. 2.1) can not be deduced from the conservation laws in the solid region. Expecting a correction due to IBM where ϕ ≥ 0, a general formulation of the tendencies is written as:

10
The RHS first term of Equation (10) is given by the conservation laws (Sect. 2.1). The ∂ ∂t ibm tendency is the correction due to the GCT in the solid region and near the immersed interface satisfying the ψ desired boundary conditions at ϕ = 0: 15 Note that ∂u ∂t ibm is taken into account in the Explicit-Runge-Kutta (ERK) temporal algorithm. The freeze of the immersed wind conditions in the ERK algorithm had also been implemented; it had shown a more unstable behaviour for large Courant number.
The integer value k l determines the cells truncated by the interface: k l = 1 (resp. k l = 2) defines the first (resp. second) layer. The calculation of these ghost layers has a computational overhead due to data exchange among processors in parallel simulations. The stencil 10 of the numerical schemes modelling the interface defines the k l value. In order to limit this overhead, a low-order version of centered explicit-in-time scheme (Sect. 2.1) is employed when ϕ > −∆. The CPU cost of the 'hybrid' advection scheme is largely compensated by the decrease of the ghosts number and parallel exchanges. Appendix B reports a comparison analysis between third-order WENO and second-order centered schemes used in the vicinity of the interface; the studied case is the inviscid flow around a circular cylinder. 15 In the classical GCT (Tseng and Ferziger, 2003) the fluid information is obtained at a mirror point (noted I, merely renamed mirror) found in the normal direction to the interface in such a way that the interface node B is equidistant to I and G. The new GCT. To overcome this problem, we define image points (noted I 1 and I 2 in Fig. 2-a, merely renamed images) having a distance to the interface that depends only on the grid spacing: GI l = l∆ + ϕ G n with l = (1; 2). Figure 2-a shows the images for one ghost. The new approach enforces a large enough value of the |I l B| distance. Figure 2-b (resp. c) illustrates the 5 classical (resp. new) GCT for several ghosts. Figure 2-b shows some mirror points associated with ghosts of the first layer in the vicinity of the interface. Figure 2-c shows the new approach ensures the image points are always located in the fluid region, irrespective of the ghost location. The definition of several images per ghost allows to build a profile of the ψ fluid information normal to the interface. ψ(I) is therefore recovered by a quadratic reconstruction f using the (B, I 1 , I 2 ) points. Two distinct calculations of f (ψ(B), ψ(I 1 ), ψ(I 2 )) noted P LI a and P LI b are tested to build the Lagrange interpolation: where L a (I) and L b (I) are the Lagrange polynomials: The accuracy of an interpolation depends on the ψ-profile. For example, a logarithmic evolution of the tangent velocity is 20 expected in LES. Otherwise when the viscous layer is modelled, a linear evolution is expected. To compare the ability of each quadratic interpolation to approach a wide variety of profiles, the recover of power laws such as ψ = ϕ 3/2 (Figure 3-a) and ψ = ϕ 1/4 ( Fig. 3-b) is studied. As it illustrates, P LI a fits better the two analytical solutions and is therefore adopted. The classical and new GCTs have been compared thoroughly and part of this analysis is deferred to the Appendix. The interpolated field of the potential flow around a single cylinder or a sphere was compared to the theoretical solution; the sensitivity of the 25 inviscid flow around the same bodies (Appendix B) to the type of GCTs had been studied. The new GCT had given the best results, especially in the symmetry preservation in the inviscid flow cases. Note that these results are also dependent to the 3D interpolations choice detailed in a following paragraph. The proposed GCT is employed in the rest of this study. The GCT implementation is divided in four main steps: the fluid information recovery, the interface basis change, the interface condition and the ghost value.
For truncated cells (at least one corner node is in the solid region), ψ I l is recovered using an inverse Distance Weighting 5 (DW) interpolation: . This formulation diverges when x i → x l and it is commonly adopted to impose The use of these interpolations was decided after comparisons with Barycentric Lagrange and Modified Distance Weighting interpolations (Franke, 1982) and tests on the α coefficient. As the boundary condition is expressed in the interface frame and the grid is staggered, the non-collocation of the u components implies firstly to interpolate three different classes of cells (with U/V /W corners, Fig. 1 ghosts, secondly to build the change of frame matrix for which the proposed GCT presents an interest during the characteriza- 15 tion of the direction tangent to the interface.
The interface basis change. Velocity vector u known in the Cartesian mesh basis at the images I l (∆n 1 = ∆ and ∆n 2 = 2∆ in Fig. 2-a) is projected in the basis of the interface (n(B), t(B), c(B)) in which the boundary conditions on each vector component are imposed. Computing the LSF gradient, the normal direction is defined. Otherwise, (t, c) are two arbitrary tangent directions. The tangent direction t is considered as the predominant tangent direction of the flow along the fluid-solid interface depending on the images values and defining the velocity vector such as u(I l ) = u n (I l )n+u t (I l )t(I l ). The cotangent 5 direction is defined such as: The (n, t, c) basis at the interface is defined by considering or not the rotation of the tangent velocity with the distance to the interface: Finally, the third component is c(B) = n ⊗ e t (B) and (inverse) projection is known.
The interface condition. Let ψ B and ∆ ∂ψ ∂n B the Dirichlet and Neumann conditions on ψ. The general formulation of the 15 boundary condition ψ(ϕ = 0) is written as a Robin condition: An interface condition depending on the characteristics of the surrounding fluid such as ψ(ϕ = 0) = F (ψ I l ; ∂ψ ∂n I l ) is a wall model. Using two (resp. three) images, simple wall models such as the constant (resp. linear) extrapolation of the ψ gradient is 25 reached by the ∂ 2 ψ ∂n 2 I l = 0 (resp. ∂ 3 ψ ∂n 3 I l = 0) computation. The consistency between the tangent component to the interface of the resolved wind and the subgrid turbulence is the subject of Section 2.3.

Cut-Cell Technique and pressure solver
First looking at the RHS of (5), the ∂(ρru * ) ∂t csl coming from the resolution of the explicit-in-time schemes near the interface and in the solid regions badly affects the ∇.u * computation (note that the GCT operates after the step projection). Therefore 15 the fictive wind of the solid region can spread errors in the fluid region during the pressure resolution. To avoid it, a correction of the pressure solver is proposed.
The elliptic problem (5) is re-written as a resolution of the linear system P.Ψ * = Q. In the standard MNH version, ∇.u * = Q is estimated using a finite difference approach. To uncouple the solid region from the fluid region our revisited version enforces a null-divergence for pure solid cells and estimates the balance of momentum fluxes by a finite volume approach for truncated 20 cells (noted Q cct ): where V = ∆ 3 is the cell volume, V f (resp. V s ) the fluid (resp. solid) part of V, S i the cell surfaces where i is the index of each surface orientation [e,w,n,s,b,f] as it illustrates in Figure 5-a.
According to the Green-Ostrogradski theorem, the u i * S i calculation is the classical way of a CCT Cut-Cell Technique (Yang 5 et al., 1997;Kim et al., 2001) to estimate the velocity divergence. A similar approach is here performed re-building the flux ∆ 2 u * for truncated and solid cells.
The ± ∆ 2 u i * calculation consists in a weighting of the out-fluxes and in-fluxes function of the fluid and cell surfaces ratio ( Fig. 5-a). Figure 5-b gives an example of the west surface (i = w, red border) where ∆ 2 u w * (j, k) is calculated using the LSF value ϕ = ϕ w and the ones of the eight adjacent nodes ϕ p (j ± 1, k ± 1). A disk of radius −1 √ π∆ is split in eight 'piece of 10 cake' segments P * p (p = [1 : 8]). A LSF linear interpolation detects or not the interface location. In presence of an interface, its distance from the studied node is 0 < δ p < −1 √ π∆. Knowing δ p , the momentum fluxes balance is formulated for a non-moving body as (p is the index of the 'piece of cake' and i the index of the cell surface): The four encountered cases correspond to a pure fluid cell ∆ 2 u i * = ∆ 2 8 8 ∑ p=1 u i * when ϕ p < 0 and ϕ i < 0 ( Fig. 6-a); a pure solid cell ∆ 2 u i * = 0 when ϕ p > 0 and ϕ i > 0 ( Fig. 6-b); two types of truncated cells depending on the fluid/solid nature of the main node for which ϕ p .ϕ i < 0 (Fig. 6-c/d). Using Eq. (23), Equations (22) are solved and lead to the RHS computation of (5).
Knowing Q cct , the reflection concerns now the P matrix to invert. The classical interface condition on the potential Ψ * is a 5 homogeneous Neumann condition ∂Ψ * ∂ϕ = 0. Using a Boundary Fitted Method (BFM), the interface condition of the moving or non-moving body (Auguste, 2010) appears only on the border of a numerical domain. Using an IBM and without any impact of this interface condition on the P-coefficients, the impermeability character of solid obstacles is not achieved. Due to the inversion of the horizontal part of P by a Fast Fourier Transform (Schumann and Sweet, 1988), the solution of calculating P cct appears to be problematic. The adopted solution consists in an iterative procedure as used in MNH for non-flat problems. The fluid domain the enforcement of the null-divergence imposed on solid cells. The resolution of the pseudo-Poisson equation where M is the number of iterations. This number is limited by a convergence criterion (compromise between incompressibility satisfaction and CPU cost). Many iterative procedures are available in MNH and originally developed for non Cartesian grids. A Richardson and a Preconditioned Conjugate-Residual algorithms had been here adapted to the obstacles immersion. The newly modified pressure solver is tested and validated in Sect. 4.1.

Consistency with the turbulence scheme
It is known that the lm lϵ → 1 is a reasonable approximation in non-homogeneous, non-isotropic turbulence such as in the nearwall region. This approximation is indeed retained in the present IBM implementation, which assumes l m = l ϵ (hereafter noted l m and called the mixing length). The Redelsperger et al. (2001) corrections near the ground is to match the similarity laws and the free-stream models constants are not activated. l m is equal to the numerical cut-off space scale sufficiently far from  The turbulent characteristics are highly affected by a surface interaction. As a consequence and for LES, the subgrid turbulence scheme (Sect. 2.3) is modified in presence of immersed obstacles on the subgrid turbulent kinetic energy equation 15 (Eq. (6)), mixing length computation and Reynolds stresses diagnosis (Eq. (7), (8) and (9)).
The Subgrid Turbulent Kinetic Energy condition. The explicit-in-time resolution of Eq.(6) claims a GCT forcing and an interface condition on the STKE e. Commonly, the STKE profile is considered parabolically in the viscous sublayer (Craft et al., 2002;Bredberg, 2000) and constant in the inertial and wake/outer layers (Kalitzin et al., 2005;Capizzano, 2011). Due to the high turbulent Reynolds number Re t ≈ O(10 4 − 10 5 ) encountered, a homogeneous Neumann condition is applied at the immersed interface. The equilibrium between production and dissipation of STKE could be discussed and controverted; this choice acts as a first stage in the IBM development.

5
The near-wall correction of the mixing length. The von Kármán limitation due to immersed walls acts through the LSF and the upper limit on the mixing length l m near interface becomes min(k z, −ϕ, ∆) with a banning of negative values in the solid region. Whatever the production of Subgrid Turbulent Kinetic Energy e (STKE) and the turbulent shear, the lower limit l m (−ϕ ≤ 0) induces a null value of the diagnosed surface fluxes. In addition, a singularity appears in the dissipative term ρ r K ϵ e 3/2 l −1 m . By a pragmatic reasoning, the singularity due to l −1 m (ϕ → 0 − ) → ∞ amounts to say that modelled length . This thickness estimate is also proportional to Eν/u * (E ≈ 9.8 is commonly employed) where the friction velocity u * is about the 15 centimeter per second. Following these estimates, the length scale due to the viscous effects z ib ν belongs to the millimetres domain in the expected atmospheric cases. Looking after a building surface and its large heterogeneity (door, windows, surface characteristics), its roughness length z ib 0 is at least in the decimeter domain and z ib 0 > z ib ν (Illustration in Fig. 7). For low Re and smooth surfaces, z ib ν > z ib 0 could be encountered. Therefore, we assume z ib = max(z ib 0 , z ib ν ) and that z ib is related to the size of smallest unresolved eddies near walls (i.e. dissipative scale). The mixing length near wall is z ib < l m < min(k z, −ϕ, ∆). 20 The turbulent fluxes correction. The ψ-gradient and the turbulent diffusion O(z ib √ e) prescribe the turbulent fluxes at the immersed interfaces (Eq. (7), (8) and (9)). As a first step in the MNH-IBM implementation, no-flux condition on the mean potential temperature is imposed bringing to a zero-value of the sensible heat flux. Writing the mean velocity field at B such as u = u t t, u t (B) is needed to recover a gradient consistent with the turbulent shear. Considering the Prandtl (1925) or Kármán 25 (1930) theories, the logarithmic profile is assumed in the vicinity of the wall according to u t (z) = u * k ln(1 + z z ib ). Considering ∆ as the limit of the resolved scales, most of the turbulent kinetic energy 1 2 (u ′2 + v ′2 + w ′2 ) is contained in the subgrid one when −ϕ < ∆ and such as K tke √ e with a constant K tke ≳ 1. This assumption is reinforced by the homogeneous Neumann condition applied on e. This approach derives from the RANS (Reynolds-Averaged Navier Stokes) approaches and the velocity friction is formulated as u * = K tke 4 √ C µ √ e where C µ is a constant evolving between 0.03 (atmospheric applications) and 0.09 30 (fluid mechanics applications). Adding a damping function for the viscous cases (low turbulent Reynolds number, Re t < 20), the tangent wall velocity at the interface is written such as: Finally the pragmatic limitation u t (B) ≤ u t (ϕ = −∆/2) operates if the STKE value is too high. The proposed dynamic wall-model evolves between the no-slip and free-slip conditions. If the subgrid turbulence is weak or if the physical problem is fully resolved, the viscous layer is well-modelled and u t (B) → 0. Otherwise for an intense subgrid turbulence or a fully unresolved problem, the shear due to the wall presence is not perceived and dut dn B → 0. The wall-model establishes an equilibrium between the production of STKE and the mean parietal friction. Note that the use of a log-law model near a singularity such 5 as sharp edges and corners could be called into question. Nevertheless, Sections 5.1 and 5.2 show LES results employing this proposition. After numerical investigations done during the single cube study, K tke 4 √ C µ ≈ 1 appears as a suitable choice.
4 Flows around a circular cylinder

Potential flow
Isolated from the rest of the code, the resolution of the pseudo-Poisson equation (5) leads to potential solutions (Sect. 2.2).

10
Theoretical ones are available for flow developed around a non-deformable obstacle such as an infinite cylinder or a sphere (Milne-Thomson, 1968;Batchelor, 2000). The two bodies are here investigated. The flow around the infinite cylinder is predominantly presented. Figure 8 illustrates the cylinder case. The fluid density is considered as constant in time and in space. The flow is initially 15 imposed as spatially homogeneous with a constant module of velocity U ∞ and parallel streamlines ( Fig. 8-a). This initialization does not respect the conservation of the momentum flux and the irrotationnal correction of the projection method goes to recover this conservation. In the same time, the impermeability of the cylinder of diameter D cyl = 2R cyl is achieved. Figure 8b shows the streamlines obtained with the MNH pressure solver modified to take into account the presence of immersed obstacles (Sect. 3.2). Defining x as the direction parallel to the initial streamlines and y as the perpendicular one, the expected 20 solution is u.
r 2 )y (single and non-confined body, (α; r) cylindrical coordinates). The numerical confinement is comment hereafter, characterized by L = L cyl /R cyl where L cyl is the distance separating the lateral domain surfaces (Fig. 8-a).
The Richardson (RICH) and the Residual Conjugate Gradient iterative (RESI) methods are tested (Sect. 3.2). Figure 9-a plots the evolution of the dimensionless residue R(k) (based on a characteristic divergence U ∞ /∆) with the iterations number k 25 and obtained with the confinement L = 16. The two algorithms converge with a weak dependence to the spatial discretizations (N = [4(red); 8(green); 16(blue); 32(purple)] nodes per R cyl ). The slope coefficient dR(k) dk (RESI) ≲ 3 dR(k) dk (RICH), so RESI demonstrates its highest velocity convergence. Even if RICH is about 20% faster per iteration than RESI, the global CPU cost of the last one is lowest for a same solver residue R(k). For this reason and due to an a priori higher radius convergence, RESI is adopted. Note that the momentum flux computed after the solver convergence at the x cyl location (Fig. 8-a)   shown here). . The confinement is defined in Fig. 9-a.
With a change of Galilean reference frame, this study corresponds to an uniform body acceleration a b in a fluid initially at rest. However a possible viscous term, the hydrodynamic force exerted on the body, is reduced to the added mass effect A is the dimensionless coefficient and m f the displaced fluid mass. A cyl theoretically equals 1 in the non-confined cylinder case (Lamb, 1932). The red curve of Figure 9-b illustrates the effect of the confinement L on A cyl for a N = 16 resolution. Unsurprisingly A cyl increases with the confinement (Brennen, 1982). The weak dependence of A cyl with L > 16 allows to consider the body as isolated for L ∼ 16. The green curve of Figure 9-b shows the impact of the space resolution for L = 16. The numerical added mass coefficient is in good agreement with the theoretical one presenting a 5 relative error of about 2% for N > 16. It induces a well-respect of the impermeability hypothesis at the immersed interface.
A similar study for a spherical body gives A sph = 1 2 + 0.4%. Figure 10 illustrates the contours of the kinetic energy around the sphere in an arbitrary symmetry plane. The green contours (numerical solutions) fit well with the red contours (theoretical solutions). Convergence study of the pressure solver is discussed in Appendix A.

10
A pure dynamic and well-documented case which naturally follows the previous ones is studied here. This physical case is the wake past a circular cylinder (non-stratified flow) at two moderate Reynolds numbers Re = (40; 140). One of the forerunners is Taneda (1956) who experimentally studied the nature of the eddies structure.

Taneda (1956) found a regular Hopf bifurcation at a critical Reynolds number
Below Re c and above Re > 5, a boundary layer separation brings to a steady recirculating region in the near-wake (Fig. 11-c). Above Re c ≈ 45, an 15 unsteady mode breaks the planar symmetry and the body wake presents an alternate vortex shedding (Fig. 11-d). The standing eddy (resp. the von Kármán street) obtained by MNH-IBM at Re = 40 (resp. 140) is visualized in Figure 11-a (resp. b) by the injection of a passive tracer on the body surface.
The standing eddies at Re < Re c are commonly described with a θ d detachment angle, l r recirculating length and (a; b) location of the vortex core (Fig. 12). The limit of the numerical domain is 10D cyl upstream the obstacle for the inlet condi-20 tion (U ∞ , the uniform incoming velocity) and lateral condition (slip condition), 15D cyl for the outlet condition allowing the  Note that the impact of the low-order (centered or WENO) modeling the advection at the immersed interface is weak for this

Flow over a surface mounted-cube
Using static pressure measurements, laser-sheet and oil-film visualizations, Martinuzzi and Tropea (1993) and Hussein and Martinuzzi (1996) generated a large data set for the study of flows around a cubic body placed in a channel ( Fig. 13-a). RANS and LES had explored in detail this physical case (Breuer et al., 1996;Shah and Ferziger, 1997;Rodi et al., 1997;Frank, 1999;Krajnovic and Davidson, 2002;Farhadi and Rahnama, 2006).

5
Physical details. A cube (H side) is placed in a channel of 2H height. The channel is sufficiently large in the span-wise direction to consider the cube as single in that direction. Turbulent flow is generated in the channel upstream the cube with a mean bulk velocity U b . Defining the dimensionless wall coordinate z + = u * .z/ν f , the stream-wise upstream velocity corresponds to a log-law for smooth walls u(z + ).u * −1 = 5.54 + 1 κ log(z + ) as described in Hussein and Martinuzzi (1996). The Reynolds number defined by the mean bulk velocity, the cube height and the molecular diffusion is Re ≈ 40000. 10 The mean flow around the cube presents a set of five recirculating regions (Fig. 13-b). Each cube surface is associated with one of these regions: the A/B vortex separations in front of the cube which spread laterally in a horseshoe D, two vortices near side-walls E, one F on the roof and main arch vortex G downstream.  (Hussein and Martinuzzi, 1996); (b) Schematic representation of the time-average vortex structure around a cube    u ′ 2 /U 2 b ≈ 2.10 −2 (Hussein and Martinuzzi, 1996) are recovered at x/H ∼ −4. We mention that the turbulence generation should deserve more details but we prefer only to concentrate our comments in the cube wake.  ; (e) MNH-IBM instantaneous visualization of the Q-criterion (Hunt et al., 1988).
Results. Figures 14-(a-b-c) show the time-average streamlines in the vertical symmetry plane of the cube obtained by MNH-IBM for the three space resolutions. The streamlines of the coarse, medium and fine resolution are respectively in red, green and blue colour. The discretization order of the fine resolution is close to that of most literature's LES except Shah and Ferziger (1997) who had used a far more precise grid near wall. The same figure obtained by the experimental investigation is given in Figure 14-d. The size of the front (resp. rear) region is characterized by the recirculating length x f /H (resp . x r /H). The  colour code corresponds to the spatial resolution (10pts/H in red; 20pts/H in green; 40pts/H in blue). Note that the U b bulk velocity was set to the unity and therefore presented variables in Figure 15 are dimensionless. In most of the sub-figures, the higher space resolution is the lower the gap with the experiment is. The top of Figure 15 brings to a similar conclusion with that of the literature's LES: the stream-wise velocity is well-recovered. The counterflow at the rear and roof of the cube near (x/H; z/H) ≈ (1; 1) informs about the existence of the bifurcation point S1 (Fig. 15-b). An underestimation appears on the 10 counterflow at (x/H; z/H) ≈ (2; 0) and is frequently observed in the literature (Fig. 15-c). Some discrepancies are found on two w-profiles ( Fig. 15-f/g). Note that Shah and Ferziger (1997) and Krajnovic and Davidson (2002) highlight the difficulty to recover it. The turbulent kinetic energy and the Reynolds stress are correctly predicted at the cube roof (

Sensitivity study. Some tests have shown a sensitivity of the solution with both the incoming turbulence and constants in
the immersed wall-model. Indeed the existence of the small vortex near the saddle node S1 (Depardon et al., 2006) turns out to be strongly dependent of the inlet turbulence (the less the turbulent intensity is the bigger the size of this vortex is); the 20 reattachment or not of the roof vortex F with the main arch G is also observed. This sensitivity follows the observations of Castro and Robins (1977) who studied the cube placed in a uniform or turbulent incoming flow. In the same way the existence of the vortex near the saddle node S2 (Frank, 1999) is dependent of the ground boundary condition. The u * −1 √ e = K tke 1/4 √ C µ ratio and the z 0 roughness length fixed in the immersed wall-model (Sect. 3.3) impact the nature of the incoming turbulent state for which the surface shear plays an important role. To give a significant example and if a null-value of K tke is applied 25 on the channel surfaces (non-slip condition), x f highly increases and the vortex at S2 appears. Otherwise K tke and z 0 do not crucially affect the nature of the dynamic in the vicinity of the cube surfaces; the pressure gradient between front and back faces governs.
To conclude this section and despite the uncertainties of the inlet condition, MNH-IBM is in a good agreement with the experiments of Martinuzzi and Tropea (1993), Hussein and Martinuzzi (1996) and other LES using a resolution of about forty 30 points per cube length. Even if the coarsest resolution looses a part of the expected physics, it maintains a suitable modelling of the largest structures of the flow. A parametric study on the inlet turbulence generation has led to fix K tke ≈ 2 in Equation (24).  Martinuzzi and Tropea (1993) data except for TKE (Hussein and Martinuzzi, 1996). The profiles are given at four longitudinal locations: (a-e-i-m)

The Mock Urban Setting Test experiment (MUST)
The MUST is an experimental campaign organized during the early Autumn 2001 in the Utah's West desert (Biltoft, 2001(Biltoft, , 2002. Its objective was to quantify the dispersion of a passive tracer (propylene) in a dry atmospheric context over a topography reproducing a near-urban canopy. The main interest lies in the similitude between this experiment and a pollution episode due to a toxic gas propagation over a city with high population densities. It provides extensive measurements of meteorological 5 variables and scalar dispersion informations.
Physical details. The near-regular array is composed of 120 containers. Figure 16 gives a photograph and a picture of the topography. The containers are equivalent in volume and shape. Their spatial dimensions are (L x , L y , L z ) =∼ (2.4, 12.2, 2.5)m and the horizontal distance between containers is O(10)m. Following the table II in Yee and Biltoft (2004), the 2681829 case   Biltoft (2001) and Yee and Biltoft (2004).
Numerical details. The externalized scheme SURFEX  models the ground friction. LSF is generated to represent the topography and IBM is used to model the containers. The smallest characteristic container length is discretized by O(10) cells. The mesh is a Cartesian and regular grid for z < 10m with a space step ∆ = 0.2m. Above ten meters, the vertical space step is released in geometric progression of 1.08 ratio. The altitude of the numerical domain top is about 8 times the container height.

5
The distance between the horizontal limit of the computational domain and the array is about 20 times the container height.
The large scale flow is forced by an open boundary condition. A mean horizontal angle of -41 • with the x direction is fixed for all altitudes; note that the low angle deviation and the turbulence observed upstream the containers in the experiment are not numerically considered. Following the experimental data given at the S tower ( Fig. 17-a, blue symbols) and assuming a log-law |u| = u * κ log(1. + z/z 0 ), a least-square regression estimates a friction velocity u * = 0.71m.s −1 and allows to build the vertical 10 profile of the mean incoming flow ( Fig. 17-a, blue line). This u * value is close to the experimental one u * exp = 0.68m.s −1 found by a sonic anemometer at the feet of the S tower. The same inlet condition is used in the numerical studies of Hanna et al. (2004), Milliez and Carissimo (2007) and Donnelly et al. (2009). Note that an additional term O(L −1 mo ) appears in their formulation but is negligible in the 2681829 case. Results. The black line (MNH-IBM) and symbols (Yee and Biltoft, 2004) of Figure 17-a (resp -b) show the impact of the 15 near-urban canopy on the vertical profile of the |u| (resp. wind angle) in the core of the array (T tower). Not surprisingly the canopy induces a global slowdown of |u| near the ground and up to z ≲ 8m. A decrease of the mean horizontal wind angle is found at the same altitudes. This deviation is related to the containers orientation which tends the flow to be aligned with the y direction. The same pattern is discussed in Milliez and Carissimo (2007). A wind acceleration and an increase of the wind angle are observed by MNH-IBM for z ≳ 8m as in LES results of Konig (2013) and Dejoan et al. (2010) on similar MUST cases. This few degrees deviation is observed by the experiment but not the acceleration. A part of this acceleration may be explained by a Venturi effect and a too closed top limit of the computational domain (numerical confinement, under investigation).  Figure 18-a (resp. Fig. 18-b) illustrates the time-averaged horizontal wind field |u| (resp. |u| at t = 200s) at the 1.6m 5 altitude. These figures highlight the fact that the incoming flow is not turbulent in the simulation. This assumed gap constitutes a perspective (Camelli et al., 2005) not directly linked to IBM. It allows also to introduce here some comments on the turbulence state. The atmospheric turbulence is dependent of the roughness length (z 0 considered as constant due to the homogeneous and flat desert over few miles upstream). The first container rows are the scene of the boundary layer transition and act as a region of strong roughness change. The turbulence observed in the urban-like canopy has two origins: the incoming turbulence 10 and that of induced by the containers presence. The contribution of both turbulence types varies in function of the altitude and distance to the first row.
The mean kinetic energy E k = 1 2 (u 2 +v 2 +w 2 ), friction velocity u * = 4 √ u ′ w ′ 2 + v ′ w ′ 2 and turbulent kinetic energy T KE = 1 2 (u ′2 + v ′2 + w ′2 ) are estimated at the T tower and indicated in Table 2. All the variables are in good agreement with the experiment at z = 4m (about two times the container height). The friction velocity at z(T ) = 4m is at least two times larger 15 than u * exp (S) = 0.68m s −1 observed at the S tower feet. That increase is the signature of the turbulence developed by the urban-like canopy. Looking after the results at z(T ) = 8m and z(T ) = 16m, the higher the altitude is the more discrepancies appear between the experimental and numerical results. The experimental measurement of the friction velocity at z = 16m is closed to the upstream value. (2004) (Yee and Biltoft, 2004) and numerical (MNH-IBM) investigations at three altitudes of the tower T.

MNH-IBM Yee and Biltoft
The Discrete Fourier Transform of the u temporal evolution (Fig. 19)   The MNH-IBM results are consistent with the experimental observations below z(T ) ≲ 10m. This well-modelling induces that the turbulence is mostly due to the container wakes upstream the T tower and not to the incoming turbulence (not modelled 5 in the simulation). The influence of the atmospheric turbulence grows with the altitude and MNH-IBM diverges with the experiment. This divergence for z(T ) ≳ 10m leads to think that the thickness of containers influence is about 4-5 times the height of the urban-like canopy.

Conclusions and perspectives
This study details the first implementation of an immersed boundary method (IBM) in the atmospheric Meso-NH (MNH) model currently based on mathematical formulations written for structured grids. The MNH-IBM aim is to explicitly model the fluid-solid interaction in the surface boundary layer developed over grounds presenting complex topographies such as cities or industrial sites.

5
A LevelSet function (Sussman et al., 1994) characterizes the geometric properties of the fluid-solid interface. Two original approaches of the (GCT) Ghost-Cell (Tseng and Ferziger, 2003) and (CCT) Cut-Cell Techniques (Udaykumar and Shyy, 1995) are implemented to correct the MNH numerical schemes. A newly proposed GCT recovers the fluid information in several image points presenting a distance to the interface independent of the ghost's distance. The CCT consists of a new finite volume approach of fluxes balance near the immersed interface. The GCT is applied to the numerical schemes based on explicit The pressure solver, adapted to the IBM and isolated from the rest of MNH, is used to model potential flows around several 15 obstacles. Compared to analytical and theoretical solutions, the numerical results demonstrate the ability of the IBM adaptation to ensure that the momentum is preserved and the continuity equation is respected. Non-dissipative flows are simulated to test the IBM forcing of the wind advection scheme (the impact of interpolations collected in Franke (1982), classical and novel GCT comparison and numerical diffusion near interface). These tests validate the proposed GCT 'three images/ghost points' using an inverse distance weighting (resp. trilinear) interpolation near (resp. far from) the interface. The modelling of the 20 advection term near the fluid-solid interface is ensured by a 2nd-order centered scheme associated with an artificial viscosity.
The artificial viscosity is calibrated after comparisons with a 3rd-order WENO scheme. With these numerical choices, MNH-IBM demonstrates its ability to model wake instabilities past a circular cylinder placed in a viscous fluid. Then Large-Eddy-Simulations of turbulent flows around bodies with sharp edges and corners are executed (i.e. a cube placed in a channel ) and a near-neutral atmospheric application over an array of containers (Biltoft, 2001)). These 25 two LES validate the proposed immersed wall-model, switching the characteristic space scale defining the turbulent Reynolds number between one obtained by either a viscous length scale (the cube case) or a roughness length (the containers case).
Future work. This study constitutes a first robust step towards a better understanding of the interactions between 'Weather and Cities' and better predictions of such interactions. The idealized character of the physical cases approached here offers some insights. One improvement would consist of a generalization of the IBM writing to terrain-following coordinates (Gal- 30 Chen and Somerville, 1975) allowing for the simulation of high-curvature bodies in the presence of non-flat ground. In the run-up of the resolution of multi-scales problems, the consistency with grid nesting (Stein et al., 2000) would be pertinent as well as the coupling to a drag model (Aumond et al., 2013). In the current paper, 'simple' bodies are investigated; the modelling of real houses and buildings with arbitrary shape in close proximity to each other is ongoing with work dedicated to a brief and intense pollutant episode due to a factory explosion in 2001 over Toulouse assuming a dry neutral case with a non-reactive gas dispersion. Such hypotheses involve a broad range of physical phenomena requiring numerical compliance with IBM including the chemical reactions, phase changes and radiative effects. Such compliance would allow the access to a large variety of atmospheric situations.
the numerical pressure and P t the theoretical one) are estimated in presence or not of an immersed cylindrical body (Fig. 20).
The space second-order of the pressure solver is recovered without IBM. The order decreases with IBM and stays consistent 5 regarding the L p=(∞;1;2) slopes. Note that an immersed square or sphere give similar results.
Agnesi hill. The irrotationnal solution around two 2D bell shape interfaces is investigated with IBM and the Boundary-Fitted Method (BFM, terrain-following coordinates). The topography is characterized by a height h a and a shape h(x) = ha 1+( ka .x ha ) 2 (k a = 4, bell 1; k a = 8, bell 2). The bells slope is here arbitrarily and respectively described as gentle or steep. Figure 21 10 shows the pressure contours obtained with IBM (left) and BFM (right) for a gentle (top) and a steep (bottom) shape. The minimal pressure value is localized at the top of each bell and goes to zero far from this location. The reference BFM and IBM simulations with the fine resolution (N = 160 nodes per h a , red colour) show a good agreement for each hill. The blue (resp. green) colour corresponds to a coarser mesh employing N/3 (resp. N/9) nodes per h a . Weak differences appear between the N and N/3 meshes for both IBM and BFM revealing a good space convergence ( Fig. 21-a/b). Numerical errors are visible 15 with IBM near the interface but the Venturi effect is well-modelled. Differences become more significant with the N/9 mesh especially with the BFM-BELL2 presenting the highest curvature value (Fig. 21-d). IBM appears less accurate than BFM when the ground presents low curvature in regard of the space resolution. Otherwise, IBM seems more pertinent than BFM to model high interfaces such as sharp edges or corners. The minimum pressure which could reach to −∞ is smoothed by IBM and allows the pressure solver not to diverge. Note that Lundquist et al. (2010Lundquist et al. ( , 2012 using compressible WRF model observe 20 similar behaviors.   Table 3. Summary of the used mean wind advection scheme (WENO = Weight-Essential-Non-Oscillatory).

Appendix B: Inviscid flow around a circular cylinder
For most atmospheric applications, the region size where the fluid molecular viscosity ν f influences the dynamic is sufficiently small to be considered as negligible (Sect. 2.3). Solving the Euler equations, the impact of the numerical diffusion could be significant especially near the fluid-solid interface. The adopted strategy with IBM is to model the advection term with a 5 low-order scheme near the interface (Sect. 3). The order decrease in IBM is not essential but allows to limit the number of ghosts in the solid region, limit the communications during a parallel computation. Indeed, the chosen implementation implies that the associated images and ghosts points have to be localized in the integration volume of each processor. The WEN3 third-order Weight-Essential-Non-Oscillatory and CEN2 second-order centered schemes are available in MNH. Far from the interface a CEN4 fourth-order centered scheme is employed. Table 3 summarizes the advection scheme nomenclature.
The vorticity equation for a 2D inviscid flow reveals no production in time. The numerical vorticity production at the 5 immersed surface of a cylindrical body is here studied initializing the simulation with the potential solution. To fit as well the potential solution, a non-trivial condition is employed on the tangent velocity ∂ 3 ut ∂n 3 = 0. Expecting a numerical vorticity sufficiently controlled to avoid the flow separation, the effect of the artificial diffusion ν art ∆u injected with CEN2 is compared to the WEN3 intrinsically diffusive behaviour. Furthermore this study had estimated the 3D interpolations impact (not detailed here).     Figure 23b corroborates the last comment presenting the vorticity contours dimensionless by D U∞ . A suitable ν art combined with CEN2 choice is also in the range of the too diffusive WEN3 results and the growth of numerical instabilities. CEN2+ν ref art 4 −1 is retained as the advection scheme of the mean wind near an immersed interface.