Area-Averaged Surface Moisture Flux over Fragmented Sea Ice: Floe Size Distribution Effects and the Associated Convection Structure within the Atmospheric Boundary Layer

: Sea ice fragmentation results in the transformation of the surface from relatively homogeneous to highly heterogeneous. Atmospheric boundary layer (ABL) rapidly responds to those changes through a range of processes which are poorly understood and not parametrized in numerical weather prediction (NWP) models. The aim of this work is to increase our understanding and develop parametrization of the ABL response to different ﬂoe size distributions (FSD). The analysis is based on the results of simulations with the Weather Research and Forecasting model. Results show that FSD determines the distribution and intensity of convection within the ABL through its inﬂuence on the atmospheric circulation. Substantial differences between various FSDs are found in the analysis of spatial arrangement and strength of ABL convection. To incorporate those sub-grid effects in the NWP models, a correction factor for the calculation of surface moisture heat ﬂux is developed. It is expressed as a function of ﬂoe size, sea ice concentration and wind speed, and enables a correction of the ﬂux computed from area-averaged quantities, as is typically done in NWP models. In general, the presented study sheds some more light on the sea ice–atmosphere interactions and provides the ﬁrst attempt to parametrize the inﬂuence of FSD on the ABL.


Introduction
Parametrization of sub-mesoscale processes related to spatially heterogeneous land/water surface remains one of the key challenges in atmosphere and climate modeling. In polar regions, sea ice surface properties vary over a wide range of scales, from centimeters to tens of kilometers, influencing the oceanic and atmospheric boundary layers through many complex interactions. The quality of weather and climate predictions depends heavily on our ability to accurately represent the effects of those interactions in numerical weather prediction (NWP) models. As most of currently used climate models have horizontal grid spacing of tens of kilometers and are therefore unable to directly reproduce processes related to smaller-scale sea ice features, their effects must be parametrized. However, the scarcity of observations of atmospheric properties over fragmented sea ice and the resulting incomplete understanding of the sub-mesoscale processes hampers the development of parametrizations. Consequently, sub-grid interactions between the surface and the atmosphere, which may account for the inaccuracies present in model predictions [1,2], are poorly represented in NWP models.
The influence of sea ice heterogeneity on the atmospheric boundary layer (ABL) has long been the subject of both modeling [3][4][5][6] and observational studies [7][8][9][10]. The general aspects of ABL's response to the presence of open water areas in the form of singular or several leads in an otherwise continuous ice cover are relatively well studied. These effects include changes of the ABL stability over and downstream of open leads [5,8], sensitivity to prevailing wind conditions [10,11], formation of steam fog and cloud plumes [12], or initiation of gravity waves [4]. Similarly, the influence of leads width on the strength of the atmospheric response has been studied as well [13,14], and it is generally assumed that for a given total surface area of leads, numerous narrow leads tend to induce stronger area-averaged turbulent exchange of heat and moisture than a few wide ones. However, recent satellite data studies indicate that such conclusion might be flawed [15] and suggest that the larger contribution from small leads to the total surface heat flux may not be associated with their efficiency, but rather number density, length and frequency of their occurrence in the Arctic. Other studies point out that not only the average lead width is important, but also their spatial arrangement and the shape of their size distribution [16]. According to modeling results [17,18], even very small open water areas might result in the development of small convective plumes which penetrate the ABL and modify its circulation and properties. The areas of updraft air motion can be recognized in the winter Arctic from the satellites due to the presence of steam fog. Steam fog occurs because of cold (−40 • C to −5 • C) air advection over water which is only slightly warmer than its freezing point. The fog is characterized by either columnar or banded structure and extends in height from 1 to 1500 meters with the liquid water content in the range of 0.01-0.5 g/m 3 [19]. This fog type can be recognized in satellite images as linear features, organized into bands of different width, interspersed with clear air areas [20], or as plume-like streaks, with a typical spacing of a few kilometers [21]. Numerical studies indicate that such features may have significant impact on atmospheric processes not only locally, but also regionally [5,18,22].
Obviously, large scale models do not directly resolve heterogeneity at scales smaller than the model grid cell. To account for the larger-scale effects of sub-grid properties of the atmosphere and the underlying surface, NWP models must rely on estimation methods and simplified parametrizations. Computation of surface turbulent heat fluxes (STHF)-the main subject of the present paper-provides a very good example of problems associated with constructing such parametrizations. Locally, STHF is a function of several variables characterizing the surface and the lower atmosphere (temperature, wind speed, etc.). Due to the nonlinear character of the relationship between those variables and spatial correlations between them, the area-averaged STHF cannot be directly computed from the area-averaged variables, available in numerical models. The necessary correction depends on the characteristics of the surface (size and spatial arrangement of "patches" of different surface types; in the case of ice-covered sea surface analyzed here: spatial distribution of sea ice and open water) and the associated local response of the lower atmosphere. Several approaches to this problem have been developed and are used in various NWP models.
In the standard, most basic approach, a single, dominant type of the surface in a grid cell is used as representative for that cell. Thus, the cell-averaged surface fluxes are calculated from cell-averaged atmospheric and surface properties, like, e.g., in the standard version of COSMO-CLM or CCLm model [23]. Such averaging of atmospheric and surface properties is termed an "aggregation" process. An important disadvantage of this approach is a possibility that the aggregation will produce biased results due to the above-mentioned nonlinear relationship between the variables involved. The problem of biased STHF estimates has long been recognized [14,[24][25][26] and several improved methods of flux calculation have been proposed. One of the most popular ones is a so-called mosaic method developed by Avissar and Pielke [27]. In this approach, the sub-grid variability is accounted for by independently coupling each surface type within a grid cell to the overlying atmosphere. Surface-specific fluxes are thus calculated from area-averaged atmospheric properties (as previously), but with the consideration of the parameters specific for each surface type. Grid cell-averaged fluxes are then obtained by summing the fluxes from different surfaces, weighted by the fractional area they cover. This method has been adopted, sometimes with slight modifications, in various modeling [26,28,29] and observational [30] studies. The mosaic approach, to some extent, includes the effects of surface heterogeneity, but due to very complex relationship between the ABL and the underlying surface, it represents only a rough estimation of surface-atmosphere interactions, completely disregarding spatial structure of those interactions. As we demonstrate in this work, in the case of ABL over fragmented sea ice it is that spatial structure that is responsible for large deviations between area-averaged fluxes computed with different methods. It is very likely that in the case of winter sea ice fragmentation, the averaging of atmospheric variables is behind the largest errors in NWP predictions of STHF. This is because differences between local ABL properties over sea ice and open water are often very large (e.g., can reach even tens of degrees in the case of air temperature).
The methods of STHF calculation described above rely solely on sea ice concentration in a grid cell and do not take into account any information about spatial arrangement of ice and water within that cell (size of ice floes, presence and orientation of leads and fractures, etc.). This limitation seems particularly serious in view of the present climate change in the polar regions, associated with a rapid transition of the sea ice cover, with decreasing extent of multiyear sea ice and expanding proportion of seasonal ice [31]. Due to the fact that seasonal, first-year sea ice is more prone to breaking and deformation due to external forcing, heterogeneous sea ice is going to appear more frequently [32]. Since the sea ice cover moderates energy transfer between ocean and atmosphere and is highly sensitive to weather patterns changes, the atmospheric response to sea ice fragmentation becomes a particularly important research topic. A specific aspect related to sea ice fragmentation that attracted a lot of scientific attention throughout the years is the floe size distribution (FSD), defined as the number of floes in different size categories in a given region, divided by the area of that region [33]. It is a key parameter in describing the state of the marginal ice zone (MIZ) and the evolution of sea ice cover. The information about floe size is essential to determine the rate of lateral sea ice melt [34], to assess the sea ice dynamics and internal stresses [35], and to resolve sea ice interactions with oceanic [36] and atmospheric boundary layers [17]. The importance of FSD drives the effort to include its representation into sea ice and climate models. Equations describing the evolution of the joint distribution of sea ice thickness and floe size were formulated and implemented in a continuum sea ice model by Horvat and Tziperman [37]. The advantage of this model, further developed by Roach et al. [38], is that it does not involve any assumptions about the form of that distribution, so that it is allowed to evolve in response to forcing from the ocean and atmosphere. Another approach, alternative to continuum models, is based on discrete-element methods (DEM), as, e.g., by Herman [39] whose joint-particle DEM allows direct simulations of floe formation due to dynamic processes. However, due to extremely high computational costs of DEMs, larger-scale sea ice simulations must be based on continuum models (coupled with atmospheric, oceanic and possibly other components). The new possibility of including FSD in those models, described above, provides a wide range of possibilities in terms of implementing new physical processes, both influencing and influenced by the FSD. On the other hand, however, this can be done only if those processes are sufficiently well understood. In the case of the response of the ABL to sea ice fragmentation, this is at present not the case. Moreover, observations that could contribute to improving our knowledge are very rare if not non-existent. Therefore, in order to implement FSD-related effects into NWP models we should first focus on better understanding of its impact on the atmospheric and oceanic boundary layers and the physical processes involved, in parallel with the development of parametrizations that would make the models less expensive computationally.
Several aspects of the ABL response to sea ice heterogeneity have been studied thoroughly in the recent numerical study of Wenta and Herman [17]. They examined the influence of small-scale sea ice surface variability (different floe size distributions or different lead widths and orientation) on domain-averaged STHF and properties of the ABL (e.g., the total moisture content). In agreement with earlier studies, mentioned above, the model predicted formation of convective plumes, rapid release of heat and moisture, strong local wind speed increase and intensification of vertical air motions. Crucially for the present analysis, at the same ice concentration, different size and spatial arrangement of ice floes resulted in changes of the domain-averaged ABL properties. This indicates that the above-described methods of computing STHF, based solely on sea ice concentration, might produce biased results, with inaccuracy dependent on the properties of FSD. Furthermore, significant correlation is observed not only between the local surface and atmospheric variables, but also within the ABL itself (i.e., between ABL properties at different heights), indicating that the transfer coefficients used in the STHF algorithms should take into account local variations of temperature and wind speed.
In this study, the influence of FSD on the ABL is studied further. The modeling results obtained by Wenta and Herman [17] are examined, with the focus on the convective structures within the ABL and the surface turbulent moisture heat flux (the second component of STHF, the sensible heat flux, will be analyzed in a subsequent study). The goal of the first part of the analysis (Section 3.1) is to explain how the spatial arrangement and size of convective structures depends on the floe size distribution, and to describe the properties of those structures found in the modeling results. (Anticipating Section 3.1, it is important to stress here that convection cells in this case do not have characteristic, regular shapes typical for Rayleigh-Benard convection, and therefore they are referred to as convective structures rather than convective cells in the rest of this paper.) The convection analysis is also used to explain the variability of domain-averaged STHF and total cloud liquid water content found in Wenta and Herman [17].
The goal of the research described in the second part of the paper (Section 3.2) is to develop a formula that would allow incorporating the FSD-related effects in the calculation of surface turbulent moisture heat flux (further referred as surface moisture heat flux, or SMHF). More specifically, the goal is to propose a correction factor, dependent on floe size, which might be used in combination with the existing, ice-concentration-based flux formulae in order to improve their accuracy over fragmented sea ice. The proposed correction takes a form of a coefficient, denoted with α, which describes the ratio between the surface moisture flux F m calculated with two different methods: α = F m,1 /F m,2 . In the first method, the effects of FSD are properly accounted for as the surface and atmospheric quantities on the sub-grid scale are taken directly from high-resolution numerical simulations. In the second method, the flux is determined with one of the algorithms described earlier, i.e., it is based on area-averaged properties of the surface and the atmosphere. Obviously, the computation of F m,2 depends on the availability of surface and atmospheric properties present in a given atmospheric model. As described in detail in Section 3.2, α can be formulated as a function of area-averaged wind speed, ice concentration and floe size, so that, for given ice and wind conditions, F m,2 can be computed in a traditional way and corrected for FSD effects by multiplying it with α.
It must be stressed that the proposed correction is based on idealized modeling results and must be validated with observations and further modeling in the future. Therefore, it can be understood as the first step towards a parametrization suitable for implementation in NWP models.

Configuration of the WRF Model
Results from the modeling study of Wenta and Herman [17] are used in the analysis in the present study. This section provides a brief description of the model configuration and initial conditions; detailed information can be found in [17].
The Advanced Research Weather Research and Forecasting (ARW WRF) model is used in an "idealized mode". In this approach, the model is initialized with a single sounding file consisting of the vertical profiles of potential temperature, water vapor mixing ratio, and the two components of wind velocity. The full air pressure and density, as well as other atmospheric variables, are initialized from that input sounding. The model domain is rectangular, with periodic boundaries in both horizontal directions, dimensions 200 × 200 grid points and horizontal resolution of 100 m, thus covering the area of 20,000 × 20,000 m. The top boundary is set at 2000 m above the surface, close to the top of the inversion layer separating the ABL from the free atmosphere in the initial sounding. At the upper model boundary, a damping zone is defined as a Rayleigh relaxation layer of 200 m depth, with a default damping coefficient of 0.002 s −1 . The vertical velocity component at the top of the model domain is set to zero. Air column is divided into 61 vertical levels with exponential thickness distribution, from ∼2 m at the surface to ∼200 m at the top. Each simulation used in this analysis was performed for 14 h, with a time step of 1 s. The results were stored every 10 minutes.
The parametrization schemes of physical processes along with diffusion options are chosen based on the WRF model validation performed in other studies under similar conditions (Table A1). A detailed description of selected settings together with the justification of their choice can be found in [17].
The model box is placed at 75 • N, 160 • E, i.e., at the location of the Surface Heat Budget of The Arctic Ocean (SHEBA) experiment [40]. The Coriolis parameter, corresponding to that latitude, is assumed constant and equal to 1.417 × 10 −4 s −1 . We assume, for simplicity, that at the start of each simulation the wind blows from west to east. For each sea ice mask (see below), the simulations are run 6 times: with wind speed equal to zero, and for 5 different wind speed profiles obtained by modifying the profile measured during SHEBA on 23rd of February 1998 (see [17] for details). The resulting vertically averaged values of wind speed for initial profiles 1-5 are as follows: 1.06 m/s, 2.1 m/s, 4.26 m/s, 7.4 m/s and 9.35 m/s, respectively.
The surface water temperature is set to freezing point, assumed equal to 271.4 K. The initial temperature within sea ice varies linearly from 271.4 K at the bottom to the air temperature at the upper ice surface. Default roughness length values are used, equal to 1 × 10 −4 m over water and 1 × 10 −3 m over ice.
In the present study, six different sea ice concentration values are considered: c = 50%, c = 60%, c = 70%, c = 80%, c = 85% and c = 90%. For each ice concentration and wind speed profile, the model is launched for a series of simulations with different spatial arrangements of ice floes. Sea ice is represented as round floes with a power-law probability distribution of their radii P(r) ∼ r −β , with an exponent β = 1.8 (values between 1.5 and 2.0 are typically observed in the marginal ice zone). Floe radii range from several tens of meters to over 4 km. The sea ice masks for WRF are obtained by first generating an ensemble of floes with a desired total ice surface area, corresponding to a prescribed c and number of floes N f ; subsequently, these floes are randomly placed within a temporary, very large rectangular domain at a very low ice concentration and without floe overlap, and a DEM model [39] is used to reduce the size of that domain and to obtain the final floe positions, from which eventually an ice concentration map is computed. In other words, the DEM is used to converge the initially loosely packed floe assemblage to obtain a desired ice concentration-which is otherwise a rather nontrivial task. A summary of all floe configurations considered is given in Wenta and Herman [17].
To summarize, in the analysis presented below we consider simulations which are initialized with five different floe-size distributions, six wind profiles and six sea ice concentrations, giving a total of 180 combinations. Different FSDs are distinguished in the first part of the analysis by N f , but in the calculation of the α coefficient a median of floe radius is calculated for every FSD in order to make the algorithm applicable in other studies (obviously, N f is domain-size dependent).

Calculation of Spatial Coverage of Convective Structures
The study of convective structures within the model domain is based on the analysis of the vertical component of the air motion. The highest values of turbulent fluxes, water vapor and cloud liquid content within the atmosphere are associated with positive (upward) values of vertical wind component, which in the present context serve as an indicator of convective motion. In the first step of the analysis, gray-colored maps of upward air motion areas are generated for every output time step of a given simulation (Figure 1a). In the next step, the program called Pixie is used to determine the size of updraft areas in every produced image. Pixie was created specifically for the presented research and is available for free on GitLab (https://gitlab.com/seadata-software/pixie). It is a Python script that applies simple methods of image recognition to distinguish the areas of upward air motion. It processes every image and counts the pixels for a given color intensity threshold from 0 (black) to 255 (white). For the presented study, a 140 threshold has been determined as the most suitable and applied for the whole analysis (Figure 1b). Based on the computed number of classified pixels and the number of pixels present in the whole model domain, the total updraft area in km 2 and their fractional coverage is computed.

Computation of Area-Averaged Surface Moisture Heat Flux
In the WRF model, the latent turbulent heat flux is obtained as a product of the moisture heat flux and the specific latent heat for vaporization of water; therefore, our study of surface moisture flux and the method of its adjustment due to sea ice fragmentation can also be applied to the latent heat flux.
The surface variables in the WRF setup described in Section 2.1 are calculated using the parametrization by Janjić [41], Janić [42]. In this high-resolution WRF application, where the planetary boundary layer schemes are turned off, the upward moisture flux F m (kg m −2 s −1 ) at the surface is given as: where ρ a (kg/m 3 ) is air density, C h (m/s) denotes the surface exchange coefficient for heat, q a (kg/kg) is water vapor mixing ratio at the lowest layer of the atmosphere, and q z 0 (kg/kg) denotes specific humidity at roughness length. The surface turbulent exchange coefficient for heat is computed using Monin-Obukhov similarity theory [43]: where k is von Kármán constant, u * (m/s) is a friction velocity, z (m) is the height of the lowest model layer, z 0m (m) is defined as a height at and below which the mean wind speed becomes zero, and z 0t (m) denotes the height below which vertical transfer is through molecular diffusion and above which mixing by air currents dominates. L denotes the Monin-Obukhov length and ψ m and ψ h are integrated similarity functions for momentum and heat, respectively. Another surface related parameter q z 0 is obtained from: with z q (m) symbolizing moisture roughness length, q s (kg/kg) specific humidity at lower boundary, and α t (m 2 /s) the thermal diffusivity of air.
In the following analysis, the whole model domain of 20,000 × 20,000 m is treated as a single grid cell of an NWP model. Sub-grid variations of the relevant variables are obtained directly from the 40,000 grid points of 100 × 100 m dimensions, which together form the modeled area. As mentioned in the introduction, the correction coefficient α is defined as the ratio between moisture heat fluxes F m,1 and F m,2 , which differ in terms of the order of averaging of surface and atmospheric properties: To obtain the value of F m,1 , the moisture heat flux from Equation (1) is calculated for every grid point of the model domain. In the next step, it is area-averaged and it thus represents the "true" area-averaged moisture flux: In the calculation of F m,2 , averaging takes place at the level of flux equations components, i.e., the values of ρ a , C h , q a , and q z 0 from every grid point are area-averaged over the model domain, and then used as input in the flux computation: This approach corresponds to the aggregation method described in the introduction. Another method of flux calculations, described as "mosaic" in the introduction, is also tested to determine whether similar patterns and relations are present. When this method is used in NWP models, the atmospheric properties are averaged over a grid cell area, but the surface quantities of q z 0 and C h are computed by one-dimensional (1D) SVAT model driven by grid-scale atmospheric quantities [26]. In the present study, the sub-grid properties of the surface are taken directly from the high-resolution model results and averaged separately over sea ice and water (e.g., C h,water , q z0,ice , etc.). Therefore, the flux is calculated from the atmospheric quantities which are averaged over the whole model domain (as in the method above), and the surface quantities which are averaged over each surface type, and weighted with ice concentration: F m,3 = (1 − c)ρ a C h,water (q a − q z 0 ,water ) + cρ a C h,ice (q a − q z 0 ,ice ).
The corresponding correction coefficient in this case, analogous to Equation (4), is: In all cases analyzed, the ratio between F m1 and F m3 is slightly smaller (by 0.1-0.3) than that between F m1 and F m2 , but overall α and γ follow the same pattern. It is found that the averaging of atmospheric properties is responsible for the largest errors in the calculation of moisture heat flux. Considering that this averaging is done in both methods, the further analysis is limited to the coefficient α, due to the straightforwardness and simplicity of computation of the term in its denominator (F m,2 ).
More than 180 values of α are obtained from simulations with 5 different FSDs, 6 initial wind speed profiles and 6 sea ice concentrations along with several additional tests with different initial conditions. As described in Section 2.1, the FSD is represented by the median floe radius in the model domain (r), and the spatial and temporal mean of wind speed U m is used to characterize wind conditions. For every combination of initial conditions, α is calculated for every output time step of a 14-h long simulation, and the median of those values is considered further. The aim of the analysis is to express α in terms of sea ice concentration, FSD and wind speed: α = α(c, r, U m ).

Horizontal and Vertical Structure of Convection
Convection found in our modeling results is driven by thermal instability, which develops when an open water area comes in contact with sufficiently cool air above. Open water areas of various sizes, scattered throughout the model domain, are locations of rapid exchange of heat and moisture. Convective plumes develop between sea ice floes and transform the atmospheric circulation in the whole ABL. Intensive turbulent exchange over open water rapidly warms the lowest atmospheric layer and produces negative pressure anomalies. The resulting pressure differences between sea ice floes and water produce a breeze-like circulation. In the weak wind conditions (wind profiles No. 1-3 along with the highest water vapor and cloud liquid water content. As in every convection, all updrafts are balanced by downdrafts, which develop just beside the uprising plumes. Therefore, convection found in our results does not take a typical cell-like form and is referred to as a convective structure, rather than a cell. In comparison to the updrafts, the downdrafts are much weaker, but occupy slightly more space. In some cases, a descending air motion occurs on the verge of the sea ice floes (Figure 3). Slightly negative values of the fluxes, along with downward vertical air motions are also found high above the surface in the areas where plumes penetrate the inversion layer, generating turbulence and, in consequence, negative entertainment flux.
Convective structures occur next to each other, creating narrow lines between sea ice floes (  (Figure 2). In the presented results the parallel alignment of those lines is slightly disturbed due to the presence of large sea ice floes. A resemblance to roll-like features, commonly forming in the shallow boundary layer [44], can be found in such strong wind situations.  In our simulations the number N f (and thus the size) of the floes controls the spatial arrangement of open water areas within the model domain. To identify the links between FSD and convection, we investigate the differences between spatial arrangement and intensity of convective structures for sea ice maps with different N f . Simulations with sea ice maps of 50, 100, 500, 1000 and 5000 floes generated for sea ice concentrations of c = 50%, c = 70% and c = 90% are considered. Similar patterns are found for different FSDs for all sea ice concentrations. The most intensive updrafts (reaching the height of 1000 m) occur in the regions of largest open water areas between sea ice floes, which are more frequent in sea ice maps with N f = 50 or N f = 100, as spaces between bigger floes, although infrequent, are generally larger. Furthermore, in those simulations the breeze-like circulation encompasses wider areas and the path of air over open water regions is longer, so that the resulting updraft is more intensive and penetrates deeper into the ABL. The linear structure of convection between the floes is also more distinctive, whereas with increasing number of sea ice floes open water areas become smaller and more scattered, sometimes covering only a few pixels between the smallest sea ice floes. The distinguishing feature of those simulations is the occurrence on convective updraft even in a very narrow areas without sea ice. However, the smaller are the individual open water areas, the less intensive is the convection as the breeze-like flow from sea ice toward the water occurs over shorter distance. Therefore, as the floe-size distribution determines the arrangement of open water areas, it also directly controls the distribution and intensity of convective updrafts in our results (Figure 4b).
In the next step, the total area covered by the convective updrafts in every simulation with varying FSD and sea ice concentration is examined. In general, the smaller the N f , the smaller is the total area of convection. This result complements the conclusions formulated above that larger number of the sea ice floes (N f = 1000 and N f = 5000) results in numerous, very small areas of open water with weak convective plumes developing over every one of them, while in simulations with N f = 50 or N f = 100 the updrafts are well developed, more intensive and form narrow lines between the sea ice floes thus covering smaller area in total (Figure 4a).
The described results provide an explanation for the differences between area-averaged values of turbulent heat fluxes and the total liquid water content between the simulations with various N f found in Wenta and Herman [17]. Stronger exchange of heat and moisture is associated with the most intensive updrafts, characterized by particularly large values of surface fluxes, liquid water and water vapor content. Such exceptional convective structures occur in simulations with smaller N f and thus larger sea ice floes. Therefore, although in total the area of convective updraft is in general larger in the simulations with N f = 1000 and N f = 5000, the area-averaged turbulent fluxes from those simulations will be smaller in comparison to values obtained for sea ice maps with N f = 50 and N f = 100 ( Figure 4). Furthermore, as there are more updrafts in the simulations with N f = 1000 and N f = 5000, there are also more (weaker) downdrafts. In summary, the described variability of convective patterns are well correlated with the variability of the area-averaged turbulent fluxes and total moisture content in different model runs observed in [17] (Figure 4c,d).

Correction Coefficient for Surface Moisture Flux
In the preceding study of Wenta and Herman [17] it was found that the basic variables (e.g., temperature, specific humidity) describing the state of the ABL over sea ice are significantly spatially correlated due to the specific three-dimensional structure of convection. The surface turbulent heat fluxes are nonlinear functions of those variables. Assuming that our model domain represents a single grid cell of an NWP model, the characteristic scales of aforementioned convective structures are smaller than the mesh size. Therefore, when only grid-scale atmospheric properties are available and local covariance between them is unknown, assessing the strength of turbulent heat exchange is biased and problematic. In this section, we develop a correction coefficient α, described in Section 2.3, for the calculation of moisture heat flux in the situations where the sub-grid variability of surface fluxes is influenced by the FSD.
As described earlier, cold atmospheric conditions present at the beginning of our simulations (251 K) and relatively warm water surface between the sea ice floes (271.4 K) generate rapid and substantial moisture exchange. Water vapor plumes develop over the areas of open water where the updraft air motions are strongest, and the local surface turbulent moisture heat flux reaches maximum. In general, relatively large values of STHF are found over open water areas, with maximum in the regions of convective updrafts, while over sea ice the STHF is close to zero or marginally negative. Therefore, as the FSD determines the spatial arrangement and intensity of convective updrafts, it also controls the distribution and intensity of STHF.
As described in Section 2.3, the correction coefficient α is a function of the mean wind speed U m , median floe size r, and sea ice concentration c. The set of those three variables defines each WRF simulation, for which a single, time-mean value of α is obtained. At first, the coefficients are examined for every sea ice concentration used in the simulations ( Figure 5). It is evident that the higher the sea ice concentration, the higher is the value of α. Hence, the lower is the sea ice concentration, the less erroneous is the averaging of surface and atmospheric properties, whereas, when sea ice dominates in the model domain (c ≥ 80%) the aggregation process will result in disregarding effects of small convective updrafts between each individual floe. Furthermore, for wind speeds higher than ∼6 m/s the value of α approaches one, indicating that the influence of floe size becomes irrelevant. It is a consequence of more intensive mixing of the air in the ABL and homogenization of its properties over sea ice floes and water. As far as the role of FSD is concerned, the values of α are slightly higher for smaller sea ice floes than for bigger floes. It seems reasonable to assume that with decreasing size of the floes, α → 1, i.e., the conditions approach those with a uniform, constant sea ice concentration everywhere in the model domain.
In general, as can be seen in Figure 5, while the significance of mean wind speed and floe-size changes with changing sea ice concentrations, the overall shape of the α(U m , r) relationship remains similar. The higher the sea ice concentration, the lower the mean wind speed and the bigger sea ice floes are present in the model domain, the higher is the value of α. Therefore, we presume that α can be expressed by a single equation which includes the influence of wind speed, sea ice concentration and floe size.
After a comprehensive analysis, the objective of which was to propose a single formula suitable for a wide range of conditions, the following function has been selected: where the coefficients a and b depend on ice concentration and d is a constant: For a given ice concentration, it exhibits the desired behavior in the limiting cases of large wind speeds (α → 1 for U m → ∞), very small ice floes (α → 1 for r → 0) and large floes (α becomes floe-size independent for r → ∞). Considering that the accuracy of each particular fit and the general shape of the fitted surface are similar for all sea ice concentrations, only three examples are presented in Figure 6. Elementary statistics, i.e., the correlation coefficient (CC) and the root mean squared error (RMSE) describing the goodness of fit are shown in Table 1. In general, very high values of CC are obtained for every ice concentration, while the values of RMSE change slightly with c. The highest RMSE is found for c = 50%, which corresponds to generally small influence of floe size on the α coefficient. Taking into account the nonlinear relationship between the analyzed variables, the obtained high values of CC must be treated with caution. In general, while a high CC is a necessary condition for a fit to be regarded as satisfactory, it is not a sufficient measure of the goodness of fit and has to be completed by other statistics and the analysis of residuals [45]. As can be seen in Figure 6, the algorithm tends to slightly overestimate the values of α in the range of higher wind speeds (where those values are close to 1), but in general, the random scatter of differences between the fitted and the original values of α, and small values of those differences, indicate that the proposed approximation produces reasonable results over the range of wind and sea ice conditions considered. Given the modeling results, the derived equation can be used to compute the coefficient α and to correct the value of the moisture flux derived from the area-averaged surface and atmospheric properties.

Discussion
We have shown in this study that the floe-size distribution plays an important role in sea ice-atmosphere interactions in the polar regions. It controls the spatial arrangement of sea ice and open water areas together with the associated convective structures in the ABL. Convective updrafts occur in every one of our simulations, but their intensity differs with varying FSD due to the changing extent and strength of breeze-like circulation. Our results are similar to the study of Esau [22], who showed that secondary circulations like the ones found in our modeling results might be of large importance in the studies of ABL convection. Furthermore, his results do not entirely agree with the assumption that narrow leads produce more intensive fluxes, claiming that the associated processes are much more complicated and still not fully understood. Furthermore, in their study of atmospheric response to various configurations of polynya, Dare and Atkinson [16] point out that the spatial arrangement of polynyas influences the magnitude and spatial distribution of vertical heat transfer in the atmosphere. The spaces between the floes in our results vary from a few meters to several kilometers, but in order to fully understand the atmospheric processes above them we also have to examine the size of adjacent floes and the extent of the associated breeze-like circulations. Therefore, our results agree with the suggestions of [16,22] that, when studying the atmospheric response to sea ice fragmentation, the effects related to the spatial arrangement and size distribution of leads and floes must be considered.
Large local differences in atmospheric properties over ice with different floe-size distributions, associated with convective structures, are disregarded when the atmospheric quantities like air temperature or humidity are averaged over the whole model domain. Furthermore, Wenta and Herman [17] pointed out that even when all other initial conditions remain identical, different results of area-averaged fluxes are obtained for the same ice concentration, but different FSDs-which further indicates that the lack of FSD information in NWP models may be at least partially responsible for the observed errors in the simulations over polar regions [1,2]. Based on our WRF modeling results, we propose a simple correction method in the form of an α coefficient that might be used in computation of the surface moisture flux from area-averaged surface and atmospheric variables. The value of α depends on the mean wind speed, sea ice concentration and the median floe radius from particular FSD. In view of the fact that models do not include FSD or the size of the floes in particular grid cell and that presented results have not been validated with observational data, the proposed coefficient is not yet applicable for the GCMs. However, it indicates that the processes associated with the FSD and atmosphere interactions can be parametrized in a simple, computationally efficient way. With the development of models like Horvat et al. [36], Roach et al. [38] such straightforward solution, which includes the effects of different spatial arrangement and sizes of the floes, may soon become applicable in GCMs.
In general, the present research complements that of Wenta and Herman [17] and provides an explanation of the effects described there. The differences in area-averaged values of turbulent fluxes, water vapor and liquid water content in simulations with different FSDs found in [17] are explained based on the fact that FSD determines the spatial arrangement and intensity convection, which in turn controls the exchange of heat and moisture. Furthermore, the problem of high variability of local atmospheric conditions due to surface heterogeneity is addressed further and the first attempt to develop a coefficient integrating the effects of FSD into the calculation of surface moisture flux is made. The development of parametrization in the form of coefficient α will be continued in the subsequent research based on field observations and modeling with the full version of WRF model. Nonetheless, the results presented here provide new insights into the complex influence of sea ice fragmentation on the atmospheric boundary layer. In our opinion, the effects of FSD on the ABL can and should be included in the regional and global weather and climate models, as it is very likely that they are responsible for some of the known inaccuracies present in models results. Furthermore, it is expected that due to the ongoing warming of the Arctic region, sea ice fragmentation in winter is going to become more prevalent, thus further increasing the importance of properly taking into account sea ice heterogeneity in numerical models. Undoubtedly, more field campaigns are needed to improve our understanding of specific processes in situations with various sizes, shapes and orientations of ice floes and open water areas.
It must be noted that present study is based on an idealized model setup and that the results have not been validated with observations. It is a serious limitation. However, as already discussed in the previous paper [17], the crucial aspects of the simulated processes in our results are realistic and have been observed in other modeling and observational studies in the polar regions (e.g., [8,13,46,47]). Furthermore, the described analysis provides an invaluable support in planning of our further research which will include observations and modeling of atmospheric response to sea ice fragmentation in the Bay of Bothnia (Baltic Sea). The aim in this follow-up study is to confirm, based on the field data that FSD does affect the ABL properties and circulation, and to determine to which degree it controls the spatial arrangement and intensity of convective structures within the ABL.
Finally, it is worth noticing that the correction method presented here, based on a coefficient defined as the ratio of the "true" and biased value of the moisture flux, is not applicable to the second component of the total turbulent heat flux, namely the sensible heat flux. Whereas the values of the latent heat flux are always positive (with very few exceptions of extremely small negative values), the sensible heat flux can be both positive and negative, and our previous results [17] showed that different estimation methods-analogous to Equations (5), (6) and (9)-might produce results with a different sign, making a correction based on their ratio meaningless. Thus, a different approach will be necessary for an analogous parametrization of the sensible heat flux.