Wind Forecasting Based on the HARMONIE
Model and Adaptive Finite Elements
Albert Oliver
1
, Eduardo Rodríguez
1
, José María Escobar
1
,
Gustavo Montero
1
, Mariano Hortal
2
, Javier Calvo
2
,
José Manuel Cascón
3
, and Rafael Montenegro
1
1
University Institute for Intelligent Systems and Numerical
Applications in Engineering (SIANI), University of Las Palmas de
Gran Canaria (ULPGC) Campus de Tafira, 35017 Las Palmas de Gran
Canaria, Spain http://www.dca.iusiani.ulpgc.es/proyecto2012-2014
2
Agencia Estatal de Meteorología (AEMET), Leonardo Prieto Castro,
8, 28040 Madrid, Spain
3
Department of Economics and Economic History, Faculty of
Economics and Business, University of Salamanca, 37007 Salamanca,
Spain
Abstract
In this paper we introduce a new methodology for wind field forecasting over
complex terrain. The main idea is to use the predictions of the HARMONIE meso-
scale model as the input data for an adaptive finite element mass-consistent wind
model. The HARMONIE results (obtained with a maximum resolution of about
1 km) are refined in a local scale (about a few metres). An interface between both
models is implemented in such a way that the initial wind field is obtained by a
suitable interpolation of the HARMONIE results. Genetic algorithms are used to
calibrate some parameters of the local wind field model in accordance to the HAR-
MONIE data. In addition, measured data are considered to improve the reliability
of the simulations. An automatic tetrahedral mesh generator, based on the mec-
cano method, is applied to adapt the discretization to complex terrains. The main
characteristic of the framework is a minimal user intervention. The final goal is to
validate our model in several realistic applications in Gran Canaria island, Spain,
with some experimental data obtained by the AEMET in their meteorological sta-
tions. The source code of the mass-consistent wind model is available on-line at
http://www.dca.iusiani.ulpgc.es/Wind3D/
1
1 Introduction
Over the last years the use of wind power to produce electric power has augmented
considerably. Wind models are tools that allow the study of several problems related
to the atmosphere, such as, the effect of wind on structures, pollutant transport Oliver
et al. [2012, 2013], fire spreading Ferragut et al. [2007], wind farm location Rodríguez
[2004], etc.
In this paper we propose a methodology for wind forecasting by coupling the HAR-
MONIE meso-scale model with a local mass-consistent wind model specially suited for
complex terrain Rodríguez et al. [2012]; similar coupling methods have been proposed
by Gasset et al. [2012] and Carvalho et al. [2013]. HARMONIE is used experimentally
at AEMET with promising results Navascués et al. [2013]. Despite the high-resolution
of the HARMONIE meso-scale model, the minimum horizontal resolution is about
1 km, which is a drawback when the micro-scale (about 1 m) is considered. For this
reason the results of the HARMONIE meso-scale model are coupled with the local
wind field model. An initial wind field is required: it is obtained by a vertical extrapo-
lation and a horizontal interpolation. The vertical extrapolation is based on a log-linear
wind profile Lalas and Ratto [1996]. Both the mass-consistent model and the interpo-
lation are defined by a set of parameters. Some of these parameters are known, but
others have to be estimated. In order to calibrate these parameters, genetic algorithms
are used Montero et al. [2005]. Algorithm 1 synthesises the main steps of the model.
Algorithm 1 Overall algorithm
1. Mesh generation with the Meccano method
2. Assimilation of HARMONIE weather meso-scale model data for its use in the
local wind field model
3. Calibration of the wind field model parameters using genetic algorithms
This paper is organised as follows. Sections 2, 3 and 4 explain in detail the main
parts of this work: the meccano mesh generation (Section 2), the HARMONIE meso-
scale model (Section 3) and the local wind field model (Section 4). Section 5 discusses
the genetic algorithms and the parameters to be estimated. Section 6 shows an experi-
ment of this method applied to Gran Canaria island.
2 Meccano Mesh Generation
The main steps of the meccano tetrahedral mesh generation algorithm are summarized
in this section. This method has been previously introduced in Montenegro et al. [2009,
2010] and Cascón et al. [2013]. The input data of the algorithm is the definition of the
solid boundary (for example a surface triangulation or CAD description) and a given
precision (corresponding to the approximation of the solid boundary). Algorithm 2
describes our mesh generation approach.
The first step of the procedure is to construct a meccano approximation by connect-
ing polyhedral pieces. The meccano and the solid must be equivalent from a topological
point of view, i.e., their surfaces must have the same genus.
2
Algorithm 2 Meccano tetrahedral mesh generation
1. Meccano: Construct a meccano, M , approximation of the solid, Ω, formed by
polyhedral pieces.
2. Mapping: Define an admissible mapping, Π, between the meccano boundary
faces, ∂M , and the solid boundary, ∂ Ω, i.e., Π : ∂ M → ∂ Ω.
3. Coarse Mesh: Construct a coarse tetrahedral mesh, T
0
(M ),of the meccano.
4. Refined Mesh: Generate a local refined tetrahedral mesh, T (M ), from T
0
(M ),
such that the surface triangulation, τ(Ω) , obtained after Π-mapping of T (M )
boundary nodes, approximates the solid boundary ∂ Ω for a given precision, ε.
5. External Node Mapping: Move the boundary nodes of T (M ) to the solid surface
according to Π.
6. Relocation and Optimization: Relocate the inner nodes of T (M ) and optimize
the resulting tetrahedral mesh by applying the simultaneous untangling and smooth-
ing procedure to obtain the final tetrahedral mesh, T (Ω), that approximates the
solid.
Figure 1: The different meccano steps
Once the meccano is assembled, we have to define an admissible one-to-one map-
ping between the boundary faces of the meccano and the boundary of the solid. If the
solid is genus-zero and its boundary is given by a triangulation, we propose in Mon-
tenegro et al. [2010] an automatic method to construct a parametrization of the solid
surface triangulation to a cube boundary. For this purpose, we first divide the solid
surface triangulation into six patches with the same topological connection as the cube
faces. Then, a discrete mapping from each surface patch to the corresponding cube
face is built using the mean value parametrization proposed in Floater [2003].
At the moment, if the genus of the surface of the solid is greater than zero, the
meccano should be defined by the user. An automatic construction of the meccano
could be difficult when the topology of the solid is complex. We also remark that a
non-optimal meccano can introduce large distortion in mesh generation. To avoid this
issue an optimization of the boundary parametrization could be included Wan et al.
[2011].
In step 3, the meccano is decomposed into a coarse tetrahedral mesh T
0
(M ) by
an appropriate subdivision of its initial polyhedral pieces. Although any tetrahedraliza-
tion algorithm could be used, we propose a partition of meccano compatible with the
Kossaczký refinement algorithm Kossaczký [1994].
This mesh is locally refined in step 4 to obtain an approximation of the solid bound-
ary within a given precision. To be more precise we have to introduce some notations.
Given a tetrahedral mesh of the meccano T (M ), we denote as τ(M ) its boundary
triangulation and τ(Ω) the surface triangulation obtained after Π-mapping of τ(M )
nodes. Note that τ
o
(Ω) is a coarse approximation of ∂ Ω. In order to improve this ap-
proximation we build a refined mesh T (M ) of T
0
(M ) such that the distance between
τ(Ω) and ∂ Ω is less than a prescribed tolerance ε. The concept of distance between
surfaces can be defined and implemented in several ways. In our case it is as follows:
3
Let T = ha,b,ci be a triangle of T (M ), where a, b and c are their vertices, and let
p
k
∈ {p
i
}
N
q
i=1
be a Gauss quadrature point of T . We define the distance, d(t), between
the triangle hΠ(a),Π(b), Π(c)i ∈ τ(Ω) and ∂ Ω as the maximum of the volumes of the
tetrahedra formed by Π(a),Π(b),Π(c) and Π(p
k
). Then, the distance between τ(Ω)
and ∂ Ω, d(τ(Ω)),∂ Ω, is the maximum of all d(T), that is
d(τ(Ω),∂Ω) = max
T ∈τ(M )
d(T ) (1)
We recall that local refinement stops when d(τ(Ω), ∂ Ω) < ε. Note that this is an
approximation of the maximum missed (or overestimated) volume per face of τ(Ω).
A more accurate approach of distance based on Hausdorff envelope can be found in
Borouchaki and Frey [2005].
Then, we construct a mesh of the solid T (Ω) by mapping the boundary nodes of
T (M ) and by relocating the inner nodes at a reasonable position. After these two
steps, the resulting mesh is generally tangled. Therefore, a simultaneous untangling
and smoothing procedure Escobar et al. [2003, 2010] is applied and a valid adaptive
tetrahedral mesh of the solid is obtained. In short, this last procedure finds the new
positions of the inner nodes of T (Ω) optimizing an objective function. Such a function
is based on a certain measurement of the quality of the local submesh N(q), formed by
the set of tetrahedra connected to the free node q. In fact, we use a suitable modification
of the objective function such that it is regular over all R
3
, Escobar et al. [2003].
An example of the different steps of this method is shown in Fig. 1.
3 HARMONIE Meso-scale Weather Model
HARMONIE (HIRLAM-ALADIN Research on Meso-scale Operational NWP in Eu-
rope) is a weather prediction model design for operational use at convective scale res-
olutions. The system has been mainly developed by Meteo-France and ALADIN Con-
sortium in collaboration with ECMWF and HIRLAM Consortium.
This model uses a 3D-Var data assimilation Fischer et al. [2005] which shares most
of the code with the ECMWF and ARPEGE models. For surface variables a statistical
interpolation algorithm is used. The Non-Hydrostatic Dynamics is based on Bubnová
et al. [1995] and the physics is adapted from Meso-NH research model.
AEMET is running HARMONIE with AROME configuration at 2.5 km horizontal
resolution since October 2011. This configuration is close to the one used operationally
at Météo-France Seity et al. [2011]. Local and extreme forecasts are improved sig-
nificantly with the HARMONIE 2.5 km model compared to coarser grid models like
HIRLAM or ECMWF Navascués et al. [2013]. The model is run 4 times per day over
2 domains (Iberian Peninsula and Canary Islands) with a forecast length of 48 hours.
4 Local Wind Field Simulation
Once the tetrahedral mesh is constructed, we consider a mass-consistent model Mon-
tero et al. [1998, 2005], Ferragut et al. [2010] to compute a wind field u in the three-
4
dimensional domain ω, with a boundary Γ = Γ
a
∪ Γ
b
, that verifies the continuity equa-
tion and the impermeability condition on the terrain Γ
a
,
∇ · u = 0 in ω
n · u = 0 on Γ
a
(2)
where n is the outward-pointing normal unit vector, being Γ
b
the boundary where the
impermeability condition is not imposed.
The model formulates a least-squares problem in the domain ω to find a wind field
u = (u,v,w), such that it is adjusted as much as possible to an interpolated wind field
u
0
= (u
0
,v
0
,w
0
). The adjusting functional for a field v = (
e
u,
e
v,
e
w) is defined as
e(v) =
1
2
Z
ω
(v − u
0
)
t
p(v − u
0
)dω (3)
where p is a 3 × 3 diagonal matrix with p
1,1
= p
2,2
= 2α
2
1
and p
3,3
= 2α
2
2
. The La-
grange multiplier technique is used to minimise the functional (3), with the restrictions
(2). Considering the Lagrange multiplier λ, the Lagrangian is defined as
l (v,λ ) = e (v) +
Z
ω
λ∇ · v dω (4)
and the solution u is obtained by finding the saddle point (u,φ) of the Lagrangian (4).
This resulting wind field verifies the Euler-Lagrange equation,
u = u
0
+ p
−1
∇φ (5)
where φ is the Lagrange multiplier. As α
1
and α
2
are constant in ω, the variational
approach results in an elliptic problem in φ, by substituting (5) in (2), that is solved by
using the finite element method.
−∇ ·
p
−1
∇φ
= ∇ · u
0
in ω (6)
−n · p
−1
∇φ = n · u
0
on Γ
a
(7)
φ = 0 on Γ
b
(8)
4.1 Construction of the Initial Field
The interpolated wind field u
0
is constructed from the HARMONIE data, specifi-
cally the values of the 10m wind, u
h
n
, given at point n of the HARMONIE grid , the
geostrophic wind u
g
. Therefore, we consider a horizontal interpolation and a vertical
extrapolation of the available measurements to construct u
0
in the whole computational
domain.
4.1.1 Horizontal Interpolation
A common technique of interpolation at a given point, placed at a height z
m
over the
terrain, is formulated as a function of the inverse of the squared distance between that
5