High-resolution numerical modeling of mesoscale island wakes and sensitivity to static topographic relief data

Recent decades have witnessed a drastic increase in the fidelity of numerical weather prediction (NWP) modeling. Currently, both research-grade and operational NWP models regularly perform simulations with horizontal grid spacings as fine as 1 km. This migration towards higher resolution potentially improves NWP model solutions by increasing the resolvability of mesoscale processes and reducing dependency on empirical physics parameterizations. However, at the same time, the accuracy of high-resolution simulations, particularly in the atmospheric boundary layer (ABL), is also sensitive to orographic forcing which can have significant variability on the same spatial scale as, or smaller than, NWP model grids. Despite this sensitivity, many highresolution atmospheric simulations do not consider uncertainty with respect to selection of static terrain height data set. In this paper, we use the Weather Research and Forecasting (WRF) model to simulate realistic cases of lower tropospheric flow over and downstream of mountainous islands using the default global 30 s United States Geographic Survey terrain height data set (GTOPO30), the Shuttle Radar Topography Mission (SRTM), and the Global Multi-resolution Terrain Elevation Data set (GMTED2010) terrain height data sets. While the differences between the SRTM-based and GMTED2010-based simulations are extremely small, the GTOPO30-based simulations differ significantly. Our results demonstrate cases where the differences between the source terrain data sets are significant enough to produce entirely different orographic wake mechanics, such as vortex shedding vs. no vortex shedding. These results are also compared to MODIS visible satellite imagery and ASCAT near-surface wind retrievals. Collectively, these results highlight the importance of utilizing accurate static orographic boundary conditions when running high-resolution mesoscale models.


Introduction
Massively parallel computing platforms now enable regional-scale numerical weather prediction (NWP) models1 to be easily integrated with fine-scale grid spacings, down to approximately 1 km horizontally.A valuable benefit of such high-resolution models is their capability to simulate orographically induced flow phenomena.Examples of such phenomena include gap winds (Mass et al., 2014), lee rotors (Ágústsson and Ólafsson, 2014), and wake vortices (Li et al., 2008).The accuracy of model simulations of orographic flows has been verified against a suite of observational data including, but not limited to, ground-based instruments e.g., lidar (Lesouëf et al., 2013), mesonets (Bieringer et al., 2013); satellite-based remote-sensing instruments e.g., synthetic aperture radar (SAR) (Miglietta et al., 2013); and airborne measurement platforms, e.g., aircraft (Gioli et al., 2014), radiosonde (Nunalee and Basu, 2014).Despite the increased resolvability, and overall fidelity, offered by finer-resolution models as it pertains to orographic flows, mesoscale NWP models are still constrained by multiple factors such as necessary physics parameterizations (Doyle et al., 2013;Draxl et al., 2014).The treatment of subgrid-scale (i.e., sub-mesoscale) processes such as turbulence, radiative transfer, moisture phase change, etc. collectively contributes to the uncertainty of model solutions (see Coiffier, 2011).At the same time, it has also been demonstrated that model uncertainty can be increased through the prescription of inaccurate, or unrepresentative, time-dependent atmospheric boundary conditions (Kumar et al., 2011;Pielke, 2013).In the past decade, advanced data assimilation techniques, coupled with improved remote-sensing capabilities, have been shown to reduce simulation uncertainty (Ancell et al., 2011;Bieringer et al., 2013) and increase forecast skill (Pu et al., 2013).While great efforts have been expended to identify sources of NWP error with respect to model configuration (i.e., physics parameterizations) and dynamic (meteorological) boundary conditions, often overlooked is the sensitivity of model solutions to static boundary conditions, namely topographic relief.
Presently, several global terrain height data sets exist which can be used by regional-scale NWP models.One of the most used surface relief data sets, named GTOPO30, was developed by the United States Geographic Survey and comprised through a synthesis of numerous international digital elevation models.GTOPO30 contains maximum spatial resolution of 30 arc seconds and is the default data set for many community models such as the Weather Research and Forecasting (WRF) model.Aside from GTOPO30 data, other satellite-derived global terrain height data sets are also available such as the Shuttle Radar Topography Mission (SRTM) (Farr et al., 2007), and the Global Multi-resolution Terrain Elevation Data (GMTED2010) (GMTED, 2011).These data sets offer higher spatial resolutions globally of 3 and 7.5 arc seconds, respectively.The construction of surface terrain height grids in NWP models from source data sets (e.g., GTOPO30, SRTM, and GMTED2010) typically involves subgrid scale averaging of the source data, grid-scale spatial interpolation during data ingestion, and/or preprocessing smoothing effects (e.g., see the WRF model Preprocessing System Documentation; NCAR, 2014).Although in many circumstances these activities are necessary, they can effectively result in under-resolved topographic relief.Under-resolved terrain height implies that the NWP model generated terrain height does not fully capture the relevant features of the natural topography described by the source data (Jiménez and Dudhia, 2012) and can result in terrain height discrepancies on the order of tens to hundreds of meters (Jiménez and Dudhia, 2013).Such discrepancies have been shown to result in significant error in simulated lowlevel wind fields (Rife and Davis, 2005;Jiménez et al., 2010;Santos-Alamillos et al., 2013).Aside from under-resolved terrain height in modeled grids, which is essentially an oversimplification of the source terrain height data, we show in this paper that uncertainty in source terrain height data sets themselves can be significant enough to result in fundamental differences in simulated orographic flow mechanics.This finding illustrates that the sensitivity of NWP models can be more complex than first-order biases recently documented by Teixeira et al. (2014).
In this paper, we simulate two realistic cases of atmospheric flow past mountainous islands; for each case, we run WRF model simulations using GTOPO30, SRTM, and GMTED2010 source terrain height data while keeping all other model configurations identical.From the results, we comment on the fundamental differences in simulated atmospheric wake patterns associated with the three terrain height fields.At the same time we compare the simulated flow features to those observed in visible satellite imagery.Our results demonstrate that selection of terrain height source data can, in some cases, be critical to successfully capturing the fundamental mechanics of mesoscale orographic wakes.

Case studies and modeling details
Two historical atmospheric events were considered in this paper, both corresponding to cases of flow past mountainous islands.Since the islands were far from any upstream surface heterogeneity, only the local terrain features associated with the islands perturbed the local winds and consequent cloud structures.For these events, the wind wake characteristics associated with each island were indicated by distinct cloud structures captured in visible satellite imagery from the Moderate Resolution Imaging Spectroradiometer (MODIS) instrument.The modeled wind wake patterns of the events were compared to one another and the differences were documented in the context of the inferred wake patterns shown in satellite imagery.
The first, and primary, case study involved the Spanish island of Gran Canaria (GC) off the west coast of northern Africa on 30 April 2007.MODIS visible satellite imagery from this day (Fig. 1, left panel) revealed a coherent pattern of dipole vortices (i.e., von Kármán vortices) being shed downstream of GC.GC has a diameter of approximately 50 km at sea level and has a peak elevation of 1948 m m.s.l.GC's SRTM-based topography is shown in Fig. 1 (right panel) for reference.
The second case study presented here involves flow past several islands which collectively comprise the Lesser Antilles (LA) in the eastern Caribbean.On 1 August 2013, MODIS visible satellite imagery of the Lesser Antilles region illustrated distinct wakes behind all of the major islands of the LA (Fig. 2 left panel).Contrary to the GC case which had a coherent vortex shedding wake regime, the LA case had weak wind wakes where the rotation behind each island was not strong enough to counter the background wind flow.Furthermore, the wakes were correlated with a reduction in cumulus cloudiness and darker sea surface color, a phenomenon investigated by Smith et al. (1997).The windward islands of the LA are generally lower than GC but are, nonetheless, predominately mountainous with peak elevations near 1 km for each island (see Table 1).The corresponding near-surface wind retrievals from the ASCAT (Advanced Scatterometer) instrument (Vogelzang et al., 2011) on the MetOp-B (Meteorological Operational) satellite are shown in Fig. 2 (right panel).The speed of the predominantly east-southeasterly winds decreased from an upstream value of 7.5-9.0 to 4.5-7.0m s −1 in the lee of the islands.The maximum speed reduction and the width/length of the wind wake were correlated with island peak elevation.The slower winds in the wake lead to a smoother sea surface and, thus, increased specular reflection and decreased diffuse reflection, while the opposite is true for the rougher sea surface areas experiencing faster winds outside the wake.

C. G. Nunalee et al.: High-resolution numerical modeling of mesoscale island wakes
Depending on the sun-satellite geometry, this difference between the relative strengths of specular vs. diffuse reflection can result in both dark and bright island wakes.When the island is farther away from the solar specular point, as in the LA case, the wake is darker than the surrounding rougher sea surface.Contrarily, when the island is close to the specular point, as in the GC case, the wake appears brighter than the surrounding ocean surface.
The numerical simulations performed in this study used the Weather Research and Forecasting (WRF) model version 3.6.1 which was initialized by ERA-Interim reanalysis data (physics configurations are shown in Table 2).The simulations used a nested four-domain configuration centered on the islands of interest.Of note, a horizontal grid spacing of 1 km was chosen in the innermost domain (d04) while the parent domains (d03-d01) used grid spacings of 3, 9, and 27 km, respectively.Additionally, in d04 the control simulations used GTOPO30 terrain height while the experimental simulations used terrain height data remapped from SRTM and GMTED2010 to 30 arcsec (i.e., SRTM30 and GMTED30, respectively).SRTM30 data were made available by http://dds.cr.usgs.gov/srtm/version2_1/SRTM30and GMTED30 data were downloaded from http://earthexplorer. usgs.gov.Domains d01, d02, and d03 used 10 min, 2 min, and 30 s GTOPO terrain height (respectively) with one-way feedback.All other modeling variables were held constant between the control simulations and experimental simulations.For each of the three terrain height fields, the default smoothing and interpolation methods were selected.That is, one pass of the built-in WRF Preprocessing System (WPS) smoother-desmoother and four-point averaging interpolation, respectively.

Gran Canaria case study
In this section, we analyze the atmospheric flow patterns downstream of GC as simulated by the WRF model with GTOPO30, SRTM30, and GMTED30 terrain fields.Before beginning the analysis, we compare the discrepancies between the three terrain height data fields in the left panels of Fig. 3.Here a southern view of the model terrain height for GC as generated by GTOPO30, SRTM30, and GMTED30 is shown.Notice that compared to GTOPO30, the SRTM30-based and GMTED30-based terrain height profiles have not only increased ruggedness but also led to a significant increase in peak terrain height of GC island of nearly 1 km.Furthermore, the differences between the SRTM30 and GMTED30 terrain profiles are largely insignificant.Figure 3 also illustrates the upstream mean potential temperature cross section in the lower troposphere on 30 April 2007.Within the potential temperature cross section, a well-mixed planetary boundary layer (PBL) can be identified by the nearly constant potential temperature in the lowest 800 m of the atmosphere.Above this layer, in the free atmosphere, a thermal capping inversion was present.Most importantly for the purposes of this paper, the increase in peak elevation of GC with the SRTM30 and GMTED30 data makes the modeled GC island penetrate into the stably stratified free atmosphere.
Given that the original GTOPO30-based elevation of GC was predominantly within the well-mixed PBL, the simulated flow around it was mostly 3-D.That is, the impinging air column was able to rise and cross the crest of the island barrier and then descend on the lee slope without significant buoyant restriction.This effect produced the unorganized wake pattern shown in the upper right panel of Fig. 3.Alternatively, with the SRTM30-based and GMTED30-based elevation, the increased topographic steepness along with the layer of stable stratification beneath the maximum height of the island caused much of the flow to split and pass around the lateral flanks of GC.This flow behavior generated coherent lee vortices (i.e., von Kámán vortices) which were shed downstream of the island, similar to what was observed in the MODIS satellite imagery shown in Fig. 1 (left panel).
In addition to invoking differences in the simulated wake pattern of GC, the GTOPO30-based, SRTM30-based, and GMTED30-based simulations also exhibited substantial variability in the wind regime very near to GC itself.In Fig. 4, an instantaneous streamwise wind speed cross section is presented for all three simulations.Of particular note is the wind speed extrema (greater than 17 m s −1 ) on the crest of GC in the GTOPO30-based simulation.This zone of high wind speed was a result of a Venturi-type effect caused by compression of the air column as it passed over the crest of the island.Alternatively, in the SRTM30-based and GMTED30-based simulations this zone of strong wind speed was not present due to the lack of significant air column compression over GC.Instead, the lateral flow around GC produced a zone of weak wind speed along the island centerline with respect to the flow direction.

Lesser Antilles case study
The second case study presented here deals with boundary layer flow impinging on the eastern slopes of the Lesser Antilles (LA) island archipelago.As can be seen in Fig. 2, the wake signatures from all of the major islands in this region persisted for up to approximately 300 km downstream.Contrary to the GC case, the wake patterns in the LA case did not contain strong enough vorticity to counter the ambient wind speed and therefore coherent wake vortices did not form.This type of wake pattern has been called a weak wake pattern by Smith et al. (1997), and forms in conditions of slower wind speed and lower island height in comparison to the vortex shedding patterns in the GC case.
In the left panels of   The differences in the depiction of Dominica's relief between the three simulations manifested in substantial differences in regards to the simulated 6 h mean surface wind speeds.The right panels in Fig. 5 show the mean surface wind fields simulated by the three model runs.Most notably, the weak wind wake associated with Dominica is nearly non-existent in the GTOPO30-based simulation while it extends hundreds of kilometers in the SRTM30based and GMTED30-based simulations.In addition, the zone of enhanced wind speed associated with funneling between Dominica and its northern neighbor of Guadeloupe is increased in the SRTM30/GMTED30-based simulations.Lastly, the unique shapes of the individual island wakes showed signs of variability between the GTOPO30-and SRTM30/GMTED30-based simulations.

Conclusions
In this work, we have simulated two realistic cases of atmospheric flow past mountainous islands using the WRF model.For each case, we explored the sensitivity of the simulated wake patterns with respect to three different terrain height source data sets (i.e., GTOPO30, SRTM, and GMTED2010).Our results show cases where differences in source terrain height corresponded to fundamental differences in simulated wake mechanics.For the GC case, the simulation which used GTOPO30 terrain height had a peak island elevation nearly 1 km lower than that in the SRTM30-based and GMTED30based model terrain.The SRTM30 and GMTED30 terrain, on the other hand, were very similar.As a consequence, the GTOPO30-based terrain did not reach the stably stratified thermal inversion above the planetary boundary layer while the SRTM30 and GMTED30 terrain extended hundreds of meters into the free atmosphere.This difference resulted in substantially less vertical vorticity downstream of GC island, along with an area of wind speed extrema on the crest of the island in the GTOPO30-based simulation.In other words, the SRTM30-based and GMTED30-based simulations produced more significant lateral flow around the island and downstream von Kármán vortices, in agreement with MODIS visible satellite imagery, while the GTOPO30-based simulation facilitated anomalous Venturi-type wind effect speed-up on the crest of the island and incoherent downstream vortices.
For the LA case, the GTOPO30-based model terrain represented the island of Dominica to be essentially flat and near sea level (i.e., 1 m m.s.l.) and consequently resulted in no surface wind wake pattern.At the same time, the SRTM30based and GMTED30-based simulations were almost identical and resulted in a weak wind wake field which extended hundreds of kilometers downstream of Dominica.The latter two results were similar to what was observed in visible satellite imagery and scatterometer surface wind retrievals.
This work explored the value of using representative terrain height source data for high-resolution mesoscale modeling activities.The results presented here indicate that the differences in simulated flow features associated with different terrain data sets is not a consequence of the terrain source spatial resolution but instead arise due to fundamental differences in the data sets.This conclusion is supported by the fact that significant differences were found despite first remapping the higher-resolution SRTM and GMTED data sets to 30s (equal to that of GTOPO30) prior to ingesting the data into the WRF model's preprocessing system.Moreover, this finding highlights the fact that considerable care should be taken when selecting orographic relief input data for simulating atmospheric flow over, around, and downstream of remote mountainous islands (e.g., Gran Canaria and Dominica).That being said, future studies should evaluate the uncertainty of global terrain data sets for other lo-cations and their representativeness for mesoscale modeling.At a basic level, this can be done by comparing the similarity, or dissimilarity, of available terrain data sets for the area of interest prior to performing numerical simulations.Furthermore, the use of parameterization methods which incorporate higher-level terrain data (e.g., the standard deviation of terrain height within a grid cell) may be able to provide improved terrain representation in simulations of island wakes (see Jiménez and Dudhia, 2012).

Figure 1 .
Figure 1.MODIS-TERRA visible satellite imagery of the Canary Islands on 30 April 2007 (left) and SRTM terrain height profile of Gran Canaria (right).

Figure 2 .
Figure 2. MODIS-TERRA true color image of the Lesser Antilles at 14:40 UTC on 1 August 2013 (left panel).Note the dark island wakes embedded in sunglint in the eastern part of the image.Corresponding ASCAT-B 6.25 km resolution near-surface winds at 13:26 UTC (right panel).For clarity, only every fourth wind vector is plotted.

Figure 3 .
Figure 3. WRF model terrain height profiles of Gran Canaria as viewed from the south corresponding to GTOPO30 (upper left), SRTM30 (middle left), and GMTED30 (lower left).Background color scheme represents ambient upstream potential temperature profile.Upper right, middle right, and lower right panels depict instantaneous relative vertical vorticity averaged vertically across the lowest 15 grid levels for the GTOPO30-based, SRTM30-based, and GMTED30-based simulations, respectively.

Figure 4 .
Figure 4. Instantaneous wind speed cross section for the GTOPO30-based (top panel), SRTM30-based (middle panel), and GMTED30-based (bottom panel) simulations at 06:00 UTC on 30 April 2007.Cross sections are oriented in the streamwise axis with inflow to the left.

Figure 5 .
Figure 5.The WRF model's GTOPO30-based terrain height (upper left), SRTM30-based terrain height (middle left), and GMTED30-based terrain height (lower left).Upper right, middle right, and lower right panels depict averaged wind speed in the boundary layer from 06:00 to 12:00 UTC on 1 August 2013 as simulated by the GTOPO30-based, SRTM30-based, and GMTED30-based runs, respectively.

Table 1 .
Peak elevations of the major islands in the Lesser Antilles archipelago.