Efficient unstructured mesh generation for marine renewable energy applications

Renewable energy is the cornerstone of preventing dangerous climate change whilst maintaining a robust energy supply. Tidal energy will arguably play a critical role in the renewable energy portfolio as it is both predictable and reliable, and can be put in place across the globe. However, installation may impact the local and regional ecology via changes in tidal dynamics, sediment transport pathways or bathymetric changes. In order to mitigate these effects, tidal energy devices need to be modelled, to predict hydrodynamic changes. Robust mesh generation is a fundamental component required for developing simulations with high accuracy. However, mesh generation for coastal domains can be an elaborate procedure. Here, we describe an approach combining mesh generators with Geographical Information Systems. We demonstrate robustness and efficiency by constructing a mesh with which to examine the potential environmental impact of a tidal turbine farm installation in the Orkney Islands. The mesh is then used with twowell-validated ocean models, to compare their flow predictions with and without a turbine array. The results demonstrate that it is possible to create an easy-to-use tool to generate high-quality meshes for combined coastal engineering, here tidal turbines, and coastal ocean


Introduction
Renewable energy generation is becoming important in stimulating economies across microetoemacro scales [1], delivering an environmentally sustainable resource [2], and has the potential to change the global energy security and interedependency landscape [3]. The potential environmental impact from renewable energy infrastructure and the cost of developing emerging technologies, leads to the requirement for numerical methods with highefidelity and accuracy in order to reliably estimate power output, optimise design and minimise environmental impact. Ocean and coastal models are increasingly used to study complex geophysical fluid dynamics and its interaction with biological and geochemical processes [4]. The potential insight from simulations has made ocean and coastal models a valuable tool in scientific research as well as engineering. Simulations are used to assess the impact of anthropogenic changes, the vulnerability of coastal and urban areas to natural hazards [5,6] and in hydrocarbon exploration and sequestration research [7,8]. In the contexts of tidal renewable energy and coastal engineering, coastal ocean models are used for estimations of power output, design, optimisation [9e11] and environmental impact assessment [12e15]. The focus of this paper is on one of the first stages in computational ocean and coastal modelling: the specification of the simulation domain and its tessellation into discrete elements, commonly referred to as mesh generation. The predictive accuracy of simulations can be significantly affected by the mesh resolution, gradation and shape of mesh elements, collectively identified as "mesh quality". Therefore, generation of highequality meshes is fundamental to ocean and coastal modelling.
The paradigm of mesh generation in coastal and ocean modelling can be broadly described as a twoestep procedure. In the first step, the domain is defined in a topologically twoedimensional space, bounded by shorelines and open boundaries [16]. A finite, twoedimensional reference surface is thus defined, often on a geodetic datum, and a mesh is generated over this reference surface. If a twoedimensional approximation, such as the deptheaveraged shallow water equations is sufficient, meshegeneration is complete. When a threeedimensional approximation is required, the second step is the projection of the surface mesh vertices at successive levels towards the ocean floor, thus creating three-dimensional elements [16]. The second step requires little user intervention, so the first step has become a synonym for mesh generation, in the ocean and coastal modelling context. Based on the regularity and structure in the mesh over the twoedimensional reference surface, a taxonomy of the mesh and relevant simulation approaches into structured and unstructured mesh methods is commonly used [17e19]. Structured meshes can be formed through a mapping of coordinate contours from an indexed space to a physical space [20,21]. The connectivity between mesh vertices is preedetermined and repeating [22]. On the contrary, in unstructured meshes, the connectivity can arbitrarily vary from vertex to vertex [23,22].
Despite the reduction of the relevant dimensions to just two, the production of quality meshes for geophysical domains can be an elaborate procedure [16,24e30]. A significant number of data sources must be combined to compose a geometrically complex domain. The geometry of geophysical domains is one of the most widely known examples of fractal geometry in nature [31,32]. Geometrical lengthescales across four orders of magnitude is typical in ocean modelling: The simulation domain size can span hundreds of kilometres, while the smallest bays can be tens of metres long. An even smaller scale may be relevant when the domain must accurately represent coastal infrastructure such as piers, pylons or embankments [26,29,30,33,34]. In addition to domain geometry, the flow typically exhibits a large range of scales, with many transient flow features, such as internal waves or jets. Therefore, a mesh must represent very complex domains with element sizes across a broad range of scales, with smaller elements in areas that require a higher fidelity, while gradation across element sizes must be smooth. Unstructured meshes are increasingly favoured in many coastal and ocean models as they tend to satisfy the above requirements with relative ease [26,36]. Alternatively, nested structured meshes embed higher resolution structured meshes into a parent coarse resolution grid [17e19,37e39]. Information from the parent grid is passed to the high-resolution grid in the form of boundary conditions, through an interpolation operation. A restriction operation, passing information from the high-resolution grid to the parent grid, is also common in ocean modelling [38e40]. Nested structured grid approaches benefit from the mesh structure as well as the desirable numerical properties common to structured mesh approaches. Nested mesh approaches thus result into efficient numerical frameworks [37,38,40,41], and are commonly used in estuarine modelling [42], as well as coastal and ocean modelling [37e41].
The aim of this paper is to present an efficient unstructured mesh generation framework suitable for coastal ocean modelling, with a particular application focus on marine renewable energy. In order to facilitate the use of multiple models, the meshing framework is not tailored to a single ocean model. Therefore, the use of multiple models in complex simulations is promoted, and comparisons between different models are facilitated. In this way, the strengths of various coastal models are highlighted, and the efficient use of their strengths is promoted. Highequality generic unstructured mesh generation packages are available, enabling routine generation of meshes. However, the interface of most mesh generators is based on Computer Aided Design and Manufacture (CADeCAM). Such interfaces have been developed for describing geometries produced by manufacturing and construction processes and do not facilitate use or manipulation of geographical data [43].
In particular, the fractal nature of ocean and coastal domains, as well as the various conventions used in geodetic coordinate reference systems are not natively expressed in CADeCAM systems [26,34,43]. Combining shoreline or bathymetry data is a typical example; one may have to combine several global, national and regional datasets as well as data from very higheresolution surveys over relatively small regions. Each dataset can also differ at least in terms of extents, resolution, coordinate reference system and vertical zeroedatums. Unlike CADeCAM, Geographical Information Systems (GIS) have been developed specifically for storing and analysing geographical information. GIS packages are robust and widely used in research as well as in operational and strategic planning contexts, where resource management, hazard mitigation and infrastructure development are a few examples. Therefore, GIS packages are ideal for supplying mesh generators with the geospatial data they require, as demonstrated by various previous projects. In Blue Kenue [44] a GIS capability was implemented as part of a complete meshing software application. Projects such as Terreno [16,24] were built on geospatial libraries, such as GDAL [45] and GMT [46], and offered a range of utilities to facilitate mesh generation. However, over recent years, with the emergence of extensible GIS packages, various extensions to and integration with these have been built, including mesh generation projects [25e28].
Another important feature of GIS systems is their capability to interact with databases, allowing concurrent data analysis and manipulation. Databases can also be used to record data origin and evolution, termed data provenance [26,47]. However, databases in GIS systems are typically not maintained in a manner pertinent to scientific research, where the primary aim of data provenance is to show reproducibility [47]. A record of data and software depended upon is usually required as evidence of reproducibility, including scientific data and software used in preparation of simulations. Therefore, the reproducibility of numerical simulations relies, at least in part, on the ability to exactly reproduce the underlying mesh, so we will show how Research Data Management (RDM) can be integrated with mesh generation in GIS. The integration of RDM and mesh generation was motivated by the increased attention on the reproducibility of scientific computation [48], perhaps as much as openesource software. Also, public research bodies are adopting policies on data and software output from publiclyefunded research to be made readily available, and provenance to be clearly identified [49e51]. In the industrial sector reproducibility, data archiving and data provenance are viewed as efficient modelling practices. Industry and governing bodies are also bound by regulatory frameworks which require public accessibility to data during the planning phase [52e55], as well as after commissioning [56,57], especially when data pertains to environmental impact of infrastructure.
Here we present the qmesh package, interfacing GIS with a meshegenerator and online data repositories. We link the abstractions offered by mesh generators and GIS packages and build tools that facilitate meshegeneration for coastal flow modelling. While the motivation for linking GIS and other geographical data manipulation tools and libraries with mesh generators is not novel [24e30, 44,58], the implementation presented here is. Unlike existing integrations of GIS and mesh generators, qmesh was principally developed as an objecteoriented software library, accessible through an Application Programming Interface (API). In qmesh GIS capability is implemented through the use of existing and robust GIS implementations as generic libraries, rather than building extensions to a particular GIS implementation, such as is described in Refs. [25e28]. qmesh is thus a library with which command line utilities and a graphical interface have also been developed. The integration of RDM functionality is another novel feature of qmesh. The broader aim of qmesh development is the creation of robust, efficient, operational and userefriendly tools for mesh generation over geophysical domains, but the particular focus here is on multiescale meshes for simulations targeted at marine renewable energy applications. Thus, the qmesh design was centred around providing the following requirements: facilitation of domain geometry and mesh element size definitions; an intuitive way of specifying boundary conditions and parameterizations; the ability to use various geodetic coordinate reference systems; promotion of the reproducibility of output and citation of data provenance and the provision of an openesource and tested package.
In Section 2 we present the qmesh package in detail, showing how the above requirements were met. An application of qmesh alongside simulation results of tidal flow around a turbine array sited in the Pentland Firth is given in Section 3, followed by discussion and conclusions in Section 4.

Overview
The user perspective of meshegeneration packages is centred around specification of two parts: domain geometry and mesh element size. Encoding domain geometry and mesh element size is a useful paradigm for describing meshes for ocean modelling [16], as it organises the necessary information in a conceptually clear way. We here follow the conventional norm in ocean modelling, discussed in Section 1, where meshes are produced in a topologically twoedimensional space [16]. Thus the shorelines and "open boundaries" represent the boundaries of the domain. The two primary data structures of GIS are used to describe linear features and field data. A vector data structure can represent points, lines and regions on a reference surface, while a raster data structure encapsulates the discrete representation of fields. The analogue to the abstraction used to drive mesh generators is clear: the domain geometry can be described with a vector data structure, while a raster can express the element size metric. Thus, the obvious route to interfacing mesh generators with GIS is to provide a translation of GIS data structures into the corresponding structures native to the mesh generator software. The data structure translation is at the heart of the qmesh package. The translation is done with little user intervention, as the user typically interacts with the GIS package and parts of the qmesh package that facilitate specification of domain geometry and mesh element size. To meet the demands around data archiving, publication and reproducibility a Research Data Management (RDM) tool is included in the qmesh package. The RDM tool facilitates the process of publishing all resources, including output such as the Gmsh mesh itself, to online citable repository services.
The GIS package chosen in this study is the QGIS package [59], the mesh generator is Gmsh [60] and the PyRDM software library [61] was used to integrate research data management [47]. The main reasons for choosing QGIS, Gmsh and PyRDM, are robustness, extensibility and permissive licences. Specifically, Gmsh is a robust mesh generator featuring a CADeCAM interface, and has been used for generating meshes in various scientific and engineering domains, including geophysical domains [60,62]. QGIS is a widely used GIS platform, with an active community of users and developers, and has been used as a user interface to mesh generation in past efforts [25,26,28]. The functionality of QGIS is available to the user as a standard GIS system with a rich graphical interface and as an objecteoriented Python module. Therefore, QGIS is a solid framework on which to develop complex applications that require GIS methods. Also, QGIS provides a framework for using such applications as extensions, via the QGIS graphical user interface. QGIS, Gmsh and PyRDM are released under the GNU General Public Licence, making possible the use of qmesh in an academic or industrial context, free of charge. Fig. 1 presents an overview of the architecture of qmesh and conveys the usual workeflow. As shown, qmesh is composed of four modules, named vector, raster, mesh and publish. The purpose of the modules vector and raster is to facilitate the definition of the domain geometry and meshesize metric and to interface qmesh to QGIS. The translation between GIS and meshegenerator data structures is performed by the mesh module, thus interfacing qmesh to Gmsh. The RDM functionality is implemented by the publish module, interfacing qmesh to online repositories and enabling identification and publication of data. A more detailed description of each module follows.

qmesh design
The vector module is used to construct a complete definition of the domain geometry in terms of domain boundaries (lines) and domain surfaces (areas). Surfaces are defined in terms of lines, so the definition of surfaces is automated by methods in the vector module such that only the boundary lines need to be supplied. Other methods allow for essential geometric operations such as checking for erroneous geometries (i.e. intersecting shorelines) or the removal of small islands and lakes, based on a threshold surface area specified by the user. The necessity of shoreline processing in ocean modelling is discussed in Ref. [24], where the Terreno project used GMT [46] to affect shoreline processing. In qmesh however, geometry processing is done primarily through the QGIS software library, also allowing use of extensive functionality built-into the GIS platform. In addition to geometry definition, methods for identifying separate parts of the domain geometry are necessary.
For example, open boundaries are associated with different boundary conditions to shorelines. The qmesh user can assign numerical identifiers to separate lines and apply different boundary conditions to separate boundaries. Numerical IDs can also be assigned to surfaces, allowing the identification of areas where different numerical treatments or parameterizations must be applied, as shown in Section 3. The QGIS library is used to store and retrieve the digital IDs as standardised feature attributes. The output of the vector module uses the ESRI shapefile [63] vector dataestructure, which also supports storage of the ID feature attributes. This way the module output as well as IDs can be visualised, assigned and edited with any GIS platform.
The aim of the raster module is to facilitate construction of raster fields that describe the desired element edge length distribution over the domain. For example, the element size might be chosen to be smaller in areas of shallow water, steep bathymetry and areas of significant variation in bathymetry slope. Therefore an optimal element size distribution is typically expressed as a function of bathymetry, its gradient and Hessian and the distance to boundaries [62,64,65]. The raster module facilitates application of various mathematical operators to be applied to raster data such as derivatives, methods for combining raster fields such as pointwise minimum and maximum operators, but also methods for calculating the distance function raster from any given vector feature. The latter is useful when specifying a mesh size gradation towards specific features in the domain: for example, the element size gradually becoming smaller as a coastline or a tidal turbine is approached. A generic method has been implemented, aimed towards the construction of element size raster fields based on the distance from a given vector feature (lines, polygons or points). This kind of operation is expressed by the arrows between the raster and vector modules inside qmesh, in Fig. 1. As with the vector module, the output of the raster module uses GIS raster data structures enabling visualisation and editing of the output via the GIS system.
The mesh module is used to translate the domain and mesh element size definitions into Gmsh data structures and, as suggested in Fig. 1 can be used to convert the mesh into a vector dataestructure. Such functionality enables mesh visualisation using QGIS, and in particular to overelay the mesh on other data. Various qualities of the mesh can thus be assessed and the workeflow can be restarted towards improving the mesh. The meshing module also allows the user to specify the coordinate reference system of the output mesh, which need not be the same as that of the domain geometry and meshemetric raster. Coordinates are reprojected to the target coordinate reference system before the data is passed to the mesh generator. The reprojection procedure uses the QGIS library; this way meshes can be obtained in all cartographic projections that QGIS supports and identifies via an EPSG code. The output mesh is twoedimensional and the EPSG specification describes the dimensions, including their units. As a particular case, the output mesh can be constructed in a threeedimensional space, where the mesh vertices lie on a sphere, using specific Gmsh functionality described in Refs. [58,62]. The vertex coordinates are specified in terms of a Cartesian reference system whose origin lies at the sphere centre, the zeaxis is the axis of rotation and the xeaxis intersects the surface of the sphere at 0 + longitude and 0 + latitude. Meshes thus constructed can be used to perform global simulations or simulations over large areas [6e8, 16,29,30,58,62].
The aim of the publish module is to facilitate provenance description and reproducibility of qmesh output. Broadly, the specific version of qmesh used to produce the mesh and all of the input data sources are stored in an online repository. In general, data provenance may seem intractable since data and software are often stored in a nonepersistent way and are not easily accessible. However, given the increasing importance of data provenance, online data repositories with efficient storage and access controls such as Zenodo and figshare, are becoming popular means of archiving and dissemination. Also, such services incorporate metaedata as means of describing hosted data and minting a unique Digital Object Identifier (DOI). The DOI is a standardised [66] citable identifier and is aimed to be assigned to digital objects, stored in a persistent way in open repositories. Therefore, DOI is a widelyeadopted identifier for digitally stored data, be that a scientific publication, the output of scientific computations or records from experiments and observations. Given the wide range of data sources that can be combined during mesh generation for realistic geophysical domains, the task of manually maintaining the provenance information of all the relevant data files can be timeeconsuming and erroreprone. As shown in Fig. 1 the publish module interfaces with the qmesh development repository, via PyRDM [47,61], to identify the exact version of qmesh used. A query is then made with the repository hosting service to establish if this version of qmesh has already been uploaded and assigned a DOI. A similar query is performed for each input data source. Each unpublished item is then uploaded and a new DOI is minted and assigned to the entire dataset. The dataset also includes citations, in the form of metaedata, of the DOI markers of already published items. The various DOI markers can be thought of as nodes of a tree, and the citations are the tree connections, a similar concept to scientific publications. Also, the output can be archived in a private repository, without a DOI, to facilitate archival of commercially sensitive information.

User interface
qmesh can be used in a graphical as well as a programmatic environment. The user interface consists of a Graphical User Interface (GUI), a Linux Terminal Command Line Interface (CLI) and a Pythonebased Application Programming Interface (API). The graphical user interface is a QGIS extension. This way qmesh can be used via the QGIS application and in combination with other QGIS functionality. The GUI has been designed to allow access to the qmesh package with ease and little knowledge of the qmesh design. The command line interface consists of a set of utilities each with welledefined input and output, making each utility a separate program. However, their purpose is to be used as diagnostic tools or, in sequences that automate operations which do not require a graphical interface. The Python-based API is an integral part of the qmesh implementation. Like the CLI it can be used to build scripts that automate series of operations, but its primary purpose is to allow the use of qmesh as a software library.

Testing and availability
A testing framework is used to ensure that code development does not affect the output. The testing framework executes various tests, classified into unitetests and regressionetests. The former compare the output of relatively small and atomic parts of the code, such as object methods, against the expected result. Regression tests are scripts that aim to emulate the typical use of qmesh, by producing the geometry and element size metric followed by mesh production of realistic domains. All tests are run in a continuous integration fashion [67,68] where a run of all tests is initiated after each new code revision is submitted to the public code repository. The qmesh source code is publicly available under the GNU General Public License (v 3). The documentation is released under the GNU Free Documentation License (v 1.3). The source code and documentation are available through [69].

Mesh generation for simulations of tidal flow with a turbine array in the Inner Sound
In this section, we demonstrate how qmesh can be used for higheresolution simulations in the context of renewable energy generation with tidal turbine farms. A mesh is constructed in a domain encompassing the Orkney and Shetland islands, as shown in Fig. 2(a and b), and this mesh is used with two wellevalidated ocean models.
The seas north of the Scottish mainland and around the Orkney islands, shown in Fig. 2(c), have been identified as having a very high potential for renewable energy generation [72,73]. A number of areas have been selected for commercial power production [70,74,75], extracting energy from waves or tidal streams. The total power output from proposed tidal stream turbine farms in the Orkney Islands is estimated at 1 GW [70,74,75], with 800 MW estimated in the Pentland Firth alone, a channel separating the Orkney islands and mainland Scotland. In this study, we focus on a tidal turbine farm site located in the Inner Sound, a narrow passage within the Pentland Firth between Gill's Bay and the Isle of Stroma shown in Fig. 2(d). The flow in the Inner Sound is characterised by strong semiediurnal tidal currents. Surveys [76] suggest a turbulent, eddying velocity field and show that local current speeds regularly exceed 3 ms À1 at various stages of the floodeebb tidal cycle. Velocity measurements at a fixed location (see Fig. 5), show that deptheaveraged current speed exceeds 4 ms À1 during spring tides. More detailed ADCP measurements of vertical velocity profiles [77] show that velocities will regularly exceed 4.5 ms À1 and may even exceed 5 ms À1 during an equinoctial spring tide. With such strong currents, the Inner Sound has been a very attractive location for pilot studies on tidal stream power generation, and the site capacity is estimated to be up to 400 MW [70,74,75]. MeyGen Ltd has successfully submitted an application towards consent for commercial power generation in the Inner Sound [78,79]. The installation of up to 86 turbines over an area at the central part of the lease area [71], producing an estimated output up to 86 MW, has received consent as a first phase of the Inner Sound plot development [80]. Fig. 2(a) and (b) show the area selected as the simulation domain in this example. As shown, the domain boundaries are formed using various data sources. Since the region of interest in this study is the Inner Sound, a relatively accurate representation of the shorelines around the Pentland Firth is required. Thus, the "full" resolution set from the GSHHG dataset [81] is used to represent the Orkney islands, as well as a part of the mainland coast. As discussed later, regions further away are simulated in order to capture the correct hydrodynamic response around the area of interest. A less detailed shoreline representation in these regions should not affect accuracy within the region of interest and will require larger elements to produce highequality multiescale grids. Therefore, the "intermediate" resolution set from the GSHHG dataset and the 0 m elevation contour extracted from the GEBCO 2014 30 arceminute resolution data [82] are used to represent mainland shorelines and the shorelines of the Shetland Islands, as shown in Fig. 2(b).

Mesh generation
To avoid simulating the flow past the continental shelf break and limit the simulation domain over the shallow area of the northewest European continental shelf, part of the 300 m depth contour from GEBCO 2014 [82] is used to define the open boundary. As seen in Fig. 2 this line closely follows the "edge" of the continental shelf. The domain is then closed by constructing smooth curved lines linking the coastline to the 300 m depth contour. The curved, yellow lines in Fig. 2(b) are defined as a linear combination of loxodromic lines; i.e. lines of constant bearing. Loxodromic lines allow specification of boundaries with accuracy, and their construction is efficient as well as simple; qmesh provides utilities for constructing such features: the user specifies the coordinates of the starting point and the angle between the meridian at the starting point and the loxodrome direction. For practical reasons qmesh utilities construct a segment of a loxodrome, from the starting point up to a given longitude or latitude. Other qmesh utilities combine two loxodromes by assigning a linearly varying weight over the The western boundary is a combination of the following loxodromic lines: Starting at 5:0 + West, 58:6 + North (a point near Loch of Strathbeg, Aberdeenshire, Scotland), at a bearing of À10 + up to longitude 60:2 + North. Starting at 3:43 + West, 60:46 + North at a bearing of À135:0 + up to latitude 5 + West.
Simulations of flow interaction with energy extraction devices for electrical power generation highlights the challenges to mesh generation in coastal modelling. The large size of the domain, relative to the turbine scale, and the long simulated periods necessary to collect results pose untenable resource requirements when an accurate geometrical description of the tidal turbines is needed. Parameterisations are thus typically used, treating the energy extraction from the flow as a subegridescale process. Individual devices can be represented as distinct areas where such a parameterisation is applied [83], or an entire tidal turbine array can be treated as a single region where a parameterisation is used to represent the energy extraction [84]. In this study, a parameterisation approximating the effect of each turbine on the flow is applied over cells contained within a 10 m radius at each turbine location. The parameterisation is discussed in Refs. [14,83] and a brief overview is given in Section 3.2. The circular markers in Fig. 2(d) show the locations of 180 tidal turbines. The boundaries of the Inner Sound tidal plot, obtained from Refs. [70,71], are also illustrated in Fig. 2(d). The coordinates of the turbines are obtained from a random number generator, subject to the constraints that all points are placed at least 15 m away from the boundaries of the plot and 100 m away from all other turbines. A 10 m radius circle was constructed around each turbine centre point, to facilitate identification of elements where turbine parameterisation should be applied. The boundaries of the Inner Sound tidal plot and the turbine circles are treated by qmesh as internal boundaries, such that mesh vertices should be placed along the boundaries of each area and the triangulation should produce edges along the internal boundaries.
Given the exclusion of areas with notable bathymetry changes from the domain in Fig. 2, the optimal mesh size, in this case, can be described in terms of proximity to the shorelines and the tidal array [65] only. Therefore, the following proximity based gradation specifications were used: Element edge length of 150 m at the GHSSG full resolution coastlines, shown with black lines in Fig. 2(b). The element edge length is increased linearly to 15 km across a distance of 1 (on an orthodrome) from the shoreline. Element edge length of 1500 m at all other shorelines. The element edge length is increased linearly to 15 km across a distance of half a degree (on an orthodrome) from the shoreline. Element edge length of 5 m inside the Inner Sound tidal array site. Outside the site, the element edge length is maintained at 5 m up to a distance of 0.01 (on an orthodrome) from the site, after that increasing linearly to 15 km across a distance of a degree. Fig. 3 shows the element edge length distribution, obtained by combining the above specifications through a pointwise minimum operator. The three panels in Fig. 3 highlight the gradations described above. The element size metric fields are also conveyed by the images of the mesh, in Fig. 4.
The Domain shown in Fig. 2 extends from 6 + West to 0:8 + East and the region of interest is centred at approximately 3:1 + West. Therefore, the UTM30 zone was selected for this study, as UTM30 was developed for regions with longitude ranging between 6 + West to 0 + . Fig. 4 shows the mesh obtained, where panel (a) shows the whole domain and panels (b) and (c) show details around the Inner Sound tidal site. The elements are coloured according to the numerical tag assigned to each area. As discussed in Section 2.2 numerical tags can be used to identify separate areas, and flexibly specify differences in the numerical treatment of the different areas in the domain. Fig. 4(d) shows a detail of the East side of the Inner Sound site, where the blue elements are located inside the circular areas utilised here to parameterise turbines, as discussed in Section 3.2. All other elements are treated identically; the different tags were here introduced for demonstration alone.
The resolution within the Inner Sound used in the present study is amongst the highest appearing in the literature to date, while the domain is relatively large. The addition of the seas around the Shetland Islands, north of the Orkney Islands, does not markedly impact on the simulation accuracy but emphasises the multiescale nature of geophysical domains and tests qmesh on assembling and meshing such domains. Similar studies have included large areas East and West of the Inner Sound, as the fast currents are closely linked to large potential differences forming around the Pentland Firth. Adcock and coeworkers discuss the issue of appropriate open boundary placement and boundary conditions in Ref. [85], where the mesh resolution increased from 20 km to 150 m in the Pentland Firth. Their domain is larger east and west of the domain in this study, but the boundary conditions are very similar to the current study. Easton and coeworkers [86e88] use a domain smaller than the one in the present study centred around the Pentland Firth, to examine the tidal dynamics, erosion and deposition patterns and estimate extractable power. The domain extends approximately 100 km East, West and North of the Inner Sound. Their grid resolution increases from 2 km at the offshore boundaries to 100 m in the Pentland Firth. In both cases, agreement with observations is very good and informed the present choice of open boundary placement and overall mesh resolution. The effect of turbine array installations is also examined in Ref. [14] where the grid resolution varies from 18 m in the Inner Sound to 20 km at the open boundaries. A detailed survey of the Inner Sound seabed topography and composition is given in Ref. [89] where the effect of the currents on the topography are also examined via highly resolved simulations, with approximately 40 m resolution in the Inner Sound.
The domain boundaries, Inner Sound plot boundaries and turbine parameterisation areas (shown in Fig. 2) are available in Ref. [90], The mesh (shown in Fig. 4) is available in Ref. [91]. The qmesh source code was uploaded separately, and is available in Ref. [92].

Model configuration
In order to show that the meshes are immediately usable by many ocean modelling packages, we present some results from simulations using the meshes shown in Fig. 4. The Fluidity and Telemac ocean modelling and fluid dynamics packages were used to solve the shallow water equations in the "deptheaveraged" formulation: where u i is the deptheaveraged velocity, h is the total water depth and n e is the effective viscosity. Z s is the free surface elevation measured from the reference level, in this case the UTM30 vertical datum. Denoting the bathymetry, shown in Fig. 2, by Z f gives h ¼ Z s À Z f . C i are the components of the Coriolis term: (2) is the seabed friction term, and S t is a momentum sink term used to parameterise the effect of the tidal turbine on the larger scales of the flow. The friction term is written as a function of depth, the friction coefficient C f and the velocity magnitude q : The value of the friction coefficient is in general a variable, depending on floor topography, roughness and fauna. Here the constant value 0.0025 was used throughout the domain, as this is a commonly used value for coastal waters and shelf seas [93,94]. However, recent studies have shown a value 0.005 to give better agreement with observations for the Pentland Firth [85,87]. Nonetheless, the meshing procedure described in Section 3 will not change with an adjustment of C f . In Telemac, the specification of the seabed friction is made in terms of a Chezy formulation, with a coefficient C, The tidal turbine parameterisation is here defined in terms of a drag force, following [83]: where A p is the horizontal area over which the turbine parameterisation is applied, r denotes density and F is the drag force of a turbine [83]: where C D is the turbine thrust coefficient and A is the representative area of the tidal turbine; for example, the area swept by the turbine fan. We have chosen the area to apply each turbine parameterisation A p to equal A in this example, so the momentum sink term can be written as: A constant thrust coefficient is used in simulations using Telemac and Fluidity, C D ¼ 0.6. It should be noted, that a more realistic representation is one where the power output of a turbine is rated and thus the thrust coefficient varies with the local velocity magnitude as discussed in Refs. [14,83]. The discretization of Equations (1) and (2) in Fluidity employs a mixed finite element pair combining a linear discontinuous function space for velocity with a quadratic continuous function space for depth. A twoetimeelevel temporal integration is used, detailed in Ref. [14]. The discretization is described in detail in Refs. [14,95e98]. The Telemac discretization also employs a mixed finite element formulation, combining a quadratic quasiebubble [99] approximation for velocity with a linear continuous Galerkin approximation for depth. Temporal integration is effected using a method of characteristics to advance the advection terms, followed by a twoelevel integration of the remaining terms. The Telemac discretization is further detailed in Refs. [34,99].
A fixed effective viscosity, n e ¼ 10 m 2 s À1 , was used in both Fluidity and Telemac simulations. No turbulence model was used, for reasons of consistency across the simulations; the differences in results due to differences in turbulence model parameterisation and implementations are thus removed.
Identical bathymetry and boundary conditions were used in both Fluidity and Telemac simulations. The bathymetry Z f was obtained from the 1 arcesecond Digimap Marine, HydroSpatial One [100] gridded bathymetry data, shown in Fig. 2(b,c,d). The bathymetry was modified in very shallow areas to ensure the sea floor was submerged throughout the simulated period; areas shallower that 10m were made 10m deep. Therefore, no wettingeandedrying algorithms were necessary to capture the exposure and submergence of tidal flats. As with turbulence parameterizations, this allows consistency to be maintained between the two simulations, by removing the differences in the treatment of wetting and drying. However, if the accurate prediction of tidal inundation is required for high predictive accuracy, the mesh generation procedure described in Section 3 would not change appreciably when tidal flats are accurately represented in the simulation.
The shoreline boundaries were treated as impermeable walls where the velocity boundary conditions did not permit flow normal to the boundary, without constraining the velocity components tangential to the boundary. At the open boundary, the tidal freeesurface elevation from the OTPS tidal dataset over the North-eWestern European shelf [101] was used to calculate the free surface height as a Dirichlet boundary condition, without constraining the velocity. As pointed out in Refs. [85,102]

Simulation results
While the focus of this study is not to perform an inedepth comparison of the two models, it was noted that the differences between Fluidity and Telemac in discretization required the adoption of different timeestep sizes. Fluidity was found to perform stably with Dt ¼ 20 s, but Telemac simulations were limited to Dt ¼ 0.2 s, as a timeestep size similar to that used in Fluidity led to unphysical oscillations in the primitive variables. Fluidity simulations were run on 480 cores and Telemac simulations were run on 300 cores of the Helen higheperformance computing resource at Imperial College London. Denoting the simulated time by t s , the duration of each run (walletime) by t w and the number of cores by N c , the Fluidity simulations run as ts twNc z0:005 while the Telemac runs were overall quicker at 0.02 of the same metric. Given the discontinuous velocity space leads to much larger linear systems and the nonelinear Picard iterations in Fluidity [14,97], this result is expected.
Five days were simulated, starting from 2 July 2011. Fluidity and Telemac were used to simulate the baseline flow, without the tidal array, as well as the flow with parameterizations of the tidal turbines. Results from only the last three days will be shown, as the two days at the beginning of the simulations are used to allow simulations to spineup and settle from transients due to nonerealistic initial conditions. The simulation start date and spin-up period were selected so that the final three days of the simulated period coincide with the observation of spring tides when the tidal currents in the Inner Sound reach peak velocities.
Data collected with an Acoustic Doppler Current Profiler (ADCP) south of the Isle of Stroma (see Fig. 2(d)) are here used to validate both models. Fig. 5 shows the temporal variation in free surface elevation and velocity components from Fluidity and Telemac simulations of the baseline flow at the ADCP location, as well as deptheaveraged velocity data from the ADCP. The grey dashed lines identify various events during a single semiediurnal tidal cycle over the 4th and 5th of July: the line labelled with (i) indicates that maximum velocity magnitude was obtained on July 4 21:50, while the free surface elevation was approaching its highest point. A high tide was seen on 5 July 2011 00:00, marked by the line labelled (ii), followed by slack flow at 01:30 indicated by line (iii). The lowest free surface elevation was nearly concurrent with ebb tide maximum velocity and was seen at 4:40, marked by line (iv). The next slack water is indicated by line (v) at 08:00. Line (vi) denotes the maximum flood velocity, concurrent with zero free surface elevation at 10:10.
The asymmetry in the free surface elevation during successive floodeebb cycles is clearly evident in the free surface elevation plot of Fig. 5. The agreement between the two models is good during ebb tides, but differences are apparent during the flood tides. Both models predict an undulation in the free surface during the flood tides, and the undulation is larger during the greater of the two semiediurnal flood tides. Fluidity results show larger undulations than Telemac. A good agreement between Fluidity and Telemac in velocity components is seen in Fig. 5. However, comparison with ADCP data demonstrates that both models over-predict peak velocity magnitudes during the flood and ebb tides. During peak flood currents a strong fluctuation, similar to the undulation in free surface elevation, is predicted by both models. ADCP data and simulation results show u 1 fluctuations during flood tides. The ADCP data also suggests that the models overepredict the peak u 1 component magnitude during peak ebb currents.
The u 2 plots in Fig. 5 show that both models are in close agreement but overepredict the u 2 component magnitude during both flood and ebb tides. During ebb tides, the plot shows that both models predict a flow slightly towards the North during the peak ebb currents at the ADCP location, while the ADCP recorded an overall flow to the West, with strong fluctuations. Overall, Fluidity gives the largest overeprediction of u 2 magnitude during flood tides, while Telemac gives a larger overeprediction during ebb tides. Fig. 6 shows contours of the free surface elevation at six instants during the semiediurnal cycle. The shown instants were identified in Fig. 5. There is a good qualitative agreement between the two models, and both models predict the expected overall behaviour in the free surface elevation: the flow is driven by the North Atlantic tides and the resonant response of the North Sea. Also, results from both models suggest that at instants of maximum flow velocity at the ADCP location there is a large difference in free surface elevation around the Pentland Firth. In particular, west of the Pentland Firth at instants (i) and (iv) offshore the Northern Highland shoreline, both models show substantial departures of the free surface elevation from the datum, while a much smaller departure from the horizontal datum is seen east of the Pentland Firth. The The velocity magnitude over the Pentland Firth area is shown in Fig. 7, for the instants (i) to (v), for both Fluidity and Telemac. The panels also show the velocity vectors, in areas where the second invariant of the deformation tensor indicates the existence of vortices [105,106]. An animated sequence of the same plots, over a period of 24 h can be found in Ref. [107]. The largest velocities are seen in the Inner Sound, between Stroma and Swona and north of Swona during both flood and ebb tides in panels (i) and (iv). During the ebb tide highevelocity magnitude is also seen near the Head of Ducansby, in panels (iv), consistent with the observation of the "Ducansby Race" during ebb tides. Fluidity results also show high velocities off the coast of South Ronaldsay while Telemac results show increased velocities, but not as high as in the Fluidity results. The overall flow direction during the flood tide (panels (i) and (ii)) is west to east between Caithness and Hoy, west of Stroma. Further downstream, the flow is deflected around the Isles of Stroma and Swona and then follows a SoutheEast direction, also forming a jet in that direction. A smaller jet emanates North of the Pentland Skerries flowing eastwards. The overall change in the flow direction in the Pentland Firth is consistent with the Kelvin wave in the North Sea. During the ebb tide, the predominant flow direction reverses.
The most prominent features revealed by the velocity vectors are vortices downstream of islands and in bays. The results from both models show that island wakes are hosts to vortices, especially during peak flows. The animated sequence in Ref. [107] indicates  Ronaldsay in panels (i) and (ii), rolling off the South-West and SoutheEast sides of the island during the flood tide. The results from both models show recirculation regions forming in the two bays West of Duncansby Head during peak flood and ebb currents, in panels (i), (ii) and (iv). South of Duncansby Head both models predict a vortex forming in the bay and then shed from the headland in panels (i) and (ii). The shed vortex appears East of Skirza Head in panels (ii). A second vortex forms inside the bay as shown in panels (ii), but is not shed into the flow as shown in Ref. [107]. Overall, the results from both models are consistent with documented locations of vortices and races in the Pentland Firth [108]. Fig. 8 shows the difference in velocity magnitude between simulations which include turbine parameterizations and simulations of the baseline flow, for both models. Negative values signify a simulations. An animated sequence can be found at [107]. reduction in velocity magnitude when the turbines are installed. Results from the first five instants identified in Fig. 5 are used to compose the panels in the two top rows, and the third row shows a detail around the Inner Sound, at the time of peak flood currents.
The results indicate a velocity magnitude decrease over the turbine farm area. The change in velocity magnitude predicted by Fluidity appears more widespread than in Telemac results. In particular, shifts in the wakes of Swona and the Skerries are larger in Fluidity results than in Telemac. As shown in Ref. [107] and noted in connection to Fig. 7, Fluidity results show larger undulation in island wakes with more frequent vortex shedding when compared to Telemac results. The addition of the tidal turbine farm to the Fluidity set-up affects the location of vortices as well as flow direction and magnitude, leading to the differences shown in Fig. 7. The details shown in panels (a) and (b) show the change in velocity magnitude due to the addition of the turbine parameterizations in the area around the tidal farm. The areas with turbine parameterisation are evident in both plots as areas of velocity magnitude reduction, as well as a wake downstream of the parameterisation area itself. Also, panels (a) and (b) show that both models show a decrease of the velocity magnitude in the wake of the turbine array, with a strong reduction west of the array.

Discussion and Conclusions
The minimisation of environmental impact from activities exploiting the coastal environment is an increasingly important aspect in the design and operation of related infrastructure. The ability to explore the effects of different designs and operating parameters makes coastal and ocean modelling an important tool in the design and optimisation of coastal infrastructure. Unstructured meshes can be advantageous in the context of multiescale ocean and coastal modelling: coarse resolutions can be used to economise computational resources, while gradating to finer resolution in the regions of interest or where accuracy, stability and parameterisation constraints impose a small element size. Large domains can be efficiently tessellated with accurate representations of complex domain boundaries, such as natural coastlines and coastal infrastructure. Also, the ability of unstructured meshes to comply with arbitrary internal boundaries is very useful in studies where parameterizations must be applied locally, within given regions of the domain. The methodology presented in this paper provides an example for how mesh generators may be linked to Geographical Information Systems in order to form an efficient and robust mesh generation framework for coastal ocean modelling. Particular emphasis is given to simulations in the context of tidal renewable power generation and the modelling of its potential environmental impact, via the modification of tidal stream patterns.
Geographical Information Systems have been developed for processing, visualisation and analysis of spatial geophysical data. The combination of GIS and mesh generation methodologies is therefore widely practised in geophysical fluid dynamics [24e28, 44,58]. The novelty of qmesh lies in its design as an objecteoriented library, with additional user interfaces and utilities, integration of Research Data Management functionality, and the applicability to multiple application codes as demonstrated here. As with similar projects, qmesh interfaces GIS and mesh generators by linking corresponding fundamental abstractions between the two. The basic abstractions are also reflected in the user interfaces of qmesh, resulting in a clear and consistent structure. The multitude of user interfaces was implemented to allow qmesh to be used in various ways: as a software library, a utility in a mnemonic command oriented terminal and a tool within the rich graphical interface of a GIS system. This way qmesh is a robust and flexible mesh generation methodology that allows effort to be focused on interpreting and combining geospatial data in a mesh generation context, rather than focusing on translating data. However, the volume of data and number of data sources can rapidly increase. Data management can become problematic as duplication and lack of source attribution can impact on the reproducibility of mesh generation as well as simulation output. Thus, RDM functionality was added to qmesh, to facilitate reproducibility and data provenance attribution of computational simulations [47,61].
A mesh was generated around the Orkney and the Shetland Islands, to show the efficiency and robustness of qmesh. This region was chosen in the current study because a number of sites have been selected for commercial power generation from waves or tidal streams. The Inner Sound site in the Pentland Firth was chosen as the area of interest in the present study. Existing measurements indicate turbulent flow [76], but also highevelocity magnitudes throughout large parts of the tidal cycle leading to commercially attractive estimates of maximum and extractable power [70,74,75,85,86,104].
The distribution of mesh resolution was chosen to be typical of a study investigating the effect of the turbine array installation at the Inner Sound. A higheresolution mesh was generated at the Inner Sound site, with moderate resolution extending around the Orkney Islands. The impact of turbine arrays is assessed via a parameterisation approximating the effect of the turbines on the larger scales of the flow. The parameterisation was applied to regions representative of individual turbines, and the mesh is compliant to the region boundaries with vertices and element edges on the region boundaries. The turbine regions were the smallest geometric scale present in the domain, as their size was representative of the turbine fan. Relatively small elements were thus needed inside and around the tidal farm. This study is therefore a typical example reflecting the advantages of unstructured multiescale meshes.
We showcase the robustness of qmesh by using the same mesh with two ocean models: Fluidity and Telemac. The parameters of Fluidity and Telemac simulations were chosen to be as close as possible, but the differences present between the two models prevent identical discretization. Thus, the results highlight the effect of differences in numerical schemes and implementation details. The potential to use the mesh generator output with multiple models and the integrated RDM functionality to support dissemination of generated meshes facilitate more extensive studies towards quantification of the uncertainty due to differences in numerics and implementation as well as the calculation of more accurate ensemble results, as practised in climate modelling and forecasting [109,110].
While a detailed validation is necessary to support conclusions on the hydrodynamic aspects, a broad comparison of the simulation results with other literature was presented, to further showcase the robustness of qmesh. The simulation results, in broad agreement with other literature [85e88, 104], show that the strong currents are closely linked to the large variations in free surface elevation around the Firth. The free surface elevation and velocity of the baseline flow from the two models were in good overall agreement. Also, results from both models show good agreement with ADCP measurements, but both models overepredict the northerly velocity component in the Inner Sound during peak ebb currents. The substantial fluctuations in ADCP results during peak ebb currents suggest that better representations of subegridescale turbulence and sea floor drag would improve accuracy.
An examination of velocity maps in the Pentland Firth shows that both models predict various transient flow features, in the form of vortices forming downstream of islands and inside bays. Fluidity results suggest that the island wakes include shed vortices, contrary to Telemac results. For that reason Fluidity results show that the introduction of the Inner Sound tidal farm has a widerespread effect than Telemac results; a small change in overall flow direction and speed caused by the tidal farm will alter the path of shed vortices, and this is reflected in the speed difference maps.
More detailed investigations, outside the realm of the present study, are aimed towards quantifying the effect of tidal turbine farms on the flow, the sedimentology and in turn the benthic environment [14,111]. Studies have shown a close link between vortices and the morphology as well as sedimentology of the sea floor [89]. A thorough comparison of coastal models should be based on comparisons with additional observational data from ADCP and tide-gauges [14,112], as well as simulation of idealised configurations [38,113,114,115]. The comparison of results with accurate laboratory or observation measurements and highequality simulation data would help identify best practices for each model and achieve high predictive accuracy. In that context, a robust mesh generation framework promotes model validation and comparison as well as the simulation of complex, realistic problems such as the quantification of the environmental impact of coastal infrastructure, the latter being more challenging from the mesh generation perspective. Thus, this study shows how robust mesh generation can be used as part of computational modelling framework facilitating the utilisation of a variety of ocean models.