High-Resolution COSMO-CLM Modeling and an Assessment of Mesoscale Features Caused by Coastal Parameters at Near-Shore Arctic Zones (Kara Sea)

Coastal Arctic regions are characterized by severe mesoscale weather events that include extreme wind speeds, and the rugged shore conditions, islands, and mountain ranges contribute to mesoscale event formation. High-resolution atmospheric modeling is a suitable tool to reproduce and estimate some of these events, and so the regional non-hydrostatic climate atmospheric model COSMO-CLM (Consortium for Small-scale Modeling developed within the framework of the international science group CLM-Community) was used to reproduce mesoscale circulation in the Arctic coast zone under various surface conditions. Mid-term experiments were run over the Arctic domain, especially over the Kara Sea region, using the downscaling approach, with ≈12 km and ≈3 km horizontal grid sizes. The best model configuration was determined using standard verification methods; however, the model run verification process raised questions over its quality and aptness based on the high level of small-scale coastline diversity and associated relief properties. Modeling case studies for high wind speeds were used to study hydrodynamic mesoscale circulation reproduction, and we found that although the model could not describe the associated wind dynamic features at all scales using ≈3 km resolution, it could simulate different scales of island wind shadow effects, tip jets, downslope winds, vortex chains, and so on, quite realistically. This initial success indicated that further research could reveal more about the detailed properties of mesoscale circulations and extreme winds by applying finer resolution modeling.


Introduction
Arctic climate and extreme weather events have attracted attention recently due to the Arctic amplification of global warming and the accompanying environmental changes. It is well known that the Arctic is the most sensitive region to global climate change, and global temperature increases have been reported to be the highest there [1][2][3]. The Arctic climate system has many complex feedbacks and teleconnections that manifest in different atmospheric circulation forms [1,[4][5][6][7], and Arctic warming exceeds the global warming signal, mainly because dynamic processes in the atmosphere provide a poleward heat advection [8].
Sea-ice area decline during the summer season also contributes to polar amplification, and although it is not the main factor [9,10], it facilitates increased extreme wind and wind wave frequencies across the Arctic Ocean region [11,12].
The Kara Sea and surrounding mainland area generally mirror these regional features. Kohnemann et al. have shown [13] maximum winter temperature increases of up to 20°C, including the surface

Model Description
The non-hydrostatic COSMO-CLM model (version 5.0) was the main tool used to simulate atmospheric dynamics in this work. This is a regional mesoscale model developed by the Consortium for Small-scale Modeling (COSMO), which includes some national weather services. The climate version of the model was developed within the framework of the international science group CLM-Community [42]. The main distinctions of this climate model from weather prediction models are modifications and extensions focused on long-term numerical experiments-such as using a deeper soil layer and taking seasonal parameter variations into account [43,44].
The COSMO-CLM model is based on primitive thermo-hydrodynamical equations describing compressible flow in a moist atmosphere. Model equations are solved in rotated latitude-longitude coordinates with a tilted pole, thus minimizing the problem of longitude convergence at the pole. The numerical scheme was implemented on the Arakawa-C grid [45]. For height coordinates, we used the terrain-following, hybrid Gal-Chen coordinate µ (σ-z system), which is an analog of the σ-coordinate from the surface, Z 0 , up to an intermediate level, Z F , and above the Z F level, it is the simple Z-coordinate [46,47]. This representation of vertical coordinates avoided problems related to surface relief inhomogeneity.
As applied to this work, the standard COSMO-CLM version 5.0 configuration included a two-time level, split-explicit, Runge-Kutta scheme, with acoustic and gravity wave splitting, together with 5th order horizontal advection numerical approximation, Smagorinsky diffusion, a two-stream radiation scheme after Ritter-Geleyn [48], and short and longwave fluxes (employing eight spectral intervals). It also included precipitation formation by bulk microphysics parameterization, including water vapor, cloud water, cloud ice, rain, and snow, with 3D transport for the precipitating phases, as well as Tiedtke mass-flux convection scheme [49] with equilibrium closure, based on moisture convergence for moist and shallow convection. Subgrid-scale turbulence was parameterized using prognostic turbulent kinetic energy closure at level 2.5, including effects from subgrid-scale condensation and thermal circulations. The model also embodied a surface layer scheme based on turbulent kinetic energy, including a laminar-turbulent roughness layer, which allowed the separation of solid surface model values based on a roughness level.
In the model, equations were split into parts associated with acoustic and gravity wave modes and those associated with slow mode atmospheric motions, in order to enhance computation efficiency. This involved evaluating the thermal buoyancy term in the equation for vertical velocity and additional stepping of the heat equation, with respect to the divergence term. The equation subset containing these terms was integrated using the Runge-Kutta scheme of the second order of precision, splitting large time steps onto smaller steps. For applications on the meso-β scale, with horizontal grid spacings of ≈10 km, the gravity wave condition became restrictive; by including gravity wave modes in the reduced set, the stability criterion was shifted to a smaller time step. This extension of the Atmosphere 2020, 11, 1062 4 of 23 mode-splitting integration scheme significantly enhanced the efficiency of model applications at the meso-β scale [50,51]. More detailed model descriptions may be found at [52].
External parameters describing some surface properties-including orography, soil properties, Leaf Area Index, Normalized Difference Vegetation Index, and so on-were obtained using the EXTPAR (external parameters) version 5.0 tool [53] from GLOBE (surface orography), MODIS (soil properties and albedo), ECOCLIMAP (forests and plants cover, root depth, land fraction, and many others), and other common datasets.
Model configurations were adopted with specific consideration of Arctic hydrometeorological conditions. In particular, the number of model vertical levels was increased to resolve surface layer processes better (50 levels total, including 10 levels in the surface layer). It was important to be able to reproduce surface wind fields more correctly, which was a major objective of this study, and increasing the vertical resolution contributed to the better reproduction of extreme wind speeds and gusts, whose calculation techniques in the model are quite simple [54,55].
The regional COSMO-CLM model has been proven suited to a broad range of tasks, including modeling high-latitude atmospheric dynamics. It has been applied to polar low dynamics experiments [37,56], for the statistical downscaling of storms over European waters [57], in a CORDEX regional downscaling project [58], and for long-term climate experiments over different areas [59][60][61], including for the Far East of Russia [62,63] and the Russian Arctic [64].

Experimental Design
Experiments were carried out using a standard nesting domain scheme; that is, we used the base domain forced by global reanalysis data as initial and boundary conditions, and the nested domain forced by modeling output over the base domain-thus reducing the horizontal grid and time step, as well as the modeling area. The location and cover of the base domain choice represent an important success determinant: firstly, they must cover the entire research area, with a margin to describe large-scale processes from the Kara Sea that influence the main meteorological features there adequately. The domain was enlarged to the N and W, considering the predominance of western motions and northern advections over the region. Moreover, since the model outputs would be used as wave model ADCIRC [65] forcing, some wave-forming patterns in the Kara Sea need to be taken into account. These were clearly defined by large-scale waves on the W boundary coming from the Barents Sea and Northern Atlantic, so we included these regions in the domain area.
ERA-Interim reanalysis data were used as initial and boundary conditions, with a horizontal grid size of 0.7 • (≈75 km) [28]; this meant that the base domain resolution had to be no less than 0.10 • to obtain the correct downscale ratio. A resolution of 0.12 • (≈13 km) was chosen, as had been used widely previously in COSMO-CLM modeling practice [66]. As for the nested domain requirements, they had to be situated deep enough toward the base domain boundaries, and they had to have a resolution up to 8 times less than that of the base domain. In due course, as settled for the model, it covered the entire Kara Sea area, spreading W, with a 0.03 • (≈2.8-3 km) grid size, as shown in Figure 1.
The model was run to cover the periods August-October 2012 and July-September 2014; in both cases, the first month was used for a "cold start" only, as it was necessary to adapt the model to the initial conditions, and results from the second and third months were used for analysis. The periods selected for our analyses were characterized by some strong storm events that were objects of interest [67]. The reason of choosing the late summer or early autumn months is that they are the most ice-free [68]. It is simplifying the analysis of surface effects on the wind speed patterns clearly, without any other side impacts. The second reason was concern with the next application of our experiments for wave simulations, for which the ice-free conditions are the most favorable in test runs [68,69]. The standard model configuration, with its nested domains, was supplemented by versions using the "spectral nudging" technique (further marked as *_sn), the reduced time step, dt (further marked as *_dt), and the enlarged nested domain (further marked as *_large). These changes were suggested by the following considerations. Studies [70][71][72] have shown that "spectral nudging" contributed to a better assimilation of large-scale meteorological fields, on account of using additional forcing data-in this case, reanalysis. This applied not only at the domain boundaries but also inside the domain, which further restricted the model from departing from real conditions. Accounting for large-scale features was based on two-dimensional Fourier decomposition of meteorological fields from the regional model and reanalysis of regional modeling with subsequent correction. "Nudging" was performed to temperature and wind fields, with 500 km set as the minimal scale for assimilated atmospheric processes from the 850 mb pressure level and above, since middle-and upper-tropospheric circulation systems generally control atmospheric dynamics.
Model time step changes influenced the Courant-Friedrichs-Levy (CFL) numerical criterion, which was crucial for strong wind reproduction: the object of interest in this work [73,74]. Standard experiments were run with a time step (dt) = 100 s and then with reduced time steps (dt) = 80 s (*_dt) over the base domain of ≈13 km, after which both experiments were downscaled to the domain of ≈3 km.
It has been reported that domain boundary changes, that is, domain enlargement or reduction, influence boundary condition information propagation and affect model adaptation to changing boundary conditions [75][76][77]. In this case, the nested domain, with a horizontal resolution of ≈3 km, was enlarged to the S on 50 model grids (*_large). This was done to account for continental properties S of the Gulf of Ob, including some considered storm conditions formed in southern streams.
A total of 14 experiments were conducted. The model was run in different configurations for the August-October 2012 and July-September 2014 periods, using the nested domains scheme. All experiments were run using shared high-performance computing resource research facilities (supercomputer "Lomonosov-2" [78]) at Lomonosov Moscow State University.
All listed experiments were verified by surface 10 m wind speeds, using all coastal and inshore meteorological station data over the Kara Sea region and surrounding areas (Figure 2), as available on www.meteo.ru [79]. Verification was carried out using the simple "nearest neighbor" method; that is, station data were compared with those of the nearest station, based on the coordinates model Interesting storm events were observed over the Kara Sea region during 9-14 October 2012 and 17-23 October 2012. These cases were related to the rapid intensification of a narrow-scaled frontal cyclone centered over the White Sea coast (see Figure S1, Supplementary Materials) and its movement through the Kara Sea ( Figure S2, Supplementary Materials). Generally, the 2012 modeling period was characterized by W and NW cyclonic activity, which was intercepted by separate anticyclones from the N and NE. S and SE anticyclonic conditions, as well as N cyclonic conditions, were more frequent during the 2014 modeling period.
The standard model configuration, with its nested domains, was supplemented by versions using the "spectral nudging" technique (further marked as *_sn), the reduced time step, dt (further marked as *_dt), and the enlarged nested domain (further marked as *_large). These changes were suggested by the following considerations. Studies [70][71][72] have shown that "spectral nudging" contributed to a better assimilation of large-scale meteorological fields, on account of using additional forcing data-in this case, reanalysis. This applied not only at the domain boundaries but also inside the domain, which further restricted the model from departing from real conditions. Accounting for large-scale features was based on two-dimensional Fourier decomposition of meteorological fields from the regional model and reanalysis of regional modeling with subsequent correction. "Nudging" was performed to temperature and wind fields, with 500 km set as the minimal scale for assimilated atmospheric processes from the 850 mb pressure level and above, since middle-and upper-tropospheric circulation systems generally control atmospheric dynamics.
Model time step changes influenced the Courant-Friedrichs-Levy (CFL) numerical criterion, which was crucial for strong wind reproduction: the object of interest in this work [73,74]. Standard experiments were run with a time step (dt) = 100 s and then with reduced time steps (dt) = 80 s (*_dt) over the base domain of ≈13 km, after which both experiments were downscaled to the domain of ≈3 km.
It has been reported that domain boundary changes, that is, domain enlargement or reduction, influence boundary condition information propagation and affect model adaptation to changing boundary conditions [75][76][77]. In this case, the nested domain, with a horizontal resolution of ≈3 km, was enlarged to the S on 50 model grids (*_large). This was done to account for continental properties S of the Gulf of Ob, including some considered storm conditions formed in southern streams.
A total of 14 experiments were conducted. The model was run in different configurations for the August-October 2012 and July-September 2014 periods, using the nested domains scheme. All experiments were run using shared high-performance computing resource research facilities (supercomputer "Lomonosov-2" [78]) at Lomonosov Moscow State University.
All listed experiments were verified by surface 10 m wind speeds, using all coastal and inshore meteorological station data over the Kara Sea region and surrounding areas (Figure 2), as available on www.meteo.ru [79]. Verification was carried out using the simple "nearest neighbor" method; that is, station data were compared with those of the nearest station, based on the coordinates model grid, and interpolation methods were not applied. Such an approach was chosen so as not to introduce noise into the model fields and to reveal the impact of surface properties on modeling results. Additional verification of some global reanalyses-namely ERA-Interim [28], NCEP-CFSv2 [32], and ERA5 [31]-was performed using the same technique to compare and estimate the "added value", or contribution, for mesoscale modeling. It should be noted that all stations are located either on the continental coast or on an island in diverse coastline structural conditions and at reliefs of various spatial scales (from km to tens of km). This allowed the effects of surface condition diversity within different sectors, scales, and grid sizes on modeling results to be assessed in detail.
Atmosphere 2020, 11, x FOR PEER REVIEW 6 of 23 grid, and interpolation methods were not applied. Such an approach was chosen so as not to introduce noise into the model fields and to reveal the impact of surface properties on modeling results. Additional verification of some global reanalyses-namely ERA-Interim [28], NCEP-CFSv2 [32], and ERA5 [31]-was performed using the same technique to compare and estimate the "added value", or contribution, for mesoscale modeling. It should be noted that all stations are located either on the continental coast or on an island in diverse coastline structural conditions and at reliefs of various spatial scales (from km to tens of km). This allowed the effects of surface condition diversity within different sectors, scales, and grid sizes on modeling results to be assessed in detail.

Results
In this section, we have reviewed the results achieved by all the listed model configurations for both experimental periods, with some statistical outcomes summarized in Tables 1 and 2.

Results
In this section, we have reviewed the results achieved by all the listed model configurations for both experimental periods, with some statistical outcomes summarized in Tables 1 and 2. The standard configuration experiment verification obtained correlation coefficients between wind speed timeseries at stations and the corresponding nearest model grids in the range 0.5-0.7 (0.6 on average) on the base domain, with a ≈13 km grid; for the 2012 period, the worst values were at Russki Island and the named Krenkel stations (≈0.45). Mean biases were quite satisfactory for most stations (less than 1 m/s, with a −0.08 average, except for Cape Bolvanskiy Nos, Antipayuta, and Malye Karmakuly stations), indicating the realistic reproduction of synoptic-scale process dynamics and variability, for two months. Wind speed root mean square errors (RMSEs) were within satisfactory ranges (2.5-3 m/s), except for the above-mentioned stations, which gave values of > 3.5-4 m/s. In particular, large errors at the Malye Karmakuly station were associated with the extreme wind speeds often observed there, which were caused by the well-known Novaya Zemlya Bora [80,81]. Mesoscale processes, highlands, and rugged shorelines contributed significantly to the formation and variability of the downslope winds.
Statistical performance was overall the same for the nested domain; however, the RMSE decreased at some stations, especially at those that had been the highest (mean RMSE was 2.25, compared to 2.84). The same tendency was indicated for biases, and it was also noted that biases and median errors could get negative values for the ≈3 km domain, which could be explained by the fact that mesoscale circulations were reproduced at some locations that were different from the actual observed sites. This problem has been reviewed in the next section. The 2014-period experiment gave rise to increased biases and median errors up to 0.5 m/s, whereas the RMSEs remained at the same level. The correlation coefficients had a large spread and were all < 0.7 (Antipayuta, Dikson Island), while RMSEs exceeded 3, and even 3.5 m/s, over almost half of the stations (4.4 m/s in Malye Karmakuly).
When we considered the experimental results achieved by the other model configurations, we first noted that the "spectral nudging" technique improved the reproduction of wind speeds in all cases and by all metrics compared to the standard configuration. The mean correlation coefficients increase up to 0.75-0.77, and RMSEs decreased to 2.1-2.2 m/s. This effect was seen across almost all stations: biases were < 1 m/s, except at Antipayuta station, while the RMSE at Malye Karmakuly station decreased to 2.5 m/s for the 2014 period and 3.6 m/s for the 2012 experimental period. The error probability distribution function analysis (see examples in Figure 3) allowed us to conclude that extreme wind speeds could be reproduced better using finer resolution, while background values remained at the same levels. All these results encouraged us to apply the "spectral nudging" technique in our further experiments.
Atmosphere 2020, 11, x FOR PEER REVIEW 8 of 23 The important role of the "spectral nudging" technique in wind speed reproduction has been illustrated using the graphical example in Figure 4, which represents the case at 09:00 GMT on 10 October 2012. A strong SE flow was observed at that time, coming from the Kara Sea and crossing the Novaya Zemlya Island ridges (see Figure S1 of the Supplementary Materials), creating wind speeds up to 25 m/s according to Malye Karmakuly station data. The synoptic conditions were determined by the sea-level pressure (SLP) dipole pattern composed of low pressure over the southern Barents Sea coast and high pressure over the northern Western Siberia. This pattern contributed to the steady strong SE flow. It can be well seen that "spectral nudging" contributed to shifting the area that experienced the maximum winds (≈20 m/s) much closer to the station location, thus creating a more realistic pattern. It was also apparent that reducing the time step did not change the pattern (Figure 4, right panel). Investigations into using finer resolution and its impact on the reproduction of wind speed and mesoscale circulation, including this example, have been considered in the next section. The important role of the "spectral nudging" technique in wind speed reproduction has been illustrated using the graphical example in Figure 4, which represents the case at 09:00 GMT on 10 October 2012. A strong SE flow was observed at that time, coming from the Kara Sea and crossing the Novaya Zemlya Island ridges (see Figure S1 of the Supplementary Materials), creating wind speeds up to 25 m/s according to Malye Karmakuly station data. The synoptic conditions were determined by the sea-level pressure (SLP) dipole pattern composed of low pressure over the southern Barents Sea coast and high pressure over the northern Western Siberia. This pattern contributed to the steady strong SE flow. It can be well seen that "spectral nudging" contributed to shifting the area that experienced the maximum winds (≈20 m/s) much closer to the station location, thus creating a more realistic pattern. It was also apparent that reducing the time step did not change the pattern (Figure 4, right panel). Investigations into using finer resolution and its impact on the reproduction of wind speed and mesoscale circulation, including this example, have been considered in the next section.
Atmosphere 2020, 11, x FOR PEER REVIEW 9 of 23 Thus, the "spectral nudging" technique was used to improve wind speed reproduction compared to the standard model configuration. Since it plays a significant role in the large-scale conditions' reproduction, the "spectral nudging" is the best choice to assimilate large-scale meteorological fields additionally inside the model domain. At the same time, there was no significant advantage obtained in using reduced time steps, with the improvements achieved being limited to just the second decimal places in the correlation coefficients and biases. Summing these standard verification results, we were able to conclude that the optimal model configuration included "spectral nudging" and time step reduction up to dt = 80 s over the base domain of ≈13 km. Experiments applying this configuration were designated as *_sn_dt.
These verification results raised a question about the comparability of model experimental quality with existing hydrometeorological datasets, such as reanalyses. Three modern, current generation, fine-resolution reanalyses were chosen for comparison-ERA-Interim (0.75°) [28], ERA5 (0.3°) [31], and NCEP-CFSv2 (0.2°) [32]-and the results have been listed in Tables 1 and 2. Reanalysis ERA5 showed the best results for our experimental periods. The best model configuration with "spectral nudging" was of a quality comparable with the reanalyses, but it was slightly less accurate according to most statistics. This result could be explained by the fact that the reanalyses represented the complete assimilation of most available data, including from satellites, while COSMO-CLM uses reanalysis data as forcing only, with no additional data assimilation. At the same time, it should be noted that the obtained scores of the best model configuration agree with estimations of regional climate modeling and reanalyses performance at the Arctic region [82][83][84]. Therefore, taking these circumstances into account, we could consider our simulation results as quite satisfactory compared to reanalyses.
However, considering the statistics for individual stations, we were able to see that those with the highest errors-such as Antipayuta, Dikson Island, and Malye Karmakuly-corresponded to the worst reanalysis statistics and were worse than the model. In particular, the RMSE for the 2012 test period at Malye Karmakuly station was 4.2 m/s, while the median error at Dikson Island station was 1.6 m/s. This meant that in some cases, the regional model caught the most intense wind speeds formed by mesoscale processes over the 2-month averaging period better. Taking into account this circumstance, we considered the model experimental results to have been satisfactory and reasonable. The above-mentioned biases could have different sources including model errors, inconsistencies of surface (e.g., land, sea, soil type), altitude, long distance between considered Thus, the "spectral nudging" technique was used to improve wind speed reproduction compared to the standard model configuration. Since it plays a significant role in the large-scale conditions' reproduction, the "spectral nudging" is the best choice to assimilate large-scale meteorological fields additionally inside the model domain. At the same time, there was no significant advantage obtained in using reduced time steps, with the improvements achieved being limited to just the second decimal places in the correlation coefficients and biases. Summing these standard verification results, we were able to conclude that the optimal model configuration included "spectral nudging" and time step reduction up to dt = 80 s over the base domain of ≈13 km. Experiments applying this configuration were designated as *_sn_dt.
These verification results raised a question about the comparability of model experimental quality with existing hydrometeorological datasets, such as reanalyses. Three modern, current generation, fine-resolution reanalyses were chosen for comparison-ERA-Interim (0.75 • ) [28], ERA5 (0.3 • ) [31], and NCEP-CFSv2 (0.2 • ) [32]-and the results have been listed in Tables 1 and 2. Reanalysis ERA5 showed the best results for our experimental periods. The best model configuration with "spectral nudging" was of a quality comparable with the reanalyses, but it was slightly less accurate according to most statistics. This result could be explained by the fact that the reanalyses represented the complete assimilation of most available data, including from satellites, while COSMO-CLM uses reanalysis data as forcing only, with no additional data assimilation. At the same time, it should be noted that the obtained scores of the best model configuration agree with estimations of regional climate modeling and reanalyses performance at the Arctic region [82][83][84]. Therefore, taking these circumstances into account, we could consider our simulation results as quite satisfactory compared to reanalyses.
However, considering the statistics for individual stations, we were able to see that those with the highest errors-such as Antipayuta, Dikson Island, and Malye Karmakuly-corresponded to the worst reanalysis statistics and were worse than the model. In particular, the RMSE for the 2012 test period at Malye Karmakuly station was 4.2 m/s, while the median error at Dikson Island station was 1.6 m/s. This meant that in some cases, the regional model caught the most intense wind speeds formed by mesoscale processes over the 2-month averaging period better. Taking into account this circumstance, we considered the model experimental results to have been satisfactory and reasonable. The above-mentioned biases could have different sources including model errors, inconsistencies of surface (e.g., land, sea, soil type), altitude, long distance between considered points, and some others. Some cases of strong wind reproduction studied in detail in the Discussion section.
Model resolution refinement led to a reduced error maxima for most stations, and narrowing of the pdf errors, while some stations regularly exhibited high errors. This could be associated with circulation features, such as extreme wind speeds, or with mesoscale surface properties-or even the basic realism of the model. The reproduction of individual case studies using the model gave better results than did reanalysis using coarser resolution. Precise extreme event reproduction, rather than descriptions of averaged parameters, is often the main advantage of mesoscale modeling applications. Our results so far have indicated that a detailed analysis of case studies within different spatial scales for diverse domain regions and experiments should be the next step in our research alongside an attempt to interpret the impact of relief on various coastal zone mesoscale effects.

Discussion
Solving the task of realistic wind speed modeling in coastal zones is important, because it is in essence the only data option where there is poor monitoring network coverage. It was clear that the ability of the model to reproduce certain classes of hydrodynamic events was determined by its scale and resolution. Some of these events will be analyzed further, in pursuit of the goal of clarifying the model limitations and its ability to reproduce certain mesoscale effects.
Model wind speed outputs were investigated during storms that occurred over the periods 9-11 October 2012 and 19-20 October 2012, as recorded by the named Popova and Dikson Island stations around Malye Karmakuly, where different obstacle overflows occur, and their effects on the wind regime could be studied. The purpose here was to compare modeling results with meteorological observations. Airflow interactions with surface objects-such as mountainous terrain, islands, coastlines, and a combination of these factors-cause certain established effects [85]. These can include internal gravity wave [86] and vortex chain formation [87,88], slope catabatic winds, cyclogenesis, and wind shadows [89,90]. These interactions can also give rise to wind speed strengthening over peaks and "corridors" formed by coastlines, horizontal vortex formation behind obstacles [86], wind speed gains from overflowing narrow capes (tip jets) [20], and so on. It should be emphasized that these effects have rarely been documented through field experiments, in numerical solutions with simplified formulations, or in laboratory research. In particular, idealized obstacle shapes and steady stream conditions were considered, which have not been observed in nature.
In this research, we applied the model across a broad spectrum of hydrodynamic effects to simulate real flows traversing obstacles. As prepared, model resolution was sufficiently coarse to not only reconstruct all the effects listed above explicitly but also to measure its integral responses to these effects. In fact, the model needed to be sized with its x, y, and z grid dimensions all at ≈10 m to reproduce internal gravity waves. This was determined by the fact that gravity wave frequency is defined by the Brunt-Väisälä frequency (N ≈0.01 s −1 ), and the wavelength is determined by the wave scale defined by either the Lyra scale or the inverse Scorer parameter where U is the leaking flow velocity. It was not possible to reproduce vortices in a numerical experiment using a sufficiently coarse resolution. However, the footprint of these structures, when repeated over several obstacle diameters, could be identified as a linearly elongated, enhanced turbulence zone in the coarse grid model. For example, it appeared that under some conditions, it could capture a vortex chain appearance representing the vortices (generated from both sides) behind the obstacle. This was manifested in the higher turbulence activity of the flow, including enhanced vorticity values.
Wind strengthening along slopes occurs because the wind is not normal to the slope, and thus, the flow is turning more in the corresponding direction. This phenomenon has also been associated with changes in Coriolis force dynamics [91], although the flow needs to transit a sufficiently long fetch along the meridionally elongated slope in order to capture this Coriolis force phenomenon. Issues concerning cyclogenesis and other synoptic-scale processes were not considered in our research for the same reason-that is, they occur beyond the scale of the events that we studied.
The Belyi Island region, where the named Popova station is situated, can be seen in Figure 5. This island, which is approximately 40 km across and flat with no uplands, is separated from the Yamal Peninsula by the 10 km-wide Malygin Strait. The region experienced strong S winds, at ≈15 m/s, on 11 October, 2012, at 1000 h GMT, and this event was established through almost the entire troposphere (see Figure S3a, Supplementary Materials). Such steady flow was caused by a combination of the deep cyclone centered over the White Sea and the Siberian high. By analyzing the *sn_dt experiment with a ≈3 km grid (Figure 5), we saw wind speed reduced over the island and observed the wind shadow phenomenon with a significant wind speed reduction (by 5 m/s and more) on the leeward (N) coast. This zone extended out over the sea for 20-30 km, and the island was framed by jets of different intensities and with widths between 10 and 20 kms. It should be noted that the observed wind speeds at 10 m height were 11-13 m/s, which were only 1-2 m/s above the modeled velocities. These results showed us that the COSMO-CLM model could capture events at the spatial scale of islands, changes in surface type (land-sea), island impacts on wind field structure, and interactions with the sea.
Atmosphere 2020, 11, x FOR PEER REVIEW 11 of 23 The Belyi Island region, where the named Popova station is situated, can be seen in Figure 5. This island, which is approximately 40 km across and flat with no uplands, is separated from the Yamal Peninsula by the 10 km-wide Malygin Strait. The region experienced strong S winds, at ≈15 m/s, on 11 October, 2012, at 1000 h GMT, and this event was established through almost the entire troposphere (see Figure S3a, Supplementary Materials). Such steady flow was caused by a combination of the deep cyclone centered over the White Sea and the Siberian high. By analyzing the *sn_dt experiment with a ≈3 km grid (Figure 5), we saw wind speed reduced over the island and observed the wind shadow phenomenon with a significant wind speed reduction (by 5 m/s and more) on the leeward (N) coast. This zone extended out over the sea for 20-30 km, and the island was framed by jets of different intensities and with widths between 10 and 20 kms. It should be noted that the observed wind speeds at 10 m height were 11-13 m/s, which were only 1-2 m/s above the modeled velocities. These results showed us that the COSMO-CLM model could capture events at the spatial scale of islands, changes in surface type (land-sea), island impacts on wind field structure, and interactions with the sea. The relative vorticity pattern ( Figure 6) was determined by elongated alternating bands of positive and negative values, which were 100 km or more long and up to 10 km wide. These bands were distinguishable by given ≈3 km grids-that is, they could be determined by at least four grid points considering the relative vorticity calculation technique. This was unrelated to the influence of the island itself, except for one feature-the manifestation of positive vorticity maxima on the E coast. The relative vorticity pattern ( Figure 6) was determined by elongated alternating bands of positive and negative values, which were 100 km or more long and up to 10 km wide. These bands were distinguishable by given ≈3 km grids-that is, they could be determined by at least four grid points considering the relative vorticity calculation technique. This was unrelated to the influence of the island itself, except for one feature-the manifestation of positive vorticity maxima on the E coast. In considering how to explain this phenomenon, we introduced natural coordinates t and n [92]. In this solution, unit vector t was oriented parallel to the horizontal velocity at each point, while unit vector n was normal to the horizontal velocity and oriented so as to be positive to the left of the flow direction. The relative vorticity in such coordinates would be defined as Here, R denotes the radius of curvature following the parallel motion. It followed that decreasing the wind speed module by transition from the background sea flow to the shadow zone should have led to an intensification of cyclonic vorticity. It was interesting that on the island W coast, the corresponding effect of an opposite sign, that is, intensification of anticyclonic vorticity, was expressed only weakly; the possible reason for this difference could have been that there was no advection of the necessary vorticity inbuilt into the oncoming flow.
The existence of elongated vortices of opposite signs alternating inside the flow could be presented as a manifestation of vortex chains-and these were not those Karman vortex streets related genetically to the flow around islands. This phenomenon reminded us of narrow zones of higher wind speeds in turbulence flows, which are sometimes referred to as "strings". Their generation and disappearance originate from dissipation bursts associated with the turbulence intermittency process, namely, irregular autofluctuations of flow with magnitudes comparable with the flow speed itself [93][94][95]. Although the genesis of this phenomenon is completely different from that considered here, the similarity of patterns manifesting in coherent structural formations has been attracting attention.
In the second example, we considered wind behavior at the E Kara Sea coast in the vicinity of the Yenisey Gulf, looking at Sibiryakov and Dikson islands (the latter having dimensions of 3-4 km). Dikson Island is separated from the Taymyr Peninsula by a narrow (1-2 km) strait ( Figure 7) and features rocky terrain (unlike Belyi Island), although this terrain does not rise above 25 m. In considering how to explain this phenomenon, we introduced natural coordinates t and n [92]. In this solution, unit vector t was oriented parallel to the horizontal velocity at each point, while unit vector n was normal to the horizontal velocity and oriented so as to be positive to the left of the flow direction. The relative vorticity in such coordinates would be defined as Here, R denotes the radius of curvature following the parallel motion. It followed that decreasing the wind speed module by transition from the background sea flow to the shadow zone should have led to an intensification of cyclonic vorticity. It was interesting that on the island W coast, the corresponding effect of an opposite sign, that is, intensification of anticyclonic vorticity, was expressed only weakly; the possible reason for this difference could have been that there was no advection of the necessary vorticity inbuilt into the oncoming flow.
The existence of elongated vortices of opposite signs alternating inside the flow could be presented as a manifestation of vortex chains-and these were not those Karman vortex streets related genetically to the flow around islands. This phenomenon reminded us of narrow zones of higher wind speeds in turbulence flows, which are sometimes referred to as "strings". Their generation and disappearance originate from dissipation bursts associated with the turbulence intermittency process, namely, irregular autofluctuations of flow with magnitudes comparable with the flow speed itself [93][94][95]. Although the genesis of this phenomenon is completely different from that considered here, the similarity of patterns manifesting in coherent structural formations has been attracting attention.
In the second example, we considered wind behavior at the E Kara Sea coast in the vicinity of the Yenisey Gulf, looking at Sibiryakov and Dikson islands (the latter having dimensions of 3-4 km). Dikson Island is separated from the Taymyr Peninsula by a narrow (1-2 km) strait ( Figure 7) and features rocky terrain (unlike Belyi Island), although this terrain does not rise above 25 m. Atmosphere 2020, 11, x FOR PEER REVIEW 13 of 23 A significant wind speed increase was seen in the Dikson stream (> 20 m/s), which had a width of up to 10 km. Its formation was determined by the declining friction force there, as well as by the rocky coastal terrain, creating a peculiar wind corridor, and we noted that the stream conserved its shape at the exit from the strait. The observed wind speeds were 21-23 m/s, which corresponded to the stream wind speed reproduced by the model. Vorticity pattern analysis (Figure 8) illustrated this vorticity band alternation of opposite signs inside the flow. Cyclonic vorticity maxima were noted again along the E coasts of Sibiryakov and Dikson islands, and they were analogous to that considered in Belyi Island, although the effect of the anticyclonic vorticity maximum appeared on the W coast of the larger Sibiryakov Island. A significant wind speed increase was seen in the Dikson stream (>20 m/s), which had a width of up to 10 km. Its formation was determined by the declining friction force there, as well as by the rocky coastal terrain, creating a peculiar wind corridor, and we noted that the stream conserved its shape at the exit from the strait. The observed wind speeds were 21-23 m/s, which corresponded to the stream wind speed reproduced by the model. Vorticity pattern analysis (Figure 8) illustrated this vorticity band alternation of opposite signs inside the flow. Cyclonic vorticity maxima were noted again along the E coasts of Sibiryakov and Dikson islands, and they were analogous to that considered in Belyi Island, although the effect of the anticyclonic vorticity maximum appeared on the W coast of the larger Sibiryakov Island. We also considered a case study of strong winds on the W coast of Novaya Zemlya Island on 10 October 2012, as observed at Malye Karmakuly station. We had noted before that the highest wind speeds over the Eurasian Arctic are observed at Novaya Zemlya [96,97]. This region is characterized by the presence of bays, bights, rugged shores, and a 700 m + mountain range standing a few km inland. Such terrain favors katabatic stock wind formation (such as the Novaya Zemlya Bora), which is strengthened over different coastal features. Theoretical approaches based on simplified models of the dynamic interaction of the incoming flow with topography, such as hydraulic and wave approaches, have been widely used to describe downslope windstorms [98]. This and other early results have been confirmed by mesoscale modeling data in some respects [80,99], where a successful attempt has been made to explain small-scale, non-modeling explicit processes by applying integral hydrodynamic parameters.
The experimental ≈3 km grid, using "spectral nudging", obtained realistic results in our study, corresponding to measurement data for Malye Karmakuly station. The maximum modeled wind speed (25 m/s) was observed over the entire length of the station's bight (Figure 9), and there was quite clear wind speed enhancement down the leeward slope.
The wind speed pattern became less homogeneous as it came ashore, reflecting the inhomogeneity of the background flow, as well as the influence of the complicated coastline configuration. It should be noted that these effects were not apparent in the ≈13 km grid experiment, showing that the model responded to more detailed surface descriptions. The wind shadow was not clearly pronounced. Ten to 20 km-wide streams were evident, as in previous examples, as well as areas of significantly lower wind speed (by 10-15 m/s), which were stretched along river valleys, with the latter wind speed reduction most likely associated with valley surface roughness. We also considered a case study of strong winds on the W coast of Novaya Zemlya Island on 10 October 2012, as observed at Malye Karmakuly station. We had noted before that the highest wind speeds over the Eurasian Arctic are observed at Novaya Zemlya [96,97]. This region is characterized by the presence of bays, bights, rugged shores, and a 700 m + mountain range standing a few km inland. Such terrain favors katabatic stock wind formation (such as the Novaya Zemlya Bora), which is strengthened over different coastal features. Theoretical approaches based on simplified models of the dynamic interaction of the incoming flow with topography, such as hydraulic and wave approaches, have been widely used to describe downslope windstorms [98]. This and other early results have been confirmed by mesoscale modeling data in some respects [80,99], where a successful attempt has been made to explain small-scale, non-modeling explicit processes by applying integral hydrodynamic parameters.
The experimental ≈3 km grid, using "spectral nudging", obtained realistic results in our study, corresponding to measurement data for Malye Karmakuly station. The maximum modeled wind speed (25 m/s) was observed over the entire length of the station's bight (Figure 9), and there was quite clear wind speed enhancement down the leeward slope. Atmosphere 2020, 11, x FOR PEER REVIEW 15 of 23 Airflow perturbation is a characteristic feature of bora and downslope winds, and the wind speed module standard deviation has usually been used to replicate the natural turbulence metric. This quantity was calculated over the period of the feature's existence (≈21 h)-that is, for a period of steady SE flow-and showed that flow turbulence reached its maximum over the Novaya Zemlya shore, as evidenced by the highest speed tip bands being immediately adjacent to the coast ( Figure 10).
The relative vorticity pattern can be seen in Figure 11. Cyclonic and anticyclonic vorticity areas were clear, with horizontal diameters of 10-30 km and with strong vorticity gradients over land corresponding to stock downslope winds. Linearly stretched vorticity zones of the same sign were quite clear again; they extended for > 100 km, were 10-15 km wide, and were probably enhanced by tip jets. Vorticity bands could also be traced in river valleys. The wind speed pattern became less homogeneous as it came ashore, reflecting the inhomogeneity of the background flow, as well as the influence of the complicated coastline configuration. It should be noted that these effects were not apparent in the ≈13 km grid experiment, showing that the model responded to more detailed surface descriptions. The wind shadow was not clearly pronounced. Ten to 20 km-wide streams were evident, as in previous examples, as well as areas of significantly lower wind speed (by 10-15 m/s), which were stretched along river valleys, with the latter wind speed reduction most likely associated with valley surface roughness.
Airflow perturbation is a characteristic feature of bora and downslope winds, and the wind speed module standard deviation has usually been used to replicate the natural turbulence metric. This quantity was calculated over the period of the feature's existence (≈21 h)-that is, for a period of steady SE flow-and showed that flow turbulence reached its maximum over the Novaya Zemlya shore, as evidenced by the highest speed tip bands being immediately adjacent to the coast ( Figure 10). shore, as evidenced by the highest speed tip bands being immediately adjacent to the coast ( Figure 10).
The relative vorticity pattern can be seen in Figure 11. Cyclonic and anticyclonic vorticity areas were clear, with horizontal diameters of 10-30 km and with strong vorticity gradients over land corresponding to stock downslope winds. Linearly stretched vorticity zones of the same sign were quite clear again; they extended for > 100 km, were 10-15 km wide, and were probably enhanced by tip jets. Vorticity bands could also be traced in river valleys. The relative vorticity pattern can be seen in Figure 11. Cyclonic and anticyclonic vorticity areas were clear, with horizontal diameters of 10-30 km and with strong vorticity gradients over land corresponding to stock downslope winds. Linearly stretched vorticity zones of the same sign were quite clear again; they extended for > 100 km, were 10-15 km wide, and were probably enhanced by tip jets. Vorticity bands could also be traced in river valleys.  It was clear to us that transition to a more detailed surface description played a key role in the model's ability to simulate actual wind patterns. The model showed itself to be capable of reproducing downslope windstorms, although the physical mechanisms responsible for this bora phenomenon could not yet be completely implemented in the model.
The mesoscale strong wind speeds modeled in the case studies were compared to global reanalysis data ERA5 and NCEP-CFSv2 for the experimental periods, with the results available in Figures S4 and S5 of the Supplementary Materials. These comparisons clearly showed that the It was clear to us that transition to a more detailed surface description played a key role in the model's ability to simulate actual wind patterns. The model showed itself to be capable of reproducing downslope windstorms, although the physical mechanisms responsible for this bora phenomenon could not yet be completely implemented in the model.
The mesoscale strong wind speeds modeled in the case studies were compared to global reanalysis data ERA5 and NCEP-CFSv2 for the experimental periods, with the results available in Figures S4  and S5 of the Supplementary Materials. These comparisons clearly showed that the differences in surface wind patterns were caused by using a coarser grid, which resulted in surface wind descriptions that were not as good as those of the reanalyses. However, as mentioned above, when case studies of extreme events were simulated by the mesoscale model, the results were significantly better in comparison with the reanalysis data.

Conclusions
The principal task of this research was to review the capability of the mesoscale model, COSMO-CLM, to reproduce sophisticated atmospheric circulation features related to the impact of surface properties on airflow. The main focus of the study was on high-latitude atmospheric dynamics, which have been closely reviewed in this paper. The investigation focused on the developing influence of obstacles, such as islands and moderately mountainous terrain, on the mesoscale processes in the atmosphere. The authors were aware that a horizontally gridded model with a scale of just a few km would not be capable of explicitly reproducing the entire spectrum of such dynamic effects, and they kept in mind that internal gravity waves were manifested in different meteorological phenomena, such as wind speed acceleration within downslope flows and circulation vortex chain formation. At the same time, we expected to detect integral effects at least, as gravity wave effects were parametrized in the model as wave resistance manifested as flow around surface obstacles.
From another point of view, well-known ideas about the physical phenomena resulting from flow around obstacles were reproduced using significantly simplified approaches. Therefore, the analysis achieved within the framework of the physically detailed COSMO-CLM model was still able to improve our event representations, despite the capability limitations mentioned above.
Analysis of model patterns showed that the mesoscale features appeared clearly on a ≈3 km grid in all the cases that were examined. In contrast, the same features either were or were not revealed by the ≈13 km grid, often appearing either indistinctly or not at all-or being translocated. This suggested that the model captured details adequately and allowed us to expect that coastal mesoscale circulations could be reproduced quite reliably. This is important for applications such as pollution calculation, coastal wind wave simulation, and so on.
The successful reproduction of mesoscale streams alternating inside a flow, that is, elongated bands of cyclonic and anticyclonic vorticity, was an important outcome from this work. We also showed that the model could reflect mesoscale circulation organization, in the form of a certain rotation vortex system (known as a hydrodynamic property of flow), and that these bands were not related to relief. The influence of obstacles met in a wind path manifests in vorticity enhancement (of the corresponding sign) on coasts, as flows wrap themselves around the obstacles. Flow organization in the form of vorticity zones called to mind the phenomena known as "intermittency zones", or "streams", which have been distinguished in turbulence flows.
Simulating flow around islands illustrated wind shadow production, whose parameters were influenced by surface roughness. These phenomena were found to extend along the flow by more than simply the size of the island, and this is an effect worthy of future research, using islands of different shapes and sizes.
Wind speed and turbulence strengthening were reproduced at the leeward base of over mountain flows, reflecting the main features of downslope wind development (the Novaya Zemlya Bora phenomenon). Moreover, although the model could not reconstruct all wind event mechanisms, those phenomena that it did simulate showed characteristics close to reality.
We noted that the high-resolution COSMO-CLM model created realistic patterns at the expense of imitation mechanisms defined by subgrid process parameterizations. Our results suggested that mesoscale modeling could be a reliable source of information on mesoscale circulations in areas such as the Arctic, where there are few weather observation platforms.
Modeling verification results should be considered critically with respect to the "nearest neighbor" technique, concerning the above-mentioned conclusion. Indeed, since station data are irreplaceable for local environmental and geographic conditions, and model simulations cannot fully replicate natural patterns, the usefulness of comparing model data with station observations becomes moot at a certain point. This implied that the further development of verification techniques, moving away from the direct grid-station point comparisons, is needed.
Model capability in reproducing some wind patterns reliably and physically, including extreme winds, will really help investigation into the physical description of these objects by applying different hydrodynamical parameters. We have shown that existing global and regional reanalysis data are unable to describe severe weather event patterns in coastal regions, although they can surpass mesoscale modeling with respect to mean climatic statistics. Reanalysis data will also facilitate the future use of regional climate models, with horizontal grid steps of just a few kilometers, to carry out adequate high-resolution wind climatological assessments, in conditions where other hydrometeorological information sources are absent. These studies could have wide applications in resolving actual tasks associated with estimating Arctic climate change. The work described in this paper will be continued in this direction by studying the quantitative parameters of mesoscale effects, assessing the influence of coastline characteristics, measuring the impact of terrain on wind properties, and so on.

Funding:
The reported study was funded by RFBR according to the research project № 18-05-60147 and Lomonosov Moscow State University project no. AAAA-A16-116032810086-4. The APC was funded by the RFBR research project № 18-05-60147.