the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
How many parameters are needed to represent polar sea ice surface patterns and heterogeneity?
Joseph Fogarty
Mitchell Bushuk
Linette Boisvert
Sea ice surface patterns encode more information than can be represented solely by the ice fraction. The aim of this paper is thus to establish the importance of using a broader set of surface characterization metrics and to identify a minimal set of such metrics that may be useful for representing sea ice in Earth system models. Largeeddy simulations of the atmospheric boundary layer over various idealized sea ice patterns, with equivalent ice fractions and average floe areas, demonstrate that the spatial organization of ice and water can play a crucial role in determining boundary layer structures. Thus, various methods used to quantify heterogeneity in categorical latticebased spatial data, such as those used in landscape ecology and Geographic Information System (GIS) studies, are employed here on a set of recently declassified highresolution sea ice surface images. It is found that, in conjunction with ice fraction, patch density (representing the fragmentation of the surface), the splitting index (representing variability in patch size), and the perimeter–area fractal dimension (representing the tortuosity of the interface) are all required to describe the twodimensional pattern exhibited by a sea ice surface. For surfaces with anisotropic patterns, the orientation of the surface relative to the mean wind is also needed. Finally, scaling laws are derived for these relevant landscape metrics, allowing for their estimation using aggregated spatial sea ice surface data at any resolution. The methods used in and the results gained from this study represent a first step toward developing further methods for quantifying variability in polar sea ice surfaces and for parameterizing mixed ice–water surfaces in coarse geophysical models.
 Article
(1916 KB)  Fulltext XML

Supplement
(1894 KB)  BibTeX
 EndNote
The polar sea ice surface, a sensitive indicator of global climate change, exhibits persistent biases in sea ice fraction and extent in coarseresolution Earth system models (ESMs) (Liu et al., 2022; Casagrande et al., 2023; Myksvoll et al., 2023). Among other factors, these biases result from the inability of ESMs to resolve finescale spatial variability in sea ice and the resulting exchanges with the ocean below (Ramudu et al., 2018) and the atmosphere aloft (Bates et al., 2006; Esau, 2007). The effect of this subgridscale sea ice variability is typically parameterized in climate models using the ice fraction (f_{i}) as the sole surface composition characteristic. Usually, either an equivalent homogeneous surface or mosaic flux aggregation is used (Elvidge et al., 2016; BouZeid et al., 2020; Elvidge et al., 2021), but both yield an average flux, weighted by the ice and water fractions, that is inaccurate as the flux does not account for the impact of surface heterogeneity on the dynamics of the lower atmosphere and the nonlinear interactions with airflow above the ice and water (de Vrese et al., 2016; Lüpkes et al., 2012). This incomplete representation of sea ice surfaces and boundary layer structures results in errors in the turbulent exchanges of heat, moisture, and momentum across polar sea ice surfaces (Nilsson et al., 2001; Bourassa et al., 2013; Taylor et al., 2018). The dynamics and secondary circulations below the first vertical grid cell level in climate models are particularly underresolved, and they have a direct impact on air–surface exchanges; it is thus imperative to understand how these features influence fluxes (Mahrt, 2000; Essery et al., 2003; de Vrese et al., 2016). These gaps in representing finescale dynamics and fluxes propagate to the projection of future changes in the Arctic climate system and the resulting surface energy budget (Persson et al., 2002; Miller et al., 2017), which may be one reason why climate model ensembles consistently underpredict Arctic sea ice sensitivity to surface temperature warming. This underprediction has persisted throughout the last three Intergovernmental Panel on Climate Change (IPCC) model development cycles (Stroeve et al., 2007; Rosenblum and Eisenman, 2016, 2017; Notz and Community, 2020). The resulting uncertainty in the ability of climate models to predict future sea ice evolution hinders effective action and decisionmaking; therefore, improving these models is imperative (Notz and Stroeve, 2018; Docquier and Koenigk, 2021).
The fringe zone that separates densely consolidated sea ice from the open ocean is known as the marginal ice zone (MIZ) – see Dumont (2022) for a review of the current state of MIZ research. In the MIZ, the size and organization of sea ice floes and water are influenced by winds, sea currents, waves, and material ice properties (Wang et al., 2016; Ren et al., 2021; Herman et al., 2021; Hwang and Wang, 2022). What makes this region unique is that nearsurface air temperatures may fall between the surface temperatures of sea ice and water, resulting in abrupt spatial transitions between stabilizing and destabilizing surface buoyancy fluxes (Lüpkes et al., 2012). Such transitions produce drastically different turbulencemeanbased nonequilibrium dynamics and timescales (as demonstrated for comparable land–water transitions by Allouche et al., 2021, 2023b), all of which affect the surface–atmosphere exchanges between air, water, and sea ice. The ice fraction (f_{i}) in the MIZ ranges from 15 % to 80 % (Strong et al., 2017); however, any region of fractured sea ice gives rise to these abrupt transitions. It is precisely in these regions where the linearweighted, averaged approaches described above will be most inadequate and where surface transitions will play a key role in the dynamics. Thus, it is important to devise better methods for quantifying the heterogeneity of a surface, characterizing its patterns, and encoding this information into coarseresolution ESMs to better represent the polar environment.
To this end, the complex geometric patterns formed by sea ice floes need to be analyzed. Larger floes will have a proportionally greater effect on surface–atmosphere fluxes, whereas smaller floes, with their more frequent transitions, will exacerbate the nonlinearity of the exchange processes. These surface–atmosphere fluxes have a large effect on the atmospheric boundary layer (ABL) overlaying the marginal ice zone (MIZABL). As a thought experiment, we can consider an ice–water surface with a very fine checkerboard pattern and a sea ice fraction (f_{i}) of 0.5; this configuration will lead to statistically homogeneous ice floes that are locally variable at the surface but effectively homogeneous with respect to the MIZABL, where turbulence rapidly mixes the floes' smallscale signatures (Brutsaert, 2005; Mahrt, 2000; BouZeid et al., 2004). However, two large patches of sea ice and water (mesoα heterogeneity; see BouZeid et al., 2020), which also have a sea ice fraction (f_{i}) of 0.5, will generate significant circulation closer to that of a sea breeze due to the abrupt transition between the two large, homogeneous surfaces (Porson et al., 2007; Crosman and Horel, 2010; Allouche et al., 2023a). The dynamics and thermodynamics in this MIZABL system, and the surface exchange therein, will thus be quite different for these two patterns, even if some key surface properties, e.g., temperature and roughness, are identical (BouZeid et al., 2007).
Given the importance of the topic and the challenges outlined above, previous work has attempted to quantify the heterogeneity of sea ice surfaces (Wenta and Herman, 2018, 2019; Michaelis et al., 2020; Horvat, 2021; Dumont, 2022), utilizing surface and meteorological properties like sea ice fraction, geostrophic velocity, lead width, and floe size distribution. Furthermore, parameterizations for flow over leads in sea ice have been developed based on noneddyresolving models (Lüpkes et al., 2008; Michaelis et al., 2021). Michaelis and Lüpkes (2022) also a conducted turbulence parameterizations (based on largeeddy simulation models) over ensembles of leads, but they used a twodimensional ice fraction geometry and a higher ice fraction. However, the smallscale patterns in the MIZ require broader and more versatile methods of heterogeneity characterization (e.g., Mandelbrot, 1967), especially as the resolution increases. In addition, the computational grid of even the highestresolution weather or climate model cannot resolve all the spatial features in the MIZ. One thus needs to consider how to represent unresolved surface characteristics in such models; these characteristics can be thought of as latticetype spatial structures, as defined by Cressie (1993). Observational data of MIZ ice patterns also have a finite resolution and are thus comparable to lattice data, allowing for the use of different metrics specifically defined for lattice surfaces and offering ways to characterize the heterogeneity patterns of these polar surfaces. In this paper, we examine approaches for this quantification that are commonly used in landscape ecology, a field that has generated a multitude of methods to study lattice spatial data (Li and Reynolds, 1994, 1995; Pickett and Cadenasso, 1995).
Studies in landscape ecology have previously sought an optimal independent group of metrics for understanding the heterogeneity of lattice surfaces. Riitters et al. (1995) used multivariate factor analysis to suggest six groups of metrics, including image texture, average patch compaction, and average patch shape. Cushman et al. (2008) used principal component analysis to propose seven broad metrics at the landscape level, including contagion, largepatch dominance, and proximity (see Table 9 in the aforementioned study). For the twodimensional, binary seaice–water surfaces considered in this study, we chose the variance inflation factor (VIF) technique to reduce these metrics to a compact set of metrics that are weakly dependent on one another, minimizing information redundancy (Miles, 2014).
The questions that will be answered in this study are as follows:

Is the sea ice fraction of a MIZ surface, combined with some measure of average floe area, sufficient to predict the behavior of the overlying MIZABL?

If not, what other surface information in a twodimensional, latticebased spatial pattern is needed to describe air–sea interactions?

How can this surface information be applied to sea ice surfaces in weather models and ESMs, considering factors such as availability of information, resolutionresampling invariance, and ease of understanding?
Section 2 details the methods employed for the idealized and realworld maps used in this study, including the steps taken to reduce multicollinearity and determine which landscape metrics, alongside sea ice fraction, provide additional information on the pattern of the sea ice surface. The largeeddy simulations used are also presented in this section, and a more complete description is found in Appendix A and B. Section 3 reports the results of the idealized sea ice surfaces in the largeeddy simulation, thereby answering the first of our questions and setting the stage for the remaining questions. Section 4 presents the results from the 2D surface analysis, thereby answering the second and third questions, with additional discussion on principal directions and climate model implications found in Sect. 5. Section 6 synthesizes the findings and outlines open questions to guide future investigations of sea ice heterogeneity.
2.1 Largeeddy simulations
Largeeddy simulations (LESs) of the MIZABL were conducted for different idealized configurations of sea ice. LESs are widely used to model heterogeneous highReynoldsnumber flows (Baidya Roy, 2002; BouZeid et al., 2004) in, for example, convective boundary layers (Courault et al., 2007; Maronga and Raasch, 2013), stable boundary layers (Huang et al., 2011), and coastlines (Allouche et al., 2023a) – see Sect. 3.6 of Stoll et al. (2020). This heterogeneous highReynoldsnumber description aptly applies to the MIZABL. Unlike direct numerical simulations (DNSs), LESs are able to attain a Reynolds number (Re) representative of the MIZABL (∼10^{7}) because the smaller turbulent eddies (smaller than the filter size, which is comparable to the numerical grid spacing in our simulations) are not explicitly resolved. However, unlike Reynoldsaveraged Navier–Stokes (RANS) approaches, which encompass all weather and climate models, LESs directly resolve and capture large turbulent eddies, the heterogeneity of the surface, advective fluxes, and largescale sea ice patterns, making them a computationally and physically appealing approach for the problem at hand. By retaining these larger structures, most of the turbulent energy and fluxes are explicitly resolved, allowing for the investigation of threedimensional flow structures that may arise over these heterogeneous surfaces.
LESs are thus used to model MIZABL flow over 10 km × 10 km patterns of idealized ice–water surfaces, modulated by a Coriolis force at a latitude (Φ) of 90° N in a horizontally periodic domain, yielding a Rossby number (Ro) of 13.7 (see Fig. 1). We note that simulations using the same Ro value, even those conducted at lower latitudes with mean wind speeds, would give similar results – see the dimensional analysis of flow over heterogeneous surfaces in Fogarty and BouZeid (2023a) and Allouche et al. (2023a). This full domain is smaller than a single grid cell in stateoftheart ESMs, underlining why simulated properties of the MIZABL in ESMs require subgridscale (SGS) parameterizations. For leads with a width of ∼1 km, a grid spacing of 10–20 m is usually chosen after grid convergence tests (Lüpkes et al., 2008; Gryschka et al., 2023). Our choice of a coarser horizontal resolution of 100 m reflects our aim of zooming out from a typical lead and examining the MIZ as a whole (which is also reflected in the horizontally periodic nature of the domain). A finer resolution in our simulations, while possibly improving the representation of turbulence and plume dynamics, may sacrifice some of the large scales and secondary circulations that arise from these heterogeneous surfaces, which we aim to capture with a large domain. The coarse resolution also does not compromise our ability to answer the first question because, as shown later, the differences in the dynamics between simulations with identical ice fractions but different surface patterns are significant and far exceed any plausible impact resulting from grid resolution. The simulations in this section are meant to demonstrate the need for surface analysis, given a constant ice fraction and an average ice floe area; thus, we do not focus on the quantitative aspects of the output. See Table 1 for all details of the simulations in this study; more information on the numerical aspects of our LES is provided in Appendix A.
The bottomboundary condition for each simulation can be thought of in terms of categorical latticebased spatial data, where each node represents either the ice class or the water class. An ice node is prescribed a surface temperature (θ_{i}) of 255 K, typical of fall and spring temperatures in the central Arctic, and a momentum roughness length of 1 mm, while a water node is prescribed a surface temperature (θ_{w}) of 271 K, roughly the freezing point of seawater, and a roughness length of 1 cm. The heat roughness length is set to 0.1 mm for the entire surface (both ice and water). There is a very large variability in roughness lengths in Arctic sea ice simulations (see Lüpkes et al., 2008; Andreas et al., 2010; Lüpkes et al., 2012; Elvidge et al., 2016; Gryschka et al., 2023), which reflects the physical differences between calm and stormy seas, as well as between new, flat ice and old ice that is deformed and refrozen. However, our choice of roughness lengths prioritized keeping constant ratios (which are the key dimensionless parameters influencing the results; Fogarty and BouZeid, 2023a; Allouche et al., 2023a) among simulations to focus on the effect of sea ice patterns. Results with different values but the same ratios will lead to almost identical conclusions. In addition, sensitivity analyses (not shown here) indicated that the roughness lengths, even when their ratios were changed, had a very minor impact on the quantitative results compared to the temperature contrast, ice fraction, and ice–water patterns. The initial potential air temperature (θ_{a,0}), a constant profile, is defined so that the areaaveraged sensible heat flux, as computed by the Monin–Obukhov surface flux parameterizations, is zero and thus lies between the ice and water surface temperatures (see Appendix B for details). The LESmodeled heat flux will, however, not be zero.
All five patterns (Pattern1 to Pattern5), displayed in Fig. 2, were simulated. They have a fixed sea ice fraction (f_{i}) of 0.46 and a mean floe area of 11.56×10^{6} m^{2}. The geostrophic wind (M_{g}=2 m s^{−1}) flows from left to right at an angle of 18° relative to the x axis in all simulations. This is expected to give a surface wind (M_{0}) value that is roughly aligned with the x axis for homogeneous neutral surfaces due to Ekman veer (Ghannam and BouZeid, 2021). These simulations are dry runs with no clouds present. The turbulence field is warmed up for about 60 h, and the statistics are then Reynoldsaveraged over an additional 12 h. A variable with an overbar denotes averaging over time, which is used as a surrogate for ensemble Reynolds averaging, and any spatial averaging over the heterogeneous domain in x and y is indicated by angle brackets.
In addition to the Rossby number and the roughness ratio discussed earlier, an important dimensionless input parameter in these simulations is the heterogeneity Richardson number, defined as
Ri_{h} encodes the competition between buoyancydriven circulations generated by the surface temperature contrast and the uniform flow resulting from the synoptic forcing (M_{g}). All simulations in this study have equivalent inputs of Ri_{h}, Ro, and ${z}_{\mathrm{0},\mathrm{w}}/{z}_{\mathrm{0},\mathrm{i}}$. Based on previous dimensional analyses and LESs of flow over heterogeneous surfaces (Omidvar et al., 2020; Allouche et al., 2023a; Fogarty and BouZeid, 2023a), matching these dimensionless inputs is required to be able to focus on the effects of surface sea ice fraction and patterns.
2.2 Ice map data
While the LES utilizes idealized surfaces to examine the influence of patterns on the MIZABL, examining what other landscape metrics might be important for surface characterization necessitates the use of realworld sea ice maps. The lattice spatial data used in the statistical analysis (see Sect. 2.3) are derived from recently declassified highresolution (1 m) Literal Image Derived Products (LIDPs) in accordance with national technical means (NTMs), as detailed in Kwok (2014). These images underwent a supervised maximumlikelihood classification algorithm, which assigned each pixel in the original LIDP to either a water class or an ice surface class (Fetterer and Untersteiner, 1998; Fetterer et al., 2008). This process converted the highresolution LIDPs into categorical latticebased spatial data, where each cell represents one of two possible surface types (ice or water).
These maps, which have a horizontal extent of up to 10 km by 10 km, comprise the dataset used to calculate the landscape metrics. Some of the images did not fully cover this full extent; thus, in order to retain the realworld sea ice geometry, we reflected the images onto areas with no data. All metric calculations and analyses were conducted on these modified surfaces. The advantage of this highresolution, largeextent dataset is that we can analyze how these metrics change with grain size. These maps are thus aggregated from a 1 m resolution to 2 m, 10 m, 50 m, 100 m, 200 m, 500 m, 1 km, and 2 km resolutions; resampling was done using the nearestneighbor method from the Python Imaging Library (PIL). These resolutions cover a range of common resolutions used in fine largeeddy simulations and numerical weather prediction (NWP) models. Due to the excessive computationalprocessing time required for the original 1 m resolution data, the highest resolution at which landscape metrics were calculated was 2 m.
One important caveat of this dataset is that the histogram illustrating sea ice fraction (f_{i}) is heavily skewed toward higher sea ice fractions (see Fig. 3). We recognize that this may lead to bias in the results; however, the analysis methods developed here are insensitive to f_{i} and can certainly be applied to other datasets with more uniform sea ice fraction distributions in the future. The high resolution provided by the present dataset remains a key factor in its adoption for the present study since it allows us to aggregate and track how landscape metrics change from actual surface states when resampled to a coarser numerical grid.
2.3 Landscape metric space reduction
FRAGSTATS, a spatialpattern analysis program, was used to calculate landscape metrics (McGarigal and Marks, 1995) based on the lattice spatial data. When given GeoTIFFformat raster lattice data, FRAGSTATS calculates patch metrics, class metrics, and landscape metrics of your choice. Patch metrics are computed for each individual patch in the landscape and are thus not relevant for the current study of sea ice surface patterns. Class metrics are computed for each patch type (class) in the landscape. In this study, this would involve calculating metrics for sea ice only and water only. Although this metric may be useful in other applications of pattern analysis, in this study, we want to look at the aggregated patterns of sea ice and water combined. Thus, only landscape metrics were calculated. Sea ice fraction (calculated as the number of ice cells divided by the total number of ice and water cells) was calculated using the PIL as it is the only metric not provided by FRAGSTATS.
This resulted in 22 landscape metrics that focus on the global patterns of the surface. Many of these landscape metrics, however, are correlated with one another; for example, patch density and mean patch size are proportional to each other. This is due to the fact that there are limited observations one can make about a surface (such as the number of patches, the area of a patch, and the proportion of edge in a patch) but an infinite number of operations that can be performed on a surface. Therefore, many of these metrics (especially at the landscape level) simply represent different methods of aggregating or statistically analyzing these observations.
While collinearity between two metrics can be easily detected through a correlation matrix, multicollinearity (where one indicator is a linear combination of two or more other indicators) is more likely in these types of datasets. It is thus possible for two or more landscape metrics to jointly define another metric. An objective and statistical way of reducing these parameters is thus needed. Here, we chose the variance inflation factor (VIF),
where R^{2} is the coefficient of determination, which is used to detect multicollinearity among these heterogeneity parameters (see Ibidoja et al., 2023). For each metric (X_{i}), when $i\in [\mathrm{1},\mathrm{\dots},\mathrm{22}]$, the VIF was calculated using a regression equation,
where α_{0} is a constant. The statsmodels
Python library (Seabold and Perktold, 2010) was used for these computations. The metric with the largest VIF was then removed from the dataset as it was not considered important in the quantification of sea ice surfaces. All VIFs were then recalculated for this “reduced” dataset, and the new metric with the highest VIF was removed. Through this process, metrics were removed one by one until all remaining metrics exhibited a VIF below a predefined cutoff, which was set as VIF<2.5. While such a low cutoff may not be necessary in certain practices of multicollinearity reduction (O'Brien, 2007), the ultimate goal of this technique is to reduce the parameter space. In other words, for climate modelers, fewer metrics in their SGS parameterizations result in more practical models.
Each of these metrics, listed in Table C1, can be categorized into one of six metric groups: Area and Edge, Shape, Core Area, Aggregation, Contrast, and Diversity. The first four metric groups are important for analyzing sea ice surfaces. Area and Edge metrics deal with the size of floes and the amount of edge they create, while Shape metrics discriminate based on patch morphologies and overall geometric complexity. Core Area metrics analyze the area within a patch beyond a specified buffer width. Aggregation metrics focus on the tendency of patches of similar types to be spatially aggregated or otherwise dispersed across the landscape.
The last two metric groups, Contrast and Diversity, are less important for the present application to sea ice. Contrast metrics refer to the magnitude of difference between adjacent patch types with respect to a given attribute – in the case of sea ice surfaces, which only correspond to two classes (ice and ocean), there is only one contrast between two categories. Thus, metrics in this group are simply represented by the contrast of surface temperature and roughness. Diversity metrics are influenced by the number of patch types present and the areaweighted distribution of these patch types. In this case, we only have two types of patches (ice and water), so the diversity is the same across all maps, and the weighted distributions of the patches are related to the ice fraction. Further information on all of these metric groups can be found in the FRAGSTATS manual (McGarigal and Marks, 1995).
The largeeddy simulation technique detailed in Sect. 2.1 was used to simulate the MIZABL across different configurations of sea ice patterns. Figure 4 displays the Reynoldsaveraged and horizontally averaged normalized vertical profiles of horizontal wind speed ($M=\langle \sqrt{{u}^{\mathrm{2}}+{v}^{\mathrm{2}}}\rangle $), normalized by the geostrophic velocity (M_{g}); potential air temperature, normalized by the initial potential air temperature (θ_{a,0}); total heat flux ($\langle \stackrel{\mathrm{\u203e}}{w\mathit{\theta}}\rangle $), normalized by M_{g}θ_{a,0}; and the total horizontal stresses ($\langle \stackrel{\mathrm{\u203e}}{uw}\rangle $ and $\langle \stackrel{\mathrm{\u203e}}{vw}\rangle $), both normalized by ${M}_{\mathrm{g}}^{\mathrm{2}}$. Note that these include both turbulent and dispersive fluxes that arise over heterogeneous surfaces as a result of the spatial correlation of the mean (timeaveraged) fields (Raupach and Shaw, 1982; Finnigan and Shaw, 2008; Li and BouZeid, 2019). The results clearly display significant differences, underlining the fact that sea ice fraction alone is not a sufficient surface metric for describing MIZABL dynamics, even in these simulations, where the mean floe size was also kept constant.
In all five simulations, the largest difference occurs between Pattern4 and Pattern5 (red and purple lines, respectively), which is expected since the geostrophic wind, and thus the nearsurface wind, flows parallel (rather than perpendicular) to the strips of ice (Willingham et al., 2014; Anderson et al., 2015; Salesky et al., 2022; Fogarty and BouZeid, 2023a). Although wind direction is not a surface property of sea ice, its orientation relative to the surface features is still an important driver that must be taken into consideration and will be discussed further. All patterns developed a lowlevel jet (LLJ), which can be seen in Fig. 4a, although the LLJs in Pattern1 and Pattern5 are weak (Tetzlaff et al., 2015; Michaelis et al., 2021). The LLJs seem to strengthen in Pattern2 and Pattern3, likely due to large swaths of ice aligned with the direction of the geostrophic wind; the unstable–stable transitions in these ice regions decouple the air from the surface friction over long, stable ice patches, allowing for lowlevel wind acceleration. This mechanism is similar to the one advanced by Blackadar (1957), which concerned the creation of a lowlevel jet via inertial oscillation over time as the ABL transitioned to a stable regime at sunset. However, in this case, the oscillation occurs in space, with columns of air advecting from hot to cold surfaces and decoupling from the surface. The secondstrongest LLJ is in Pattern3 and is likely caused by the one large ice floe. It persists across all wind directions, and the same can be said for Pattern2. The strongest LLJ, found in Pattern4, is likely reinforced by secondary circulations (not shown), consisting of streamwisealigned rolls driven by the lateral contrast in surface temperature, and by the fact that the small strips of ocean between the ice floes are not wide enough to interrupt the circulations. However, unlike Pattern2 and Pattern3, Pattern4 is highly anisotropic, meaning that the geostrophic wind direction is more important in this case, as evidenced by the fact that Pattern5 has no LLJ despite being a 90° rotated version of Pattern4.
Each simulation showed an increase in potential temperature of 2.1 %–2.2 %; this warming is consistent with the large ocean fraction in each of the patterns. Despite initializing the air temperature to produce zero fluxes according to the Monin–Obukhov flux models, which assume equilibrium between the air and water above each surface grid cell, the upward fluxes over warm water are larger than the downward fluxes over cooler ice. This is due to the effect of advection, which perturbs the equilibrium.
Major differences are also observed in the total streamwise and crossstream stresses, displayed in Fig. 4d and e, respectively. Again, Pattern4 and Pattern5 exhibit the greatest differences from one another due to the geostrophic wind direction. All simulations show similar negative streamwise stress in the surface layer, but at higher altitudes in the MIZABL, the differences between simulations are greater. Above the LLJ, some of the stresses turn positive, implying an upward transfer of momentum from the LLJ, which explains the differences that may be observed below the blending height (Wood and Mason, 1991; Mahrt, 2000; Brunsell et al., 2011). The crossstream components of the total momentum flux, shown in Fig. 4e, are all quite distinct, indicating significant differences in the Ekman rotation of wind and stress with height.
To further explain the differences seen in these patterns, we consider the decomposition of the total flux into its dispersive and turbulent contributions. Vertical turbulent heat flux is denoted as $\stackrel{\mathrm{\u203e}}{{w}^{\prime}{\mathit{\theta}}^{\prime}}$, while vertical turbulent streamwise and crossstream stresses are denoted as $\stackrel{\mathrm{\u203e}}{{u}^{\prime}{w}^{\prime}}$ and $\stackrel{\mathrm{\u203e}}{{v}^{\prime}{w}^{\prime}}$, respectively. Dispersive fluxes arise from a timeaveraged but spatially variable mean flow (Raupach and Shaw, 1982; Li and BouZeid, 2019). Since our Reynolds averaging is performed over time, we can spatially decompose any Reynoldsaveraged variable. For example, the vertical velocity can be decomposed as follows: $\stackrel{\mathrm{\u203e}}{w}=\langle \stackrel{\mathrm{\u203e}}{w}\rangle +{\stackrel{\mathrm{\u203e}}{w}}^{\prime \prime}$, where the brackets represent spatial averaging (as defined in Sect. 2.1) and the double prime represents variations in the mean planar fields in space. We then calculate the local dispersive fluxes, using vertical heat flux as an example, as follows:
We then spatially average the entirety of Eq. (5) over the horizontal plane to obtain
The middle two terms in Eq. (5) are zero due to spatial averaging, since $\langle {\stackrel{\mathrm{\u203e}}{w}}^{\prime \prime}\rangle =\mathrm{0}$ and $\langle {\stackrel{\mathrm{\u203e}}{\mathit{\theta}}}^{\prime \prime}\rangle =\mathrm{0}$ by definition; therefore, these terms have no impact on spatially averaged surface–atmosphere exchanges. Furthermore, $\langle \stackrel{\mathrm{\u203e}}{w}\rangle $ is assumed to be very small (unless strong, largescale subsidence or uplift is present). In our LES, $\langle \stackrel{\mathrm{\u203e}}{w}\rangle $ must be zero since there cannot be accumulation or depletion of mass below a given horizontal plane in a periodic domain with an incompressible flow. Thus, the first term on the righthand side of Eq. (5) is also negligible, leaving only one term remaining in Eq. (6). This term, denoting dispersive flux, is of particular interest as it represents the coherent spatial correlation of vertical velocity and potential temperature in regions with consistent secondary structures, such as consistently warm updrafts, cool downdrafts, or streamwise rolls.
Figure 5 shows the horizontally averaged vertical profiles of the total, turbulent, and dispersive fluxes of each pattern with respect to heat flux, streamwise momentum flux, and crossstream momentum flux, allowing us to decompose and analyze the total fluxes shown in Fig. 4c–e. For example, it is clearer now that the crossstream components of the total momentum flux (Fig. 4e) are distinct from one another due to these dispersive fluxes (dotted green lines). Over these heterogeneous surfaces, dispersive crossstream stress dominates over its turbulent counterpart in all patterns except Pattern5 (see panels c, f, i, and l in Fig. 5). The magnitudes of these dispersive fluxes are not equal when the heterogeneous surfaces differ from one another. Thus, it is not the ice fraction or average floe area but the surface pattern itself that leads to these differences in total crossstream flux. The streamwise stresses, on the other hand, seem to balance the dispersive and turbulent stresses, which are not always of the same sign (panels e, h, k, and n in Fig. 5). Thus, the dispersive components that directly result from the secondary motions imprinted by the surface pattern on the atmosphere are also critical here. These secondary circulations can be seen in the twodimensional cross sections in the supporting information.
The total heat flux in all these simulations linearly decreases with height, as dictated by the LES setup. As seen in Fig. 4c, the variations are not as impressive as those concerning the momentum fluxes, but they can still result in a difference of up to 30 %, especially near the surface. However, analyzing the left column of Fig. 5 (i.e., panels a, d, g, j, and m) shows that the relative importance and profiles of the dispersive and turbulent heat fluxes exhibit more significant differences. The various surface patterns seem to lead to differences in the dispersive fluxes, but in all cases, these are balanced out by the turbulent fluxes. Nevertheless, in all of the figures, there is strong variability near the surface in the dispersiveflux and turbulentflux profiles. This is partially due to the shallow internal boundary layers that are created by the mean flow and the secondary circulations over the floes, and these differing secondary circulations arise from differences in the surface patterns. However, one should also note that the first few points are strongly influenced by the transition from the wall model to the SGS model in representations of resolved turbulence – a persistent challenge in LESs (Piomelli and Balaras, 2002; Brasseur and Wei, 2010). Thus, the quantitative details of the results from such regions should be interpreted with care.
Lastly, analyzing the ratios between dispersive and total atmospheric vertical fluxes (Table 2) provides insight into the differences between simulations and allows for comparisons with other surface types. For example, in Pattern1 and Pattern2, $\left\overline{v}\overline{w}\right/\left\stackrel{\mathrm{\u203e}}{vw}\right>\mathrm{1}$, indicating that the dispersive and turbulent fluxes along the crossstream axis are in opposite directions, as seen in the green and orange profiles in Fig. 5c and f. The ratios in Pattern4 and Pattern5 are also closer to unity than the other values, which, again, highlights the differences in secondary circulations due to the sea ice pattern. These values have been observed before – for example, over urban or forest canopies (Moltchanov et al., 2015; Boudreault et al., 2017; Li and BouZeid, 2019).
Overall, these LES results indisputably indicate that ice–water patterns hold key information on how the MIZABL interacts with the underlying surface; thus, the rest of the paper is dedicated to characterizing these patterns.
Now that we have established the need for surface characteristics beyond sea ice fraction, we aim to examine the indicators that can be used for this purpose. For each of the nine resolutions considered, 44 observed sea ice images were analyzed. When conducting this analysis at the 2 m resolution, four metrics (including sea ice fraction) remained after the VIF elimination process and did not exhibit multicollinearity with one another: sea ice fraction (f_{i}), patch density (PD), the splitting index (SPLIT), and the perimeter–area fractal dimension (PAFRAC). These remaining metrics were then categorized into the different metric groups defined in Sect. 2.3: sea ice fraction is an Area and Edge metric, PD and SPLIT are both Aggregation metrics, and PAFRAC is a Shape metric. No metrics remained in the Diversity, Core Area, or Contrast metric groups. We expected no metrics to remain in the Contrast metric group since an edge can only exhibit one contrast for a seaice–water surface – however, this could change if the analysis were conducted with continuously variable surface temperatures or roughness lengths. We also did not expect any metrics from the Diversity group since many metrics, such as evenness or Simpson's diversity index, are functions of sea ice fraction. The absence of metrics from the Core Area group was not predicted, but it is likely related to the fact that the Shape group and the Area and Edge group can collectively represent the Core Area characteristics.
The first aggregation metric, PD (Riitters et al., 1995; Šímová and Gdulová, 2012), with a VIF of 1.9, represents an areanormalized number of patches and is described as follows:
where PD represents patch density, A_{t} represents the total area of the surface, and n represents the total number of distinct patches (of either sea ice or water). As the PD of a sea ice surface increases, one would expect to find more instances of ice–water edges and thus more regions of stabletounstable and unstabletostable stratification transitions. We also note that reducing PD increases the average time a parcel spends over the stable (or unstable) surface, which affects how the parcel adjusts to the transition to a new stability regime. Patch density may also work in tandem with geostrophic wind direction, as discussed in Sect. 3, since geostrophic wind flowing in one direction may encounter more ice–wateredge transitions than in another direction (see Pattern4 and Pattern5, for example).
The second metric, SPLIT, with a VIF of 1.9, was first described in Jaeger (2000). It represents the inverse probability that two randomly chosen points on a map will be in the same patch, and the corresponding equation is expressed as follows:
where a_{i} is the area of patch i and the index (i) iterates over all patches. SPLIT has a value of n if there is only one patch or if all patches are of an equal size. However, generally, SPLIT has a value greater than n, with lower values indicating a larger variance in patch sizes. Since we already use PD, the new information provided by SPLIT is specifically about the variance in patch sizes.
The only Shape metric, PAFRAC, with a VIF of 2.1, is obtained by regressing each patch's perimeter (P_{i}) against its area (A_{i}) on a log–log plot, resulting in the following:
where k is a constant and PAFRAC is the perimeter–area fractal dimension, which measures the tortuosity, or jaggedness, of an ice–water interface, similar to any fractal dimension (Mandelbrot, 1982).
Thus, in addition to sea ice fraction (f_{i}), the three metrics that would be useful in describing a sea ice surface are SPLIT, PD, and PAFRAC. Table 3 details the values of these metrics for each of the sea ice patterns simulated in Sect. 3. We observe that SPLIT is invariant to the shape of the floes as long as the areas of the floes are equivalent (i.e., no variance), which is why SPLIT remains consistent for all patterns except Pattern3. This is where PAFRAC proves useful as it provides three different values among the five ice–water patterns.
For the real ice maps obtained from Arctic images, the ice fraction varies from 0.19 to 0.99, PD varies from 5.4 to 42.5, SPLIT varies from 1.23 to 4.07, and PAFRAC varies from 1.338 to 1.726 at a 2 m map resolution. In many cases, however, numerical simulations also require the resampling of highresolution surfaces by increasing the grain (pixel) size. For example, sea ice maps from reconnaissance satellites may have a resolution of up to 1 m, but this is computationally impractical for numerical weather models. Largeeddy simulations of the ABL can have resolutions down to 50 m, while NWP models have resolutions of 2 to 10 km. Therefore, even with highresolution data, the aggregation and resampling of these surface patterns are inevitable in modeling. Furthermore, considering the operational use of these metrics, regularly updating these values would likely involve multiple satellite products with differing resolutions; thus, metrics that can be extrapolated/interpolated across different grid cell sizes would allow for consistent computation when standardized to a single weather model grid cell.
Therefore, it is useful to examine how these chosen metrics vary as an image is aggregated to a resolution applicable to numerical weather models (or other numerical models, such as LESs); an appealing metric would be one that is invariant to resolution changes. Sea ice fraction (f_{i}) is a good example of an invariant indicator, as shown in Fig. 6a, where it is calculated for all images and then averaged over the corresponding resolution. A secondbest case would involve a metric that exhibits a clear scaling law relative to the resolution, such as the PD depicted in Fig. 6b. In this case, the PD for the “real” 2 m resolution surface can be extrapolated to higher resolutions based on powerlaw scaling, expressed as
where m is the metric, Δ is the map resolution, and k and D_{q} are the scaling coefficients.
Some metrics, such as SPLIT and PAFRAC, seem to exhibit nearinvariant behavior after a certain jump in the resolution. For example, when starting from the 10 m resolution, SPLIT stays fairly constant as the resolution decreases. There is also variation in PAFRAC as the resolution decreases from 10 m. This is consistent with results from previous studies as some landscape metrics exhibit large errors when surfaces are aggregated to lower resolutions (Moody and Woodcock, 1994, 1995). Given the nearinvariant scaling of three of the four metrics and the predictable powerlaw scaling of the fourth, we can proceed with these four metrics since they are usable (or translatable) across scales.
Thus far, we have identified four surface pattern indicators that characterize the MIZ surface: sea ice versus water concentration (f_{i}), patch density (PD), variance in the sizes of patches (SPLIT), and the tortuosity of patch edges (PAFRAC). However reexamining Pattern4 and Pattern5 in Fig. 4 raises the question of why these maps show the largest differences in their respective MIZABLs. Although their geometric patterns are the same (and thus the four metrics are identical for the two configurations), the difference in the geostrophic wind direction results in large differences in the surface–atmosphere interactions. This suggests that another important attribute relates to how the surface patterns are oriented relative to the wind. If the surface is isotropic, the wind angle should be irrelevant, but most water–seaice patterns in the MIZ display a significant degree of anisotropy (Feltham, 2008). Therefore, quantifying the impact of surface orientation and including it in the set of metrics obtained from the VIF analysis, as discussed in Sect. 4, may provide additional information, allowing modelers to parameterize MIZABL dynamics in global climate models.
We observe that in Pattern4 and Pattern5, the difference in geostrophic direction is related to the directionality of sea ice organization. In Pattern4, the wind blows consistently over an infinitely repeating pattern of sea ice and water at regular intervals. In Pattern5, the wind blows over much longer strips of ice and water, even though some seaice–water transitions are present. Any other oblique flow is thus in between these two parallel and perpendicular regimes. We characterize the differences between these regimes by examining the variance in the surface exposed to the wind. In other words, Pattern4 exhibits high variance since the surface wind flows over a maximum of eight ice–water transitions within one domain length, whereas Pattern5 exhibits low variance since the wind flows over a maximum of two seaice–water transitions. This then raises the question of how to determine a principal direction for a more complex surface.
We therefore attempted to characterize this anisotropy by computing the direction of the eigenvector (the eigendirection) for the surface with the least amount of variance, resulting in the fewest ice–water transitions possible. This was done using the scikitlearn
Python package, specifically via principal component analysis (PCA) using the sklearn.decomposition.PCA
class (Pedregosa et al., 2011). This method performs an eigendecomposition of the covariance matrix of the sea ice map, yielding two orthogonal eigenvectors. The principal eigendirection points in the direction of minimal variance, as indicated by the longer of the two arrows in Fig. 7. In other words, the principle or most coherent mode that best explains the pattern is aligned with the direction of least variability. The secondary eigendirection is, by definition, orthogonal to the principal eigenvector. Some of these eigendirections are intuitive as one can imagine trying to select a geostrophic wind direction that passes over a minimal number of ice–water edges. The maps in Figs. 7a and 7g provide two such examples. However, the map in Fig. 7f, for example, is a bit less intuitive – visual inspection may suggest that the principle eigendirection should be aligned from the lowerright corner of the map to the upperleft corner, but the results indicate a less obvious orientation.
It is hypothesized that, for a fixed sea ice fraction (f_{i}), geostrophic wind flowing in the principal direction with minimal variance will behave more like Pattern5, while the perpendicular angle to said principal direction (the secondary direction) will behave more like Pattern4 – similar to the parallel and perpendicular cases in the simulations conducted in Fogarty and BouZeid (2023a). Nevertheless, further LESs are needed to elucidate the exact impact of this relative orientation and the other parameters we identified on the MIZABL.
However, some of these maps exhibit a higher degree of anisotropy than others, such as Fig. 7c. To measure the degree of anisotropy in these maps, one can look at the percentage of variance (POV), defined as
where λ_{i} is an eigenvalue in an ndimensional matrix. The POV for a twodimensional surface (n=2) thus describes the amount of variance that can be explained (or reconstructed) by that eigenvector alone, resulting in the following relationship: $\text{POV}\left({\mathit{\lambda}}_{\mathrm{0}}\right)+\text{POV}\left({\mathit{\lambda}}_{\mathrm{1}}\right)=\mathrm{1}$. In theory, a sea ice map with a high POV(λ_{0}) value, and thus a low POV(λ_{1}) value, would be anisotropic since the secondary eigendirection would account for much more of the variance than the principal eigendirection, and the surface would thus have a preferential direction of variability (one would expect Fig. 7g to have a high POV(λ_{0}) value). Conversely, a map with a low POV(λ_{0}) value would be fairly isotropic.
By definition, POV(λ_{0})≥0.5 (where $\text{POV}\left({\mathit{\lambda}}_{\mathrm{0}}\right)=\text{POV}\left({\mathit{\lambda}}_{\mathrm{1}}\right)=\mathrm{0.5}$ indicates a truly isotropic surface) since POV(λ_{0}) represents the POV for the principal eigendirection. However, most of the ratios in these examples lie in the range of $\mathrm{0.50}<\text{POV}\left({\mathit{\lambda}}_{\mathrm{0}}\right)<\mathrm{0.60}$, showing little variability among these maps. Thus, in this small subset of sea ice maps, examining only POV(λ_{0}) would not provide information on how much influence the principal eigendirection has.
Although stability over an ice or waterdominated surface depends on many factors, such as wind direction and potential air temperature, in the cases where potential air temperature falls between the surface temperatures of ice and water, the ice fraction of a sea ice surface can be a fair indicator of the behavior of the MIZABL. However, this is only the case when the ice fraction approaches 0.0 (all ocean), leading to an unstable atmosphere, or 1.0 (all ice), leading to a stable atmosphere. On the other hand, when f_{i} falls between these limits (in other words, when the surface flow alternates between very stable and very unstable), the ice fraction alone is insufficient to predict the dynamics and thermodynamics of the MIZABL. Largeeddy simulations conducted for five different sea ice surfaces, as detailed in Fig. 2, have shown that surfaces with the same ice fraction, number of floes, and mean floe size can result in very distinct atmospheric dynamics. Differences were observed in horizontal wind speed (Fig. 4a) and total surface stresses (Figs. 4d–e and 5).
While Fig. 5 shows moderate differences in the total heat flux, which, in this case, is constrained by the simulation setup, more significant differences are observed in the dispersive and turbulent fluxes that make up this total flux (see also Table 2). The total, turbulent, and dispersive fluxes of the streamwise and crossstream momenta are even more sensitive to surface patterns. These dispersive fluxes have been shown to drive many of the differences and are thus nonnegligible in climate models (Margairaz et al., 2020; Fogarty and BouZeid, 2023a; Lu et al., 2023).
To understand what other information can be obtained from a twodimensional, binary lattice surface, we examined 44 spatial metrics traditionally used in the field of landscape ecology since knowing the cover fraction (ice fraction) and the number and median area of the floes is not enough to able to fully describe the ice–water–atmosphere physics. These 44 additional spatial metrics were applied to LIDPs of realworld satellite sea ice imagery to determine which metrics were important, and the variance inflation factor was used to detect and remove any multicollinearity in this dataset. The remaining metric set included the following metrics: ice fraction, patch density (representing the number of sea ice floes and, thus, their mean size in a given area), the splitting index (representing variance in floe sizes), and the perimeter–area fractal dimension (representing edge tortuosity). We also proposed using the surface eigendirection relative to the mean wind direction to characterize the influence of surface anisotropy and its interaction with the wind direction.
The resulting set of five metrics, including eigendirection, is useful for describing a twodimensional surface. However, based on the VIF analysis, it also represents a minimal set of indicators needed to describe such a surface since these metrics contain distinct and important information. Nevertheless, the development of practical parameterizations for sea ice and the MIZABL will ultimately need to include additional considerations, including the ease of obtaining these parameters for modeling applications, the computing time needed to calculate these surface metrics dynamically versus that needed to resolve surface features when running an ESM, and the availability of easiertocompute surrogate metrics.
The first step in answering this broad question is to investigate to what degree these other metrics affect the MIZABL compared to the firstorder effect of ice fraction on the MIZABL. In other words, given a specific ice fraction, how will changing any of the metrics in the resulting set affect the overlying MIZABL? While this study addressed this question for idealized surfaces and established the relevance of these parameters under certain conditions, determining which parameters will be critical for real ice maps – and understanding their frequency and impact – requires additional simulations, and a followup to this study is underway (see Fogarty et al., 2024). More turbulenceresolving numerical simulations of real ice surfaces are thus needed.
Another crucial step in answering this question involves figuring out how one would go about creating an accurate parameterization based on available external grid cell variables; the resources needed to answer this question may be extensive, especially considering that, in this age of machine learning, highresolution synthetic satellite imagery is generated more frequently (see AuBoehm et al., 2024). Even more largeeddy simulations over real sea ice maps, beyond what has been discussed here, are imperative for answering this question. Lastly, this open question also requires us to look at how to incorporate the resulting metrics and eigendirections into these climate models. This includes examining (i) how the geostrophic wind (at certain principal directions) interacts with the resulting metric set, (ii) other possible metrics of anisotropy, and (iii) how a sea ice model in an ESM can provide the data needed to capture the heterogeneity of the sea ice surface, as well as other questions that remain unanswered at this time.
In this study, incompressible filtered Navier–Stokes equations (along with a Boussinesq approximation for the mean state) and a heat budget are solved for a horizontally periodic flow, where a variable with a tilde represents a quantity filtered via the numerical grid spacing (Δ). The equations are expressed as follows:
The equations above use the Einstein summation rule, where i represents the free index and j represents the repeated index. Moreover, u_{i} represents the velocity vector; x_{i} represents the position vector; p represents modified pressure (see BouZeid et al., 2005, for details); θ represents the potential temperature;θ_{r} and $\widehat{\mathit{\theta}}$ represent the Boussinesq reference (the planar mean in our calculations) and the perturbation from said reference for potential temperature, respectively; ρ_{r} represents the reference mean density corresponding to θ_{r}; and F_{i} represents the main flowdriving force (a synoptic pressure gradient). The Coriolis force is represented by the third term on the righthand side of Eq. (A2), where f_{c} is the Coriolis parameter and ϵ_{ij3} represents the LeviCivita symbol. Buoyancy is represented by the fourth term on the righthand side of Eq. (A2), where δ_{ij} represents the Kronecker delta.
An overbar denotes averaging over time, which is used as a surrogate for ensemble Reynolds averaging, while spatial averaging over the heterogeneous domain (in both x and y) is denoted by angle brackets. The subgridscale stress (${\mathit{\tau}}_{ij}=\stackrel{\mathrm{\u0303}}{{u}_{i}{u}_{j}}{\stackrel{\mathrm{\u0303}}{u}}_{i}{\stackrel{\mathrm{\u0303}}{u}}_{j}$) and buoyancy flux (${q}_{j}=\stackrel{\mathrm{\u0303}}{{u}_{j}\mathit{\theta}}{\stackrel{\mathrm{\u0303}}{u}}_{j}\stackrel{\mathrm{\u0303}}{\mathit{\theta}}$), which result from the filtering, are modeled using a scaledependent Lagrangian dynamic model (BouZeid et al., 2005) with a constant subgridscale Prandtl number (Pr) of 0.4. As noted before, the numerical grid serves as the inherent filter of the model, but any explicit filtering needed to compute the dynamic Smagorinsky constant (c_{s}) is done at scales of 2Δ and 4Δ (2Δ for the local wall model); for these computations, a sharp spectral cutoff filter is used. This model was validated by BouZeid et al. (2005) for boundary layer flows over both homogeneous and heterogeneous terrain by reproducing experimental velocity and stress profiles obtained by Bradley (1968) following a change in surface roughness. It was further validated for urban flows (Tseng et al., 2006; Li et al., 2016), as well as for both stable and unstable boundary layers (Kleissl et al., 2006; Kumar et al., 2006; Huang and BouZeid, 2013). Therefore, the ability of this model to successfully capture the impacts of stability and spatial transitions in surface properties was not tested further in this paper.
The LES employs boundary conditions that are periodic in the horizontal direction, with zero vertical velocity at the top and bottom of the domain, as well as a stressfree top lid (${\partial}_{z}{u}_{i}=\mathrm{0}$, where $i=\mathrm{1},\mathrm{2}$) with zero heat flux. These conditions mimic a very strong top inversion and are adequate for our setup since the top of the domain is not stably stratified, eliminating the need for a sponge layer that prevents wave reflection. This allows for the surface characteristics to be isolated from z_{i} and the inversion strength. The indices ($i=\mathrm{1},\mathrm{2},\mathrm{3}$) represent the x, y, and z directions, oriented along the streamwise, crossstream, and vertical directions, respectively. At the bottom of the domain, the surface stress and heat flux are computed by a wall model based on a local lawofthewall formulation (BouZeid et al., 2005), incorporating a Monin–Obukhov buoyancy correction. Numerically, a pseudospectral approach is employed in the horizontal direction, and an explicit secondorder centereddifference scheme is used in the vertical direction. Time advancement is carried out using the fully explicit secondorder Adams–Bashforth scheme. Dealiasing of the convective terms is performed using the “$\mathrm{3}/\mathrm{2}$ rule” (Orszag, 1971). Pressure is computed from a Poisson equation obtained by taking the divergence of the momentum equation and applying the incompressibility assumption.
The initial potential temperature was chosen so that the mean heat flux over the entire domain was zero; in other words, the heat flux going into the ice was equal in magnitude to the heat flux coming from the water, based on the area fraction of the domain (in this case, the ice fraction). Thus, the initial potential air temperature of the largeeddy simulation was chosen to ensure that the icefractionweighted heat flux over the ice (f_{i}H_{i}) was equivalent to the waterfractionweighted heat flux over the ocean (f_{w}H_{w}), resulting in the following equation:
where ${f}_{\mathrm{i}}+{f}_{\mathrm{w}}=\mathrm{1}$. Monin–Obukhov flux profile relations were used to express the surface fluxes as follows:
where θ_{a} represents the bulk potential air temperature; θ_{i} represents the temperature of the ice surface; θ_{w} represents the temperature of the ocean surface; κ≈0.4 represents the von Kármán constant; ${u}_{*,\mathrm{i}}$ and ${u}_{*,\mathrm{w}}$ represent the friction velocities of the ice and water surfaces, respectively; ρ represents the density of air; c_{p} represents the specific heat of air; z represents the height near the surface (taken at z=50 m); z_{0h,i} and z_{0h,w} represent the scalar roughness lengths of the ice and water surfaces, respectively; L_{i} and L_{w} represent the Obukhov lengths over the ice and water surfaces, respectively; and Ψ_{s} and Ψ_{u} represent the stable and unstable correction functions, respectively, as described in Brutsaert (2005). These Obukhov lengths are defined as
where the values of H_{i} and H_{w} in these equations are taken as firstorder estimates obtained by rearranging Eqs. (B2) and (B3) without the stability functions (i.e., assuming a neutral atmosphere). H_{i} and H_{w} are given as
This allows one to derive a function by substituting Eqs. (B2) and (B3) into Eq. (B1), resulting in a function that can be solved for θ_{a} via a numerical rootfinding approach. Of course, this is only a first guess for initializing the LES, which will then dynamically create its potential air temperature field during the warmup period. As noted in the main text, the actual domainaveraged heat flux in the LES will not be zero.
Table C1 lists the landscape metrics used in the VIF analysis conducted in Sect. 4. For more information on each individual metric (apart from ice fraction), consult the FRAGSTATS manual (McGarigal and Marks, 1995).
A dataset containing the simulation results for the five patterns and the FRAGSTATS output for the sea ice maps are publicly available at https://doi.org/10.34770/5x2y5485 (Fogarty and BouZeid, 2023b). FRAGSTATS is publicly available for download at https://doi.org/10.2737/PNWGTR351 (McGarigal and Marks, 1995), and the preprocessed sea ice maps are available at https://doi.org/10.7265/N5PK0D32 (Fetterer et al., 2008).
The supplement related to this article is available online at: https://doi.org/10.5194/tc1843352024supplement.
JF: conceptualization, methodology, formal analysis, investigation, writing (original draft), and visualization. EBZ: methodology, investigation, formal analysis, writing (review and editing), supervision, and funding acquisition. MB: writing (review and editing) and validation. LB: writing (review and editing) and validation. All authors reviewed the results and approved the final version of the paper.
The contact author has declared that none of the authors has any competing interests.
The statements, findings, conclusions, and recommendations are those of the authors and do not necessarily reflect the views of the National Oceanic and Atmospheric Administration or the US Department of Commerce.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
We would like to acknowledge the support provided by NCAR's Computational and Information Systems Laboratory for highperformance computing on Cheyenne (Computational and Information Systems Laboratory, 2019), sponsored by the National Science Foundation as part of projects UPRI0007 and UPRI0021.
This research has been supported by the National Science Foundation (grant no. AGS 2128345) and the US Department of Commerce (grant no. NA18OAR4320123).
This paper was edited by Michel Tsamados and reviewed by Ian Brooks and Christof Lüpkes.
Allouche, M., Katul, G. G., Fuentes, J. D., and BouZeid, E.: Probability law of turbulent kinetic energy in the atmospheric surface layer, Phys. Rev. Fluids, 6, 074601, https://doi.org/10.1103/PhysRevFluids.6.074601, 2021. a
Allouche, M., BouZeid, E., and Iipponen, J.: The Influence of Synoptic Wind on LandSea Breezes, Q. J. Roy. Meteor. Soc., 149, 3198–3219, https://doi.org/10.1002/qj.4552, 2023a. a, b, c, d, e
Allouche, M., BouZeid, E., and Iipponen, J.: Unsteady LandSea Breeze Circulations in the Presence of a Synoptic Pressure Forcing, ESS Open Archive [preprint], https://doi.org/10.22541/essoar.170542134.41279506/v1, 2023b. a
Anderson, W., Barros, J. M., Christensen, K. T., and Awasthi, A.: Numerical and experimental study of mechanisms responsible for turbulent secondary flows in boundary layer flows over spanwise heterogeneous roughness, J. Fluid Mech., 768, 316–347, https://doi.org/10.1017/jfm.2015.91, 2015. a
Andreas, E. L., Horst, T. W., Grachev, A. A., Persson, P. O. G., Fairall, C. W., Guest, P. S., and Jordan, R. E.: Parametrizing turbulent exchange over summer sea ice and the marginal ice zone, Q. J. Roy. Meteor. Soc., 136, 927–943, https://doi.org/10.1002/qj.618, 2010. a
AuBoehm, C., Tsamados, M., Manescu, P., and Takao, S.: ARISGAN: Extreme SuperResolution of Arctic Surface Imagery using Generative Adversarial Networks, Front. Remote Sens. [preprint], 5, 1417417, https://doi.org/10.3389/frsen.2024.1417417, 2024. a
Baidya Roy, S.: Impact of land use/land cover change on regional hydrometeorology in Amazonia, J. Geophys. Res., 107, 8037, https://doi.org/10.1029/2000JD000266, 2002. a
Bates, N. R., Moran, S. B., Hansell, D. A., and Mathis, J. T.: An increasing CO_{2} sink in the Arctic Ocean due to sea–ice loss, Geophys. Res. Lett., 33, L23609, https://doi.org/10.1029/2006gl027028, 2006. a
Blackadar, A. K.: Boundary Layer Wind Maxima and Their Significance for the Growth of Nocturnal Inversions, B. Am. Meteorol. Soc., 38, 283–290, https://doi.org/10.1175/1520047738.5.283, 1957. a
Boudreault, L.É., Dupont, S., Bechmann, A., and Dellwik, E.: How Forest Inhomogeneities Affect the Edge Flow, Bound.Lay. Meteorol., 162, 375–400, https://doi.org/10.1007/s1054601602025, 2017. a
Bourassa, M. A., Gille, S. T., Bitz, C., Carlson, D., Cerovecki, I., Clayson, C. A., Cronin, M. F., Drennan, W. M., Fairall, C. W., Hoffman, R. N., Magnusdottir, G., Pinker, R. T., Renfrew, I. A., Serreze, M., Speer, K., Talley, L. D., and Wick, G. A.: HighLatitude Ocean and Sea Ice Surface Fluxes: Challenges for Climate Research, B. Am. Meteorol. Soc., 94, 403–423, https://doi.org/10.1175/BAMSD1100244.1, 2013. a
BouZeid, E., Meneveau, C., and Parlange, M. B.: Largeeddy simulation of neutral atmospheric boundary layer flow over heterogeneous surfaces: Blending height and effective surface roughness: LES OF NEUTRAL ATMOSPHERIC BOUNDARY LAYER FLOW, Water Resour. Res., 40, W02505, https://doi.org/10.1029/2003WR002475, 2004. a, b
BouZeid, E., Meneveau, C., and Parlange, M.: A scaledependent Lagrangian dynamic model for large eddy simulation of complex turbulent flows, Phys. Fluids, 17, 025105, https://doi.org/10.1063/1.1839152, 2005. a, b, c, d
BouZeid, E., Parlange, M. B., and Meneveau, C.: On the Parameterization of Surface Roughness at Regional Scales, J. Atmos. Sci., 64, 216–227, https://doi.org/10.1175/JAS3826.1, 2007. a
BouZeid, E., Anderson, W., Katul, G. G., and Mahrt, L.: The Persistent Challenge of Surface Heterogeneity in BoundaryLayer Meteorology: A Review, Bound.Lay. Meteorol, 177, 227–245, https://doi.org/10.1007/s10546020005518, 2020. a, b
Bradley, E. F.: A micrometeorological study of velocity profiles and surface drag in the region modified by a change in surface roughness, Q. J. Roy. Meteor. Soc., 94, 361–379, https://doi.org/10.1002/qj.49709440111, 1968. a
Brasseur, J. G. and Wei, T.: Designing largeeddy simulation of the turbulent boundary layer to capture lawofthewall scaling, Phys. Fluids, 22, 1–21, https://doi.org/10.1063/1.3319073, 2010. a
Brunsell, N. A., Mechem, D. B., and Anderson, M. C.: Surface heterogeneity impacts on boundary layer dynamics via energy balance partitioning, Atmos. Chem. Phys., 11, 3403–3416, https://doi.org/10.5194/acp1134032011, 2011. a
Brutsaert, W.: Hydrology: An Introduction, Cambridge University Press, https://doi.org/10.1017/CBO9780511808470, 2005. a, b
Casagrande, F., Stachelski, L., and De Souza, R. B.: Assessment of Antarctic sea ice area and concentration in Coupled Model Intercomparison Project Phase 5 and Phase 6 models, Int. J. Climatol., 43, 1314–1332, https://doi.org/10.1002/joc.7916, 2023. a
Computational and Information Systems Laboratory: Cheyenne: HPE/SGI ICE XA System (University Community Computing), Boulder, CO, National Center for Atmospheric Research, https://doi.org/10.5065/D6RX99HX, 2019. a
Courault, D., Drobinski, P., Brunet, Y., Lacarrere, P., and Talbot, C.: Impact of surface heterogeneity on a buoyancydriven convective boundary layer in light winds, Bound.Lay. Meteorol., 124, 383–403, https://doi.org/10.1007/s105460079172y, 2007. a
Cressie, N.: Statistics for Spatial Data  Revised Edition, Wiley Series in Probability and Statistics, John Wiley & Sons, Inc., ISBN 9781119115151, https://doi.org/10.1002/9781119115151, 1993. a
Crosman, E. T. and Horel, J. D.: Sea and lake breezes: A review of numerical studies, Bound.Lay. Meteorol., 137, 1–29, 2010. a
Cushman, S. A., McGarigal, K., and Neel, M. C.: Parsimony in landscape metrics: Strength, universality, and consistency, Ecol. Indic., 8, 691–703, https://doi.org/10.1016/j.ecolind.2007.12.002, 2008. a
de Vrese, P., Schulz, J.P., and Hagemann, S.: On the Representation of Heterogeneity in LandSurface–Atmosphere Coupling, Bound.Lay. Meteorol., 160, 157–183, https://doi.org/10.1007/s1054601601331, 2016. a, b
Docquier, D. and Koenigk, T.: Observationbased selection of climate models projects Arctic icefree summers around 2035, Commun. Earth Environ., 2, 144, https://doi.org/10.1038/s43247021002147, 2021. a
Dumont, D.: Marginal ice zone dynamics: history, definitions and research perspectives, Philos. T. Roy. Soc. A, 380, 20210253, https://doi.org/10.1098/rsta.2021.0253, 2022. a, b
Elvidge, A. D., Renfrew, I. A., Weiss, A. I., Brooks, I. M., LachlanCope, T. A., and King, J. C.: Observations of surface momentum exchange over the marginal ice zone and recommendations for its parametrisation, Atmos. Chem. Phys., 16, 1545–1563, https://doi.org/10.5194/acp1615452016, 2016. a, b
Elvidge, A. D., Renfrew, I. A., Brooks, I. M., Srivastava, P., Yelland, M. J., and Prytherch, J.: Surface Heat and Moisture Exchange in the Marginal Ice Zone: Observations and a New Parameterization Scheme for Weather and Climate Models, J. Geophys. Res.Atmos., 126, e2021JD034827, https://doi.org/10.1029/2021JD034827, 2021. a
Esau, I. N.: Amplification of turbulent exchange over wide Arctic leads: Large–eddy simulation study, J. Geophys. Res.Atmos., 112, D08109, https://doi.org/10.1029/2006jd007225, 2007. a
Essery, R. L. H., Best, M. J., Betts, R. A., Cox, P. M., and Taylor, C. M.: Explicit Representation of Subgrid Heterogeneity in a GCM Land Surface Scheme, J. Hydrometeorol., 4, 530–543, https://doi.org/10.1175/15257541(2003)004<0530:EROSHI>2.0.CO;2, 2003. a
Feltham, D. L.: Sea Ice Rheology, Annu. Rev. Fluid Mech., 40, 91–112, https://doi.org/10.1146/annurev.fluid.40.111406.102151, 2008. a
Fetterer, F. and Untersteiner, N.: Observations of melt ponds on Arctic sea ice, J. Geophys. Res., 103, 24821–24835, https://doi.org/10.1029/98JC02034, 1998. a
Fetterer, F., Wilds, S., and Sloan, J.: Arctic Sea Ice Melt Pond Statistics and Maps, 1999–2001, Version 1, National Snow & Ice Data Center [data set], https://doi.org/10.7265/N5PK0D32, 2008. a, b, c
Finnigan, J. J. and Shaw, R. H.: Doubleaveraging methodology and its application to turbulent flow in and above vegetation canopies, Acta Geophysica, 56, 534–561, https://doi.org/10.2478/s116000080034x, 2008. a
Fogarty, J. and BouZeid, E.: The Atmospheric Boundary Layer Above the Marginal Ice Zone: Scaling, Surface Fluxes, and Secondary Circulations, Bound.Lay. Meteorol., 189, 53–76, https://doi.org/10.1007/s1054602300825x, 2023a. a, b, c, d, e, f
Fogarty, J. and BouZeid, E.: LargeEddy Simulation and Statistical Metric Results for Patterned Sea Ice Surfaces, Princeton University [data set], https://doi.org/10.34770/5x2y5485, 2023b. a
Fogarty, J., BouZeid, E., Bushuk, M., Calaf, M., Allouche, M., and Ghannam, K.: Numerical Simulations of SatelliteSensed Surface Maps in the Marginal Ice Zone, https://doi.org/10.22541/essoar.172251979.90440727/v1, 2024. a
Ghannam, K. and BouZeid, E.: Baroclinicity and directional shear explain departures from the logarithmic wind profile, Q. J. Roy. Meteor. Soc., 147, 443–464, https://doi.org/10.1002/qj.3927, 2021. a
Gryschka, M., Gryanik, V. M., Lüpkes, C., Mostafa, Z., Sühring, M., Witha, B., and Raasch, S.: Turbulent Heat Exchange Over Polar Leads Revisited: A Large Eddy Simulation Study, J. Geophys. Res.Atmos., 128, e2022JD038236, https://doi.org/10.1029/2022JD038236, 2023. a, b
Herman, A., Wenta, M., and Cheng, S.: Sizes and Shapes of Sea Ice Floes Broken by Waves–A Case Study From the East Antarctic Coast, Front. Earth Sci., 9, 655977, https://doi.org/10.3389/feart.2021.655977, 2021. a
Horvat, C.: Marginal ice zone fraction benchmarks sea ice and climate model skill, Nat. Commun., 12, 2221, https://doi.org/10.1038/s41467021220047, 2021. a
Huang, H.Y., Margulis, S. A., Chu, C. R., and Tsai, H.C.: Investigation of the impacts of vegetation distribution and evaporative cooling on synthetic urban daytime climate using a coupled LESLSM model, Hydrol. Process., 25, 1574–1586, https://doi.org/10.1002/hyp.7919, 2011. a
Huang, J. and BouZeid, E.: Turbulence and Vertical Fluxes in the Stable Atmospheric Boundary Layer. Part I: A LargeEddy Simulation Study, J. Atmos. Sci., 70, 1513–1527, https://doi.org/10.1175/JASD120167.1, 2013. a
Hwang, B. and Wang, Y.: Multiscale satellite observations of Arctic sea ice: new insight into the life cycle of the floe size distribution, Philos. T. Roy. Soc. A, 380, 20210259, https://doi.org/10.1098/rsta.2021.0259, 2022. a
Ibidoja, O. J., Shan, F. P., Sulaiman, J., and Ali, M. K. M.: Detecting heterogeneity parameters and hybrid models for precision farming, J. Big Data, 10, 130, https://doi.org/10.1186/s40537023008108, 2023. a
Jaeger, J. A.: Landscape division, splitting index, and effective mesh size: new measures of landscape fragmentation, Landscape Ecol., 15, 115–130, https://doi.org/10.1023/A:1008129329289, 2000. a
Kleissl, J., Kumar, V., Meneveau, C., and Parlange, M. B.: Numerical study of dynamic Smagorinsky models in largeeddy simulation of the atmospheric boundary layer: Validation in stable and unstable conditions, Water Resour. Res., 42, W06D10, https://doi.org/10.1029/2005WR004685, 2006. a
Kumar, V., Kleissl, J., Meneveau, C., and Parlange, M. B.: Largeeddy simulation of a diurnal cycle of the atmospheric boundary layer: Atmospheric stability and scaling issues: LES OF A DIURNAL CYCLE OF THE ABL, Water Resour. Res., 42, https://doi.org/10.1029/2005WR004651, 2006. a
Kwok, R.: Declassified highresolution visible imagery for Arctic sea ice investigations: An overview, Remote Sens. Environ., 142, 44–56, https://doi.org/10.1016/j.rse.2013.11.015, 2014. a
Li, H. and Reynolds, J. F.: A Simulation Experiment to Quantify Spatial Heterogeneity in Categorical Maps, Ecology, 75, 2446, https://doi.org/10.2307/1940898, 1994. a
Li, H. and Reynolds, J. F.: On Definition and Quantification of Heterogeneity, Oikos, 73, 280, https://doi.org/10.2307/3545921, 1995. a
Li, Q. and BouZeid, E.: Contrasts between momentum and scalar transport over very rough surfaces, J. Fluid Mech., 880, 32–58, https://doi.org/10.1017/jfm.2019.687, 2019. a, b, c
Li, Q., BouZeid, E., Anderson, W., Grimmond, S., and Hultmark, M.: Quality and reliability of LES of convective scalar transfer at high Reynolds numbers, Int. J. Heat Mass Trans., 102, 959–970, https://doi.org/10.1016/j.ijheatmasstransfer.2016.06.093, 2016. a
Liu, C., Yang, Y., Liao, X., Cao, N., Liu, J., Ou, N., Allan, R. P., Jin, L., Chen, N., and Zheng, R.: Discrepancies in Simulated Ocean Net Surface Heat Fluxes over the North Atlantic, Adv. Atmos. Sci., 39, 1941–1955, https://doi.org/10.1007/s0037602213607, 2022. a
Lu, J., Nazarian, N., Hart, M. A., Krayenhoff, E. S., and Martilli, A.: Representing the effects of building height variability on urban canopy flow, Q. J. Roy. Meteor. Soc., 150, 46–67, https://doi.org/10.1002/qj.4584, 2023. a
Lüpkes, C., Gryanik, V. M., Witha, B., Gryschka, M., Raasch, S., and Gollnik, T.: Modeling convection over arctic leads with LES and a noneddyresolving microscale model, J. Geophys. Res.Oceans, 113, C09028, https://doi.org/10.1029/2007JC004099, 2008. a, b, c
Lüpkes, C., Gryanik, V. M., Hartmann, J., and Andreas, E. L.: A parametrization, based on sea ice morphology, of the neutral atmospheric drag coefficients for weather prediction and climate models, J. Geophys. Res.Atmos., 117, D13112, https://doi.org/10.1029/2012JD017630, 2012. a, b, c
Mahrt, L.: Surface heterogeneity and vertical structure of the boundary layer, Bound.Lay. Meteorol., 96, 33–62, 2000. a, b, c
Mandelbrot, B.: How Long Is the Coast of Britain? Statistical SelfSimilarity and Fractional Dimension, Science, 156, 636–638, https://doi.org/10.1126/science.156.3775.636, 1967. a
Mandelbrot, B. B.: The fractal geometry of nature, Freeman, San Francisco, CA, https://cds.cern.ch/record/98509 (last access: 5 December 2023), 1982. a
Margairaz, F., Pardyjak, E. R., and Calaf, M.: Surface Thermal Heterogeneities and the Atmospheric Boundary Layer: The Relevance of Dispersive Fluxes, Bound.Lay. Meteorol., 175, 369–395, https://doi.org/10.1007/s1054602000509w, 2020. a
Maronga, B. and Raasch, S.: LargeEddy Simulations of Surface Heterogeneity Effects on the Convective Boundary Layer During the LITFASS2003 Experiment, Bound.Lay. Meteorol., 146, 17–44, https://doi.org/10.1007/s105460129748z, 2013. a
McGarigal, K. and Marks, B. J.: FRAGSTATS: spatial pattern analysis program for quantifying landscape structure., Tech. Rep. PNWGTR351, U.S. Department of Agriculture, Forest Service, Pacific Northwest Research Station, Portland, OR [data set], https://doi.org/10.2737/PNWGTR351, 1995. a, b, c, d
Michaelis, J. and Lüpkes, C.: The Impact of Lead Patterns on Mean Profiles of Wind, Temperature, and Turbulent Fluxes in the Atmospheric Boundary Layer over Sea Ice, Atmosphere, 13, 148, https://doi.org/10.3390/atmos13010148, 2022. a
Michaelis, J., Lüpkes, C., Zhou, X., Gryschka, M., and Gryanik, V. M.: Influence of Lead width on the Turbulent Flow Over Sea Ice Leads: Modeling and Parametrization, J. Geophys. Res.Atmos., 125, e2019JD031996, https://doi.org/10.1029/2019JD031996, 2020. a
Michaelis, J., Lüpkes, C., Schmitt, A. U., and Hartmann, J.: Modelling and parametrization of the convective flow over leads in sea ice and comparison with airborne observations, Q. J. Roy. Meteor. Soc., 147, 914–943, https://doi.org/10.1002/qj.3953, 2021. a, b
Miles, J.: Tolerance and Variance Inflation Factor, John Wiley & Sons, Ltd, ISBN 9781118445112, https://doi.org/10.1002/9781118445112.stat06593, 2014. a
Miller, N. B., Shupe, M. D., Cox, C. J., Noone, D., Persson, P. O. G., and Steffen, K.: Surface energy budget responses to radiative forcing at Summit, Greenland, The Cryosphere, 11, 497–516, https://doi.org/10.5194/tc114972017, 2017. a
Moltchanov, S., BohbotRaviv, Y., Duman, T., and Shavit, U.: Canopy edge flow: A momentum balance analysis, Water Resour. Res., 51, 2081–2095, https://doi.org/10.1002/2014WR015397, 2015. a
Momen, M. and BouZeid, E.: LargeEddy Simulations and DampedOscillator Models of the Unsteady Ekman Boundary Layer, J. Atmos. Sci., 73, 25–40, https://doi.org/10.1175/JASD150038.1, 2016. a
Moody, A. and Woodcock, C.: Scaledependent errors in the estimation of landcover proportions: Implications for global landcover datasets, Photogramm. Eng. Remote Sens., 60, 585–594, 1994. a
Moody, A. and Woodcock, C. E.: The influence of scale and the spatial characteristics of landscapes on landcover mapping using remote sensing, Landscape Ecol., 10, 363–379, https://doi.org/10.1007/BF00130213, 1995. a
Myksvoll, M. S., Britt Sandø, A., Tjiputra, J., Samuelsen, A., Çaǧlar Yumruktepe, V., Li, C., Mousing, E. A., Bettencourt, J. P., and Ottersen, G.: Key physical processes and their model representation for projecting climate impacts on subarctic Atlantic net primary production: A synthesis, Prog. Oceanogr., 217, 103084, https://doi.org/10.1016/j.pocean.2023.103084, 2023. a
Nilsson, E. D., Rannik, É., Swietlicki, E., Leck, C., Aalto, P. P., Zhou, J., and Norman, M.: Turbulent aerosol fluxes over the Arctic Ocean: 2. Winddriven sources from the sea, J. Geophys. Res.Atmos., 106, 32139–32154, https://doi.org/10.1029/2000JD900747, 2001. a
Notz, D. and Community, S.: Arctic Sea Ice in CMIP6, Geophys. Res. Lett., 47, e2019GL086749, https://doi.org/10.1029/2019GL086749, 2020. a
Notz, D. and Stroeve, J.: The Trajectory Towards a Seasonally IceFree Arctic Ocean, Current Climate Change Reports, 4, 407–416, https://doi.org/10.1007/s4064101801132, 2018. a
Omidvar, H., BouZeid, E., Li, Q., Mellado, J.P., and Klein, P.: Plume or bubble? Mixedconvection flow regimes and cityscale circulations, J. Fluid Mech., 897, A5, https://doi.org/10.1017/jfm.2020.360, 2020. a
Orszag, S. A.: On the Elimination of Aliasing in FiniteDifference Schemes by Filtering HighWavenumber Components, J. Atmos. Sci., 28, 1074, https://doi.org/10.1175/15200469(1971)028<1074:OTEOAI>2.0.CO;2, 1971. a
O'Brien, R. M.: A Caution Regarding Rules of Thumb for Variance Inflation Factors, Quality & Quantity, 41, 673–690, https://doi.org/10.1007/s1113500690186, 2007. a
Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., and Duchesnay, E.: Scikitlearn: Machine Learning in Python, J. Mach. Learn. Res., 12, 2825–2830, 2011. a
Persson, P. O. G., Fairall, C. W., Andreas, E. L., Guest, P. S., and Perovich, D. K.: Measurements near the Atmospheric Surface Flux Group tower at SHEBA: Nearsurface conditions and surface energy budget, J. Geophys. Res.Oceans, 107, SHE 211–SHE 2135, https://doi.org/10.1029/2000JC000705, 2002. a
Pickett, S. T. A. and Cadenasso, M. L.: Landscape Ecology: Spatial Heterogeneity in Ecological Systems, Science, 269, 331–334, https://doi.org/10.1126/science.269.5222.331, 1995. a
Piomelli, U. and Balaras, E.: WallLayer Models for LargeEddy Simulations, Annu. Rev. Fluid Mech., 34, 349–374, https://doi.org/10.1146/annurev.fluid.34.082901.144919, 2002. a
Porson, A., Steyn, D. G., and Schayes, G.: Seabreeze scaling from numerical model simulations, part II: Interaction between the sea breeze and slope flows, Bound.Lay. Meteorol., 122, 31–41, 2007. a
Ramudu, E., Gelderloos, R., Yang, D., Meneveau, C., and Gnanadesikan, A.: Large Eddy Simulation of Heat Entrainment Under Arctic Sea Ice, J. Geophys. Res.Oceans, 123, 287–304, https://doi.org/10.1002/2017jc013267, 2018. a
Raupach, M. R. and Shaw, R. H.: Averaging Procedures for Flow Within Vegetation Canopies, Bound.Lay. Meteorol., 22, 79–90, https://doi.org/10.1007/BF00128057, 1982. a, b
Ren, H., Zhang, C., and Zhao, X.: Numerical simulations on the fracture of a sea ice floe induced by waves, Appl. Ocean Res., 108, 102527, https://doi.org/10.1016/j.apor.2021.102527, 2021. a
Riitters, K. H., O'Neill, R. V., Hunsaker, C. T., Wickham, J. D., Yankee, D. H., Timmins, S. P., Jones, K. B., and Jackson, B. L.: A factor analysis of landscape pattern and structure metrics, Landscape Ecol., 10, 23–39, https://doi.org/10.1007/BF00158551, 1995. a, b
Rosenblum, E. and Eisenman, I.: Faster Arctic Sea Ice Retreat in CMIP5 than in CMIP3 due to Volcanoes, J. Climate, 29, 9179–9188, https://doi.org/10.1175/JCLID160391.1, 2016. a
Rosenblum, E. and Eisenman, I.: Sea Ice Trends in Climate Models Only Accurate in Runs with Biased Global Warming, J. Climate, 30, 6265–6278, https://doi.org/10.1175/JCLID160455.1, 2017. a
Salesky, S. T., Calaf, M., and Anderson, W.: Unstable turbulent channel flow response to spanwiseheterogeneous heat fluxes: Prandtl's secondary flow of the third kind, J. Fluid Mech., 934, A46, https://doi.org/10.1017/jfm.2022.15, 2022. a
Seabold, S. and Perktold, J.: statsmodels: Econometric and statistical modeling with python, in: 9th Python in Science Conference, 28 June–23 July 2010, Austin Texas, 2010. a
Šímová, P. and Gdulová, K.: Landscape indices behavior: A review of scale effects, Appl. Geogr., 34, 385–394, https://doi.org/10.1016/j.apgeog.2012.01.003, 2012. a
Stoll, R., Gibbs, J. A., Salesky, S. T., Anderson, W., and Calaf, M.: LargeEddy Simulation of the Atmospheric Boundary Layer, Bound.Lay. Meteorol., 177, 541–581, https://doi.org/10.1007/s10546020005563, 2020. a
Stroeve, J., Holland, M. M., Meier, W., Scambos, T., and Serreze, M.: Arctic sea ice decline: Faster than forecast, Geophys. Res. Lett., 34, L09501, https://doi.org/10.1029/2007GL029703, 2007. a
Strong, C., Foster, D., Cherkaev, E., Eisenman, I., and Golden, K. M.: On the Definition of Marginal Ice Zone Width, J. Atmos. Ocean. Tech., 34, 1565–1584, https://doi.org/10.1175/JTECHD160171.1, 2017. a
Taylor, P. C., Hegyi, B. M., Boeke, R. C., and Boisvert, L. N.: On the Increasing Importance of AirSea Exchanges in a Thawing Arctic: A Review, Atmosphere, 9, 41, https://doi.org/10.3390/atmos9020041, 2018. a
Tetzlaff, A., Lüpkes, C., and Hartmann, J.: Aircraftbased observations of atmospheric boundarylayer modification over Arctic leads, Q. J. Roy. Meteor. Soc., 141, 2839–2856, https://doi.org/10.1002/qj.2568, 2015. a
Tseng, Y.H., Meneveau, C., and Parlange, M. B.: Modeling Flow around Bluff Bodies and Predicting Urban Dispersion Using Large Eddy Simulation, Environ. Sci. Technol., 40, 2653–2662, https://doi.org/10.1021/es051708m, 2006. a
Wang, Y., Holt, B., Erick Rogers, W., Thomson, J., and Shen, H. H.: Wind and wave influences on sea ice floe size and leads in the Beaufort and Chukchi Seas during the summer–fall transition 2014, J. Geophys. Res.Oceans, 121, 1502–1525, https://doi.org/10.1002/2015JC011349, 2016. a
Wenta, M. and Herman, A.: The influence of the spatial distribution of leads and ice floes on the atmospheric boundary layer over fragmented sea ice, Ann. Glaciol., 59, 213–230, https://doi.org/10.1017/aog.2018.15, 2018. a
Wenta, M. and Herman, A.: AreaAveraged Surface Moisture Flux over Fragmented Sea Ice: Floe Size Distribution Effects and the Associated Convection Structure within the Atmospheric Boundary Layer, Atmosphere, 10, 654, https://doi.org/10.3390/atmos10110654, 2019. a
Willingham, D., Anderson, W., Christensen, K. T., and Barros, J. M.: Turbulent boundary layer flow over transverse aerodynamic roughness transitions: Induced mixing and flow characterization, Phys. Fluids, 26, 025111, https://doi.org/10.1063/1.4864105, 2014. a
Wood, N. and Mason, P.: The influence of static stability on the effective roughness lengths for momentum and heat transfer, Q. J. Roy. Meteor. Soc., 117, 1025–1056, https://doi.org/10.1002/qj.49711750108, 1991. a
 Abstract
 Introduction
 Methods and data
 Results: the MIZABL across idealized configurations
 Results: statistical analysis
 Principal direction of the sea ice patterns
 Conclusions
 Appendix A: Details of the largeeddy simulation
 Appendix B: Temperature initialization
 Appendix C: Landscape quantification metrics
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement
 Abstract
 Introduction
 Methods and data
 Results: the MIZABL across idealized configurations
 Results: statistical analysis
 Principal direction of the sea ice patterns
 Conclusions
 Appendix A: Details of the largeeddy simulation
 Appendix B: Temperature initialization
 Appendix C: Landscape quantification metrics
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement