Jeffrey Mirocha

Pierre-Antoine Joulin
email: pierre-antoine.joulin@ifpen.fr

Maria Laura Mayol

Valéry Masson

Frédéric Blondel

Quentin Rodier

Marie Cathelain

Christine Lac

The Actuator Line Method in the Meteorological LES Model Meso-NH to Analyze the Horns Rev 1 Wind Farm Photo Case

Keywords: wind turbine wake, wind farm wake, meteorological interaction, actuator line method, large-eddy simulation, new MEXICO experiments, horns Rev

The offshore wind market is expected to see massive growth worldwide in the next decade. Recent research in this sector has led to a new generation of wind turbines characterized by larger rotors and higher capacity factors. Furthermore, the trend is to develop larger and denser wind farms which leads to potential deeper interactions with the troposphere. In order to provide an accurate modeling of the processes between the wind turbines, the surface and the atmospheric boundary layer, a new numerical tool has been developed. The ALM has been implemented into the opensource non-hydrostatic mesoscale atmospheric model Meso-NH used in the Large-Eddy Simulation (LES) framework. First, the tool is validated by comparing simulated results with data acquired during the New MEXICO experiments. In particular, good correlation are obtained for normal and tangential loads along the blades and the axial PIV traverse of instantaneous wind velocity components. Once the tool is validated, its ability to reproduce the Horns Rev 1 photo case is evaluated, since it is an emblematic case of wind farm-atmosphere interaction. Simulation output are post-processed using a render algorithm, resulting in synthetic images of cloud radiance. These images show great similarity with the photographs taken. All the results highlight the capacity of this new tool to represent interactions between wind farms and the lower atmosphere with a high level of details.

INTRODUCTION

Considerable deployments of offshore wind farms can be expected in the near future (IEA, 2019). After many years of development, the swept areas of modern, large rotors have increased substantially, resulting in wind plants influencing deeper portions of the atmospheric boundary layer. In fact, turbulent wakes produced by wind turbines can significantly affect the flow dynamics within wind farms and downstream areas. Operating offshore farms have already shown losses on electricity production [START_REF] Barthelmie | Modelling and measuring flow and wind turbine wakes in large wind farms offshore[END_REF] and impacts on the local meteorology [START_REF] Hasager | Wind farm wake: the horns rev photo case[END_REF].

To investigate these interactions, a better understanding of processes occurring between the atmospheric boundary layer, the surface and the wind turbines is necessary. This investigation can be done with numerical tools, among which Navier-Stokes based Computational Fluid Dynamics (CFD) models are frequently used. Several approaches can be implemented to solve the governing flow equations, such as the Reynolds averaged Navier Stokes (RANS), the Large-Eddy Simulation (LES), or the Direct Numerical Simulation (DNS). The RANS method computes mean flow and models an important part of turbulence aspects: it is not the proper tool to explore small scale wake eddies nor atmospheric boundary layer eddies. The DNS method is more precise but computationally too expensive as it resolves all turbulent scales. Because the LES approach consists in solving the large scales and modeling the smaller ones, it can be seen as an intermediate approximation and a good compromise.

In order to introduce wind turbines into a CFD software, different approaches can be used. One of those approaches is to simulate a fully-resolved blade, but this technique is highly resource-consuming for a whole wind farm for the moment. Introduced by [START_REF] Sørensen | Numerical modeling of wind turbine wakes[END_REF], the Actuator Line Method (ALM) represents the three blades and their motion with a body force approach. Compared to the Actuator Disc Model from [START_REF] Rankine | On the mechanical principles of the action of propellers[END_REF] and [START_REF] Froude | On the part played in propulsion by differences of fluid pressure[END_REF], the ALM allows a better description of the wake behind the wind turbines, by reproducing the helical shape of the unsteady wake and capturing tip and root vortices as mentioned by [START_REF] Mikkelsen | Actuator Disc Methods Applied to Wind Turbines[END_REF], [START_REF] Ivanell | Numerical computations of wind turbine wakes[END_REF][START_REF] Troldborg | Actuator Line Modeling of Wind Turbine Wakes[END_REF]. This technique has already been coupled and validated with different LES solvers to simulate wind turbine and wind farms flows [START_REF] Porté-Agel | Largeeddy simulation of atmospheric boundary layer flow through wind turbines and wind farms[END_REF][START_REF] Churchfield | A large-eddy simulation of wind-plant aerodynamics[END_REF]. It has also been implemented into the Weather Research and Forecasting model (WRF) by [START_REF] Marjanovic | Implementation of a generalized actuator line model for wind turbine parameterization in the weather research and forecasting model[END_REF].

The present paper introduces a new numerical tool allowing the study of wind farms impacts on local meteorology. This is achieved by coupling the ALM with the meteorological model Meso-NH [START_REF] Lac | Overview of the meso-nh model version 5.4 and its applications[END_REF]. This model can simulate the atmospheric boundary layer using an LES framework, taking into account atmospheric phenomena such as buoyancy, thermal exchanges, humidity, and even cloud formation and precipitation.

The study is organized as follows: the numerical methods are described and validated in section 2. Section 3 introduces the Horns Rev 1 photo case, provides the simulation details and presents the results. Concluding remarks are provided in section 4.

METHODS

Meso-NH Model

Meso-NH is a non-hydrostatic mesoscale atmospheric model [START_REF] Lac | Overview of the meso-nh model version 5.4 and its applications[END_REF] developed by the Centre National de Recherches Météorologiques (CNRM -Météo France/CNRS) and the Laboratories d'Aérologie (LA -UPS/CNRS). It can simulate atmospheric flows from meso to micro-scales. The governing equations are based on the conservation equations of mass, momentum, humidity and the thermodynamic equation derived from the conservation of entropy. A description of the dynamical core can be found in [START_REF] Lafore | The meso-nh atmospheric simulation system. part i: adiabatic formulation and control simulations[END_REF]. The model is anelastic and based on the pseudo-compressible system of [START_REF] Durran | Improving the anelastic approximation[END_REF], filtering elastic effects from sound waves.

The numerical model can be used with a Large-Eddy Simulation framework. The second-order turbulent fluxes are computed through an 1.5-order closure proposed by Redelsperger andSommeria (1981, 1986) and discussed in detail by [START_REF] Cuxart | A turbulence scheme allowing for mesoscale and large-eddy simulations[END_REF]. The closure scheme is based on a prognostic equation of the subgrid turbulent kinetic energy and a diagnostic mixing length. The mixing length of [START_REF] Deardorff | Stratocumulus-capped mixed layers derived from a threedimensional model[END_REF] was chosen for this study, taking into account the mesh size and the stratification of the atmospheric boundary layer.

The approach is Eulerian and the domain is discretized using a staggered Arakawa C-grid [START_REF] Mesinger | Numerical methods used in atmospheric models[END_REF]. Wind advection is treated using a 4th-order centered spatial scheme combined with the explicit Runge Kutta scheme for time integration. Computations are horizontally parallelized [START_REF] Jabouille | Parallelization of the french meteorological mesoscale model mésonh[END_REF], and the vertical coordinate can be stretched to increase the resolution in the region of interest.

Actuator Line Method

The ALM consists in modeling each blade by one line divided into blade element points [START_REF] Sørensen | Numerical modeling of wind turbine wakes[END_REF]. These points are applying aerodynamic forces into the flow. Each point carries a two-dimensional (2D) airfoil as shown in Figure 1, where c denotes the chord and -→ U rel the relative wind velocity. The relative velocity is a combination of the wind speed ---→ U wind and the opposite of the tangential blade element speed -→ U , as seen by the airfoil. For a blade element point, the force -→ dF is the sum of the drag force -→ dF D acting in the direction of -→ U rel (along the unit vector -→ e D ), and the lift force -→ dF L , perpendicular to -→ U rel (along the unit vector -→ e L ).

The aerodynamic forces -→ dF are evaluated using the lift coefficient C L and the drag coefficient C D :

-→ dF = 1 2 ρc -→ U rel 2 (C L -→ e L + C D -→ e D )dr, (1)
where ρ and dr are respectively, the air density and the length of the blade element. Both the lift and drag coefficients (C L and C D ) depend on the airfoil characteristics, and are evaluated through a cubic spline interpolation of the tabulated data.

Coupled System

Meso-NH and the Actuator Line Method are exchanging information with each other. At each time step, the position, the velocity and the orientation matrix of blade element points are evaluated. Then, the ALM can use air properties (as ρ and ---→ U wind ) from the three-dimensional (3D) fields of Meso-NH at the current time step in order to compute the aerodynamic force -→ dF (Equation 1). The wind velocity ---→ U wind is obtained using a tri-linear interpolation of the eight neighboring cell faces, surrounding the blade element point. Then, the aerodynamic force is projected in the global frame and added to the momentum equation of Meso-NH:

m ∂ ---→ U wind ∂t = -→ Adv + -→ Cor + -→ Pres + -------→ Turb( ---→ U wind ) - -→ dF, (2)
where -→ Adv is the advection term, -→ Cor the Coriolis term, -→ Pres the pressure gradient term, and --→ Turb the turbulence term. To avoid numerical instabilities, a commonly used regularization kernel is used [START_REF] Sørensen | Numerical modeling of wind turbine wakes[END_REF][START_REF] Mikkelsen | Actuator Disc Methods Applied to Wind Turbines[END_REF]. Further studies also give guidelines for optimal volume force distribution regarding the kernel shape [START_REF] Jha | Guidelines for volume force distributions within actuator line modeling of wind turbines on large-eddy simulation-type grids[END_REF][START_REF] Martínez-Tossas | A highly resolved large-eddy simulation of a wind turbine using an actuator line model with optimal body force projection[END_REF][START_REF] Nathan | Near-wake analysis of actuator line method immersed in turbulent flow using large-eddy simulations[END_REF]. Nevertheless, in this study, the forces are distributed smoothly over nearby flux points by computing the spatially averaged force of the two closest mass points. It shows similar results, avoiding instabilities and the computational cost of a convolution product. Otherwise, the tip loss correction of [START_REF] Glauert | Airplane propellers[END_REF] is used to alleviate the over predicted loads at the blade tip region, as proposed by [START_REF] Mikkelsen | Actuator Disc Methods Applied to Wind Turbines[END_REF]. One can note that this correction should only be used with the Actuator Disc Model to correct for finite number of blades: an improved load distribution will be investigated in future work, as advised by [START_REF] Churchfield | An advanced actuator line method for wind energy applications and beyond[END_REF].

To go further in computational cost savings, a time-splitting method was introduced, also known as the Actuator Sector Method [START_REF] Storey | An actuator sector method for efficient transient wind turbine simulation[END_REF]. The Courant-Friedrichs-Lewy (CFL) criterion of Meso-NH imposes a time step t MNH . Nevertheless, the ALM often requires a smaller time step in order to ensure that a blade element point will not skip a mesh cell during this time step (otherwise the information of the blade passage will be missed at this cell). The criterion for the Actuator Line Method time step t ALM can thus be expressed as:

t ALM < min( x, y, z) r tip , (3)
where x, y, and z are the cell dimensions, the angular speed of the rotor and r tip the tip radius of the blade. Because the meteorological model does not need such a small time step, a t MNH respecting the CFL criterion -and not moreis preserved. In parallel, the ALM algorithm is called N split times over this duration in order to respect the ALM time step criterion (Equation 3) and to reduce by N split -1 the number of iterations from the LES solver. One can note that during a Meso-NH time step, the wind field is frozen: the same velocity vectors are sent to the blades at each sub-time step to compute the forces. In the end, for each blade element point, the applied force is the aerodynamic force evaluated by the ALM divided by N split to respect the momentum budget over the time increment t MNH .

Validation

The New MEXICO Experiments

To assess the accuracy of the results obtained with this coupled system, a validation study based on the New MEXICO (Model Experiments in Controlled Conditions) experiments is carried out [START_REF] Boorsma | New MEXICO experiment: Preliminary Overview With Initial Validation[END_REF]. These experiments have been performed in the German Dutch Wind Tunnel, using the 9.5 × 9.5 m open section of the Large Scale Low Speed Facility.

During these experiments, several flow cases were investigated [START_REF] Boorsma | New MEXICO experiment: Preliminary Overview With Initial Validation[END_REF]. The validation has been done for the three MexNext cases, but this study only focuses on the 10 m/s case corresponding to a high tip speed ratio (TSR ≈ 10). Indeed, in the work presented in section 3, the wind turbines were operating at a same TSR or higher.

The wind tunnel wind turbine is a three-bladed 4.5 m diameter rotor. The twisted and tapered blades are based on DU91-W2-250, RISØ-A1-21, and NACA64-418 airfoils from root to tip. The wind turbine is rotating at a constant rotational velocity = 44.55 rad/s and a constant blade pitch angle β = -2.3 • . The inflow is laminar and no initial turbulence is introduced. The normal and tangential loads are [START_REF] Boorsma | New MEXICO experiment: Preliminary Overview With Initial Validation[END_REF], CASTOR and SOWFA simulations [START_REF] Blondel | Improving a bem yaw model based on newmexico experimental data and vortex/cfd simulations[END_REF].

available at five positions thanks to the integration of pressure sensors measurements.

Numerical Set-Up

The dimension of the domain is about L x = 40 m × L y = 30 m × L z = 30 m, with a x = y = z = 0.1 m resolution, leading to N x = 400 × N y = 300 × N z = 300 uniformly spaced grid points. The discretization of 45 cells across the rotor diameter is thus consistent with the recommendations of [START_REF] Jha | Guidelines for volume force distributions within actuator line modeling of wind turbines on large-eddy simulation-type grids[END_REF]. For the Actuator Line Method, 42 elements are used to describe each blade. The hub is placed in the domain at (15, 15, 15 m). The wind turbine is placed at a high vertical level with the aim of avoiding the effect of the tiny boundary layer developing at the ground. The lateral and upper boundaries are open conditions [START_REF] Bougeault | The meso-nh Atmospheric Simulation System: Scientific Documentation[END_REF]. The upstream boundary generates a constant wind flow. The nacelle and the tower are not taken into account. As mentioned in section 2.3, the average method for the force distribution and the tip loss correction are used. A time step about t = 0.0005 s is chosen such that the blade tip does not pass through more than one grid cell. The time-splitting method will be investigated afterwards.

The wind speed, the potential temperature and the pressure are set as inputs according to the experimental data : respectively u = 10.05 m/s, θ = 293.65 K, p = 1013.98 hPa [START_REF] Boorsma | New MEXICO experiment: Preliminary Overview With Initial Validation[END_REF]. The density at the hub height, computed by Meso-NH regarding the aforementioned data, corresponds to the measured one (i.e., ρ h = 1.197 kg/m 3 ). The wind velocity profile is constant and the thermal conditions are neutral.

Validation Results

The numerical results are compared to the experimental data. The loads along the blade span (shown in Figure 2) and axial Particle Image Velocimetry (PIV) traverses (shown in Figure 3) are used. The obtained results are also compared to CASTOR and SOWFA simulations [START_REF] Blondel | Improving a bem yaw model based on newmexico experimental data and vortex/cfd simulations[END_REF]. CASTOR is a free wake, vortex filament, lifting-line solver developed at IFPEN, and SOWFA is a library based on the CFD solver OpenFOAM developed at the National Renewable Energy Laboratory [START_REF] Churchfield | Overview of the Simulator for Wind Farm Application (sowfa)[END_REF].

Figure 2 shows the time-averaged normal and tangential loads along a blade after 7.552 s of simulation (i.e., more than 50 rotor rotations). The average loads are evaluated during the last rotation of the rotor (i.e., t average = 0.141 s). The computed normal load F N fits well with the comparative data. However, the computed tangential load F T shows an overestimation for the second half of the blade. As explained in [START_REF] Boorsma | New MEXICO experiment: Preliminary Overview With Initial Validation[END_REF], the experimental tangential load is more difficult to obtain and is often overpredicted by numerical models. In fact, important uncertainties of tangential loads evaluation from pressure taps measurements are mentioned. The results can also be affected by the adopted assumptions, as the consideration of a 2D flow over airfoils, the use of non Reynolds-dependent airfoils, or the absence of the nacelle. Nevertheless, these results are in quite good agreement with both experimental and numerical data for the most part.

In Figure 3, the axial PIV traverse of instantaneous wind velocity components (u, v, and w) after 7.552 s of simulation are shown. These measurements are taken along an axial line ranging from x = -4.5 m to x = +5.9 m for a hub considered at position x = 0 m. This line crosses the rotor plane at radius r = 1.5 m in the 9 o'clock plane when looking from upstream, and when a blade is at azimuth 0 • (upper vertical position). Numerically, it is not evident to reach this azimuth with a finite time step: at 7.552 s, the blade is at 2.1 • .

The upper plot of Figure 3 shows the evolution of u component. A very good correlation between the Meso-NH results and the experiment can be observed. Below, both v and w components show good agreements in the induction area upstream the wind turbine, but do not fully capture the oscillating wake. It might be due to the mesh resolution of 10 cm, which is not enough to take into account smaller vortices. Furthermore, an underestimation of w can be noticed downstream the rotor, as for the others numerical results [START_REF] Boorsma | New MEXICO experiment: Preliminary Overview With Initial Validation[END_REF]. In any case, even if the nacelle and the tower are not modeled, and in spite of the quite low resolution, the results fit quite well with the experimental data.

Time-Splitting Results

In order to verify the conformity of the time-splitting method, another numerical study is carried out. The loads along the blades and the axial traverses are compared to the experimental data for different time-splitting factors N split . The numerical set-up is the same as the one presented in 2.4.2, only the time step t MNH is modified to obtain the desired time-splitting, as presented in Table 1. The duration about 7.552 s is conserved, but not for case D. In fact, this time is not reachable with a time step of 0.006 s.

Because of that, an additional rotation is done for case D (i.e., 7.692 s), in order to reach an azimuth close to 0 • . Figure 4 shows a slight overestimation of the calculated loads when N split increases. Nevertheless, the loads evaluation seems correct. The curve for case D unsplit shows what would happen without the time-splitting method, and assess its functionality. A deterioration of the wake deficit can be observed in Figure 5 for the different cases. In fact, an important value of N split leads to higher amplitude of oscillations. It might be due to the deformation of the numerical modeling of the blade. When the time-splitting is used, the modeled blade can be seen as a larger blade applying aerodynamic forces on a wider sector with less intensity.

Table 1 also compares the computational time, evaluated during the last rotation of the wind turbine. It should be mentioned that the writing time of the output files is also taken into account. In any case, the gain obtained is directly linked to the meteorological solver, since unnecessary calculation steps have been avoided. The difference between D and D unsplit about 3.6% is mainly due to the sub-calls of the ALM algorithm. According to the results, a factor of N split = 7 seems to be a good compromise: it is enough to save computational cost and obtain a representative wake.

APPLICATION TO THE HORNS REV 1 PHOTO CASE

Description

The Horns Rev 1 wind farm is well known for being the first large scale offshore wind farm in the world. It is located in the North Sea, west of Denmark and consists of 80 Vestas V80 wind turbines with 70 m hub height, 80 m rotor diameter and rated power of 2.0 MW. The parallelepiped shaped layout of the wind farm is arranged in 10 vertical and 8 horizontal rows, with a minimal distance of 7 diameters between the closest turbines.

In the morning of February 12, 2008 at around 10:10 UTC, the photographer Christian Steiness took two photos of the Horns Rev 1 wind farm from a helicopter (Figures 12A,13A). These two photographs showed a physical evidence of the interaction between the wind farm and the lowest levels of the atmosphere, illustrating wind turbine-induced cloud formation downwind the wind farm.

A detailed description of the meteorological conditions at the time the photos were taken and the process that would have lead to the clouds formation can be found in [START_REF] Hasager | Wind farm wake: the horns rev photo case[END_REF]. This rare phenomenon is a consequence of the combination of very specific conditions, involving a layer of cold and humid air above a warmer sea surface that re-condensates to fog in the turbines wakes. The torque exerted by the flow on the wind turbine blades induces rotation in the wake, causing upward movements of humid warm air from near the sea surface and downward movements of relatively cooler and dryer air. The rotor-induced updraft causes air parcel to rise, as the pressure decreases with height, the air parcel expands and its temperature decreases. Therefore, the relative humidity of the lifted parcel increases. At some height, the parcel becomes saturated and a subsequent condensation of water vapor leads to cloud formation. In summary, the large scale structure of the wake fog is produced by the rotational pattern of the spiraling bands.

Despite the fact that Horns Rev I is one of the most studied wind farms in the world, with three meteorological masts installed on-site, one of the main difficulties to simulate the photo case is the need of very precise information to reproduce the specific conditions that led to the wake fog formation. [START_REF] Hasager | Wind farm wake: the horns rev photo case[END_REF] conducted an exhaustive compilation of much of the available information, including satellite data, radiosounding profiles and surface observations. The mentioned on-site data come from the meteorological tower (M6) installed in the park which provides measurements of temperature at 16, 64, and -3 m (sea temperature). This met mast also gives measurements of wind velocity at three different heights (30, 40, and 70 m) and wind direction at 28 and 68 m. To have a more complete description of the vertical structure of the atmosphere for the region and moment of interest, the present work also used reanalysis data from the Modern-Era Retrospective analysis for Research and Applications (MERRA) [START_REF] Rienecker | Merra: Nasa's modern-era retrospective analysis for research and applications[END_REF].

The wake fog simulation is performed using Meso-NH through two steps. The objective of the first step is to establish the precursor meteorological conditions without the wind turbines. The second one includes the wind farm that leads to the wake clouds development.

Precursor Simulation

In order to start the spin-up simulation, various vertical profiles of the initial state of the atmosphere must be set. In Figure 6, the vertical profiles (horizontally and 10-min averaged) of wind speed, temperatures, and vapor mixing ratio, used to initialize the whole domain (blue dashed lines) are shown. The profiles after 2 h of simulation (blue solid lines) are also plotted. The potential temperature θ is defined as:

θ = T (P 0 /P) 0.286 , ( 4
)
where P is the air pressure and P 0 a reference pressure set to 1,000 hPa. The initial conditions were deduced from the met mast (M6) and the reanalysis data. Because of its position, the met mast M6 is considered to be not affected by the wake, nor by the fog development. For the wind, several uniform conditions (i.e., intensity and direction) have been tested to initialize the atmospheric state. Indeed, during the simulation, the atmospheric flow will tend to an equilibrium between the geostrophic wind forcing, the surface friction and the Coriolis effect. The chosen value for the initial profile corresponds to the geostrophic wind and allow to obtain the measured intensity and direction provided in [START_REF] Hasager | Wind farm wake: the horns rev photo case[END_REF]. As it can be seen in Figure 6, the 10 min averaged wind profile after 3 h of simulation shows good agreement with the met mast data at the moment before the photo, 10:00 UTC. At this time, by taking into account the reanalysis data, a boundary layer height of 400 m is determined. A strong inversion of potential temperature is established at this level. In the aim of maintaining the sea smoke under the rotor level and because it is the beginning of the day, a surface layer inside the boundary layer with a weaker inversion must be set. After a sensitivity analysis, a surface layer with a height of 90 m is introduced with a second and weaker inversion of θ = -0.4 K. It was shown that for any value greater than the selected one, no cloud would develop: the air mass would become too warm in the upper part of the rotor to condense. As it is shown in Figure 6, both surface and residual layers are neutral at the starting point of the simulation, and later, after 3 h of simulation, they become slightly stable but still distinguishable. It can also be noticed that the simulated temperature is in good agreement with the observed temperature from the mast near the surface. The initial condition for the sea temperature has been obtained from the met mast data at -3 m. The reported value at 10:00 UTC was 277.85 K with a ± 0.7 K precision. After a sensitivity study, this value was set to 277.35 K. A warmer sea would lead to a more unstable air layer, with greater associated thermal turbulence which would lead to a cloudy development that would not correspond to the layered sea smoke that can be observed in the photo (Figures 12A,13A). The sea level pressure is about 103890 Pa. As there are no in-situ data for the vapor mixing ratio r v , the MERRA reanalysis data have been used to define its order of magnitude. In Figure 6, a typical shape of r v profile can be observed, with higher values in the atmospheric boundary layer, and lower ones in the free atmosphere.

The precursor simulation is configured with a domain of 7,500 × 25,000 × 837.5 m, a horizontal resolution of 5 m and a time step of 1 s. The vertical levels are defined with a vertical resolution of 5 m below 200 m and then a 20% stretching is applied until a resolution of 20 m is reached. The boundary conditions are cyclic. The long-wave and short-wave radiation are computed following [START_REF] Fouquart | Computations of solar heating of the earth's atmosphere-a new parameterization[END_REF][START_REF] Mlawer | Radiative transfer for inhomogeneous atmospheres: Rrtm, a validated correlated-k model for the longwave[END_REF], respectively. To deal with the micro-physical phenomena that enable clouds formation, the single moment 6-class scheme ICE3 is selected (see [START_REF] Pinty | A mixed-phase cloud parameterization for use in mesoscale non-hydrostatic model: simulations of a squall line and of orographic precipitations[END_REF]. Finally, the sea surface is represented by using the Sea Flux Scheme detailed in [START_REF] Belamari | Report on uncertainty estimates of an optimal bulk formulation for surface turbulent fluxes[END_REF].

Wind Farm Simulation

Set-Up

Once the precursor 3D state of the atmosphere is established, the Horns Rev 1 wind farm is placed in the domain and a second simulation is performed for a period of 32 min. The domain configuration is the same as the previous one with 1,500 points in the -→ x direction, 5,000 points in the -→ y direction and 75 levels in the vertical one with 16 levels within the rotor area. As shown in Figure 7, this configuration allows to have a distance of 5 km upstream the wind farm, and more than 15 km downstream to observe the evolution of the wake. On each side, a distance of 1,250 m (i.e., 15 rotor diameters) separates the wind farm from the lateral boundaries. The layout is based on the data from [START_REF] Hansen | Presentation of horns rev offshore wind farm and vestas v80 wind turbine[END_REF].

The prevailing wind direction is 180 • at the hub height (i.e., parallel to -→ x ) and the yaw angle of the wind turbines is set to zero. The Vestas V80 is a three-bladed rotor. The twisted and tapered blades are based on FFA-W3-301, FFA-W3-241, FFA-W3-221, NACA-63-221, and NACA-63-218 airfoils from root to tip, according to [START_REF] Hansen | Presentation of horns rev offshore wind farm and vestas v80 wind turbine[END_REF]. The different airfoil properties can be found in [START_REF] Fuglsang | Wind tunnel tests of the FFA-W3-241[END_REF], [START_REF] Abbott | Theory of Wing Sections, Including a Summary of Airfoil Data[END_REF] and [START_REF] Bertagnolio | Wind Turbine Airfoil Catalogue[END_REF]. The SCADA data of the whole Horns Rev 1 wind farm is reported in [START_REF] Hasager | Wind farm wake: the horns rev photo case[END_REF]. The rotational velocities of the wind turbines recorded at 10:00 UTC are used. The pitch angle is set to 0 • because the turbines operated near the cut-in speed. For the Actuator Line Method, the blades are divided into 42 element points and the tip loss correction is used [START_REF] Glauert | Airplane propellers[END_REF]. The time step is set to 0.25 s, using a time-splitting factor of 7.

Results

After 32 min of simulation with the wind farm, several variables are studied to analyze the wind farm interaction with the local meteorology.

The relative humidity horizontal cross-sections, at 30 m and 110 m above the sea are presented in Figure 8. On the two figures, a quite high value of relative humidity (97-98%) can be observed around the wind farm. Differences can be noticed under and over the wind farm wake. At the 30 m level, the relative humidity values are mostly lower than the surroundings. An opposite behavior is shown at the 110 m level, with higher values of RH. This configuration is the consequence of the upward and downward displacements of air masses produced by the rotational pattern of the operating rotors, and can be explained by analyzing the vapor mixing ratio and the potential temperature. A moisture transport is clearly observable within the region of the wake, in Figure 9, where an area of lower moisture content in relation to its environment for low levels and higher moisture content for high levels can be observed. The wake of the wind turbines causes a vertical homogenization of the vapor mixing ratio. This homogenization is also noticeable in Figure 10, where the potential temperature cross-sections are plotted: the potentially cooler air mass under the farm is mixed with the potentially warmer air mass above it, in the wake of the turbines. The combination generates within the wake a potentially warmer air mass with less vapor below the farm and a potentially cooler air mass with more vapor above.

From these figures (Figures 8-10), it is also possible to identify the Ekman spiral effect as the wake is slightly tilted to the west in the 30 m plot and, on the contrary, deflected toward the east in the level of 110 m.

In Figure 11, the vertical cross-section of the cloud mixing ratio (r c ) through the fourth row of the wind farm is shown. The development of the wake cloud above the hub height along the row is observed. As in the photo, behind each first wind turbines, a constant cloud about 2 or 3 diameters long can be observed at the top of the rotor. Downstream these clouds, at the tip, a vortex shedding is also distinguishable on the photos, and reproduced by the numerical model. Thereafter, eddies periodically interact with downstream wind turbines and become bigger. After a few downstream wind turbines, the chaotic structure of the cloud becomes more established. An underrepresentation of the clouds can be noticed: the sea smoke recognizable on the photo is not reproduced, even if high values of relative humidity are reached.

Figures 12B, 13B are synthetic images of atmospheric fields (temperature, pressure, liquid and water vapor), simulated by the open source htrdr Monte Carlo code [START_REF] Villefranque | Path-tracing monte carlo libraries for 3d radiative transfer in cloudy atmospheres[END_REF]. To render these realistic cloud scenes, a virtual camera is placed approximatively at the positions of the camera of Christian Steiness when he took the photos (Figures 12A,13A). From the comparison between Figures 12A,B, the capacity of the Meso-NH/ALM tool to reproduce the Horns Rev I photo case phenomenon can be observed. The cloud structures presented in the figures are very similar, with the wakes widening downstream and the spiral pattern noticeable in each row.

Between the photo taken from the southeast (Figure 13A) and the post-processing render (Figure 13B) there is a great agreement. As the first turbine of the second row was not working at the time of the photo, it was also stopped in the simulation. This explains the absence of cloud wake downstream of that turbine in both the photo and the render.

CONCLUSION

As wind power production is increasing all over the world and the wind turbines are getting larger and taller over the years, the study of their interaction with the lowest levels of the atmosphere becomes relevant. To address this challenge, a new numerical tool was developed to explore the impacts of wind farms on the local meteorology. This new tool is a coupled model between the atmospheric mesoscale model Meso-NH used in the LES framework and the Actuator Line Method. In order to validate the tool, comparisons against New MEXICO measurements and other numerical models have been made. The predicted loads along the blades and the axial wind intensity at hub height showed good correlation with the measurements for the high tip speed ratio (TSR) value. The time-splitting method has also been validated, allowing to save computation time and obtain a representative wake.

With the goal of evaluating the capability of this new tool to reproduce a real case of interaction between a wind farm and the atmospheric boundary layer, a simulation of an idealized case of the Horns Rev I Photo case was performed. First, a spin-up simulation was carried out to reproduce the existing weather conditions before the photo, based on in situ met-mast measurements and MERRA reanalysis data. The resulting atmospheric 3D state showed good agreement with the measurements. This atmospheric field was then used as initial underestimated by the model, especially at the lower levels. It may be improved by taking into account aerosols such as the sea salts. However, the mechanism that gave rise to the wake cloud is well-described by the model, with rotating rotors carrying potentially cooler and wetter air from lower levels; and potentially warmer and drier air from higher ones. The condensation that occurred in the upward movement leads to the development of the cloud.

It should be underlined that the phenomenon showed by the Horns Rev 1 photo case can only appear under a specific set of atmospheric conditions. In fact, the sensitivity studies carried out without the wind farm (not presented) showed that small changes in temperatures would have led to a global fog formation or no cloud at all.

Based on the results of this study, the capability of this new tool to simulate the wake dynamics and to represent their impact and interaction with the atmospheric quantities was shown. Furthermore, this work demonstrates more generally the ability of coupled weather and wind plant simulation frameworks to capture their complex physical interactions.

Moreover, the atmospheric boundary layer is typically shallower over the sea, while offshore wind turbines are getting larger. It will then become possible for blades to cross the top of the boundary layer, which is often characterized by an abrupt temperature inversion and strong wind shear. At that point, simulations of new generations of wind farms will be conducted with the aim of investigating those interactions.

FIGURE 1 |

 1 FIGURE 1 | 2D blade element velocities and forces.

FIGURE 2 |

 2 FIGURE 2 | Comparison of computed normal F N and tangential F T forces per unit length along the blades for the 10 m/s case with the New MEXICO experiments[START_REF] Boorsma | New MEXICO experiment: Preliminary Overview With Initial Validation[END_REF], CASTOR and SOWFA simulations[START_REF] Blondel | Improving a bem yaw model based on newmexico experimental data and vortex/cfd simulations[END_REF].

FIGURE 3 |

 3 FIGURE 3 | Comparison of axial PIV traverses for the case 10 m/s, with the New MEXICO experiments[START_REF] Boorsma | New MEXICO experiment: Preliminary Overview With Initial Validation[END_REF], CASTOR and SOWFA simulations[START_REF] Blondel | Improving a bem yaw model based on newmexico experimental data and vortex/cfd simulations[END_REF].

FIGURE 4 |

 4 FIGURE 4 | Comparison of predicted normal F N and tangential F T forces per unit length along the blades for the 10 m/s case with the New MEXICO experiment (Boorsma and Schepers, 2014) for different time-splitting factors.

FIGURE 5 |

 5 FIGURE 5 | Comparison of axial PIV traverses for the case 10 m/s with the New MEXICO experiment (Boorsma and Schepers, 2014) for different time-splitting factors.

FIGURE 6 |

 6 FIGURE 6 | Wind speed, temperature and vapor mixing ration profiles. M6 Data from Hasager et al. (2013).

FIGURE 7 |

 7 FIGURE 7 | Layout of the Horns Rev 1 offshore wind farm. The blue bullets (•) indicate the location of wind turbines. The black crosses (+) indicate the met masts position. The thick black box denotes the domain boundaries.

FIGURE 8 |FIGURE 9 |

 89 FIGURE 8 | Horizontal cross-sections at 30 m (A) and 110 m (B) above the sea of the relative humidity RH. The black plain lines (-) delimits the locations of wind turbines.

FIGURE 10 |FIGURE 11 |

 1011 FIGURE 10 | Horizontal cross-section at 30 m (A) and 110 m (B) above the sea of the potential temperature θ. The red plain lines (-) indicate the locations of wind turbines.

FIGURE 12 |

 12 FIGURE 12 | (A) Photograph of the Horns Rev 1 wind farm (12/02/2008 around 10:10 UTC) from south (Hasager et al., 2013). Courtesy Vattenfall. Photographer: Christian Steiness. (B) Post-processing (see Villefranque et al., 2019) of simulated cloud mixing ratio. Wind turbines are not represented.

FIGURE 13

 13 FIGURE 13 | (A) Photograph of the Horns Rev 1 wind farm (12/02/2008 around 10:10 UTC) from southeast (Hasager et al., 2013). Courtesy Vattenfall. Photographer: Christian Steiness. (B) Post-processing (see Villefranque et al., 2019) of simulated cloud mixing ratio. Wind turbines are not represented.

TABLE 1 |

 1 Time step values used to obtain different time-splitting factors.

	Case		A	B	C	D	E	D unsplit
	t MNH	[s]	0.0005 0.002 0.004 0.006 0.008	0.006
	N split	[-]	1	3	5	7	9	1
	Computational time [%]	100	40.4	33.3	27.4	27.0	23.9

Frontiers in Earth Science | www.frontiersin.org

January 2020 | Volume 7 | Article 350

ACKNOWLEDGMENTS

New MEXICO data have been supplied by the consortium which carried out the EU FP5 project MEXICO: Model rotor Experiments In Controlled conditions. We wish to thank the IEA Wind Task 31 WakeBench which provided the description of the Horns Rev 1 wind farm and of the Vestas80 wind turbine model. We would like to acknowledge the Wake Conference 2019 for publishing our preliminary results under the IOP Proceedings Licence (https://iopscience.iop.org/article/10.1088/1742-6596/ 1256/1/012019). We also acknowledge Vattenfall to let us publish the photographs (https://creativecommons.org/licenses/ by-nd/4.0/).

FUNDING

MM received funding from the Make Our Planet Great Again (MOPGA) initiative for a research stay at CNRM.

AUTHOR CONTRIBUTIONS

P-AJ developed the numerical coupling, performed the numerical experiments, and led the analysis of the study. MM gathered all the meteorological data available for the Hors Rev 1 Photo case, launched the firsts Horns Rev simulations and participated to the redaction of the manuscript. VM, FB, QR, MC, and CL substantially contributed to the development and the manipulation of the different models. All authors contributed on the comprehension of the study, as well as on the revision and correction of the manuscript. Agreement to be accountable for all aspects of the work in ensuring that.

Conflict of Interest:

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2020 Joulin, Mayol, Masson, Blondel, Rodier, Cathelain and Lac. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.


Alexis Mathian

Suzanne Mouries-Martin

Karim Dorgham

Hans Yssel

Zahir Amoura

Response to: 'Circulating interferon-α measured with a highly sensitive assay as a biomarker

Of note, the area under the ROC curve value of 0.81 for IFNαconcentration as a discriminator between active versus inactive juvenile IIM determined by Melki and colleagues is similar to the value of 0.83 that we had previously reported in SLE. In Melki and colleagues' study, the optimal cutoff threshold of IFNαconcentration permitting the ability to distinguish active from inactive juvenile IIM was 42 fg/ml, which is not only substantially lower than the value observed in our study (266 fg/ml) as being optimal to distinguish active from inactive SLE, but also lower than the value (136 fg/ml) calculated as the threshold below which IFNαvalues are within the range of those in healthy controls. The discrepancy between the 2 studies might be explained by the different types of antibodies used in the digital IFNα ELISAs, the different types of disease studied, as well as by age and demographic differences between the patient populations. It is therefore clear that the implementation of the digital IFNα ELISA to meet clinical grade standards still requires further investigation. Moreover, the determination of the specific roles of type I IFNs versus type II and type III IFNs in the pathogenesis of various autoimmune diseases will be facilitated by the use of this highly sensitive digital ELISA technology that enables the measurement of cytokines at physiologic concentrations.

The prognostic value of elevated levels of IFNs in the prediction of a disease flare also remains to be studied. It is therefore important that this technology be developed for the measurement of other types of IFNs, particularly the IFNβ and IFNλ subtypes. In conclusion, the translation of digital ELISA technology into the clinical laboratory setting will enhance the potential for routine measurement of IFNα in the monitoring of several "nonmonogenic" interferonopathies.

  Response to: 'Circulating interferon-α measured with a highly sensitive assay as a biomarker for juvenile inflammatory myositis activity.' by Melki et al Alexis Mathian, MD, PhD, Suzanne Mouries-Martin, MD, MSc, Karim Dorgham, PhD, Hans Yssel, PhD, Zahir Amoura, MD, MSc Sorbonne Université, AP-HP Hôpital Pitié-SalpêtrièreFrench National Referral Center for SystemicLupus Erythematosus, Antiphospholipid AntibodySyndrome and Other Autoimmune DisordersService de Médecine Interne 2, Institut E3M, INSERM UMRS, CIMI-Paris, Paris, France We thank Dr. Melki and colleagues for their interest in our work on IFNα assessment in SLE. In their correspondence, they show that IFNα protein measurement by digital ELISA is also correlated with disease activity in other IFN-dependent autoimmune diseases, and in this case, juvenile IIMs. Since IFNs (especially IFNα) are involved in the pathogenesis of several autoimmune diseases such as SLE, myositis, Sjögren's syndrome, and systemic sclerosis (1), monitoring their presence in the circulation of patients might help physicians to better evaluate disease activity and/or predict future flares. However, until recently, clinicians lacked a robust and standardized method for the measurement of IFNs in serum because IFN concentrations are usually very low and often not detectable with the use of classic immunoassays. Tools for indirect measurement of IFNα activity through an IFN gene expression score have a very low availability. These scores are not only notoriously difficult to standardize but are also less sensitive to change than direct measurement of IFNα. The advent of the ultrasensitive Simoa digital ELISA technology now makes it possible to reliably measure very low concentrations of IFNα (2). The studies by Melki et al and ourselves show that direct measurement of IFNα might improve monitoring of clinical activity in myositis and SLE, respectively.

  1. Muskardin TL, Niewold TB. Type I interferon in rheumatic diseases. Nat Rev Rheumatol 2018;14:214-28 2. Rodero MP, Decalf J, Bondet V, Hunt D, Rice GI, Werneke S,et al. Detection of interferon α protein reveals differential levels and cellular sources in disease. J Exp Med 2017;214:1547-55)