Large eddy simulation model for wind-driven sea circulation in coastal areas

In the present paper a state-of-the-art large eddy simulation model (LES-COAST), suited for the analysis of water circulation and mixing in closed or semi-closed areas, is presented and applied to the study of the hydrodynamic characteristics of the Muggia bay, the industrial harbor of the city of Trieste, Italy. The model solves the non-hydrostatic, unsteady Navier–Stokes equations, under the Boussinesq approximation for temperature and salinity buoyancy effects, using a novel, two-eddy viscosity Smagorinsky model for the closure of the subgrid-scale momentum fluxes. The model employs: a simple and effective technique to take into account wind-stress inhomogeneity related to the blocking effect of emerged structures, which, in turn, can drive localscale, short-term pollutant dispersion; a new nesting procedure to reconstruct instantaneous, turbulent velocity components, temperature and salinity at the open boundaries of the domain using data coming from large-scale circulation models (LCM). Validation tests have shown that the model reproduces field measurement satisfactorily. The analysis of water circulation and mixing in the Muggia bay has been carried out under three typical breeze conditions. Water circulation has been shown to behave as in typical semi-closed basins, with an upper layer moving along the wind direction (apart from the anti-cyclonic veering associated with the Coriolis force) and a bottom layer, thicker and slower than the upper one, moving along the opposite direction. The study has shown that water vertical mixing in the bay is inhibited by a large level of stable stratification, mainly associated with vertical variation in salinity and, to a minor extent, with temperature variation along the water column. More intense mixing, quantified by sub-critical values of the gradient Richardson number, is present in near-coastal regions where upwelling/downwelling phenomena occur. The analysis of instantaneous fields has detected the presence of large crosssectional eddies spanning the whole water column and contributing to vertical mixing, associated with the presence of sub-surface horizontal turbulent structures. Analysis of water renewal within the bay shows that, under the typical breeze regimes considered in the study, the residence time of water in the bay is of the order of a few days. Finally, vertical eddy viscosity has been calculated and shown to vary by a couple of orders of magnitude along the water column, with larger values near the bottom surface where density stratification is smaller.


Introduction
Coastal basins are in general shallow and characterized by complex geometry arising from rapid varying bathymetry, coastline and anthropic structures.Such features may produce three-dimensionality, breaking waves and along-shore currents.In semi-closed basins, wind shear stress is the main forcing term directly driving water surface layers and triggering mechanical turbulence production.In such conditions, interaction with the coastline develops downwelling/upwelling along vertical planes and the inversion of the mean velocity field in the bottom layer of the water column with respect to the superficial one.The resulting additional shear enhances further mechanical turbulence production.Further, in coastal regions, variations in temperature and salinity along the water column give rise to buoyancy-driven mixing.
Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.
The aforementioned features make modeling of coastal hydrodynamics quite challenging for general circulation ocean models, which, on the other hand, have been conceived and calibrated for analysis of meso, or even larger, scales.Further, coastal hydrostatic or three-dimensional models developed in the hydraulic engineering community, which make use of basic turbulence parametrization, are often not suited for detailed and accurate analysis of internal mixing in coastal areas.For example, it has been recognized that the hydrostatic approximation, adopted by most models, cannot be applied successfully for the coastal environment where upwelling/downwelling phenomena play an important role; moreover, most coastal models adopt the Reynolds-averaged (RANS) approach, which is not designed to reproduce the interplay between physical mechanisms occurring in coastal regions (see Burchard et al., 2008).
On the other hand, the increased computational capability available over the years has made possible new modeling approaches to the problem.In particular, large eddy simulation (LES) has recently emerged as a powerful and promising methodology to afford real-scale coastal problems characterized by the occurrence of complex physical processes in complex geometries (Burchard et al., 2008).In LES the large scales of motion, typically threedimensional and anisotropic, are directly solved through a 3-D unsteady simulation, whereas the spatial scales of motion whose dimensions are smaller than the grid spacing (sub-grid scales, SGS), more dissipative and isotropic, are parametrized through the use of a proper model.
Large eddy simulations for ocean applications were first carried out by Skyllingstad and Denbo (1995).The authors considered simple geometry (a Cartesian box) to study the dynamics of plumes under convective conditions.Other studies have made use of LES for the analysis of mixing in the open ocean (see, among others, Wang et al., 1996).A review of the use of LES for marine application is in Scotti (2010).
As briefly described above, in coastal regions complications arise, and properly designed numerical models must be used.A three-dimensional coastal basin exhibits horizontal and vertical length scales very different from each other.The horizontal one is of the order of a few kilometers, whereas the vertical one usually lies within the range of 10 to 100 m.Turbulence itself displays large horizontal structures that are ruled mostly by the Earth's rotation and by the free surface and solid boundaries, giving rise to strongly anisotropic flow features.Unfortunately, in LES, eddy-viscosity SGS models are mostly designed for isotropic or nearly isotropic grids, making prohibitive mesh requirements for real-scale coastal studies (for a discussion see Scotti et al., 1993).In highresolution LES of real-scale coastal hydrodynamics, grid resolution in the horizontal direction can be pushed up to 10to-1 m, whereas the vertical resolution may be of the order 0.5 m, introducing large anisotropy to the grid topology itself and, hence, to the implicit filter widths.When isotropic SGS models are used in the presence of large cell anisotropy, they overestimate the contribution of the unresolved scales of motion, producing over-dissipation, which strongly affects the accuracy of the simulation.Recently, Roman et al. (2009c) proposed an SGS model properly suited for coastal applications, thus able to take into account the physics and the geometric complexities usually encountered in real-case applications.It implements a modified anisotropic Smagorinsky model (ASM) that takes advantage of the two-eddy-viscosity idea borrowed from geophysical fluid dynamics concepts.The model, applied to the analysis of water mixing and renewal in Barcelona harbor (see Galea et al., 2012), has been shown to be able to reproduce correctly the velocity of the wind-induced current as well as the complex flow patterns developing within the harbor.The results corroborate the findings of Ramachandran et al. (2013) about the ability of LES used in conjunction with the ASM model to reproduce accurately flow features in coastal hydrodynamics.
However, additional issues are still to be addressed for the application of LES methodology to real-scale coastal studies.The assessment of proper conditions at the boundaries of the computational domain has a relevant impact.From one side, the suitable assignment of momentum and heat fluxes at the free surface is of crucial importance for the accurate reproduction of the wind-driven circulation at a coastal scale.Further, since the near-shore circulation is always affected by mesoscale circulations, proper procedures are required to nest the high-resolution LES simulation within large-scale circulation models.
To summarize, analysis of water dynamics in areas characterized by interplay of different physical processes and geometrical complexities, like harbors and lakes, may take advantage of novel state-of-the-art, properly designed, numerical tools able to reproduce such features.
The scope of the present paper is to show a novel, stateof-the-art, numerical model suited for analysis of mixing and water renewal in coastal semi-closed regions (i.e., harbors) and lakes.The model solves the Boussinesq form of 3-D nonhydrostatic, primitive-variable, filtered Navier-Stokes equations, considering temperature and salinity effects on buoyancy in the water column as well as dispersion of passive scalars (i.e., pollutants).The unresolved subgrid scales deriving from the filtering operations are parametrized through a novel two-eddy viscosity Smagorinsky model.The model employs: (a) a simple and effective procedure to take into account wind-stress inhomogeneity at a local scale, related to topographic effects in the low atmosphere or, among other things, to the presence of ships or jetties; (b) a novel nesting methodology combining interpolation procedure of the flow variable from the LCM domain to the LES one, together with a synthetic generation of turbulence within the LES domain, to reproduce realistic dynamics within the area of interest.Finally, the model is applied to the analysis of water circulation and mixing in the Muggia bay, a site of great interest in Italy from the water quality point of view, due to the presence of potentially dangerous industrial plants.
The paper is organized as follows: in Sect. 2 the numerical model (LES-COAST) is presented along with a detailed discussion of the boundary conditions at the free-surface and at the open sections.In Sect. 3 the model is firstly validated against field data, and successively water circulation and mixing in the Muggia bay are discussed.Finally, concluding remarks are given in Sect. 4.

Governing equations
LES-COAST solves the Boussinesq form of the filtered, 3-D, non-hydrostatic Navier-Stokes equations together with the transport equations for temperature and salinity.Transport of pollutants treated through an advection-diffusion equation for their concentrations (or as dispersed Lagrangian particles) is also implemented in the model.The governing equations read as: ∂T ∂t + ∂u j T ∂x j = k T ∂ 2 T ∂x j ∂x j − ∂λ T j ∂x j (3) ∂S ∂t + ∂u j S ∂x j = k S ∂ 2 S ∂x j ∂x j − ∂λ S j ∂x j (4) where u i , p, T , and S are, respectively, the velocity component in the x i direction, the hydrodynamic kinematic pressure, temperature (in Kelvin) and salinity (in PSU); g is the gravitational acceleration, i the ith component of the Earth rotation vector, and ρ and ρ 0 are density anomaly and its reference value, respectively.ν, k T and k S are kinematic viscosity, and thermal and salinity diffusivity, respectively.Since the density anomaly due to salinity and temperature variations is much smaller than the bulk density of the water, in coastal applications the buoyancy force can be accounted for by the Boussinesq approximation through the equation of state: where ρ tot is the total density, T 0 and S 0 are reference values giving the density ρ 0 , and β T and β S are expansion coefficients for temperature and salinity, respectively.C n in Eq. ( 5) refers to the concentration of the nth pollutant and k Cn is its kinematic diffusivity.Over-bar refers to the filtering operation.LES-COAST uses an implicit filter, meaning that filtering is implicitly performed in the physical space through a top-hat filter function represented by the cell size.The latter term on the right-hand side of Eqs. ( 2)-( 5) represents the subgrid contributions to momentum and scalar fluxes, respectively.The sub-grid scale momentum fluxes are parametrized by the two-eddy-viscosity anisotropic Smagorinsky model (ASM) described in Roman et al. (2010) and briefly summarized hereafter.
In the standard Smagorinsky model, SGS stresses are expressed as: where S ij is the resolved strain rate tensor and ν t is SGS eddy viscosity.The Smagorinsky model is isotropic and eddy viscosity is evaluated as the product of a length scale C , proportional to the grid size, and a velocity scale C |S|, where C is a constant and |S| is the contraction of the resolved strain rate tensor.In moderate-to-strong anisotropic grids it is difficult to define a single length scale.In particular, setting = ( x y z ) 1/3 leads to an excessive overestimation of the eddy viscosity in all directions.Thus, in Roman et al. (2010) a two eddy-viscosity model was developed, one for the horizontal direction and another for the vertical one, respectively.The two eddy viscosities are defined as: where L h = ( 2 x + 2 z ) 1/2 and L v = y are proper length scales for horizontal and vertical directions, respectively.The velocity scales in the horizontal and vertical directions are obtained through the following contractions of the strain rate tensor: The model contains two empirical constants, C h and C v , respectively, that need calibration.This has been carried out simulating the well-recognized, turbulent-plane Poiseuille flow (also referred to as plane channel flow) with increasing grid anisotropy, comparing first-and second-order statistics with reference experimental and numerical data (details are in Roman et al., 2010).In Fig. 1 we show the optimal coefficients obtained in the calibration process.
Reynolds analogy is assumed to hold for scalars, setting P r sgs = Sc sgs = 0.8.against field data and water circulation and mixing in the Muggia bay is discussed.Finally, concluding remarks are given in Section 4.

Governing equations
LES-COAST solves the Boussinesq form of the filtered, 3D non-hydrostatic, Navier-Stokes equations together with the transport equations for temperature and salinity, T and S respectively.Transport of pollutants treated through an advection-diffusion equation for their concentrations (or as dispersed Lagrangian particles) is also implemented in the model.The governing equations read as: where u i , p, T , and S are respectively the velocity component in the x i direction, the hydrodynamic kinematic pressure, temperature (in Kelvin) and salinity (in PSU); g is the gravitational acceleration, Ω i the i-th component of the Earth rotation vector, ρ and ρ 0 are density anomaly and its reference value respectively.ν, k T and k S are kinematic viscosity, thermal and salinity diffusivity respectively.Since the 200 density anomaly due to salinity and temperature variations is much smaller than the bulk density of the water, in coastal applications the buoyancy force can be accounted for by the Boussinesq approximation through the equation of state: where S ij is the resolved strain rate tensor and ν t is SGS eddy viscosity.

225
The Smagorinsky model is isotropic and eddy viscosity is evaluated as product of a length-scale C∆, proportional to the grid size, and a velocity scale C∆|S|, where C is a constant and |S| is the contraction of the resolved strain rate tensor.In moderate-to-strong anisotropic grids it is dif-230 ficult to define a single lengths-cale.In particular, setting ∆ = (∆ x ∆ y ∆ z ) 1/3 leads to an excessive overestimation of the eddy viscosity in all directions.Thus, in Roman et al. (2010) a two eddy-viscosity model was developed, one for the horizontal directions and another for the vertical one re-235 spectively.The two eddy viscosity are defined as: where scales for horizontal and vertical direction respectively.The velocity scales in the horizontal and vertical directions are The complex geometry of the coastal hydrodynamic problem is treated by a combination of a body-fitted curvilinear grid (Zang et al., 1994) and the immersed boundary method (IBM) recently developed by Roman et al. (2009b).
In IBM, the governing equations are solved over a regular Cartesian or curvilinear, structured computational domain; internal obstacles are reproduced through proper modification of the governing equations at these special locations.Among others, we use the direct forcing approach, meaning that a body force is added to the momentum equations to mimic the presence of solid boundaries (see Fadlun et al., 2000 and the successive paper by Balaras ( 2004) for a description of the methodology in Cartesian coordinates).In Roman et al. (2009b) the method was improved and extended to curvilinear-coordinate formulation of the Navier-Stokes equations.The advantage of IBM over unstructured mesh solvers is that particular care is not required to adapt the computational mesh to the domain shape, and complex geometry can be treated in an easy and efficient way.Our strategy is to adopt curvilinear coordinates to shape the computational domain over the physical one to minimize the number of inactive computational cells.Immersed boundaries are used to reproduce geometric complexities such as coastlines and anthropic structures such as docks, jetties and breakwaters.

Boundary conditions over solid boundaries
In wall-bounded turbulence, LES can be performed either by solving the near-wall viscous sublayer or through parametrization of the near-wall effects on the inertial region of turbulence.The computational cost of LES resolving the near-wall viscous layer (wall-resolving LES) is proportional to Re 2.5 (see Piomelli, 2008), limiting such a class of simulations to low values of Re.Moreover, in wall-resolving LES, it is not clear how to consider wall roughness explicitly.The alternative strategy, consisting of skipping the solution of the near-wall layer through parametrization of the wall shear stress, allows the application of the methodology to full-scale values of Re, also including the effect of wall roughness on the dynamics of the boundary layer.This is the strategy adopted in LES-COAST.Specifically, we use two different approaches, depending on whether the solid boundary is body-fitted or reproduced by using immersed boundaries.
For body-fitted solid boundaries the logarithmic law of the wall is used at the first node off the wall: where u + (1) denotes the tangential velocity at the first grid point off the wall, made non-dimensional with the friction velocity u τ , k = 0.41 is the von Karman constant, y + (1) is the distance from the wall of the first computational mesh point, scaled with the viscous length scale ν/u τ , and B is a coefficient that also includes roughness effects.This equation is solved iteratively to determine the friction velocity used to compute the wall shear stress, which is employed as a boundary condition.This technique can hardly be applied to solid walls reproduced using immersed boundaries, because in the general cases, cell faces do not coincide with immersed boundaries.In order to overcome this issue, a novel approach has been proposed by Roman et al. (2009a).It consists of a two-step procedure: first the velocity at the first off-the-wall node (with respect to the immersed boundary) is calculated through Eq. ( 13) using the velocity field from the interior; second, a RANS-like eddy viscosity is set at the interface between the fluid region and the solid one as ν t = ku τ y, where y is the distance from the surface of the immersed boundary and the first fluid node.Details of the method are in Roman et al. (2009a).

Boundary conditions at the free surface
In LES-COAST, wind forcing over the free surface is taken into account by means of the formula proposed in Wu (1982), in which the induced stress at the sea surface, τ , is computed from the wind velocity 10 m above the mean sea level U 10 as: where C 10 represents the drag coefficient.In the literature it is established that the kind of forcing at the free surface may play a role in the dynamics of the upper layer.Although in the literature better methods can be found, we use a zeromean random fluctuation with 20 % variance, added to the reference τ , to give a more realistic forcing action and, at the same time, to facilitate the generation of turbulence in the upper layers.Since randomization acts over the shear-stress components, it implicitly acts on the stress angle as well.This method has proven to be reliable and able to trigger turbulence in a very efficient way.In particular, by analysis of the instantaneous velocity field, we found sub-surface coherent structures typically developing in wall-bounded (or interfacial) turbulence.
Close to obstacles such as docks, ships and breakwaters, the wind strays locally from the mean path and gives rise to complex recirculating patterns.Such local phenomena can play an important role in the surface layer dynamics at harbor scale.A detailed reproduction of such effects would come from a fully-coupled simulation of the low-atmosphere dynamics together with water basin dynamics, with integrated boundary conditions able to conserve momentum and heat fluxes at the interface.Such an integrated model is not available at the moment and, on the other hand, it would require a computational cost at least as large as twice that required for the marine basin simulation.Rather, we employ a simplified, computationally inexpensive, engineering approach that is currently in use in numerical simulations of short-range pollutant dispersion in the low atmosphere (Scire et al., 2000).Specifically, LES-COAST incorporates local effects into spatial distribution of the wind shear stress, taking advantage of some relevant literature formulations concerning studies on turbulent flow patterns around three-dimensional obstacles on a plane surface.Following Hosker (1984) and Hanna et al. (1982), we denote with H and W , respectively, the characteristic height of the obstacle and its own length perpendicular to the wind direction, and with L its characteristic width in the cross-wind direction (see the upper right panel of Fig. 2).The relationships between the aspect ratios of the obstacle and the extension of the recirculation region are discussed in detail in Hosker (1984).The authors suggest: for the leeward recirculation length L lw , where  Hanna et al. (1982) and reads: regardless of other parameters.In the model the influence of the obstacles is considered through a linear damping of the wind stress along the recirculation pattern, from the reference value at distance L lw or L ww from the obstacle to zero at the solid boundary.In case of time-varying wind conditions, the model calculates successive damping-factor maps.A linear interpolation procedure computes the updated values in the intermediate steps between two consecutive maps.2. This allows to generate a divergence-free fluctuating field which can be superimposed to the incoming mean velocity field.In order to generate fluctuating velocity components representative of the turbu-410 lence level in the flow field, the intensity of the body force is automatically adjusted based on the mean velocity profiles, U RAN Si1 , and the T KE RAN S content coming from the interpolation procedure.The transport equation of T KE equation reads:

Boundary conditions at open boundaries
where, with ., we denote an averaging operation.In Equation 19 T R represent transport terms, Φ is the dissipation rate of T KE, the second term on the right hand side is mechanical production of T KE and the third one is synthetic 420 production of turbulence.
The production terms extract energy from the mean flow to feed the fluctuating counterpart.At the inlet, the mean velocity component is provided from the LCM and the production term is zero.Therefore its role is taken by the synthetic pro-  in a discrete form as: where, at the first time step, the velocity fluctuations are computed from the mean velocity profile as being c i a coefficient and R a suitable zero mean colored noise function.After the first time step, the body-force gives rise to velocity fluctuations in the buffer sub-domain, the shear production term starts to extract energy from the mean 440 flow and the forcing term must adapt itself to the flow conditions.This is accomplished by the proper definition in Equation 20 of the velocity scale, which is based on the actual fluctuation level in the flow: an increased value of u ′ i corresponds to a decreased value of b ′ i and vice-versa.The velocity 445 fluctuations can be easily computed because the mean flow is known at the boundary.
In order to provide a smooth transition from the buffer subdomain to the reference one, the forcing term is exponentially damped to a target value b T , which is zero, at the end of 450 the buffer region, subtracting to the function b Validation tests have been performed for the turbulent plane-channel flow case, at Re τ = u τ δ/ν = 20000 (with u τ = τ w /ρ the friction velocity, τ w the wall shear stress, 455 and δ the half height of the channel).It corresponds to a bulk Reynolds number of the order of 10 6 , (Pope, 2000, pp.279) , if any) at the inflow/outflow (i/o) boundaries of the LES domain.Since LCM simulations usually have larger time steps than LES, linear interpolation between two successive LCM time steps data is performed for the intermediate LES time steps; second, we place a buffer sub-domain between the i/o boundaries and the interior LES domain, in which turbulence is triggered by means of a divergence-free synthetic, zero-mean, fluctuating body force b i , added to the right-hand side of momentum Eq. ( 2).This allows one to generate a divergence-free fluctuating field that can be superimposed on the incoming mean velocity field.In order to generate fluctuating velocity components representative of the turbulence level in the flow field, the intensity of the body force is automatically adjusted based on the mean velocity profiles, U RANS i 1 , and the TKE RANS content coming from the interpolation procedure.The transport equation of TKE reads: where, with ., we denote an averaging operation.In Eq. ( 19) TR represents transport terms, is the dissipation rate of TKE, the second term on the right-hand side is mechanical production of TKE and the third one is synthetic production of turbulence.The production terms extract energy from the mean flow to feed the fluctuating counterpart.At the inlet, the mean velocity component is provided by the LCM and the production term is zero.Therefore its role is taken by the synthetic production term b i u i .
At the first time step, in the buffer sub-domain, velocity fluctuations are absent, therefore the left-hand side of Eq. ( 19) reduces to the target value for the turbulent kinetic energy (TKE RANS ) divided by the time step.All terms of the rhs are zero but b i u i .The body force can be constructed in a discrete form as:  2. This allows to generate a divergence-free fluctuating field which can be superimposed to the incoming mean velocity field.In order to generate fluctuating velocity components representative of the turbulence level in the flow field, the intensity of the body force is automatically adjusted based on the mean velocity profiles, U RAN S i 1 , and the T KE RAN S content coming from the interpolation procedure.The transport equation of T KE equation reads: where, with ., we denote an averaging operation.In Equation 19 T R represent transport terms, Φ is the dissipation rate of T KE, the second term on the right hand side is mechanical production of T KE and the third one is synthetic production of turbulence.
The production terms extract energy from the mean flow to feed the fluctuating counterpart.At the inlet, the mean velocity component is provided from the LCM and the production term is zero.Therefore its role is taken by the synthetic production term b ′ i u ′ i .At the first time step, in the buffer sub-domain, velocity fluctuations are absent, therefore the left hand side of Equation 19 reduces to the target value for the turbulent kinetic energy (T KE RAN S ) divided by the time step.All terms of the rhs are zero but b ′ i u ′ i .The body-force can be constructed 1 Here with the index RAN S we denote data coming from LCM simulations based on the solution of the Reynolds Averaged Navier-Stokes equations  in a discrete form as: where, at the first time step, the velocity fluctuations are computed from the mean velocity profile as being c i a coefficient and R a suitable zero mean colored noise function.After the first time step, the body-force gives rise to velocity fluctuations in the buffer sub-domain, the shear production term starts to extract energy from the mean In order to provide a smooth transition from the buffer subdomain to the reference one, the forcing term is exponentially damped to a target value b T , which is zero, at the end of 450 the buffer region, subtracting to the function b Validation tests have been performed for the turbulent plane-channel flow case, at Re τ = u τ δ/ν = 20000 (with u τ = τ w /ρ the friction velocity, τ w the wall shear stress, and δ the half height of the channel).It corresponds to a bulk Reynolds number of the order of 10 6 , (Pope, 2000, pp.279) , where, at the first time step, the velocity fluctuations are computed from the mean velocity profile as c i being a coefficient and R a suitable zero mean colored noise function.After the first time step, the body force gives rise to velocity fluctuations in the buffer sub-domain, the shear production term starts to extract energy from the mean flow and the forcing term must adapt itself to the flow conditions.This is accomplished by the proper definition in Eq. ( 20) of the velocity scale, which is based on the actual fluctuation level in the flow: an increased value of u i corresponds to a decreased value of b i and vice-versa.The velocity fluctuations can be easily computed because the mean flow is known at the boundary.
In order to provide a smooth transition from the buffer sub-domain to the reference one, the forcing term is exponentially damped to a target value b T , which is zero, at the end of the buffer region, subtracting to the function b Validation tests have been performed for the turbulent plane-channel flow case, at Re τ = u τ δ/ν = 20000 (with u τ = √ τ w /ρ the friction velocity, τ w the wall shear stress, and δ the half height of the channel).It corresponds to a bulk Reynolds number of the order of 10 6 (Pope, 2000, p. 279), typical of environmental flows.Such a configuration is well suited for testing the nesting methodology, since there is no external instability source present in the channel except for the turbulence generated in the buffer sub-domain (see Keating et al., 2004                A schematic representation of the case study is given in Fig. 3.The RANS inflow is applied at the left side of the buffer sub-domain, colored in gray, in which turbulent fluctuations are triggered.The indexes refer to different crosssections in which time-averaged velocity profiles are monitored along with second-order turbulence statistics.In the example depicted in Fig. 3, the buffer sub-domain spans from section 1 (A) to section 16 (H), whereas the reference subdomain begins at section 17 (I).pical of environmental flows.Such a configuration is well ited to test the nesting methodology, since there is no exrnal instability source present in the channel except for the rbulence generated in the buffer sub-domain (see Keating        typical of environmental flows.Such a configuration is well suited to test the nesting methodology, since there is no external instability source present in the channel except for the turbulence generated in the buffer sub-domain (see Keating et al. (2004) for a discussion).The results obtained with the nesting procedure are compared with those of a simulation of the same plane channel flow, carried out with an imposed driving streamwise pressure gradient and with periodic boundary conditions at the inlet/outlet sections (as in Cabot and Moin, 1999).In the latter case, turbulence is fully developed along the entire domain.Circulation and mixing in the bay of Muggia 520 Here, LES-COAST is applied to study water circulation and mixing in the bay of Muggia, the harbor area of the city of Trieste, Italy (see fig. 9).The bay develops along the estwest axis for about 4 km and communicates with the Gulf of Trieste through the western section, about 3 km wide; there, 525 the bay is delimited by three breakwaters that, from one side protect the bay from western sea storms, and, from the other side, inhibit water circulation.
The bathymetry displays a shallow water basin, with a maximum depth of approximately 20 m in the western re-530 gion close to the breakwaters, while in the eastern portion the depth decreases to about 5-10 meters.The northern side of the coast is characterized, from west to east, by docks of the harbor area of Trieste and by two wharfs belonging to an oil-terminal international company (SIOT-TAL S.p.A.).On 535 the north-eastern side, an artificial channel with the runoff  Stravisi (1990).
Different lengths of the buffer sub-domain were considered, while the length of the reference sub-domain and grid resolution were kept the same in all cases; hereafter we show the results of the channel with dimensions 3πδ × 2δ × 2 3 δ in the streamwise, wall-normal and spanwise directions, respectively, with 49×33×33 grid points along the three directions.The mesh size employed herein is typical of LES where a wall-layer approach is used (see Piomelli (2008) for a discussion).Here we use a standard Smagorinsky SGS model with C = 0.085.
To gain a better comprehension of the role of the body force in the present nesting procedure, the plane channel flow was first run without the forcing term.The mean velocity profiles at different sections in the streamwise direction are shown in Fig. 4. In all figures of the present section, black solid lines represent profiles at different streamwise cross sections, whereas dashed red lines represent reference values from the periodic channel flow simulation.
Without synthetic generation of turbulence (b i = 0), the flow in the channel tends to laminarize, the Reynolds stresses tend to vanish and the velocity profile tends to be parabolic.As discussed above, this behavior occurs because turbulent production in Eq. ( 19) is absent, and thus there are no terms able to extract energy from the mean flow to sustain turbulence.
Conversely, the synthetic generation of turbulence by means of the body-force approach allows matching of the reference values.Simulations with b i = 0 have been performed.Among the different lengths of the buffer sub-domain tested, we show in Figs. 5, 6, 7 and 8 the result for a buffer as long as 16 x, x being the cell dimension in streamwise direction.The plots show the mean velocity profile as well as the rms values of the streamwise, wall-normal and spanwise velocities, respectively.Only cross sections in the reference subdomain are shown.
The analysis of the plots clearly shows the effectiveness of the proposed methodology in triggering the correct level of turbulence within the reference sub-domain.The method supplies accurate results, also reducing the length of the buffer sub-domain (not shown), allowing one to reduce substantially the number of grid cells employed in the simulation.

Circulation and mixing in the bay of Muggia
Here, LES-COAST is applied to study water circulation and mixing in the bay of Muggia, the harbor area of the city of Trieste, Italy (see Fig. 9).The bay develops along the eastwest axis for about 4 km and communicates with the Gulf of Trieste through the western section, about 3 km wide; there, the bay is delimited by three breakwaters that from one side protect the bay from western sea storms, and, from the other side, inhibit water circulation.
The bathymetry displays a shallow water basin, with a maximum depth of approximately 20 m in the western region close to the breakwaters, while in the eastern portion the depth decreases to about 5-10 m.The northern side of the coast is characterized, from west to east, by docks of the harbor area of Trieste and by two wharfs belonging to an international oil-terminal company (SIOT-TAL S.p.A.).On the northeastern side, an artificial channel with the runoff of a creek (Torrente Rosanda) is present.The eastern side has a coastline characterized by the presence of an old, dismissed refinery; a small river runoff (Rio Ospo) is present at the southeastern part of the bay.On the southern side, the city of Muggia develops along the shoreline, with rocky beaches and small touristic marinas.In the simulations we do not consider the Rio Ospo and Torrente Rosandra runoff directly, being intermittent over the years and not well determined with respect to their own flow rates.However, as described in the following, their effect on the dynamics of the bay is indirectly accounted for through imposition of temperature and salinity profiles measured in the bay.

Numerical setup
The anthropic structures, represented in the model by IBM, are the three aforementioned breakwaters, the docks, the coastline and the two wharfs where, typically, oil tankers are moored.The numerical grid consists of 1024 × 512 cells on the horizontal plane and 24 cells along the vertical direction.up to 2 m in the complex region around the wharfs; along the vertical direction, grid spacing ranges from 1 meter in the deepest regions of the bay to 0.20 m close to the coastline.The Muggia bay is characterized by a diurnal/nocturnal breeze regime often interrupted by events of strong and cold wind blowing from the first quadrant (Bora).
Such events can persist from hours to days.Regarding the characterization of the wind breezes in the Gulf of Trieste, its complete overview is given in the studies of Stravisi (1990) and Stravisi et al. (1980) through a statistical analysis of field data retrieved in the last decades.
In this study we consider the values of wind velocity and direction probed every second at a meteorological station, located 28 m a.s.l. on the western wharf of the oil terminal.Three typical wind scenarios are identified and considered for the present study, namely the scenarios of the Ponente (west), Maestrale (northwest) and Bora (northeast) winds, each one characterized by its own direction and intensity.In particular, the blowing angle of the Bora wind in the Muggia bay is found to veer in the range 70 • -83 • , considering only the period with more intense wind conditions.The parameters used in the simulations are summarized in Table 1, where σ is the standard deviation.A cautionary principle has been applied in order to identify worst-case scenarios, with respect to successive water-quality analysis.For this reason for the Ponente and Maestrale scenarios, the values corresponding to the 90 • percentile of the sample have been considered, whereas the 70 • percentile has been considered for the Bora case.Wind-stress non-homogeneity due to the presence of ships and anthropic structures is considered through the simplified approach described in Sect.2.2.The characteristic lengths used in the simulations are reported in Table 2.
Data from the MITgcm simulation of the northern Adriatic sea including the Gulf of Trieste are used to provide suitable inflow/outflow conditions on the western side (and over the small fluid portion of the northern side) of the numerical domain.The available data are the daily averaged values of temperature, salinity and velocity fields.Data are processed by using the technique described and validated in Sect.2.3.From the MITgcm simulations database, three characteristic days have been chosen through a comparison with meteorological data provided in the database available at http: //www.dst.univ.trieste.it/OM/mens_TS/OM_mens.html.For Table 2. Characteristic lengths of the anthropic structures as adopted in the simulations.ter circulation and mixing in the bay based on the numerical results. 635

Validation with field data
In spite of the importance of the bay with respect to water quality in the Gulf of Trieste, due to the presence of potentially dangerous industrial plants, only few field campaigns have been carried over the years to quantify water circulation 640 and mixing.Specifically, velocity profiles have been measured under different wind conditions by Stravisi (1990).The author measured the vertical profile of module and direction of the water horizontal velocity at specific locations depicted in fig.9.They are placed close to the breakwaters, and be-645 tween the breakwaters and the coastline, at the mouths of the bay.Among the available data, for the comparison we extracted those relative to the summer period, with meteorological conditions comparable with those of the numerical simulations.The stations labels, positions, dates and preva-650 lent wind conditions are reported in Table 3. Figs.11a-f, show the comparison between numerical and field data velocity profiles at stations S1-S49 of fig. 9.The module and the direction of the horizontal velocity are respectively on the left and right panel.The general agreement 655 is fairly good, with velocity module somewhat underestimated and velocity angle better reproduced by the simulations.These discrepancies may be related to different causes, among the others, the presence of tides and seiches not reproduced by the model, some underestimation of the meso-scale 660 current obtained by the LCM close to the breakwaters and some differences between the wind conditions relative to the field data and the values used in our investigation.each wind scenario analyzed, we extract MITgcm data relative to the corresponding wind conditions.
Salinity and temperature profiles are initialized based on real data.Specifically, three separate experimental campaigns were carried out, respectively, during September 2011, March 2012 and June 2012; in these periods, measurements of temperature and salinity vertical profiles were taken below the western wharf of SIOT-TAL S.p.A.In all periods, the water column displays stable stratification, mainly because of low salinity in the upper layers.This is due to the combination of small water exchange with the open sea, caused by the breakwaters, and fresh water runoff coming from the two torrents in the eastern part of the bay.Although there are little differences in the stratification conditions along the year, here we consider the salinity and temperature profiles of June 2012.The variation in picnocline is found 6 m below the free surface, corresponding to the steep variations in both temperature and salinity.In Fig. 10 the field data are shown with symbols, while the curves obtained by a best-fit method are represented by solid and dot-dashed lines for temperature and salinity, respectively.The fitted data have been used to initialize the temperature and salinity fields into the domain.During June 2012, the temperature at the surface was 22.62 • C, and decreased by a rate of −0.16 • C per meter up to a depth of 6 m, then the rate changed to −0.082 • C per meter; at the free surface, salinity was 34.73 PSU, and increased by 0.12 PSU per meter up to the picnocline; beyond that, the rate changed to 0.035 PSU m −1 .For each wind scenario of Table 1, simulations were initially carried out until the achievement of a statistical steady state.After that, statistics were accumulated in time.Here we first show the results of a validation process against field data, and successively we will discuss the water circulation and mixing in the bay based on the numerical results.

Validation with field data
In spite of the importance of the bay with respect to water quality in the Gulf of Trieste, due to the presence of potentially dangerous industrial plants, only a few field campaigns have been carried over the years to quantify water circulation and mixing.Specifically, velocity profiles have been measured under different wind conditions by Stravisi (1990).The author measured the vertical profile of the module and direction of the water horizontal velocity at specific locations depicted in Fig. 9.They are placed close to the breakwaters, and between the breakwaters and the coastline, at the mouths of the bay.Among the available data, for the comparison we extracted those relative to the summer period, with meteorological conditions comparable with those of the numerical simulations.The station labels, positions, dates and prevalent wind conditions are reported in Table 3.
Figure 11a-f shows the comparison between numerical and field data velocity profiles at stations S1-S49 of Fig. 9.The module and the direction of the horizontal velocity are, respectively, on the left and right panel.The general agreement is fairly good, with velocity module somewhat underestimated and velocity angle better reproduced by the simulations.These discrepancies may be related to different causes, among others the presence of tides and seiches not reproduced by the model, some underestimation of the meso-scale current obtained by the LCM close to the breakwaters, and some differences between the wind conditions relative to the field data and the values used in our investigation.ter circulation and mixing in the bay based on the numerical results. 635

Validation with field data
In spite of the importance of the bay with respect to water quality in the Gulf of Trieste, due to the presence of potentially dangerous industrial plants, only few field campaigns have been carried over the years to quantify water circulation 640 and mixing.Specifically, velocity profiles have been measured under different wind conditions by Stravisi (1990).The author measured the vertical profile of module and direction of the water horizontal velocity at specific locations depicted in fig.9.They are placed close to the breakwaters, and be-645 tween the breakwaters and the coastline, at the mouths of the bay.Among the available data, for the comparison we extracted those relative to the summer period, with meteorological conditions comparable with those of the numerical simulations.The stations labels, positions, dates and preva-650 lent wind conditions are reported in Table 3. Figs.11a-f, show the comparison between numerical and field data velocity profiles at stations S1-S49 of fig. 9.The module and the direction of the horizontal velocity are respectively on the left and right panel.The general agreement 655 is fairly good, with velocity module somewhat underestimated and velocity angle better reproduced by the simulations.These discrepancies may be related to different causes, among the others, the presence of tides and seiches not reproduced by the model, some underestimation of the meso-scale 660 current obtained by the LCM close to the breakwaters and some differences between the wind conditions relative to the field data and the values used in our investigation.We carry out additional qualitative comparison, comparing steady-state streamlines obtained in simulations with trajectories of floating drifters, launched during the three field campaigns of September 2011, March 2012 and June 2012.The drifters have a tronco-conical shape, with a maximum diameter of 0.30 m and height of 0.20 m.Each drifter is equipped with three wings to keep the trajectory stable during the motion and to avoid rotation around its own vertical axis.It is equipped with a GPS satellite system, a GSM modem, a   away along the wind direction.The velocity vectors obtained in the simulation give full explanation of the observed behavior.In fact, the presence of a cyclonic eddy at the lee-715 ward side of the ship, obtained in the simulation, explains the initial trajectory of the drifter along a direction parallel to the ship.Successively the drifter enters the main southeast current moving along a weakly curved trajectory.The curvature of the trajectory is explained by the velocity vec-720 tors of the numerical simulations, showing that a south-east surface jet flowing across the gap between the wharfs and the coastline encroaches with the main eastern current caus- recording memory and a battery.Each trajectory has been remotely sampled every 2 min by the GPS system and transmitted via GSM to the operator, with an overall error smaller than one meter.The drifters were released in the proximity of the wharfs close to the ships moored over there, and left in the water for a time interval of up to 15 h.Often the drifters landed either at the shore or at the booms around the wharfs.
The qualitative comparisons between the trajectories of the drifters and numerical streamlines, represented through the velocity vector field, are discussed in the following.Among the many available trajectories, representative ones are considered for each wind scenario, in which wind conditions were found to be comparable with those considered in the simulations.In particular, we consider field data where wind is nearly constant for a period long enough.In the Figs.12-14 the trajectory of the drifter is depicted with green dots, the red spot is the release point of the drifter, and the black vectors represent the surface velocity obtained in the simulations.The vectors are skipped by 10 in both horizontal directions for the sake of visualization.The background pictures are satellite views taken from Google Maps, and the violet objects are IBM representations of oil tankers present at the wharfs.In the compass at the bottom right of the figures, the black arrow shows the wind direction used in the simulation, whereas the green area represents the range of variation in the wind direction during the field measurement.
In Fig. 12 the comparison for a Ponente scenario is depicted.The vector field shows a large cyclonic eddy between the two wharfs, and another one at the leeward side of the ship moored at the eastern wharf.This behavior is expected mainly due to the presence of the ship and, to a minor extent, due to the presence of oil barriers placed under the wharfs, submerged by about 0.50 m.On the other hand, the surface remotely sampled every 2 minutes by the GPS system and transmitted via GSM to the operator, with an overall error 675 smaller than one meter.The drifters have been released in proximity of the wharfs close to the ships moored over there, and left in the water for a time interval up to 15 hours.Often the drifters landed either at the shore or at the booms around the wharfs.

680
The qualitative comparisons between the trajectories of the drifters and numerical streamlines, represented though velocity vector field, are discussed in the following.Among the many available trajectories, representative ones are considered for each wind scenario, in which wind conditions were 685 found to be comparable with those considered in the simulations.In particular we consider field data where wind is nearly constant for a period long enough.In the figs.12-14 the trajectory of the drifter is depicted with green dots, the red spot is the release point of the drifter, and the black 690 vectors represent the surface velocity obtained in the simulations.The vectors are skipped by 10 in both horizontal direction for sake of visualization.The background pictures are satellite views taken from Google Maps, the violet objects are IBM representation of oil-tankers present at the wharfs.

695
In the compass at the bottom right of the figures, we show the wind direction used in the simulation by the black arrow, whereas the green area represents the range of variation of the wind direction during the field measurement.
In fig.
12 the comparison for a Ponente scenario is de-700 picted.The vector field shows a large cyclonic eddy between the two wharfs, and another one at the leeward side of the ship moored at the eastern wharf.This behavior is expected mainly due to the presence of the ship and, at a minor extent, due to the presence of oil barriers placed under the wharfs, 705 submerged by about 0.50 m.On the other hand, the surface current is directed east in the region between the mooring systems and the coastline, where oil barrier are absent.The drifter was released by the leeward side of an oil-tanker berthed at the eastern wharf, under wind direction fluctuat-710 ing between 317 • and 260 • .After the release, the drifter initially moved parallel to the ship, and successively it ran away along the wind direction.The velocity vectors obtained in the simulation give full explanation of the observed behavior.In fact, the presence of a cyclonic eddy at the lee-  ing a curvature of the streamlines (here represented through velocity vectors).The results of the numerical simulation ap-725 pear thus consistent with the drifter trajectory.Such accurate reproduction of the surface current arises from the high resolution employed on one side, and, on the other side, because of the reproduction of wind stress in-homogeneity adopted in the present work.

730
In fig.13 we show results for the Maestrale scenario.The field data refer to a drifter released under (north-east) Bora wind condition later on turned onto north-west (Maestrale) wind.The trajectory between the release point and the point in which the wind veers is not shown in the figure.After the 735 change of direction, the wind angle fluctuated between 202 • and 314 • .In the simulation the wind was considered with a direction of 315 • .The qualitative comparison shows that the trajectory of the drifter is aligned with the velocity vectors of the Maestrale scenario.In particular the drifter enters into the 740 cyclonic eddy between the two wharfs, and eventually lands at the boom at the base of the eastern one.
Finally in fig.14   current is directed east in the region between the mooring systems and the coastline, where oil barriers are absent.The drifter was released by the leeward side of an oil tanker berthed at the eastern wharf, under wind directions fluctuating between 317 • and 260 • .After the release, the drifter initially moved parallel to the ship, and successively it ran away along the wind direction.The velocity vectors obtained in the simulation give a full explanation of the observed behavior.In fact, the presence of a cyclonic eddy on the leeward side of the ship, obtained in the simulation, explains the initial trajectory of the drifter along a direction parallel to the ship.Successively the drifter enters the main southeast current, moving along a weakly curved trajectory.The curvature of the trajectory is explained by the velocity vectors www.nonlin-processes-geophys.net/20/1095/2013/ Nonlin.Processes Geophys., 20, 1095-1112, 2013 of the numerical simulation, showing that a southeast surface jet flowing across the gap between the wharfs and the coastline encroaches on the main eastern current, causing a curvature of the streamlines (here represented through velocity vectors).The results of the numerical simulation thus appear consistent with the drifter trajectory.Such accurate reproduction of the surface current arises from the high resolution employed on the one hand, and, on the other, because of the reproduction of wind stress inhomogeneity adopted in the present work.
In Fig. 13 we show results for the Maestrale scenario.The field data refer to a drifter released under the (northeast) Bora wind condition, later on turned into a northwest (Maestrale) wind.The trajectory between the release point and the point at which the wind veers is not shown in the figure.After the change of direction, the wind angle fluctuated between 202 • and 314 • .In the simulation the wind was considered with a direction of 315 • .The qualitative comparison shows that the trajectory of the drifter is aligned with the velocity vectors of the Maestrale scenario.In particular, the drifter enters into the cyclonic eddy between the two wharfs, and eventually lands at the boom at the base of the eastern one.
Finally, in Fig. 14 we show the surface current predicted under Bora regime, and its comparison with the trajectory of the drifter released under moderate Bora conditions.The wind angle is within the range of 70 • -95.7 • , whereas a mean direction of 83.7 • is set in the simulation.Apart from an initial mismatch between numerical and field data, the trajectory of the drifter is perfectly aligned with the mean streamlines of the numerical simulation.

Water circulation and mixing
Figure 15a-c shows the horizontal velocity vectors, averaged over a period of three hours, at three different planes, for the Maestrale scenario.
The figures clearly show the effect of the breakwaters, inhibiting water exchange between the Gulf of Trieste and the bay.At the surface level (Fig. 15a), outside the bay the current is pushed northeast in the northern region, and eastward in the southern part of the bay.Water is driven in the bay through all the inlets between the breakwaters but the northern one, through which surface water leaves the bay moving northward.Downstream from the breakwaters, within the bay, the surface current is dominated by the presence of upwind obstacles developing complex flow patterns.Welldetectable cyclonic eddies develop at the southern edges of the southern and mid-breakwaters, explaining the vertical velocity profile shown in Fig. 11a.Moving downstream within the bay, the surface current moves southeast at an angle of approximately 150 • , deviated by about 15 • with respect to the wind direction, due to the effect of the Coriolis force (for a discussion see Zikanov et al., 2003).In the southern part of the bay the surface current moves parallel to the coastline.Five meters below the free surface, see Fig. 15b, water from the Gulf of Trieste moves inside the bay through the southern channel, whereas the other inlets constitute outflow sections.The eddies observed at the free surface maintain their structure along the water column; in the downwind region close to the breakwaters the current moves north, parallel to the obstacles, almost opposite to the wind direction, leaving the bay through the three northern inlets.Within the bay the main current develops in a direction nearly opposite with respect to the surface one.Ten meters below the free surface, see Fig. 15c, the same pattern as at 5 m below the sea surface can be recognized.
Complex patterns develop in the northeast, close to docks and to the wharfs.A fast current through the inlet between the wharfs and the coastline can be noticed, changing direction going down along the water column.
Figure 16a, b and c show the horizontal vector velocity for the Ponente condition.Because of the wind angle, the currents at the mouths of the bay are stronger with respect to the Maestrale condition.In the core of the bay the main currents move along a direction of 118 • , deviating by 28 • with respect to the wind direction.
Five meters below the free surface, see Fig. 16b, the current flows within the bay through the inlets, apart at the northern mouth and close to the southern coast.A cyclonic eddy develops downward from the southern edge of the southern breakwaters and the currents do not exhibit a prevalent direction in the core of the bay.Conversely, in the eastern region, a large cyclonic eddy confined by the coastline is present, driving water from the southern region north along the eastern coastline.
Ten meters below the free surface, see Fig. 16c, the current moves northeastward, almost opposite to the wind direction, and leaves the bay through all the mouths except the southern one, where a prevalent direction is not well recognized due to the presence of a large cyclonic eddy developing along the whole inlet.The flow pattern described here lends justification to the velocity profile shown in Fig. 11b.
Finally, Fig. 17a-c shows the horizontal vector velocity for the Bora condition.
The surface current is almost aligned with the wind direction (see Fig. 17a) and leaves the bay through all mouths.In the core of the bay the currents move along a direction of 297 • with a deviation of 33 • north with respect to the wind.
Five meters below the free surface, see Fig. 17b, the current tends to reverse mainly at the southern mouth and in the eastern part of the bay.An along-shore current develops upward from the breakwaters, driving the current north, and a large anti-cyclonic eddy develops in the southeastern part of the bay, driving relatively high-speed water in the eastern and northeastern regions of the bay.
Ten meters below the free surface (see Fig. 17c), a current opposite to the wind direction is well established in the southern part of the bay, while several anti-cyclonic and cyclonic eddies characterize the northern part.Comparing, for example, the flow patterns obtained in the simulations with the velocity profile shown in Fig. 11f, at the free surface the current leaves the bay at an angle of 270 • , whereas from 2 to 6 m below the surface the current rotates to 90 • when entering the bay.At the deeper layer 10 m below the free surface, the flow is completely reversed and moves inside the bay.
The analysis of field data of Fig. 11d and Fig. 11e shows that the current tends to be parallel to the breakwaters moving in and out the bay, as described by our simulations.
Figure 18 shows a vertical section of the bay with a contour plot of the mean velocity in the west-east direction, together with contour lines of the density anomaly.The data have been averaged over a period of 3 h.In the figure, the Bora, Maestrale and Ponente conditions are shown from top to bottom.All cases exhibit flow-reversal along the water column, with velocity roughly aligned with the wind direction at the surface and moving along the opposite direction at deeper layers.Structure such as jetties, breakwaters and oil tankers have a blocking effect on wind action determining local upwelling and downwelling phenomena observable from the inclined contour lines of the density anomaly.
The current along the water column can be summarized as follows: under the Bora condition the surface layer is driven west, out of the bay, whereas the bottom layer moves east inside the bay.The upper layer has a mean velocity of the order of approximately 0.04 m s −1 , and is thinner and faster than the bottom layer that has a velocity of the order of 0.02 m s −1 , with larger values in the southern part of the bay, close to the coast.In the remaining cases the situation is similar, but the two layers flow in the opposite direction with respect to the Bora case.
The estimation of the mean velocity field under the three representative conditions allows evaluation of water renewal within the bay.Water renewal is accomplished by computation of the residence time R t , i.e., the average time required to flush the corresponding sub-domain.By definition the average residence or flushing time, R t , is the time required to replace a fluid volume in a coastal area, under the assumption of steady-state inflow/outflow conditions (Geyer et al., 2000) and a constant volume of the system R t = V /q, see Fischer et al. (1979), where V is the capacity of the system to hold the fluid and q is the flow rate through the system.Here we analyze two regions, illustrated in Fig. 9: the outer zone, delimited by the dashed line, and the inner zone, delimited by a dashed-dot-dot line.The computed values are summarized in Table 4.
For the Bora scenario the residence time for the inner zone is 1.28 days, and about 2 days for the outer zone; for the Ponente case the values are about 1 and 2.5 days for the inner and outer zones, respectively; in the Maestrale scenario it appears that the inner zone is flushed in 1.5 days approximately, while the outer zone requires 27 h only.(Geyer et al. (2000)) and a constant volume of the system R t = V /q, see Fischer et al. (1979), where V is the capacity of the system to hold the fluid and q is the flow rate through the system.ticed that these structures are not generated by the Stokes drift, not considered in the present simulation; rather they are generated by the free surface stress, giving rise to turbulent large-scale sub-surface coherent structures aligned with the direction of the surface stress(not shown), typical of wall-890 bounded or interface turbulent flows.These vertical eddies contribute to the generation of the vertical Reynolds shearstress u ′ v ′ 2 + w ′ v ′ 2 responsible of vertical mixing in the water column, which, on the other hand, can be quantified through the vertical eddy viscosity.

895
The vertical eddy viscosity is calculated as in Coleman et al. (1990):   The analysis of the instantaneous velocity field shows very interesting features with regard to mixing properties within the bay.In Fig. 19 the instantaneous velocity field is depicted along a south-north cross section for the three conditions analyzed.Counter-rotating eddies are well evidenced, spanning along the whole cross-sectional area.It has to be noticed that these structures are not generated by the Stokes drift, not considered in the present simulation; rather, they are generated by the free surface stress, giving rise to turbulent large-scale sub-surface coherent structures aligned with the direction of the surface stress (not shown), typical of wallbounded or interface turbulent flows.These vertical eddies contribute to the generation of the vertical Reynolds shear stress u v 2 + w v 2 responsible for vertical mixing in the water column, which, on the other hand, can be quantified through the vertical eddy viscosity.the corresponding sub-domains.By definition the averesidence or flushing time, R t , is the time required to rea fluid volume in a coastal area.Under the assumptions ady-state inflow/outflow condition (Geyer et al. (2000)) a constant volume of the system R t = V /q, see Fischer . (1979), where V is the capacity of the system to hold uid and q is the flow rate through the system.Here we ze two regions, illustrated in fig. 9  bay.In fig.19 the instantaneous velocity field is depicted along a south-north cross section for the three conditions analyzed.Counter-rotating eddies are well evidenced spanning along the whole cross-sectional area.It has to be no-885 ticed that these structures are not generated by the Stokes drift, not considered in the present simulation; rather they are generated by the free surface stress, giving rise to turbulent large-scale sub-surface coherent structures aligned with the direction of the surface stress(not shown), typical of wall-890 bounded or interface turbulent flows.These vertical eddies contribute to the generation of the vertical Reynolds shearstress u ′ v ′ 2 + w ′ v ′ 2 responsible of vertical mixing in the water column, which, on the other hand, can be quantified through the vertical eddy viscosity.

895
The vertical eddy viscosity is calculated as in Coleman et al. (1990): in which we consider the sum of the resolved part and the SGS part.For the three cases analyzed, fig.20 contains a 900 contour plot of the vertical eddy viscosity, together with contour lines of the density anomaly along a vertical section of the bay.The data has been averaged for a period of 3 hours.
In the panels from top to bottom the Bora, Maestrale and The vertical eddy viscosity is calculated as in Coleman et al. (1990): in which we consider the sum of the resolved part and the SGS part.For the three cases analyzed, Fig. 20 contains a contour plot of the vertical eddy viscosity, together with contour lines of the density anomaly along a vertical section of the bay.The data has been averaged for a period of 3 h.In the panels from top to bottom, the Bora, Maestrale and Ponente cases are shown, respectively.In all cases, it can be observed that larger values are present close to the walls, where mechanical production of turbulence is higher.Moreover, eddy viscosity values are higher in the bottom part of the water column, below the picnocline mainly associated with a relatively high Reynolds shear stress compared with the vertical mean shear.Above the picnocline the situation is the opposite: the strong density stratification inhibits the vertical Reynolds shear stress in a region where the mean shear is relatively high.The values obtained in our simulations 0.01-0.001m 2 s −1 are consistent with canonical literature data; however, our results clearly show that this quantity is far from being constant along the water column, spanning over two orders of magnitude from the free surface to the bottom.the eastern part, determine a strong mixing, also due to the 935 fact that the depth in the eastern part of the bay is lower with a stronger interaction between the bottom and upper boundary layer.Lower values can be observed in the upper layer close to the coastline, close to anthropic structures and around bottom bumps.Generally at these locations the vertical shear is enhanced and produces turbulent mixing working against the stable effect of density stratification.
The Ponente wind condition exhibits the smallest values for the Richardson number, as observed in Fig. 22, where horizontal sections at three different levels are shown.In this case, due to the wind angle and the bay configurations, a strong along-shore current develops in the southern part of the bay up to the eastern part, and determines a strong mixing, also due to the fact that the depth in the eastern part of the bay is lower, with a stronger interaction between the bottom and upper boundary layers.

Conclusions
In the present paper we present a state-of-the-art numerical model (LES-COAST) to simulate water renewal and mixing in closed or semi-closed regions.The model is suited to studying pollution dispersion in harbor areas or lakes.The model is unsteady and three-dimensional and solves the where mechanical production of turbulence is higher.Moreover eddy viscosity values are higher in the bottom part of the water column, below the picnocline mainly associated to a relatively high Reynolds shear-stress compared to the vertical mean shear.Above the picnocline the situation is opposite: the strong density stratification inhibits the vertical Reynolds shear-stress in a region where the mean shear is relatively high.The values obtained in our simulations 0.01 − 0.001m 2 /s are consistent with canonical literature data, however our results clearly show that this quantity is far from being constant along the water column, spanning over two orders of magnitude from the free surface to the bottom.Finally we show the vertical distribution of the gradient Richardson number Ri g (Ri g = N 2 /(dU/dy) 2 with N = (−g/ρ 0 ∂ρ/∂y) 1/2 the Brunt-Väisälä frequency) along a vertical plane for the three cases considered, see fig. 21.The water column appears strongly stratified with values larger than the critical value of 0.25 over most of the fluid column.
Lower values can be observed in the upper layer close to coastline, close to anthropic structures and around bottom bumps.Generally at these locations the vertical shear is enhanced and produce turbulent mixing working against the stable effect of density stratification.
Ponente wind condition exhibits the smallest values for the Richardson number, as observed in fig.22, where horizontal sections at three different levels are shown.In this case, due to the wind angle and the bay configurations, a strong alongshore current develops in the southern part of the bay up to  the eastern part, determine a strong mixing, also due to the 935 fact that the depth in the eastern part of the bay is lower with a stronger interaction between the bottom and upper boundary layer.where mechanical production of turbulence is higher.Moreover eddy viscosity values are higher in the bottom part of the water column, below the picnocline mainly associated to a relatively high Reynolds shear-stress compared to the   the eastern part, determine a strong mixing, also due to the 935 fact that the depth in the eastern part of the bay is lower with a stronger interaction between the bottom and upper boundary layer.Boussinesq form of the governing equations to take into account vertical stratification associated with temperature and salinity variation along the water column.Turbulent mixing is accomplished using large eddy simulation with a novel subgrid-scale, two-eddy-viscosity Smagorinsky model.Geometric complexity is treated through a combination of curvilinear coordinates and immersed boundaries.The model contains: a special treatment of free-surface wind stress, to take into account inhomogeneity deriving from the blocking effect caused by the presence of natural or anthropic structures; a new nesting procedure to assimilate data from LCM and to generate the proper level of turbulent kinetic energy in the flow field by synthetic turbulence production.
The methodology is employed here for the analysis of wind-driven circulation, water renewal and vertical mixing in the Muggia bay, the industrial harbor of the city of Trieste (Italy), for three typical wind scenarios.Validation was carried out in two different ways: first, we compared velocity profiles obtained in the simulation with corresponding profiles measured in field campaigns; second, we compared the trajectories of a floating drifter released into the Muggia bay with velocity streamlines obtained in our simulation under corresponding conditions.Overall the comparison is satisfactory.Specifically, the model reproduces well the vertical profile of the angle of the horizontal velocity component at the locations of the field measurements, whereas it slightly underestimates the module of the horizontal velocity.Among the possible causes of disagreement, it has been mainly attributed to the fact that the wind velocity used in our simulations was smaller than that of the field data.Further, the comparison between drifter trajectories and numerical data shows the ability of the model to reproduce local phenomena such as small-scale eddies associated with the presence of geometric complexities that can have a large impact on small-scale pollutant dispersion.
The simulations reproduce the inversion of the horizontal velocity direction along the vertical one; in the three cases examined, the upper layer is thinner and has a velocity of the order of 0.04 m s −1 , whereas the bottom layer is thicker, with a velocity of the order of 0.02 m s −1 .Further, the analysis shows the presence of instantaneous large-scale, crosssectional structures, deriving from the sub-surface elongated turbulent structures, which typically occur in wall-bounded turbulence.The calculation of vertical eddy viscosity shows that it may vary by two orders of magnitude along the water column, and larger values are observed in the bottom layer of fluid where density stratification is weaker and Reynolds shear stress is large compared with the vertical mean shear.The spatial distribution of the gradient Richardson number is similar in the three cases examined, with values larger than 0.25 in most of the water body, thus showing the scarce level of vertical mixing present in the bay.Locally, values of Ri g < 0.25 are found close to the coast, where upwelling/downwelling phenomena increase the vertical mean shear and enhance turbulent mixing.

Fig. 1 .
Fig. 1.Variation of optimal values constants C h and Cv with the grid anisotropy.
) 205 where ρ T is the total density, T 0 and S 0 are reference values giving the density ρ 0 , β T and β S are expansion coefficients for temperature and salinity respectively.C n in Equation 5 refers to the concentration of the n-th pollutant and k Cn is its kinematic diffusivity.Over-bar refers to the filtering opera-210 tion.LES-COAST uses implicit filter, meaning that filtering is implicitly performed in the physical space through a tophat filter function represented by the cell size.The latter term on the right hand side of Equations 2-5 represents the subgrid contributions to momentum and scalar fluxes respec-215 tively.The sub-grid scale momentum fluxes are parametrized by the two-eddy-viscosity Anisotropic Smagorinsky Model (ASM) described in Roman et al. (2010) and briefly summarized hereafter.In the standard Smagorinsky model, SGS stresses are ex-220 pressed as:

Fig. 1 .
Fig. 1.Variation in optimal value constants, C h and C v , with the grid anisotropy.

Fig. 2 .
Fig. 2. Schematic view of the recirculation zones of the wind flow around obstacles.Figure reproduced with modifications from Kuehn and Coleman (2005).
time step, in the buffer sub-domain, velocity fluctuations are absent, therefore the left hand side of Equation 19 reduces to the target value for the turbulent kinetic energy (T KE RAN S ) divided by the time step.All terms of 430 the rhs are zero but b ′ i u ′ i .The body-force can be constructed 1 Here with the index RAN S we denote data coming from LCM simulations based on the solution of the Reynolds Averaged Navier-Stokes equations

Fig. 4 .
Fig. 4. Comparison of the mean velocity profile for nested channel flow simulation without body-force.Red dashed line: reference results from the periodic channel simulation; black lines: mean velocity profile at different sections of the channel with the buffer region.

Fig. 3 .
Fig. 3. Sketch of the channel flow case with the buffer sub-domain.The sections used in the comparison are denoted with indexes.
Fig. 3. Sketch of the channel flow case with the buffer sub-domain.With indexes are denoted the sections used in the comparison.

Fig. 4 .
Fig. 4. Comparison of the mean velocity profile for nested channel flow simulation without body-force.Red dashed line: reference results from the periodic channel simulation; black lines: mean velocity profile at different sections of the channel with the buffer region.

440
flow and the forcing term must adapt itself to the flow conditions.This is accomplished by the proper definition in Equation 20 of the velocity scale, which is based on the actual fluctuation level in the flow: an increased value of u ′ i corresponds to a decreased value of b ′ i and vice-versa.The velocity 445 fluctuations can be easily computed because the mean flow is known at the boundary.

Fig. 4 .
Fig. 4. Comparison of the mean velocity profile for nested channel flow simulation without body force.Red dashed line: reference result from the periodic channel simulation; black lines: mean velocity profiles at different sections of the channel with the buffer region.
for a discussion).et al.: LES model for wind driven sea circulation in coastal areas 7

Fig. 5 .
Fig. 5. Comparison of the mean velocity profile for nested channel flow simulation with body-force.Red dashed line: reference results from the periodic channel simulation; black lines: mean velocity profile at different sections of the channel with the buffer region.

Fig. 6 .
Fig. 6.Comparison of the streamwise rms profile.Red dashed line: reference results from the periodic channel simulation; black lines: stream-wise rms profile at different sections of the channel with the buffer region.

Fig. 7 .
Fig. 7. Comparison of the wall normal rms profile.Red dashed line: reference results from the periodic channel simulation; black lines: wall-normal rms profile at different sections of the channel with the buffer region.

Fig. 8 .
Fig. 8.Comparison of the spanwise rms profile.Red dashed line: reference results from the periodic channel simulation; black lines: span-wise rms profile at different sections of the channel with the buffer region.

Fig. 5 .
Fig. 5. Comparison of the mean velocity profile for nested channel flow simulation with body force.Red dashed line: reference result from the periodic channel simulation; black lines: mean velocity profiles at different sections of the channel with the buffer region.

Fig. 5 .
Fig. 5. Comparison of the mean velocity profile for nested channel flow simulation with body-force.Red dashed line: reference results from the periodic channel simulation; black lines: mean velocity profile at different sections of the channel with the buffer region.

Fig. 6 .
Fig. 6.Comparison of the streamwise rms profile.Red dashed line: reference results from the periodic channel simulation; black lines: stream-wise rms profile at different sections of the channel with the buffer region.

Fig. 6 .
Fig. 6.Comparison of the streamwise rms profile.Red dashed line: reference result from the periodic channel simulation; black lines: streamwise rms profiles at different sections of the channel with the buffer region.
. Comparison of the mean velocity profile for nested channel w simulation with body-force.Red dashed line: reference results m the periodic channel simulation; black lines: mean velocity ofile at different sections of the channel with the buffer region.

Fig. 8 .
Fig. 8.Comparison of the spanwise rms profile.Red dashed line: reference results from the periodic channel simulation; black lines: span-wise rms profile at different sections of the channel with the buffer region.

Fig. 7 .
Fig. 7. Comparison of the wall-normal rms profile.Red dashed line: reference result from the periodic channel simulation; black lines: wall-normal rms profiles at different sections of the channel with the buffer region.

Fig. 5 .
Fig. 5. Comparison of the mean velocity profile for nested channel flow simulation with body-force.Red dashed line: reference results from the periodic channel simulation; black lines: mean velocity profile at different sections of the channel with the buffer region.

Fig. 6 .
Fig. 6.Comparison of the streamwise rms profile.Red dashed line: reference results from the periodic channel simulation; black lines: stream-wise rms profile at different sections of the channel with the buffer region.

Fig. 7 .
Fig. 7. Comparison of the wall normal rms profile.Red dashed line: reference results from the periodic channel simulation; black lines: wall-normal rms profile at different sections of the channel with the buffer region.

Fig. 8 .
Fig. 8.Comparison of the spanwise rms profile.Red dashed line: reference results from the periodic channel simulation; black lines: span-wise rms profile at different sections of the channel with the buffer region.

Fig. 8 .
Fig. 8.Comparison of the spanwise rms profile.Red dashed line: reference result from the periodic channel simulation; black lines: spanwise rms profiles at different sections of the channel with the buffer region.

Fig. 9 .
Fig. 9. Top view of the numerical domain.Two regions are defined for sake of computation of water-residence time, as discussed in Section 3.3: the Outer zone, delimited by the dashed line; the Inner zone, delimited by the dashed dot dot line.The black spots indicate the measurement points in the field campaign ofStravisi (1990).

Fig. 9 .
Fig. 9. Top view of the numerical domain.Two regions are defined for the sake of computation of water residence time, as discussed in Sect.3.3: the outer zone, delimited by the dashed line; the inner zone, delimited by the dashed-dot-dot line.The black spots indicate the measurement points in the field campaign ofStravisi (1990).
The resulting spatial resolution ranges, on the horizontal, from 10 m on the western boundary outside the breakwaters

Fig. 12 .
Fig. 12. Qualitative comparison of the numerical flow field, black vectors, respect to the drifter trajectory, green dots, for Ponente scenario.

Fig. 13 .
Fig. 13.Qualitative comparison of the numerical flow field, black vectors, respect to the drifter trajectory, green dots, for Maestrale scenario.

Fig. 12 .
Fig. 12. Qualitative comparison of the numerical flow field, black vectors, with respect to the drifter trajectory, green dots, for the Ponente scenario.

715Fig. 13 .
Fig. 13.Qualitative comparison of the numerical flow field, black vectors, with respect to the drifter trajectory, green dots, for the Maestrale scenario.
Fig. 14.Qualitative comparison of the numerical flow field, black vectors, respect to the drifter trajectory, green dots, for Bora scenario.

Fig. 14 .
Fig. 14.Qualitative comparison of the numerical flow field, black vectors, with respect to the drifter trajectory, green dots, for the Bora scenario.

Fig. 15 .
Fig. 15.Horizontal velocity vector field for Maestrale wind simulation averaged over a period of three hour, for clarity purpose only one vector is shown every 15 points in x and z. a) Sea surface; b) 5 meters below the sea surface; c) 10 meters below the free surface.

Fig. 15 .
Fig. 15.Horizontal velocity vector field for Maestrale wind simulation averaged over a period of three hours; for clarity purposes, only one vector is shown every 15 points in x and z.(a) Sea surface; (b) 5 m below the sea surface; (c) 10 m below the free surface.

Fig. 16 .Fig. 16 .
Fig. 16.Horizontal velocity vector field for Ponente wind simulation averaged over a period of three hour, for clarity purpose only one vector is shown every 15 points in x and z. a) Sea surface; b) 5 meters below the sea surface; c) 10 meters below the free surface.

Fig. 16 .Fig. 17 .
Fig. 16.Horizontal velocity vector field for Ponente wind simulation averaged over a period of three hour, for clarity purpose only one vector is shown every 15 points in x and z. a) Sea surface; b) 5 meters below the sea surface; c) 10 meters below the free surface.

Fig. 17 .Fig. 18 .
Fig. 17.Horizontal velocity vector field for Bora wind simulation averaged over a period of three hours; for clarity purposes, only one vector is shown every 15 points in x and z.(a) Sea surface; (b) 5 m below the sea surface; (c) 10 m below the free surface.

Fig. 19 .
Fig. 19.Contour plots of the west-east (u) instantaneous velocity component in the vertical cut indicated with black line in the small picture on the right side, for the three scenarios considered: a) Bora; b) Maestrale; c) Ponente.The black lines indicate instantaneous stream-tracers.
in which we consider the sum of the resolved part and the SGS part.For the three cases analyzed, fig.20 contains a 900 contour plot of the vertical eddy viscosity, together with contour lines of the density anomaly along a vertical section of the bay.The data has been averaged for a period of 3 hours.In the panels from top to bottom the Bora, Maestrale and Ponente cases are shown respectively.In all cases, it can be 905 observed that larger values are present close to the walls,

Fig. 18 .
Fig. 18.Contour plots of the west-east mean velocity component, averaged on a time of 3 h, in the vertical cut indicated in the small picture on the right.Solid lines represent the density anomaly.From top to bottom the plots refer to the Bora, Maestrale and Ponente scenarios.The vertical direction is magnified by a factor 100.
etronio et al.: LES model for wind driven sea circulation in coastal areas 15 18.Contour plots of the west-east mean velocity component, ged on a time of 3 hours, in the vertical cut indicated in the picture on the right side.Solid lines represent the density aly.From top to bottom the plots refer to Bora, Maestrale and nte scenarios.The vertical direction is magnified by a factor 4. Residence time of the two bay zones considered.

Fig. 19 .
Fig. 19.Contour plots of the west-east (u) instantaneous velocity component in the vertical cut indicated with black line in the small picture on the right side, for the three scenarios considered: a) Bora; b) Maestrale; c) Ponente.The black lines indicate instantaneous stream-tracers.

Fig. 19 .
Fig. 19.Contour plots of the west-east (u) instantaneous velocity component in the vertical cut indicated with a black line in the small picture on the right, for the three scenarios considered; from top to bottom: Bora, Maestrale and Ponente.The black lines indicate instantaneous stream tracers.

Fig. 20 .
Fig. 20.Contour plot for the vertical eddy viscosity, averaged on a time of 3 hours, in the vertical cut indicated in the small picture on the right side.Solid lines represent the density anomaly.From top to bottom the plots refer to Bora, Maestrale and Ponente scenarios.The vertical direction is magnified by a factor 100.

Fig. 21 .
Fig. 21.Contour plot for the gradient Richardson number, averaged on a time of 3 hours, in the vertical cut indicated in the small picture on the right side.From top to bottom the plots refer to Bora, Maestrale and Ponente scenarios.The vertical direction is magnified by a factor 100.

Fig. 22 .
Fig. 22. Contour plot for the gradient Richardson number, averaged on a time of 3 hours, at three horizontal levels.From top to bottom the plots refer to Bora, Maestrale and Ponente scenarios.

Fig. 20 .
Fig. 20.Contour plot of the vertical eddy viscosity, averaged on a time of 3 h, in the vertical cut indicated in the small picture on the right side.Solid lines represent the density anomaly.From top to bottom the plots refer to the Bora, Maestrale and Ponente scenarios.The vertical direction is magnified by a factor 100.
Fig. 20.Contour plot for the vertical eddy viscosity, averaged on a time of 3 hours, in the vertical cut indicated in the small picture on the right side.Solid lines represent the density anomaly.From top to bottom the plots refer to Bora, Maestrale and Ponente scenarios.The vertical direction is magnified by a factor 100.

Fig. 21 .
Fig. 21.Contour plot for the gradient Richardson number, averaged on a time of 3 hours, in the vertical cut indicated in the small picture on the right side.From top to bottom the plots refer to Bora, Maestrale and Ponente scenarios.The vertical direction is magnified by a factor 100.

Fig. 22 .
Fig. 22. Contour plot for the gradient Richardson number, averaged on a time of 3 hours, at three horizontal levels.From top to bottom the plots refer to Bora, Maestrale and Ponente scenarios.

Fig. 21 .
Fig. 21.Contour plot of the gradient Richardson number, averaged on a time of 1 h, in the vertical cut indicated in the small picture on the right side.From top to bottom the plots refer to the Bora, Maestrale and Ponente scenarios.The vertical direction is magnified by a factor 100.
Fig. 20.Contour plot for the vertical eddy viscosity, averaged on a time of 3 hours, in the vertical cut indicated in the small picture on the right side.Solid lines represent the density anomaly.From top to bottom the plots refer to Bora, Maestrale and Ponente scenarios.The vertical direction is magnified by a factor 100.

910
vertical mean shear.Above the picnocline the situation is opposite: the strong density stratification inhibits the vertical Reynolds shear-stress in a region where the mean shear is relatively high.The values obtained in our simulations 0.01 − 0.001m 2 /s are consistent with canonical literature 915 data, however our results clearly show that this quantity is far from being constant along the water column, spanning over two orders of magnitude from the free surface to the bottom.Finally we show the vertical distribution of the gradient Richardson number Ri g (Ri g = N 2 /(dU/dy) 2 with N = 920 (−g/ρ 0 ∂ρ/∂y) 1/2 the Brunt-Väisälä frequency) along a vertical plane for the three cases considered, see fig. 21.The water column appears strongly stratified with values larger than the critical value of 0.25 over most of the fluid column.Lower values can be observed in the upper layer close 925 to coastline, close to anthropic structures and around bottom bumps.Generally at these locations the vertical shear is enhanced and produce turbulent mixing working against the stable effect of density stratification.Ponente wind condition exhibits the smallest values for the 930 Richardson number, as observed in fig.22, where horizontal sections at three different levels are shown.In this case, due to the wind angle and the bay configurations, a strong alongshore current develops in the southern part of the bay up to

Fig. 21 .
Fig. 21.Contour plot for the gradient Richardson number, averaged on a time of 3 hours, in the vertical cut indicated in the small picture on the right side.From top to bottom the plots refer to Bora, Maestrale and Ponente scenarios.The vertical direction is magnified by a factor 100.

Fig. 22 .
Fig. 22. Contour plot for the gradient Richardson number, averaged on a time of 3 hours, at three horizontal levels.From top to bottom the plots refer to Bora, Maestrale and Ponente scenarios.

Fig. 22 .
Fig. 22. Contour plot for the gradient Richardson number, averaged on a time of 1 h, at three horizontal levels.From top to bottom the plots refer to the Bora, Maestrale and Ponente scenarios.

www.nonlin-processes-geophys.net/20/1095/2013/ Nonlin. Processes Geophys., 20, 1095-1112, 2013 Fig
. 3. Sketch of the channel flow case with the buffer sub-domain.With indexes are denoted the sections used in the comparison.ifany) at the inflow/outflow (i/o) boundaries of the LES domain.Since LCM simulations usually have larger time steps 400 than LES, linear interpolation between two successive LCM time steps data is performed for the intermediate LES time steps; second, we place a buffer sub-domain between the i/o boundaries and the interior LES domain, in which turbulence is triggered by means of a divergence-free synthetic, Comparison of the wall normal rms profile.Red dashed reference results from the periodic channel simulation; black wall-normal rms profile at different sections of the channel wi buffer region.
Comparison of the spanwise rms profile.Red dashed reference results from the periodic channel simulation; black span-wise rms profile at different sections of the channel wi buffer region.A schematic representation of the case study is giv fig.3.The RANS inflow is applied at the left side o buffer sub-domain, colored in gray, in which turbulent 465 tuations are triggered.The indexes refer to different c sections in which time averaged velocity profiles are Comparison of the streamwise rms profile.Red dashed line: ference results from the periodic channel simulation; black lines: eam-wise rms profile at different sections of the channel with the ffer region.
Comparison of the wall normal rms profile.Red dashed line: reference results from the periodic channel simulation; black lines: wall-normal rms profile at different sections of the channel with the buffer region.

Table 1 .
Wind statistics from data recorded from 12-16 September 2011

Table 1 .
Wind statistics from data recorded from 12 to 16 September 2011.

Table 4 .
Residence time of the two bay zones considered.

Table 4 .
For the Bora scenario the residence time for the Inner zone is 1.28 days, and about 2 days for the Outer 875 zone; for the Ponente case the values are about 1 and 2.5 day for the Inner and Outer zones respectively; in the Maestrale scenario it appears that the Inner zone is flushed in 1.5 days approximately, while the Outer zone requires 27 hours only.The analysis of instantaneous velocity field shows very in- 880teresting features with regard mixing properties within the

Table 4 .
Residence time of the two bay zones considered.