Sensitivity of Mesoscale Model Forecasts of Icing Conditions to Land Use: Dynamic and Diabatic Impacts

: In-cloud ice mass accretion on wind turbines is a common challenge faced by energy companies operating in cold climates. On-shore wind farms in Scandinavia are often located in regions near patches of forest, the heterogeneity length scales of which are often less than the resolution of many numerical weather prediction (NWP) models. The representation of these forests–including the cloud water response to surface roughness and albedo effects related to them–must therefore be parameterized in NWP models used as meteorological input in ice prediction systems, resulting in an uncertainty that is poorly understood and to present date not quantiﬁed. The sensitivity of ice accretion forecasts to the subgrid representation of forests is examined in this study. A single column version of the HARMONIE-AROME 3D NWP model is used to determine the sensitivity of the forecast of ice accretion on wind turbines to the subgrid forest fraction. Single column simulations of a variety of icing cases at a location in northern Sweden were examined in order to investigate the impact of vegetation cover on ice accretion in varying levels of solar insulation and wind magnitudes. In mid-winter cases, the wind speed response to surface roughness was the primary driver of the vegetation effect on ice accretion. In early season cases, the cloud water response to surface albedo effects plays a secondary role in the impact of in-cloud ice accretion, with the wind response to surface roughness remaining the primary driver for the surface vegetation impact on icing. Two different surface boundary layer (SBL) forest canopy subgrid parameterizations were tested in this study that feature different methods for calculating near-surface proﬁles of wind, temperature, and moisture, with the ice mass accretion again following the wind response to surface vegetation between both of these schemes.


Introduction
In the past several decades, commercial wind power has grown from humble beginnings in small research wind farms to an important energy resource worldwide. Utilities, independent power producers, and investors have especially taken advantage of regions in cold climates such as those in Scandinavia, where sparse human population and high air density combined with terrain-induced wind flows make for an excellent region for wind farm development. According to [1], 69 GW of wind energy were located in cold climates in 2012, with 10 GW/Year being built from 2013-2017.
A challenging aspect of operating wind farms in cold climate regions is that of ice accretion on the wind turbine blades. Accretion of ice affects the aerodynamic efficiency of the turbine, yielding lower power output. Also, increased loads and vibration can decrease the life of turbine components [2] and ice throw from wind turbine blades can be hazardous to plant personnel. The forecast of ice accretion is therefore crucial for plant operators and dispatchers who must plan for wind operations at energy utilities and independent power producers. As a result, ice prediction methods have been developed in the form of short term forecast models [3] and regional ice climate atlases [4]. Although much work has been done in the past decade to develop ice prediction products, this forecast problem continues to be a complex challenge [5]. Icing models, such as that described in [6], have been developed and tested in lab settings, but rely upon prognostic meteorological data input from numerical weather prediction (NWP) models in order to produce predictions of ice load. The ice load forecasts are sensitive to changes in low level liquid water content (LWC), wind magnitude, temperature and humidity [6].
On-shore wind farms in Scandinavia are often located within regions of complex terrain and near heterogeneous patches of boreal forest. In these regions, the typical icing event occurs due to low level clouds containing supercooled liquid water intersecting terrain, impacting the turbine structure [3]. [5] discusses the impact of errors due to localization that can occur in regions with subgrid complex topography and mitigating these errors with an ensemble method. Land cover representation and the modelling of surface-atmosphere fluxes in areas of heterogeneous land surface cover has been a hotly contested topic for decades, with papers focusing on the sensitivity of the forecast of temperature and wind to changes in surface albedo [7] and roughness [8]. Several studies have focused on the impacts of land surface type on the timing and evolution of low level clouds and fog [8][9][10].
The aim of this study is to test the sensitivity of an operational mesoscale NWP model forecast of wind turbine icing events to changes in subgrid representation of vegetation fraction. A single-column model is used in this analysis, which uses the same physics and parameterization packages as its 3D cousin, but integrates in time only vertical fluxes in a vertical column covering one horizontal grid box. While they do not allow the horizontal advection of upstream large scale meteorological variables into and out of the atmospheric column, single column models allow the user to change certain model parameters while holding large scale conditions constant. This has several advantages: first, as subgrid conditions are only integrated in one column the computational expense is cheap relative to 3D models, allowing for one to run several experiments; second, as large scale conditions are held constant, analysis of subgrid responses to changes in parameters is less unwieldy. As a result, we were able to run several simulations with the same initial conditions, only changing the vegetation cover in each. The output from the experiments was then used as input into an icing model. From this, we could analyze the sensitivity of the turbine icing forecast to vegetation cover. In addition, we were able to test the impact of the icing forecast to vegetation fraction when two different surface boundary layer (SBL) parameterization schemes were used: a single-layer scheme and a multilayer scheme with explicit levels within the forest canopy. Section 2 of this paper will describe the models used as well as the experimental design. Section 3 will present the results of the experiments in terms of both the icing model and its input variables, as well as show the response in two specific cases. Section 4 will discuss the results and Section 5 will conclude.

Method
The model used in this study-Modèle Unifié Simple Colonne (MUSC) [11]-is a single column version of the HARMONIE-AROME (cy40h.1.1.1) configuration of the ALADIN-HIRLAM 3D NWP system [12]. HARMONIE-AROME was developed from the AROME mesocale model, itself born as part of the ALADIN (Aire Limitée Adaptation Dynamique Développement International) -HIRLAM (HIgh Resolution Limited Area Model) consortium. ALADIN-HIRLAM is a collaboration of 26 countries in Europe and North Africa with the aim to improve short-range weather prediction in this region. HARMONIE-AROME itself is a convection-permitting non-hydrostatic limited area model used in operational forecast offices in Sweden, Norway, Denmark, Ireland, The Netherlands, Finland, Spain, Iceland, Estonia, and Lithuania. MUSC contains the same physics and parameterization packages as HARMONIE-AROME with the exception that horizontal advection being ignored in the 3 of 19 simulation; although MUSC contains an added capability to implement large scale forcing based on a horizontal thermal gradient. The result is a computationally cheap research tool that can be used to understand vertical turbulent exchanges or test the performance of new parameterization schemes that may be incorporated in new versions of 3D NWP models [13]. MUSC was configured using similar settings to the operational version of HARMONIE-AROME used at SMHI. This configuration uses a 2.5 km horizontal resolution with 65 vertical levels.

Surface Model
The surface model used in MUSC is version 7.3 of the Surface Externalisee (SURFEX) externalized surface scheme [14]. SURFEX is a package of physical parameterizations that calculate the subgrid energy exchange between the surface and atmosphere. When SURFEX is coupled with an atmospheric model, forcing data from each gridpoint of the model atmosphere is sent to a series of models that calculate fluxes of momentum, sensible heat, latent heat, gas-, and aerosol feedbacks for each of four tiles representing a different surface type: sea, lake, nature (consisting of vegetation and soil), and urban areas [14]. In addition to the fraction of surface type tiles, SURFEX has a capability that can further divide the nature tile into up to 12 so-called patches. An average of each individual turbulent flux is then calculated such that it is weighted by the fraction of each surface type and the resulting flux is sent back to the atmosphere, along with information about surface roughness, albedo, emissivity, and surface temperature. In SURFEX version 7.3, the surface type fractions are determined from ECOCLIMAP II [15], which combines satellite and land maps to create a database of surface descriptions at 1-km resolution. The surface description from ECOCLIMAP II, upper air temperature, specific humidity, horizontal wind components, pressure, total precipitation, aerosol concentration, gas concentration, downwelling short-and longwave radiation from the atmospheric model are sent to SURFEX at each model time step.
Available within the SURFEX parameterization package is a 1D multilayer surface boundary layer (SBL) scheme developed by [16] (see schematic in Figure 1b). When the [16] canopy scheme is not activated, SURFEX uses a single layer scheme coupled to the lowest atmospheric model level using a variation of Monin-Obukov similarity to calculate wind, temperature and humidity profiles in the SBL (see schematic in Figure 1a).  [16] proposed that the single layer scheme, while accurate in flat, homogeneous landscapes, failed to estimate mass and momentum fluxes within tree canopies as the lowest atmospheric model level is assumed to be above the vegetation. Multilayer schemes, with explicit atmospheric model levels within the tree canopies, provide the capability to calculate more accurate atmospheric profiles within the tree canopy, but are computationally expensive when used in operational forecast models. [16] proposed a compromise multilayer SBL scheme that solves the atmospheric governing equations by utilizing the subgrid turbulence parameterization scheme combined with a large scale atmospheric term and canopy drag terms calculated from the leaf area index (LAI) of the forest. The turbulence scheme used in MUSC and HARMONIE-AROME is called HARMONIE with RACMO Turbulence (HARATU) [12]. HARATU computes a prognostic TKE equation on half-levels using separate diagnostic length scales for heat and momentum [12], providing the basis for flux computations within SURFEX. The large scale forcing term is the tendency of wind, temperature, and humidity at the lowest atmospheric level. Below the lowest atmospheric level, the aforementioned variables and TKE are calculations taken from the turbulence parameterization scheme used in the NWP system. The aim of the [16] canopy SBL scheme is to provide a computationally cheap solution to improving forecast performance in forested areas. It should be noted that in their analysis of the performance of HARMONIE-AROME (cy40h.1.1.1), [12] states that operational experiences found that the [16] canopy SBL scheme yields temperatures that are too low in stable conditions with weak winds, which is contrary to what was found in [16].

Icing Model
Output from each simulation was run through an icing model based on the one described in [6] and adapted to wind turbine use in [17]. This is the same icing model used by [5]. The model is based upon the equation: where M is the mass of the accumulated ice in kg, t is time, w is the liquid water content, V is the wind speed, D is the diameter of the cylinder, and α 1 , α 2 and α 3 are the collision, sticking, and accretion efficiencies, respectively. L is an ice loss term that combines a series of functions to take ice shedding, melting, sublimation, and wind erosion into account [5]. Equation 1 uses as input from the NWP output wind speed, temperature, pressure, cloud water, ice, rain, graupel, snow, and specific humidity. Collision efficiency α 1 -or the fraction of supercooled particles that make contact with the structure-is determined by the concentration of the particles, size distribution of the particles, wind magnitude, and size of the structure. The size distribution of the supercooled water droplets is estimated using the median volume diameter (MVD) calculation described in [18]. As the specific MVD is unavailable as NWP output, it was estimated using the liquid water content output and assuming that the number concentration of liquid water particles is constant at 100cm −1 following [5]. The sticking efficiency was determined from the work conducted by [19], using α 2 = 1/V 0.75 . The icing model includes functions that calculate ice mass due to different phases of precipitation, including supercooled water, graupel, rain, and snow. The model also includes functions that calculate the approximate decrease of ice mass due to wind erosion, sublimation, and shedding off the turbine blades. The input variables from the NWP model are interpolated within the model to the turbine hub height and the ice mass is calculated at that level.

Experimental Setup
In order to test the sensitivity of the model forecast to the vegetation fraction at the model's lower boundary, MUSC was initialized from a point in northern Sweden using 3D HARMONIE-AROME forecasts of six different low level in-cloud icing events in 2013 and 2014. The aim here was to examine the dynamic and diabatic response of the model to changes of vegetation fraction under a range of differing meteorological conditions and incoming solar radiation profiles. This way, we could examine the response of the icing model to vegetation fraction changes when forced with the NWP output. From this information, one can analyze the icing model input variables to find the contribution of diabatic and dynamic components to the icing model response and the authors could look at the impact of vegetation fraction on the evolution of these events.
The initial conditions were chosen from relatively well-forecasted icing events described in the work done by [5]. Also, the cases chosen represent a wide range of top of atmosphere (TOA) downwelling shortwave (SW) radiation profiles from both early and mid-season icing events in order to test the response due to changes in albedo. The cases chosen for the study are listed in Table 1 For each case, the initial conditions remained the same, with the only adjusted variable being the vegetation percentage. For every case MUSC was initialized from the +03 hour output of the 3D model run, spun up for 3 hours, and integrated for 48 hours, with 60 second timesteps. The simulation was run for 48 hours in order to analyze the consistency of any diabatic impacts over more than one diurnal cycle. Table 2 lists the fractions of the different ground cover types for each experiment. The surface cover from the raw 3D model output featured 86.3% vegetation cover, mostly including coniferous forest native to northern Sweden. In each successive simulation, 10% of bare ground (ECOCLIMAP-II abbreviation SFX.COVER 538) was included in the total surface cover fraction while the vegetation was decreased proportionally. The bare land representation in this study assumes a homogeneous patch of soil with an average raw visible albedo of 0.23 and an infrared albedo of 0.46, according to the ECOCLIMAP-II look-up tables. This bare ground representation was an idealized choice for the experiment and does not necessarily depict the soil conditions at the location the simulation was initialized from. The resulting grid cell albedo for each experiment as a function of vegetation cover is shown in Figure 2. Note that increasing the forest and other vegetation cover decreases the overall albedo in the gridcell. Also note that the albedo values in Figure 2 are from the ECOCLIMAP-II lookup tables and do not include the high albedo related to snow cover, which is present in each of the simulations. As a result, the actual relative difference in albedo from the bare ground experiment is larger than shown. In order to test the sensitivity of the experiments to the type of SBL scheme, each simulation was run with both the single and multilayer SBL schemes described in Section 2.1. In order to determine the subgrid roughness length z 0 to be coupled with the atmosphere, the SURFEX code weights each fraction of cover type i within each grid cell, summing n cover elements using the following: where z 0i is the roughness length for each cover type, Γ i is the weight for each vegetation fraction, and z j is the averaging level. The average subgrid z 0 is then determined using where Γ tot is the sum of all surface cover weights. When the multilayer canopy scheme is used, z j is the height of the lowest canopy level. Otherwise, z j is the height of the lowest atmospheric level. More information about subgrid averaging of surface properties in SURFEX can be found in [14]. When the single layer SBL scheme is used, the roughness length in the SURFEX output represents z 0 aggregated over all surface cover types within the grid cell and is an accurate representation of the surface effects felt by the model atmosphere at its lowest level. In contrast, since the multilayer SBL scheme uses canopy drag terms rather than roughness length to determine frictional effects related to forests, z 0 is only aggregated for bare ground and low vegetation when this SBL scheme is used. Therefore, the SURFEX output for grid cell average z 0 when using the multilayer SBL scheme appears as a low value relative to that of the single layer scheme. As a result, an effective roughness z 0e f f was calculated for both the single layer and multilayer SBL scheme experiments assuming a logarithmic wind profile [20] from the surface to the lowest model level, in which k=0.4 is the von Karman constant, u * is the friction velocity, and u(12) is the wind speed at the lowest model level of 12 m averaged for all cases for each respective SBL scheme. As friction velocity was not explicitly included in the model output, this was calculated from the surface momentum flux using the relationship defined in [20]: where |u w s | is the surface momentum flux calculated in SURFEX and averaged over all cases for each respective SBL scheme. Figure 3 shows z 0e f f as a function of vegetation fraction when using the single layer and multilayer SBL scheme, respectively. Note the difference between z 0e f f as a function of vegetation between the two SBL schemes. While the formulation for surface roughness described in Equation 4 assumes a logarithmic profile-which is not always the case-it is assumed to be a reasonable average approximation for depicting the relative difference between the two SBL schemes. That said, the effective roughness length for the single layer SBL scheme increases steadily with increasing Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 7 July 2020 doi:10.20944/preprints202007.0131.v1 vegetation fraction. This is in contrast to z 0e f f for the multilayer SBL scheme, which increases sharply from 10 to 40% vegetation fraction and only slightly increasing with increasing vegetation fraction above that. The multilayer canopy approach-which includes a forest drag formulation-enables blockage of the lower canopy levels as the amount of tall vegetation increases, something that has previously been shown to lower the roughness length [21]. The impact of the SBL formulation on the icing forecast will be examined in the next sections.

Impacts on Icing Forecast
All simulations were run through the icing model with both the single layer and multilayer SBL schemes used. These results are shown in Figure 4, in which the maximum ice mass is shown for the 48 hour period in which MUSC was run. With the default single layer SBL scheme (Figure 4a), the average icing mass response generally decreases as a function of increased vegetation fraction. This relative decrease is nearly linear, with a change of about -10% per 0.20 increase in vegetation fraction. The winter cases (red, black and blue markers) feature a generally decreasing trend in ice mass as vegetation increases; following the mean trend. The 24 December 2013 case (red markers) follows the mean trend closely and will be further examined in Section 3.3.2. The October cases were a bit more variable, with strong nonlinearity connected to the evolution of cloud cover. Of note is the 18 October 2014 case (cyan markers), which experienced a local maximum in icing mass near 0.60 vegetation fraction and a larger relative impact overall as vegetation fraction is increased. With the single layer SBL scheme, the total ice mass in the case initialized from 18 October 2014 decreased by nearly 70% from the case with bare land to the case with .86 vegetation fraction. The 18 October 2014 case will be further examined in Section 3.3.1.
With the multilayer SBL scheme activated, the percent difference was 25% less than the bare ground control case for vegetation fractions up to 0.35. The midwinter cases display a systematic decrease as vegetation fraction increases. The October cases vary between 0 and 50% less ice mass from bare ground. As vegetation fraction is increased above 0.40, however, the average ice mass impact remains between 25 and 30% less than bare ground. The 24 December 2013 case follows the average trend quite well, with a nearly flat response in ice mass as vegetation fraction is increased above 0.40. As with the experiments with the single layer SBL scheme, an outlier could be seen in the 18 October 2014 simulation, which experienced a near-linear decreasing trend in ice mass as vegetation fraction increased to 0.50, with a well-defined local maximum around 0.80 vegetation fraction. Both the 24 December 2013 and 18 October 2014 cases will be examined in further detail in Section 3.3.

Assessment of Icing Model Input Variables
Recall that the icing model described in Section 2.2 uses input from a NWP model to describe the flux of all phases of subfreezing water through the wind turbine sweep. As the efficiency of the accretion of ice on the wind turbine is dependent upon the content and velocity of the water molecules-as well as air temperature-this section will show the impact of vegetation on the input variables from the MUSC simulations. Here we show the impact of vegetation fraction relative to the bare ground control case for the NWP wind speed, temperature, and liquid water content at 117 meters. The impacts of vegetation cover on rain, snow, ice, and graupel were not included here, as the values of these were negligible throughout each simulation.
In the experiments using the single layer SBL scheme, the mean wind magnitude impact follows the trend in roughness length (see Figure 3a) as vegetation is increased (Figure 5a). With increasing  (Figure 5b) follows a similar pattern seen in the icing model trends for those simulations, with a nearly 35% decrease in wind speeds up to 0.50 vegetation fraction and a flatter response at vegetation fractions higher than this. The mid winter cases (red, black, and blue circles) feature a nearly flat wind magnitude response with vegetation fractions above 0.50, but the slope of the October cases (magenta, cyan, and green circles) is steeper. Note that the winter cases follow the average trend quite closely, while the October cases display relatively more variation.
The average impact of vegetation fraction on temperature ( Figure 6) is similar for both SBL schemes, with a general increasing trend to between 1 and 1.5 degree K. The spread amongst the cases examined is large, however. The winter cases experienced a very small temperature impact, owing to low levels of solar insulation. The temperature impact for the October cases as vegetation is increased varies from 1 K in the 23 October 2014 case to nearly 4 K in the 22 October 2014 case.
In both the single layer and multilayer SBL experiments, cloud water response (Figure 7) is not as systematic as the other variables examined, with the average response less than ± 25%. The cloud water response to vegetation percentage in the winter cases is variable around 0%, with the case initialized on 24

Cases
Several cases were analyzed in order to test the sensitivity of simulated events featuring varying diabatic and large scale influences to the represented forest fraction. To provide further insight into the physical processes of the icing events and the modelling thereof, two cases will be examined in the next sections. The case initialized from the morning of 18 October 2014 featured relatively weak synoptic dynamics under surface high pressure. The 24 December 2013 case featured stronger wind magnitudes owing to a more dynamic synoptic environment. The temporal evolution of meteorological variables for these cases will be examined to highlight the diurnal cycle under different surface albedos and roughness values. The experiments examined in the next section used the single layer SBL scheme. The simulations using the multilayer SBL scheme are not shown here as the evolution of the cases using that configuration were similar, with slightly different magnitudes of meteorological variables as discussed in the previous section.

18 October 2014 case
On the morning of 18 October 2014, a deep surface low pressure system was present in the north Atlantic with an occluded front stretching from southern Sweden into central Europe, where a weak low pressure system was present. In contrast, a region of high pressure was present over the northern Baltic Sea, with weak pressure gradients owing to light winds at the target area in northern Sweden. Figure 8 shows the 900 mb cloud fraction, wind field and temperature at 0900 UTC from ECWMF ERA5 reanalysis data. Surface temperatures were below freezing across the entire region on that day with a subsidence inversion present at 400 meters. Winds were westerly across much of Sweden, with stronger south-southwesterly winds converging near the northern parts of the border with Norway. Low clouds developed over the next 24 hours along the line of convergence, signalling the beginning of a near surface in-cloud icing event at the point of interest symbolized by a black dot on the aforementioned chart. It was from this point of interest that MUSC was initialized. The cloud water evolution for the 18 October 2014 case is greatly impacted by changing the vegetation cover in the single-column experiments. As vegetation cover is increased, the activation of turbine level cloud water is delayed into the early morning hours, with higher values occurring close to the ground (see Figure 9). Focusing on the elevation of 117 meters, Figure 11b shows a progressive trend towards a later onset of cloud water as vegetation increases.
The cloud deck develops soon after the single column model is initialized for the experiments with lower vegetation cover percentages and persists through the daytime hours, but is not present in the higher vegetation cases until later in the simulation. This is likely due to a higher sensible heat flux owing to a greater absorption of downwelling radiation with the lower albedo of higher forest fractions. The temperature at turbine level reflects this (Figure 11a), as temperatures in the experiments with higher vegetation fractions peaked at up to 4 degrees K warmer than the bare ground case on 1500 UTC on 18 October 2014. As wind magnitude ( Figure 10) at 117 meters was similar amongst all experiments in the first hours of the simulation, it is likely that the warmer temperatures alone played a significant role in precluding low level cloud formation for the experiments with higher vegetation fractions. Conversely, under light winds and clear conditions, temperatures at 117 meters are allowed to drop more dramatically overnight in the higher vegetation cases examples after 0000 UTC. With the presence of weak low level winds and very low temperatures during the early morning hours, the formation of nocturnal radiation fog is implied in the simulations with vegetation percentages above 65% coverage. In all cases, these low level clouds lift through the turbine sweep and into the upper boundary layer throughout the simulation. In the cases with higher vegetation fraction, the clouds lift dramatically through the morning hours on 19 October 2014 as low level temperatures rebound. It is noted that during the time that the supercooled clouds lift through the turbine sweep, the wind speeds increase and remain constant through the rest of the simulation. While winds were relatively light in this case, this result is significant as icing is impacted by wind speed as seen in the analysis in the previous sections. The evolution and magnitude of the liquid water content at turbine level ( Figure 11b) varies dramatically as a function of vegetation fraction, with a maximum of 1.9x10 −3 for the simulation using 69.1% vegetation fraction as the clouds lift through turbine level. Recall that the icing response in the 18 October case relative to bare ground was nearly flat at -50% from bare ground between vegetation fractions of 0.50 and 0.70.

24 December 2013 case
On the morning 24 December 2013, a surface low pressure system was present off the coast of Scotland, with an associated occlusion and warm front stretching eastward into southern Norway and Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 7 July 2020 doi:10.20944/preprints202007.0131.v1 12     Sweden. An area of surface low pressure was present in far northern Sweden near the Arctic Circle. Figure 12 shows the regional wind field and cloud cover at 0900 UTC from ECMWF ERA5 reanalysis data. Temperatures near a developing warm front over Baltic Sea were below freezing, but warming close to the freezing line toward the east. Relatively strong northwesterly winds associated with low pressure to the northeast were present in Northern Sweden, with a weak deformation zone near the Norwegian border. Near the line of convergence, a north-south oriented line of clouds formed. A low level icing event occurred at the point of interest over the next few days. Time-height sections of the cloud water content are seen in Figure 13 The wind speed at 117 meters is seen in Figure 14 for each simulation. The initial wind speeds vary from nearly 12 m/s in the bare ground simulations to 6.5 m/s in the simulation with the highest vegetation fraction. The wind speeds for each simulation gradually decrease during the first 24 to 30 hours before the boundary layer decouples, with wind speeds dropping steeply. This decoupling occurs earlier as a function of vegetation fraction. The low level temperatures (Figure 15a) cool to condensation as the boundary layer decouples and fog forms. The cloud water at 117 meters peaks between .08 and .1 g/kg for all vegetation percentages. However, the timing of the onset of the supercooled water varies by about 6 hours, with the onset occurring earlier in the higher vegetation case (Figure 15b). Simulations with higher vegetation fractions relative to the bare ground case led to an earlier development of low level cloud cover. The clouds gradually lift during the day as the relatively lower albedo leads to an increase in upward sensible heat flux in response. In the bare ground cases, clouds formed later, lifting in a similar manner to the higher vegetation cases. Low incoming solar radiation levels led to weaker turbulent mixing than the 18 October 2014 case, allowing fog to form in all cases, even during the daytime hours.
In the 24 December 2013 case, there is less downwelling shortwave radiation due to the low December sun angle at the latitude of the test location and a smaller difference in initial temperature between the experiments (Figure 15a). However, the 117 m wind speed ( Figure 14) is higher in the 24 December 2014 case than in the 18 October 2014 simulation and the sensitivity to changes in vegetation is greater in this case. Recall that the icing response follows the trend in the relative decrease in wind speed from the bare ground case. Also note that the cases with vegetation fractions above 60.3% feature a slightly higher spike in cloud water than the other experiments. Despite the increase in maximum liquid water content in higher vegetation fractions, the ice mass still decreased in this case with the mean wind response.

Discussion
The experiments show that vegetation representation has an impact on both the total accretion and evolution of ice formation during simulated wind farm icing events. The icing model is most sensitive to the impact of the effective roughness on wind magnitude at wind turbine hub height. The sensitivity of the icing model to wind speed is further highlighted when the SBL parameterization is changed to a scheme that aims to include tree canopy effects, in which the icing response also appears to have a strong correlation to surface roughness. With the [16] SBL scheme used, the sensitivity of simulated ice accretion to vegetation cover seems to decrease at higher forest fractions. The sensitivity of simulated ice mass accumulation remains high at lower vegetation percentages, however. This suggests that the simulation of low level icing events may be impacted by not only the uncertainty in vegetation type, but by the surface model, which feature tuning parameters to accurately depict boundary layer processes [16]. Compared with the single-layer SBL scheme, the multilayer SBL scheme gives an effective roughness length response to changing vegetation cover that suggests increasing homogeneity with increased forest canopy cover, which follows the response expected from theory (see [21]). In Section 2.1, [12] was referenced, who discussed that forecasted low level temperatures were caused to be much too low in stable conditions when coupling the HARMONIE-AROME model to SURFEX with the multilayer scheme turned on. This suggests that while the multilayer SBL scheme used in this study gives a surface roughness response that follows observations, proper tuning is required for the best results. Recent configurations of the HARMONIE-AROME model have moved away from the multilayer canopy scheme in favour of calculating surface fluxes over two separate patches of nature types within a grid cell before determining the weighted average flux over all other surface types. [22] suggested decreasing the height of the lowest atmospheric model level from 10 to 5 meters in the related AROME model to explicitly resolve atmospheric layers within the tree canopy; they found that this precluded the need for the [16] multilayer scheme. Diabatic effects due to surface albedo differences which impact the diurnal evolution of cloud cover, also impacted the icing cases Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 7 July 2020 doi:10.20944/preprints202007.0131.v1 12     that were examined. This is evident in the 18 October 2014 case, which featured greater solar insulation than the other mid-winter cases. In this case, the impact of vegetation on the total ice mass was not as systematic as the mid-winter cases, with the cloud evolution impacted by the convective response to changes in surface albedo. This shows that, while the dynamic effects due to surface roughness remain dominant, the albedo effect has a highly relevant impact on the timing and evolution of the icing event. Recall that in Section 2.3 the model surface is covered by snow, which enhanced the albedo of the bare ground contribution to each experiment. In a test without snow cover (not shown), the impact of vegetation on the cloud evolution was reduced, suggesting a decreased heat flux response with vegetation changes. This is not surprising, as the impact of the representation of snow cover in the vicinity of forests on the surface energy balance has been examined for decades (see [7]). As the icing model seems to be most sensitive to changes in wind speed with increasing forest fraction, wind production can be impacted by the uncertainty of land cover representation and SBL parameterization in addition to the icing forecast. The ECOCLIMAP-II package [23] used in SURFEX is a 1-km resolution land surface database created from several remote sensing datasets of higher resolution. While land cover representation datasets of this resolution have been successfully used in global climate models and NWP models of relatively coarse resolution, some further consideration may be needed in the future as icing and wind production forecasts use NWP input from finer scale models. As one moves to higher resolution, one also must consider the heterogeneity, or patchiness, of the forest represented in the model forecast. The land cover represented in the experiments described in this study was assumed to be homogeneously distributed inside the grid cell. Several large eddy simulation (LES) studies (for example, [8,13,[24][25][26]) have explored the varying impacts of land surface heterogeneties on boundary layer cloud processes. As forests surrounding onshore wind farms in Scandinavia often feature patchy heterogeneities due to nearby forestry operations and clearings for turbine equipment, and extension of this work is to examine how this patchiness of vegetation can impact wind turbine icing processes. With a properly-tuned SBL scheme, future mesoscale models may take advantage of the latest generations of surface description maps, such as the objective roughness approach (ORA) described in [27] that uses airborne lidar scans of tree height and density to create detailed snapshots of surface roughness.

Conclusions
A single column model was used to test the sensitivity of ice accretion predictions at a high latitude location to subgrid land use representation from the meteorological input into the icing model. The single column model demonstrates that the icing on wind power turbines is sensitive to the forest fraction in a gridcell. In the mid-winter cases examined in which the downwelling solar radiation was relatively low, the icing forecast seems to be most sensitive to the wind response due to surface roughness changes as forest fraction is increased.
In the early season icing case examined with increased diabatic forcing in the form of increased solar insulation, the trend of total predicted ice mass was less systematic than the mid-season cases, suggesting a sensitivity to the evolution of the in-cloud icing event. In this case, the cloud formation in the single column model was more transient and this impact can be seen in the ice mass predictions as vegetation increases. However, the general trend still appears sensitive to the wind response when surface roughness is increased.
When simulating the cases using a more sophisticated SBL scheme in the single column model with a different calculation of surface roughness, one can see a similar response to wind magnitude in the icing prediction. While this suggests that uncertainty in the represented subgrid forest fraction can lead to errors in the wind turbine icing prediction, the surface layer model parameterization can also lead to uncertainty.