A computational fluid dynamics model for the small-scale dynamics of wave, ice floe and interstitial grease ice interaction

The marginal ice zone is a highly dynamical region where sea ice and ocean waves interact. Large-scale sea ice models only compute domain-averaged responses. As the majority of the marginal ice zone consists of mobile ice floes surrounded by grease ice, finer-scale modelling is needed to resolve variations of its mechanical properties, wave-induced pressure gradients and drag forces acting on the ice floes. A novel computational fluid dynamics approach is presented, that considers the heterogeneous sea ice material composition and accounts for the wave-ice interaction dynamics. Results show, after comparing three realistic sea ice layouts with similar concentration and floe diameter, that the discrepancy between the domain-averaged temporal stress and strain rate evolutions increases for decreasing wave period. Furthermore, strain rate and viscosity are mostly affected by the variability of ice floe shape and diameter.


Introduction
The region where ocean processes affect sea ice, known as marginal ice zone (MIZ; [46]), is a highly complex system [27,3] particularly during the ice formation and melt season. Most contemporary dynamic-thermodynamic sea ice models used for predicting global climate are phenomenological large-scale models (order of 100 km [12]) in which internal stresses are related to the strain rate via a viscous-plastic rheology [17,21,47]. These models adopt a smeared approach in which an area with heterogeneous ice characteristics (fractures, leads, open water and ice type) is modelled as a single homogeneous material with averaged properties, and in which ice-ocean interactions are parameterised [34]. The dynamics is governed by the momentum balance equation and two continuity equations controlling ice thickness and concentration, accounting for deformation and growth-related effects [22]. However, sea ice deformation cannot be exclusively described by a viscous-plastic rheology as sea ice drift on the scale of less than 10 km is only accurate to a certain extent and fails to reproduce sea ice deformation at finer scales [36,10].
The MIZ is regulated by both large-and small-scale processes. Large-scale processes, on a spatial scale of 100 km and a temporal scale of several days, are the Coriolis force and large-scale ocean circulation [31,30]. Fine-scale processes, on the scale of hours and less than 10 kilometers, are the variability in the sea ice cover and stress acting at the air-ocean-ice interface [8,25]. Among small-scale processes in the Antarctic MIZ, waves play a crucial role in the formation of sea ice [40,4]. Furthermore, the propagation of waves from open ocean into the MIZ and wave-ice interaction are linked to wave scattering and dissipation through the momentum transfer to ice floes, which mainly depend on sea ice characteristics like ice floe geometry and floe size distribution [45]. In particular, at the ice edge, close to the open ocean, high energy incident waves create complicated sea states resulting in violent ice floe motion and collision [42,40]. During the Antarctic winter, when temperatures are low, sea ice initially forms as small circular floes [39,2,1] (1 − 10m in diameter, smaller than the characteristic wavelength), known as pancakes, surrounded by grease ice under the action of ocean waves penetrating deep into the ice covered ocean [23,43].
Dynamic-thermodynamic sea ice models on a floe-scale are scarce [10], in particular with respect to the characteristic pancake ice floes found in the winter MIZ. Most of existing fine-scale numerical models are based on a molecular dynamics schemes based on Hertzian collision dynamics [39,19,14] that consider sea ice physics and dynamics at a floe level [13]. The Discrete-Element bonded-particle model by Herman et al. [16] highlights the significance of skin drag on wave attenuation and floe collision dynamics due to prolonged collision and reduced restitution coefficient. The model uses an heuristic contact detection algorithm and includes simplified overwash (water overflowing the ice [7,41,32]), and elastic and inelastic contact force contributions, with the latter linked to the restitution coefficient; tangential friction is not accounted for. Damsgaard et al. [9] uses a discrete element framework for the approximation of Lagrangian sea ice dynamics also at the floe level. The modelling of fragmented sea ice in the MIZ is improved by using exact solutions for mechanical nonlinearities with arbitrary sea ice concentrations. Rabatel et al. [35] describes the dynamics of an assembly of rigid ice floes using an event-driven algorithm [29]. Their approach focuses on collisions of individual ice floes of both, arbitrary size and shape derived from satellite images from the Arctic. Ice classified as "new ice" of recently frozen sea water that is not solid ice yet, such as grease ice, has not been considered in these models. Moreover, with the exception of Rabatel et al. [35], all other models have in common the limitation of using highly simplified ice floe shapes. This paper introduces a novel computational framework, developed using the computational fluid dynamics (CFD) software OpenFOAM , for the small-scale wave-ice interaction dynamics. The governing equations for the small-scale model are introduced in Section 2 and the numerical implementation in Section 3. We focus on the pancake and grease ice rheology variables over short time periods (¡ 1 minute), during which thermodynamic effects on ice thickness and concentration are negligible and ice floe ridging can be disregarded. We demonstrate the suitability of the proposed computational framework for modelling the elastic collision dynamics of freefloating ice floes embedded in interstitial grease ice under the action of ocean waves by providing a convergence analysis with regards to the domain grid refinement and the sea ice layout in Section 4. We show that the dynamics depends on (i) ice characteristics, (ii) sea ice rheology and (iii) external forcing. With these findings at hand, an initial investigation of realistic ice floe distributions is undertaken in Section 5, focusing on the inter-dependency of floe size, wave forcing and grease ice viscosity. Final remarks and conclusions are presented in Section 6.
2 Small-scale model

Momentum equation
The sea ice dynamics is governed by the momentum balance equation where U represents the sea ice velocity vector and m the mass of ice per area where ρi is the ice density and h the ice thickness. The internal reaction forces are represented by the Cauchy stress tensor, σ = {σgrease, σ f loe }, depending on the ice constituent. In-plane wind and ocean current stress vectors applied to the ice, i.e. the external forcing, are represented by τa and τo, respectively. The in-plane stress due to the waves, τw (also an external forcing), is derived from the linear wave theory [18] and consists of two components: where τ sd is the viscous component representing the skin drag acting on the ice floe basal plane and τ f k the Froude-Krylov force due to wave pressure field acting on the ice floe circumference.
The skin drag is given as where ρw and θw represent the water density and the ice-ocean turning angle, respectively, k a unit normal vector to the surface of the ice and the ice-ocean drag coefficient, Cw, is composed of different values for pancake ice floes and grease ice, which can be found in Table 1. The orbital velocity of the water, Uw, for monochromatic waves is defined as where the x is taken along the main direction of wave propagation and z in the vertical direction, such that the y-velocity component can be assumed zero. In Equations (5) and (7) a, ω and k represent the wave amplitude, wave frequency and wave number, respectively. The wave frequency and wave number, ω = 2π/T and k = 2π/λ, are computed from the wave period T and the wavelength λ using the deep water dispersion relation ω 2 = gk.
The Froude-Krylov force, τ f k , accounts for the horizontal surge force due to the wave-induced pressure acting at the interface between ice floe and water [15] where hw represents the submerged ice floe thickness portion, calculated assuming that 90% of both ice types is submerged. The unit vector n acts normal to the circumference of the ice floes directed outwards. The wave-induced pressure, p, is written as where the water density and gravitational acceleration are denoted as ρw and g, respectively.
The form drag acting on the ice floe circumference due velocity differences of floes and surrounding grease ice is implicitly included by the continuum model comprising both ice constituents enforcing velocity continuity.
Ocean current drag, τo, at the basal plane of the ice due to relative velocity of the ocean current and the ice is given by whereas the wind drag, τa, on the apical plane by The variables Ua and Uo represent the velocities of wind and ocean boundary layers, respectively. The air density, ice-air drag coefficient, and wind turning angle are indicated by ρa, Ca and θa, respectively.

Sea ice rheology
The rheology in the model is based on the material characteristics of each constituent, grease ice and pancake floes, respectively, and each is represented by its own rheology. Grease ice behaves like a fluid whereas ice floes have a solid-like behaviour. The ice thickness of both is assumed spatially constant, since ice thickness is negligibly small compared to the domain size in lateral direction. As mentioned before, since only small time periods are considered in the model, ice thickness is assumed unaffected by thermodynamic and rafting-related growth and does not vary in time.
The sea ice strain rate tensor,˙ , is used in both grease ice and ice floe rheology, and is written in terms of the sea ice velocity vector, U , as˙

Grease ice rheology
The fluid-like viscous-plastic behaviour of grease ice is modelled applying a viscous-plastic flow rheology similar as proposed by e.g. Thorndike and Rothrock [44] or Hibler [17] where P represents the internal grease ice strength that controls its compressibility, ζ and η are strain ratedependent grease ice viscosities, namely the spherical and deviatoric contributions, respectively. The main difference to the above mentioned models is the definition of the ice strength parameter as where P * represents an empirical constant excluding ice growth and concentration effects. Both strain ratedependent viscosities are coupled via where the ratio between the in-plane principle axes of the elliptical yield curve is given by parameter, e, and the effective strain rate, ∆, is given by The Cartesian components of the strain rate tensor are denoted by˙ 11,˙ 22 and˙ 12. As the strain rate approaches zero, the viscosity tends to infinity. To address this, a lower limit of the effective strain rate is enforced, i.e. ∆ = 2 · 10 −7 s −1 [24].
Equation (12) can be, after substitution of Equation (11), written in terms of the velocity vector, U , as

Ice floe rheology
The solid-like material response of ice floes exhibits relatively small deformations and its constitutive behaviour floes is therefore described using the generalised Hooke's law where µ and λ represent the Lamé constants. Note as Hooke's law is a linear function of the strain, , Equation (11) is substituted in the linearised equation to write it in terms of the velocity vector, U , resulting in where n represents the time index discretised with time step ∆t [20].

Numerical implementation
The framework described in Section 2 is implemented in the Computational Fluid Dynamics (CFD) software OpenFOAM, which is based on the finite volume method (FVM) where a continuum body is approximated by a discrete model [11]. The domain is subdivided, by means of a mesh, into a finite number of small control volumes (cells). The volume of fluid (VOF) method is used in OpenFOAM to distinguish between two immiscible and incompressible fluids based on the non-dimensional parameter α [37,38]. The VOF method is complemented with an interface compression scheme [37] that enforces a sharp interface between dominantly solid-like pancake ice floes and viscous fluid-like interstitial grease ice. Cells containing exclusively pancake ice or grease ice have α values of 1 and 0, respectively, and all intermediate values (0 < α < 1) have no physical interpretation but are needed in the FVM which does not allow for discontinuities. The region of 0 < α < 1 should be numerically constrained to the thin interface between the two ice materials. We note that as both ice floes and grease ice are modelled in a continuum fashion, the floe-floe and the floe-grease ice interactions are naturally accounted for in terms of their normal and tangential force components.
Ice floes and interstitial grease ice are discretised in OpenFOAM with finite volume elements using the leastsquares gradient scheme and the Gauss linear divergence scheme for the spatial integration. The Euler implicit method is used for temporal discretisation.
The explicit influence of the grease ice viscosity on wave attenuation is commonly not accounted for, because dissipation per wavelength is small (< 0.5%) [4], i.e. the wave amplitude remains almost constant within a small domain. The domain boundary is numerically approximated by horizontal zero-gradient boundary conditions. As the wave forcing is imposed, the boundary condition introduces an error but only to cells directly adjacent to the boundary itself but does not affect the velocity field in the inner domain. To overcome this inconsistency the domain is enlarged to exclude undesired boundary effects for the inner domain. All terms in the momentum balance equation, Equation (1), are normalized by the respective height of pancake and grease ice, resulting in stress values in terms of kg s −2 . Table 1 shows the used simulation parameters. The ice-ocean drag coefficients are based on values found in [26,28,6].
The convergence analysis with regard to grid size and domain size is conducted to find the optimal grid size to be used in the discretisation of the FVM (in Section 4.1) and the critical ratio of floe diameter to domain size as linked to the imposed wave forcing (in Section 4.2). Simplified ice layouts are used consisting of randomly distributed disk-shaped pancake ice floes of constant diameter surrounded by grease ice with a temporally and spatially varying viscosity ν ≈ 0.04m 2 s −1 . In the grid size convergence analysis the domain is exposed to a harmonic propagating wave with period T = 18s (λ = 506m) and amplitude a = 4.8m, typically encountered in the Antarctic MIZ [3]. In the domain size convergence analysis three different wave forcing scenarios are considered with periods T = 6s, T = 12s and T = 20s (corresponding to λ = 56m, λ = 225m and λ = 625m) to cover the majority of occurring wave conditions and amplitudes a = 0.5m, a = 2.1m and a = 6m, respectively. The wave parameters are chosen to maintain a constant wave steepness (ak) across the considered cases corresponding to storm waves propagating into the MIZ to highlight their effect on heterogeneous sea ice conditions [4,5]. The time step ∆t = 0.01s for the Euler implicit method provides numerically stable simulations.
To study the high-resolution mechanical response of sea ice due to wave-ice interaction and ice composition in terms of ice type as well as floe shape and diameter in Section 5, three realistic 100 × 100m 2 -sea ice layouts are extracted from in-situ image and video material recorded by stereo cameras of the dynamics of sea ice in the Antarctic MIZ [2]. The images are collected in close vicinity to each other and, therefore, in homogeneous sea ice conditions with only slightly differing sea ice properties in terms of ice floe concentration and median ice floe caliper diameter. The different layouts also help to verify the robustness of the numerical framework with respect to natural variability of the ice floe layout in the MIZ.

.1 Grid size convergence analysis
To accurately resolve the interface of pancake floes and the surrounding grease ice, discretisations featuring a higher number of cells provide a sharper boundary as indicated by the white colour in the pancake-grease ice layouts depicted in Figure 1 at t = 0s. For each discretisation refinement level, the domain-averaged stress and strain rate magnitudes, σmag and˙ mag, as well as the bulk viscosity, ζ, are computed for comparison.
A 100 × 100m 2 -inner domain size with a 41% ice floe area fraction is embedded in a 300 × 300m 2 outer domain to exclude undesired boundary effects for the inner domain, as shown in Figure 1 by the black rectangle. The ice floes have a constant diameter of 20m. Multiple simulations are carried out, each time with an increasing number of cells in the domain ranging from 2500 to 160,000 cells.
(1) (2) (3) Figure 1: Three different grid refinement levels of a 300 × 300m 2 -outer domain of randomly distributed disk-shaped 20m-diameter pancake floes each with 100 × 100m 2 -inner domains of following discretisations: (1) 2500 cells with a 2 × 2m 2 cell size, (2) 10,000 cells with a 1 × 1m 2 cell size and (3) 160,000 cells with a 0.25 × 0.25m 2 cell size. Red and blue represent ice floes and grease ice, respectively. The numerical interface is indicated by the white colour. Figure 2 shows the domain-averaged stress, viscosity and strain rate results for the four different grid refinement levels within the 100 × 100m 2 -inner domain. The blue and red curves, representing 90,000 and 160,000cell-inner domains, respectively, show similar values indicating convergence. It can be concluded that a total of 90,000 cells for the inner domain, corresponding to a grid size of 0.33 × 0.33m 2 , is sufficient and will be used for all numerical simulations in the subsequent sections.

Domain size convergence analysis
The threshold for convergence with respect to the minimum required domain size is found by calculating the domain-averaged stress and strain magnitudes, σmag and˙ mag, as well as bulk viscosity, ζ, for increasing inner domain sizes. Once the average is unaffected by changes in the domain size, the transition from small-to largescale modelling is identified, and the link to a phenomenological model using homogeneous material properties can be established.
A total domain size of 3600 × 3600m 2 is considered, from which results are extracted for smaller inner domains, ranging from 50 × 50m 2 to 3200 × 3200m 2 . This ensures a boundary zone of 400m around the largest inner domain to exclude unwanted boundary effects due to the zero-gradient boundary conditions. Figure 3 shows the smallest four inner domain sizes ranging from 50 × 50m 2 to 400 × 400m 2 . With 20m-diameter ice floes the ice floe concentration from the largest domain size, 3600 × 3600m 2 , to the smallest domain size, 50 × 50m 2 , ranges from 40% to 33%, due to a randomly distributed sea ice layout. Figure 4 shows the overall domain-averaged sea ice rheology variables (stress, strain rate and bulk viscosity) in both ice floes and grease ice for different domain sizes and three different wave periods. Clearly, smaller wave periods reduce the required domain size for convergence. For the largest wave period of T = 20s, the threshold is identified at an inner domain size of 1600 × 1600m 2 , where the blue and red curves in Figures 4(c), (f) and (i), which represent an inner domain size of 1600 × 1600m 2 and 3200 × 3200m 2 , respectively, converge. In contrast, for the smallest wave periods T = 6s, the minimum required inner domain size is only 400 × 400m 2 .

Analysis of the sea ice layout and rheology
In this section the sea ice rheology variables are separately analysed for floes and grease ice. Section 4.2 indicates that any inner domain size smaller than 400 × 400m 2 results in a different mechanical response for any wave period larger than T = 6s. An inner domain size of 100 × 100m 2 is therefore expected to produce temporal and spatial fluctuations of the strain rate distribution depending on the ratio of median pancake floe diameter and wavelength, D p,median /λ, and thus, warrant detailed small-scale modelling. Three realistic layouts and one idealised sea ice configuration comparable in terms of sea ice concentration are illustrated in Figure 5. The ice floe caliper diameter (Dx and Dy; in x -and y-direction), standard deviation (SDx and SDy; in x -and y-direction), ice floe concentrations, wave characteristics and viscosity are summarized in Table 2 for all four layouts.
Layout differences in ice floe concentration, 57.9% (±3.2%), and median ice floe caliper diameter with 10.8m (±2.2m) in x -direction and 8.7m (±1.7m) in y-direction, respectively, are small. The distribution of the ice floe caliper diameters has been mapped for all realistic sea ice layouts through the use of a boxplot, as shown in Figure 6. Layout 2 has the smallest interquartile range. This is also reflected in the standard deviation, which show the lowest values in both x -and y-direction for layout 2. This indicates that it features the largest portion of medium-size floes and is the most homogeneous in size. Layout 3, on the other hand, exhibits the largest spread of floe sizes, in particular with regards to large floes.
As layout 2 is the most homogeneous, its mechanical response is expected to be closest to the idealised sea ice layout 4 and is therefore chosen to study the error introduced by completely disregarding floe shape and variations of floe diameter. In order to specifically focus on the sensitivity regarding those two sea ice (1) (2) (3) (4) Figure 5: Three realistic sea ice layouts (1-3) and one idealised with disk-shaped floes (4), each with a 100 × 100m 2 inner domain embedded in a 300 × 300m 2 outer domain.
characteristics and study their effect on stress, strain rate and viscosity variables in both the grease ice and ice floe rheology, the ice floe concentration is chosen to be the same and the mean ice floe caliper diameter in layout 4 has been derived from layout 2, such as that the average area per ice floe, A f loe = 293m 2 , is identical. All sea ice layouts in this analysis are subjected to three different wave forcings via Equation 3 using wave periods ranging between T = 8 − 16s and for constant wave steepness ak = 0.06 (by varying the wave amplitude). The domain-averaged grease ice viscosity value is ν ≈ 0.04m 2 s −1 in agreement with literature values [46,33,48]. The viscosity of the grease ice rheology, Equation (16), is strain rate-dependent via the ice strength parameter, P , which in turn depends on the empirical constant, P * . The value of P * is chosen such as that the domain-averaged viscosity provides a close match to the predefined values above. For each sea ice layout, three simulations of 30s are carried out according to the three aforementioned wave periods to elucidate the impact of wave length on the mechanical response.  Figure 7 shows the spatially-averaged stress and strain rate magnitudes in both grease ice and pancake ice floes for wave periods T = 8s, T = 12s and T = 16s with ν ≈ 0.04m 2 s −1 . In Figure 7(a-c) the spatially-averaged ice floe stress magnitude evolution over time is shown for all sea ice layouts. The average stress in ice floes is similar for all wave periods due to a prescribed wave steepness, ak = 0.06. An increasing wave period results in an increasing stress amplitude and a decreasing stress frequency. The discrepancy between the stress curves of the four considered layouts increases for a decreasing wave period signifying the mounting influence of the detailed heterogeneous sea dynamics description. Clearly, the realistic sea ice layout 1 exhibits a smaller floe stress magnitude compared to layouts 2 and 3, in particular for the smallest wave period, T = 8s, as illustrated in Figure 7(a). Additionally, it can be observed that the floe stress differences between layouts change with the wave period and that the grease ice stress fluctuations, shown in Figure 7(d-f), in particular for T = 8s, are considerably lower than seen for the floe stress. With regards to influence of floe shape and diameter variations, the floe stress and strain rate response of the idealized sea ice composition (layout 4) is distinctly different from the realistic layout 2, both being underestimated by the idealization of ice floe geometry. The average discrepancy in the ice floe stress curves between layout 2 and 4 is approximately 7%.
As mentioned previously, Figure 7(d-f) show the spatially-averaged grease ice stress for sea ice layouts 1-4. Results of the three realistic sea ice layouts (layout 1-3) are very similar for all wave periods, due to a comparable concentration of grease ice in all layouts. On the other hand, the uniformity of floe diameter and shape results in an average discrepancy in the grease ice stress between layout 2 and the idealized one (layout 4) of about 7%.
Lastly, Figure 7(g-l) show the spatially-averaged grease ice strain rate and bulk viscosity for all sea ice layouts. As for the flow stress, the discrepancy between the strain rate evolution in time increases for decreasing wave period looking at the realistic sea ice layouts and a distinctly smaller strain rate magnitude is exhibited for layout 1. Comparing the grease ice strain rate curves of layout 2 and 4, we clearly see that the average strain rate is lower for layout 4, as it was the case for floe and grease ice stress. The average discrepancy in the grease ice strain rate between layout 2 and 4 is substantial with approximately 103%. From Equation (14) we know that the relation between the grease ice strain rate and viscosity is inversely proportional. An increase in grease ice viscosity, results in a decreasing strain rate and vice versa. As these variables are directly related, the average discrepancy in grease ice viscosity between layout 2 and 4 is also quite high, with a value of roughly 42%. (a)

Conclusions
In this paper, we presented a new numerical approach to model the small-scale interaction of waves, pancake ice floes and interstitial grease ice. For this, the material behaviour of ice floes and grease ice is separately described considering for both characteristic material properties. Ice floes are modelled as elastic solids with high stiffness using generalised Hooke's law. Grease ice behaves as a viscous fluid obeying a viscous-plastic flow rheology. This two-phasic description is in contrast to the commonly used smeared phenomenological approaches in which a large heterogeneous sea ice area is homogenised as one isotropic, continuous material with averaged quantities.
To study the influence of wave forcing on the mechanical response of sea ice, a linearised harmonic propagating wave has been imposed on the domain via the Froude-Krylov force emulating the wave pressure field acting on the ice floe circumference and the skin drag acting on the ice floe basal plane. The form drag on the ice floe due to grease ice is implicitly accounted for by the continuum approach comprising both ice constituents.
To justify small-scale modelling a domain-size threshold was identified as smaller than 400 × 400m 2 where temporal and spatial fluctuations of the sea ice rheology variables (stress, strain rate and viscosity) become significant for all considered wave periods T = 6 sto20 s with 20m-diameter ice floes. This wave period-dependent threshold marks the transition from small-to large-scale modelling where a phenomenological model with homogenized material properties can be utilized.
Simulating three realistic sea ice layouts of 100 × 100m 2 extracted from in-situ images of the Antarctic MIZ, the robustness of the approach has been demonstrated. Furthermore, it was found that mechanical sea ice response tends to become independent of the detailed distribution of floes for wave periods larger than T = 16s in case of homogeneous sea ice conditions with similar ice concentrations and median floe caliper diameters. Moreover, the discrepancy between stress and strain rate curves increases for smaller wave period. An increasing ratio of ice floe diameter to wavelength results in a reduced stress and strain rate response in the considered realistic sea ice layouts. Based on the observation of significantly larger differences in floe stress between layouts compared to that of grease ice, it seems that the Froude-Krylov force is the main source for those observations considering that it solely acts on the ice floe circumference. As such, the dependency on the ratio of floe diameter to wave length becomes significant. It might be also plausible that the relatively large amount of larger floe diameters in layout 1 reduces the wave-floe interaction and thus, floe collision due to larger inertia. Both hypotheses, however, warrant further in-depth investigations.
Lastly, the influence of floe shape and diameter variations was shown to be significant considering additionally an idealized sea ice layout with disk-shaped floes with identical ice concentration and average area per ice floe. The domain-averaged strain rate magnitude and bulk viscosity are mostly affected by the detailed geometrical floe properties, with an average discrepancy of 103% and 42%, respectively.
In summary, results of this paper demonstrate the importance of detailed small-scale modelling on the subkilometer scale to resolve the mechanical response of a heterogeneous sea ice cover due to the complex interaction of waves, floes and grease ice. If large-scale regional models are to be informed by the actual material behaviour as originated on smaller scale, averaged material parameters need to be obtained from homogenization and up-scaling procedures. These rely on the identification of the minimum domain size threshold where the actual kinematics and material composition of the problem must be addressed in detail and the averaged quantities are statistically representative for larger domains.