A New Method for Modeling Effects of Surface Ice on Waves

: Accurate prediction of ocean surface wave attenuation in polar marginal ice zones remains a challenge. In this article, an alternative approach to the problem is introduced, in which the ice layer is represented with a modiﬁed version of the vegetation damping parameterization in a phase-resolved wave model. The new representation is evaluated by comparison to theory and measured data under varied wave and ice conditions. Model-estimated proﬁles of RMS water velocity and Reynolds stress under ice layers with different drag coefﬁcients are found to be qualitatively comparable to a range of nondimensional proﬁles computed using viscous layer theory. Modeled proﬁles appear somewhat vertically “stretched” relative to theoretical results, and in this respect, they more closely resemble measurements obtained during a recent wave–ice laboratory experiment. Estimated values of the wave attenuation coefﬁcient and wavenumber in ice from the adapted model align well with theory and with a range of lab and ﬁeld datasets. Several additional model ice parameters are available to facilitate a more nuanced representation of surface ice effects and will be investigated further in an upcoming companion study.


Introduction
As polar regions become increasingly open to navigation, there is an expanding need for accurate prediction of wave propagation in water covered by sea ice.Ships navigating the growing marginal ice zones (MIZs) will require additional guidance to avoid hazards, access resources, conduct research, and ensure compliance with international maritime law.Under stormy conditions, large waves can penetrate far into a MIZ, turning heavy ice plates into battering rams that threaten vessels, structures, and lives.There has been considerable progress toward improving our predictive capability and our understanding of the complex physical processes of wave-ice interaction [1].Nevertheless, at present, it is still not possible to reliably and consistently predict wave transmission and attenuation rates for the highly variable ice distributions that are often found in a MIZ.
Over the past decade or more, a wide range of field experiments have been conducted to measure the propagation of waves through fields of ice floes.Large-scale expeditions to the Arctic and Antarctic have compiled substantial datasets for a range of different sea and ice environments (e.g., [2][3][4]).As with all field efforts, participants and sponsors must of course accept and work with the specific conditions that they encounter during the measurement period; a field experiment cannot generally be configured to match a model simulation or produce a specific combination of waves and ice.The challenging and variable MIZ environment can also make it difficult or impossible to measure and/or parameterize some quantities that are important to model forecasts, such as effective water and ice viscosities, properties of boundary layers at the water-ice interface, and irregularities in the structure or distribution of the ice.
Smaller-scale wave-ice experiments in laboratories offer a more configurable and economical alternative, and many have addressed questions inaccessible to field research (e.g., [5][6][7]).Of course, though this approach allows for much more focused testing of specific hypotheses, it also has shortcomings due to the approximate nature of the simulated environment, which can include unrealistic wave conditions, inconsistent scaling of ice dimensions and/or material properties, and interference/reflection from test region boundaries.However, working in the laboratory does allow analyses to focus on specific physical questions using selected ice types and/or wave frequencies, to measure the wave flow field and wave-ice boundary layer much more precisely, and to directly test model forecasts of wave dispersion and attenuation.
There has also been substantial progress in recent years on the theory and modeling of wave-ice interaction in MIZs, for which a variety of different approaches have been used (e.g., see reviews by [8][9][10]).Phase-resolved wave models based on the Navier-Stokes equations (e.g., [11,12]), nonlinear potential flow (e.g., [13,14]), or the Boussinesq equations (e.g., [15][16][17]) capture details of individual waves and flow fields but are generally limited to shorter simulation periods and often shallower water.Phase-averaged [18] or spectral [9] approaches generally sacrifice smaller-timescale information in exchange for longer timespans in larger domains.Sea ice representations range from bonded discrete elements [19] to elastic plates [20] to viscoelastic layers [21], each of which offers a different level of approximation.Models of wave propagation in ice generally include one or more of four primary processes: scattering by floes, viscoelastic damping, energy loss due to ice flexure, and energy loss due to friction [10].Models that simulate the ice response to waves tend to focus either on ice fracturing or on the evolution of the floe distribution [9].Each approach offers a specific perspective on different aspects of the same questions.
Though it has been common practice in the literature to express wave attenuation in ice as a function of frequency alone, it has become clear that the specific characteristics of the ice itself also play an important role [22].As noted above, ice conditions in a MIZ are often highly variable, including a mix of grease ice, frazil, and pancakes along with irregularly shaped floes of widely varying dimensions.To remain computationally practical, a model generally cannot account for all of these variations but must instead incorporate effective parameterizations in which their combined effects are represented by a small number of tunable quantities (i.e., drag coefficient, effective viscosity, etc.).Extensive analysis of available lab and field data, supplemented by selective model simulations, is necessary to identify and calibrate these parameters.
Among the questions remaining about wave attenuation in MIZs are several relating to the fluid boundary layer that is generated beneath the ice by the propagating waves.Although the wave-attenuating effects of this under-ice wave boundary layer can be significant [23,24], its detailed characteristics have not been fully examined.Zhao and Shen [21] used a three-layer viscous model to estimate the attenuating effects of such a boundary layer and qualitatively demonstrated that it increased wave damping compared to the effects of the ice alone.Rabault et al. [25] and Orzech et al. [6] have had some initial success in mapping the characteristics of under-ice boundary layers in laboratory experiments, but their results are limited to a few specific test configurations, and more comprehensive measurements are required to fully quantify the properties of this layer and its contribution to wave attenuation.
This article is one of a series of related studies investigating the small-scale interactions between ocean waves and surface ice floes (also including [6,12,19,22,26]).The primary goal of the present study is to introduce and demonstrate a new method for modeling the velocity fields and attenuation of ocean waves propagating through surface ice at smaller scales.This somewhat heuristic approach has been implemented through modifications to NHWAVE [27,28], an existing non-hydrostatic wave model.The model's original representation of wave damping due to vegetation has been adapted to recast it as a configurable ice zone that affects wave fluxes only in a user-defined narrow region near the water surface.This reformulation relies on several broad physical similarities between the damping of waves in a marginal ice zone and wave damping by vegetation.In particular, both an ice layer and a vegetation region are permeable to some degree, allowing water to flow through them while also subjecting it to drag forces and diffusive effects.The size, spatial distribution, and effective roughness of ice and vegetation all vary to a greater or lesser degree, and each may be configured in this model representation.NHWAVE's original formulation for vegetation includes several additional adjustable parameters that can also be adapted for use with surface ice; these will be discussed in greater detail below.It is anticipated that this more realistic representation of the surface ice layer will help us to gain a better understanding of physical mechanisms related to wave attenuation and the wave-generated boundary layer beneath the ice.This analysis will examine wave-generated flow fields under surface ice and determine wave attenuation rates for several combinations of wave and ice conditions.Model results will be compared with theoretical predictions and available measurements from laboratory and field experiments.The adapted model will be applied to map out and examine the under-ice boundary layer at high resolution for several combinations of wave and ice conditions.These results are intended to extend and supplement the findings of related investigations by Orzech et al. [6] and Yu [26], and specific results from these two studies will be revisited below for comparison and evaluation of model output.
The wave model, its vegetation parameterization, and the adaptations for surface ice are detailed in Section 2. Model results are presented and compared with theory and measurements in Section 3. Additional discussion and conclusions are provided in Sections 4 and 5, respectively.

Original Model
The model utilized for this study, NHWAVE, is a nonhydrostatic, finite volume, phaseresolving wave model developed at the University of Delaware [27].Shortly after its initial release, the original model was updated to include the dissipative effects of "rigid" vegetation canopies [29] and, more recently, "flexible" vegetation [30].
NHWAVE solves the incompressible 3D Navier-Stokes equations in a domain with Cartesian horizontal dimensions (x, y) and a vertical dimension represented by a surface/bed following coordinate (σ).The governing equations are written in conservative form as where Q is the vector of conserved variables, Φ(Q) contains the flux variables, and ∇ is the total spatial derivative, given by In (1), Q and Φ are both functions of water depth h, free surface η, and velocities u, v, and ω (the contravariant of the vertical velocity in the σ dimension).The combined source term S includes separate components representing the effects of pressure gradients, turbulent mixing, bottom slope, and vegetation.The vegetation source term, S v , includes user-adjustable parameters for a number of quantifiable physical characteristics that affect the flow field.
The wave model's parameterizations for the damping effects of vegetation are configured to represent plants that grow up from the seabed to a height less than or equal to that of the water surface.Though users can specify either "rigid" or "flexible" vegetation configuration in the original model, the present analysis will utilize only a modified form of "rigid" vegetation to represent surface ice.
Turbulence production in waves due to effects of rigid vegetation was implemented in NHWAVE by Ma et al. [29] as a nonlinear k − ε model, which was calibrated and validated using data from open channel flow and random wave propagation in vegetation fields.The authors found that the drag produced by the k − ε model could significantly alter turbulent flow structure but had a weaker effect on wave attenuation.In the Cartesian reference frame, the model's continuity and momentum equations for the flow field in the presence of vegetation are written as in which x i,j and u i,j represent the coordinates and velocity components, respectively.Water density is represented by ρ, pressure by p, and gravity by g i (nonzero in vertical dimension only, and the shear matrix τ ij includes a dispersive flux term along with turbulent and viscous stresses.The viscous and form drag of the vegetation are combined into a single term, f di , which is modeled as where C D is a generalized bulk drag coefficient, and λ = b v N is vegetation density, with b v representing the stem diameter, and N equal to the number of stems per unit cross-sectional area.The inertial or "virtual mass" force, f vmi , is non-negligible for oscillatory wave flow and is determined from in which C M is the virtual mass coefficient.The primary vegetation parameters that are configurable in the model's input file include C D , b v , N, and C M .Additional parameters that may be modified include two coefficients related to turbulent drag in the k-ε formulation.

Modifications to Represent Surface Ice Layer
For this study, the rigid vegetation parameterization in NHWAVE is rewritten to limit the imposition of the drag and shear effects to a user-defined region near the water surface.Instead of at the bed, the wave energy is now damped only above a specified ice depth, within a horizontally fixed surface layer.The damping zone is comparable to a semi-permeable mesh, allowing the water to pass through it but subjecting the internal and nearby flow to shear and form drag that vary depending on specified parameter settings.Below this ice region, the water is subjected only to viscous effects produced by its own (user-specified) kinematic viscosity.The ice layer is configured to track the surface motion, moving up and down with wave crests and troughs and maintaining the user-specified ice thickness at all locations (for further details, see [12]).
To limit the scope of this initial analysis of model performance, only the drag coefficient C D from ( 6) is allowed to vary to represent different types of ice layers in the tests to follow.The number of stems per unit area is fixed at N = 8000 m −2 , the stem width is held constant at b v = 0.01 m, and the virtual mass coefficient C M is set to 1 in the ice (as the inertial force of the ice on the waves, f vmi , is assumed to be significant).The turbulent drag coefficients for the ice layer are set to the default values provided by [29]: C f k = 1.0 and C f = 1.33.

Model Test Configurations
To evaluate the performance of the new ice representation, the adapted NHWAVE model is applied to estimate wave and flow-field characteristics in a series of different test configurations, including both smaller-and larger-scale simulations.Subsurface profiles of RMS velocity and Reynolds stress are computed at selected ice-covered locations for different wave and ice types in smaller-scale "lab" tests.Following this, wave attenuation rates and wavenumber values are calculated for ice-covered regions at both smaller and larger scales.Model output is compared to theoretical predictions, laboratory measurements, and field data.
All model simulations are limited to two dimensions and have a similar general configuration (Figure 1).An internal wavemaker generates waves near the left end of a tank that is made deep enough to avoid bottom shoaling effects.The waves propagate several wavelengths in open water before they encounter a surface ice layer.After propagating through the ice for several more wavelengths, the waves are sampled at two vertical transects to determine velocity profiles, attenuation rates, and wavenumbers.The measurement region is set near the horizontal center of the ice layer in all cases (to avoid numerical effects from wave interaction with either ice edge).After passing the test range, the waves then propagate several wavelengths farther through the ice and enter open water again before finally reaching the far end of the tank.Waves are fully damped with numerical sponges as they approach either the left or the right boundary.
To evaluate the performance of the new ice representation, the adapted NHWAVE model is applied to estimate wave and flow-field characteristics in a series of different test configurations, including both smaller-and larger-scale simulations.Subsurface profiles of RMS velocity and Reynolds stress are computed at selected ice-covered locations for different wave and ice types in smaller-scale "lab" tests.Following this, wave attenuation rates and wavenumber values are calculated for ice-covered regions at both smaller and larger scales.Model output is compared to theoretical predictions, laboratory measurements, and field data.
All model simulations are limited to two dimensions and have a similar general configuration (Figure 1).An internal wavemaker generates waves near the left end of a tank that is made deep enough to avoid bottom shoaling effects.The waves propagate several wavelengths in open water before they encounter a surface ice layer.After propagating through the ice for several more wavelengths, the waves are sampled at two vertical transects to determine velocity profiles, attenuation rates, and wavenumbers.The measurement region is set near the horizontal center of the ice layer in all cases (to avoid numerical effects from wave interaction with either ice edge).After passing the test range, the waves then propagate several wavelengths farther through the ice and enter open water again before finally reaching the far end of the tank.Waves are fully damped with numerical sponges as they approach either the left or the right boundary.

Results
A total of eight different wave-ice simulations are included in this initial model evaluation (Table 1).This combination of configurations is used to test how the modeled ice layer affects both the flow field and the propagation characteristics of several different wave types at multiple scales and resolutions.For all tests, as noted earlier, only the drag coefficient, D C , is adjusted for the ice, and the other ice parameters are fixed at the values provided in Section 2.2.Six "lab-scale" model simulations (L1-L6) are completed to facilitate an examination of velocity fields near the ice boundary.The simulations feature two different monochromatic wave types (col. 2 and 3) and three different ice drag settings (col.6).For all of these tests, modeled ice thickness (col.4) is fixed at 1 cm, and tank depth (col.5) is set to 1.5 m.These six scenarios include a combination of high resolution (for resolving waves, boundary layers and ice thickness) and moderate to large horizontal extent (to avoid numerical effects near the wavemaker and ice edges).The horizontal and vertical grid-spacing of The distance between wavemaker and the ice edge, ∆X wi , is 6 m for lab-scale waves and 600 m for field-scale waves (i.e., roughly 2.5 to 4 wavelengths, depending on the case).The length of the test range is four times the open-water wavelength, λ = gT 2 /2π, in each test case.The total length of the ice is at least twelve wavelengths.

Results
A total of eight different wave-ice simulations are included in this initial model evaluation (Table 1).This combination of configurations is used to test how the modeled ice layer affects both the flow field and the propagation characteristics of several different wave types at multiple scales and resolutions.For all tests, as noted earlier, only the drag coefficient, C D , is adjusted for the ice, and the other ice parameters are fixed at the values provided in Section 2.2.Six "lab-scale" model simulations (L1-L6) are completed to facilitate an examination of velocity fields near the ice boundary.The simulations feature two different monochromatic wave types (col. 2 and 3) and three different ice drag settings (col.6).For all of these tests, modeled ice thickness (col.4) is fixed at 1 cm, and tank depth (col.5) is set to 1.5 m.These six scenarios include a combination of high resolution (for resolving waves, boundary layers and ice thickness) and moderate to large horizontal extent (to avoid numerical effects near the wavemaker and ice edges).The horizontal and vertical grid-spacing of these smaller-scale simulations are dx = 2 cm and dz = 2 mm, respectively, the tank length, X tank , is 26.2 m, and output is sampled at ∆t = 0.1 s.
Two additional "field-scale" simulations (F1-F2) feature different wave types in a much larger domain (X tank = 2.6 − 5.2 km) with significantly lower resolution (dx = 0.5 − 2 m, dz = 0.5 − 1 m, and ∆t = 0.5 s).Modeled ice thickness is fixed at 50 cm, and water depth is set to 50 m (F1) and 100 m (F2) as required to maintain a deep-water wave environment.These lower-resolution tests are less computationally expensive than the lab-scale simulations, but they do not provide detailed information about the flow field near the ice boundary.
In Section 3.1, examples of modeled flow fields are presented for two of the lab-scale cases to illustrate the variations in the ice layer's effect on fluid velocities.In Sections 3.2 and 3.3, model-estimated velocity and Reynolds stress profiles for lab-scale cases L1-L6 are directly compared to theoretical predictions and laboratory measurements.The lower resolution cases F1-F2 are not sufficiently resolved to capture the details of the flow field near the ice, so they are excluded from these results.In Sections 3.4 and 3.5, both the lab-scale and the field-scale model estimates of wave attenuation and wavenumber in ice are compared with theory and measurements from the lab and field.As noted earlier, this modeling study is part of a series of investigations of wave-ice interactions and the water-ice boundary layer.Some results from the theoretical analysis of [26] and laboratory experiments by [6] are reproduced below for comparison with model output.

Flow Field Effects
Sample velocity fields extracted from a near-surface, ice-covered location in model cases L1 and L6 are presented in Figure 2. Vorticity is elevated in the vicinity of the ice layer in both panels, as the velocity damping generally has reduced flow in the direction of wave propagation.However, vorticity is considerably lower for case L1 (top panel), in which shorter waves (T = 1 s, H = 2 cm) propagate through a 1 cm ice layer with a relatively low drag coefficient, (C D = 0.1).For case L6 (lower panel), the ice layer has a significantly larger drag coefficient, (C D = 1.0), and imposes a visible negative correction on horizontal velocities near the surface, affecting both positive and negative components for the larger waves (T = 1.34 s, H = 6 cm).This correction is too small to be seen in the flow field for case L1.There is little or no damping of the vertical velocity components in either case.
It is worth noting that, because the surface ice boundary in these simulations is moving vertically with the waves, it produces a different effect on the fluid than would be experienced along a relatively fixed bottom boundary at the seabed (for example, under shoaling waves near the beach).This surface boundary moves higher under wave crests and lower under wave troughs.Over several wave periods, the net shear and drag effects of the ice on the water are also spread out vertically over the water column, to an extent determined by the amplitude of the propagating waves.This spreading effect plays a role in how we define the "thickness" of the ice layer and is addressed further in the following section.

RMS Velocity Profiles
Vertical profiles of RMS velocity from the adapted wave-ice model are computed for cases L1-L6 along a transect at a location corresponding to X a in Figure 1.Profiles are obtained after each simulation has reached a steady state (i.e., wave height at location X a remains approximately constant).In each case, along-tank RMS velocities are computed as where u 1 (z, t) is horizontal velocity at depth z and time t, ∆t is the sample rate, and T 10 is equivalent to the time for 10 wave periods 1 .Near-surface sections of the profiles for the six lab-scale wave configurations are presented in Figure 3 (solid lines), along with corresponding profiles that would be obtained in open water for each wave type (dashed lines).A cyan-shaded region is included near the top of each panel to indicate the approximate depth range that was primarily or wholly occupied by ice during each wave period (arbitrarily set to wave amplitude plus half of ice thickness in each case).It is worth noting that, because the surface ice boundary in these simulations is moving vertically with the waves, it produces a different effect on the fluid than would be experienced along a relatively fixed bottom boundary at the seabed (for example, under shoaling waves near the beach).This surface boundary moves higher under wave crests and lower under wave troughs.Over several wave periods, the net shear and drag effects of the ice on the water are also spread out vertically over the water column, to an extent determined by the amplitude of the propagating waves.This spreading effect plays a role in how we define the "thickness" of the ice layer and is addressed further in the following section.

RMS Velocity Profiles
Vertical profiles of RMS velocity from the adapted wave-ice model are computed for cases L1-L6 along a transect at a location corresponding to a X in Figure 1.Profiles are obtained after each simulation has reached a steady state (i.e., wave height at location a X remains approximately constant).In each case, along-tank RMS velocities are computed as where 1 ( , ) u z t is horizontal velocity at depth z and time t, t  is the sample rate, and 10 T is equivalent to the time for 10 wave periods 1 .Near-surface sections of the profiles for the As expected, the magnitudes of the RMS velocities are significantly greater for the larger waves in cases L4-L6.The divergence of the modeled profiles (solid lines) from the open-water profile (dashed line) in the upper several centimeters of the water column provides evidence of a wave-generated boundary layer immediately below the ice interface.Due to the eddy generation and drag effects of the semi-permeable ice cover, this type of boundary layer is expected to be significantly larger in magnitude and vertical extent than that seen in open water when ice is not present [26], and such an under-ice boundary layer can often be turbulent [23].The increased velocity variance is expected to be greatest near the water-ice interface and taper with increasing depth.In these model results, the boundary layer region extends to roughly three times greater depth for the L4-L6 waves (H = 6 cm) than it does for the L1-L3 waves (H = 2 cm), implying that the vertical depth to which the boundary layer affects RMS velocities is dependent on the wave amplitude.The magnitudes of the boundary layer velocities also increase with increasing values of the ice drag coefficient, C D , although this effect appears somewhat more muted for the larger coefficient values (i.e., between L2-L3 or L5-L6).
depth range that was primarily or wholly occupied by ice during each wave period (arbitrarily set to wave amplitude plus half of ice thickness in each case).As expected, the magnitudes of the RMS velocities are significantly greater for the larger waves in cases L4-L6.The divergence of the modeled profiles (solid lines) from the open-water profile (dashed line) in the upper several centimeters of the water column provides evidence of a wave-generated boundary layer immediately below the ice interface.Due to the eddy generation and drag effects of the semi-permeable ice cover, this type of boundary layer is expected to be significantly larger in magnitude and vertical extent than that seen in open water when ice is not present [26], and such an under-ice boundary layer can often be turbulent [23].The increased velocity variance is expected to be greatest near the water-ice interface and taper with increasing depth.In these model results, the boundary layer region extends to roughly three times greater depth for the L4-L6 waves (H = 6 cm) than it does for the L1-L3 waves (H = 2 cm), implying that the vertical depth to which the boundary layer affects RMS velocities is dependent on the wave amplitude.The magnitudes of the boundary layer velocities also increase with increasing values of the ice drag coefficient, D C , although this effect appears somewhat more muted for the larger coeffi- cient values (i.e., between L2-L3 or L5-L6).
As detailed field data are not presently available, results from Yu [26] and Orzech et al. [6] are instead used here to establish "ground truth" values for the magnitude and vertical extent of RMS velocity profiles and boundary layers that would be expected under these conditions.A comprehensive, detailed picture of the velocity field properties expected in this near-surface region is provided by the authors of [26], who completed a theoretical wave-ice analysis in which the surface ice was represented as a continuous viscous layer.Direct measurements of these under-ice velocity fields were also obtained in lab tests conducted by the authors of [6], who employed a particle imaging velocimetry (PIV) system to record subsurface water velocities as waves propagated through ice floes in a refrigerated laboratory wave flume.As detailed field data are not presently available, results from Yu [26] and Orzech et al. [6] are instead used here to establish "ground truth" values for the magnitude and vertical extent of RMS velocity profiles and boundary layers that would be expected under these conditions.A comprehensive, detailed picture of the velocity field properties expected in this near-surface region is provided by the authors of [26], who completed a theoretical wave-ice analysis in which the surface ice was represented as a continuous viscous layer.Direct measurements of these under-ice velocity fields were also obtained in lab tests conducted by the authors of [6], who employed a particle imaging velocimetry (PIV) system to record subsurface water velocities as waves propagated through ice floes in a refrigerated laboratory wave flume.
Sample velocity profiles from both the theoretical analysis and the laboratory experiment are replotted in Figure 4 below for comparison with the model results for cases L1-L6.As the theoretical profiles were originally calculated for a non-dimensional domain, both theoretical and laboratory results are presented in a non-dimensional form to allow them to be directly compared.Following the conventions of [26], depths in each panel are normalized by the nominal ice thickness, h, and velocities are normalized by the linear open-water velocity at the water-ice interface.Depths are also adjusted so that z/h = 0 corresponds to the normalized depth of the interface.Both the theoretical and lab results of Figure 4 appear qualitatively similar to the dimensional profiles presented for the model in Figure 3.In all of the profiles (solid lines), the boundary layer is visible as a divergence of RMS velocity from the (dashed) open-water profile in the region near the water-ice interface.In the left panel of Figure 4, the three different theoretical curves demonstrate how the vertical extent of each boundary layer is affected by both the relative viscosity of the water and the relative thickness of the ice layer.
in Figure 4, the plotted velocities for cases L1-L6 are now normalized by the theoretical open-water velocity at the interface.Depths are normalized by a representative ice region thickness, ', h which is set equal to the wave amplitude plus half the modeled ice thick- ness, as in Figure 3. Normalized depths are also adjusted so that zero depth now corresponds to the bottom of the representative ice region, as with the theoretical and experimental profiles.To more directly compare modeled RMS velocity profiles with the theory and lab results, the model profiles of Figure 3 are replotted in a normalized form in Figure 5.As in Figure 4, the plotted velocities for cases L1-L6 are now normalized by the theoretical open-water velocity at the interface.Depths are normalized by a representative ice region thickness, h , which is set equal to the wave amplitude plus half the modeled ice thickness, as in Figure 3. Normalized depths are also adjusted so that zero depth now corresponds to the bottom of the representative ice region, as with the theoretical and experimental profiles.
Examining the first two theoretical cases in Figure 4 (with viscosity ratios ν i /ν w = 1 × 10 4 and 1 × 10 2 , respectively), it is apparent that when effective water viscosity ν w increases relative to ice viscosity ν i , the vertical extent of the boundary layer also increases.Comparing the second and third theoretical cases (both with ν i /ν w = 1 × 10 2 , but with the latter case featuring a 3× thicker ice layer), it can be seen that an increase in the ice layer thickness h (relative to wavelength) can produce a further increase in the boundary layer's relative vertical extent.
The RMS velocity profile obtained from lab measurements by [6] (Figure 4, right panel) contrasts significantly with the theoretical profiles.This profile was measured during a simulation in which small, monochromatic waves (T = 1 s, H = 4 cm) propagated through a 3.8 cm ice layer.The measured boundary layer exhibits a considerably greater relative vertical extent than those predicted by viscous layer theory.Based only on the results described for the theoretical profiles, it is logical to conclude that the velocity variance of the water immediately beneath the ice was significantly enhanced by turbulence in the experiment, and that the experimental ice thickness was unusually large compared to the waves (both of which were likely true).However, it seems doubtful that these characteristics of the lab simulation would be sufficient in themselves to extend the measured velocity profile so much deeper than the theoretical profiles.Examining the first two theoretical cases in Figure 4 (with viscosity ratios iw  = 1 × 10 4 and 1 × 10 2 , respectively), it is apparent that when effective water viscosity w  in- creases relative to ice viscosity i  , the vertical extent of the boundary layer also increases.
Comparing the second and third theoretical cases (both with iw  = 1 × 10 2 , but with the latter case featuring a 3× thicker ice layer), it can be seen that an increase in the ice layer thickness h (relative to wavelength) can produce a further increase in the boundary layer's relative vertical extent.
The RMS velocity profile obtained from lab measurements by [6] (Figure 4, right panel) contrasts significantly with the theoretical profiles.This profile was measured during a simulation in which small, monochromatic waves (T = 1 s, H = 4 cm) propagated through a 3.8 cm ice layer.The measured boundary layer exhibits a considerably greater relative vertical extent than those predicted by viscous layer theory.Based only on the results described for the theoretical profiles, it is logical to conclude that the velocity variance of the water immediately beneath the ice was significantly enhanced by turbulence in the experiment, and that the experimental ice thickness was unusually large compared to the waves (both of which were likely true).However, it seems doubtful that these characteristics of the lab simulation would be sufficient in themselves to extend the measured velocity profile so much deeper than the theoretical profiles.
In normalized form, the results for the six model cases (Figure 5) are surprisingly consistent in both magnitude and the vertical extent of the boundary layer region.The normalized velocity magnitudes in the boundary layer are still affected by the size of the ice drag coefficient, D C , but they do not appear to be strongly affected by the specific wave period.Employing the wave amplitude in the normalization of the water depth appears to have largely eliminated the variation in the vertical extent of the boundary layer that was seen for the two different wave configurations in Figure 3.The model profiles are qualitatively comparable to both the theoretical and the measured RMS profiles of Figure 4.However, though the vertical extent of the theoretical boundary region is only a fraction of the ice thickness, for both the measured and modeled boundary layers, this normalized In normalized form, the results for the six model cases (Figure 5) are surprisingly consistent in both magnitude and the vertical extent of the boundary layer region.The normalized velocity magnitudes in the boundary layer are still affected by the size of the ice drag coefficient, C D , but they do not appear to be strongly affected by the specific wave period.Employing the wave amplitude in the normalization of the water depth appears to have largely eliminated the variation in the vertical extent of the boundary layer that was seen for the two different wave configurations in Figure 3.The model profiles are qualitatively comparable to both the theoretical and the measured RMS profiles of Figure 4.However, though the vertical extent of the theoretical boundary region is only a fraction of the ice thickness, for both the measured and modeled boundary layers, this normalized extent is roughly of order one, i.e., the boundary layer depth is roughly the same as the nominal ice thickness value.
To explain this difference, we must consider that, unlike the theoretical analysis, the water-ice interfaces in the lab and model environments are both vertically oscillating, so that the boundary layers themselves are moving up and down several centimeters as each wave passes.The region occupied by ice in the lab and the model is determined by the amplitude of the waves, which carry the ice higher in their crests and deeper in their troughs.Within this expanded ice region, locations near mean sea level have ice in them at all times, whereas other neighboring depths are only occupied by ice for part of a wave cycle.As described in Section 3.1, this oscillation effectively spreads out the measured boundary layer region over a larger range, roughly proportional to the wave amplitude.When averaged over multiple wave periods, this vertical displacement of the existing boundary layer affects RMS velocities deeper in the water column in a manner that is not considered in the theoretical analysis.
In addition, both the modeled and lab ice include some degree of permeability, whereas the theoretical approach of Yu [26] enforces an impermeable, fixed boundary at the mean water-ice interface.The flux of water into and out of the ice layer in the lab and the model augments the ice-generated turbulence, an additional effect that is not considered by a viscous layer representation.As wave height and/or ice thickness increase, this turbulent contribution grows in magnitude and is advected deeper into the water column below the ice.

Reynolds Stresses
The wave-induced Reynolds stress, τ uv , is computed from the time-averaged product of the wave-induced velocities, uv , and represents the mean momentum flux that is generated by the wave fluctuations.For rotational wave motion, uv = 0, but when this motion becomes irrotational (as in the vicinity of surface ice), the Reynolds stress becomes non-zero.Several examples of Reynolds stress profiles computed from viscous layer theory by Yu [26] and a single profile measured in the lab by Orzech et al. [6] are reproduced in Figure 6    Within the surface ice layer (z/h > 0), the values of both theoretical and experimental Reynolds stresses are positive, which is consistent with the expectation that the phase of the v oscillation is leading that of the u oscillation by less than 90 • in the viscous fluid [26].As we move downward and away from the wave-ice interface inside the wave boundary layer (z/h ≤ 0), the stress first decreases rapidly in all profiles, then reaches a negative peak (τ uv,min < 0), before finally returning to zero at greater depths.As seen earlier with the RMS velocity profiles, the experimental profile again appears vertically elongated relative to the theoretical profiles, likely due to both different conditions and differing treatments of the ice layer.The experimental profile also has a much larger negative peak and almost no positive Reynolds stress values.
Profiles of normalized Reynolds stresses are presented in Figure 7 for datasets from cases L1-L6 of the adapted wave-ice model, using the same normalization methods as above.These model profiles, each computed over 10 wave periods, again exhibit behavior and magnitudes that are qualitatively similar to the theoretical results.In all but one of the six model profiles, stress values consistently decrease in the region immediately below the water-ice interface (z/h ≤ 0), usually reaching a negative minimum value before returning to zero at greater depths.The vertical extent of the boundary layer and negative peak region is somewhat greater in the model results than in Reynolds stress profiles obtained from viscous layer theory, and in this respect, the modeled profiles are once again more like the experimental profile.In magnitude, however, the range of the modeled stress profiles is this time much closer to that of the theoretical profiles, with consistently positive values in the shallowest region close to the ice interface and only a slightly negative minimum stress value 2 .

Wave Attenuation
Wave attenuation is estimated through an analysis of surface elevation time series at two selected ice-covered locations in the model domain.The model is allowed to reach a steady state at both locations (i.e., height remaining approximately constant) before the time series data are extracted.The wave amplitude attenuation rate, k i , is calculated over a central portion of the ice region via the formulation where ∆X test is the distance between the two measured locations (four times the open-water wavelength), and H RMS,j is the RMS wave height at location X j , with j = a or b (see Figure 1).In cases L1-L3 from Table 1, ∆X test = 6.24 m, and in cases L4-L6, ∆X test = 11.2 m.In the cases F1 and F2, the values of ∆X test are 224.5 m and 400 m, respectively.The RMS wave heights are determined from the standard deviation, σ, of the surface time series at the two locations as where η is free surface elevations, η is their mean, and M is number of timesteps averaged.Attenuation coefficients obtained for the eight ice and wave configurations are summarized in Table 2.These results are consistent with expectations for changing values of either the ice drag coefficient or the wave properties.As C D is increased, the wave attenuation rate k i also increases for a given wave type.When ice drag is held constant, the shorter period waves are consistently more attenuated than the longer waves.

Case
T (s) As with the velocity and stress profiles of the preceding sections, it is again helpful to express wave attenuation rates in terms of non-dimensional dynamic parameters, so that data from different scales may be directly compared.When the effective ice thickness, h , is used as a normalizing quantity for these data (again following the conventions of [26]), the non-dimensional wave frequency, ω, can be computed as where ω = 2π/T is the open-water wave frequency, and g = 9.83 m/s is taken as the gravitational acceleration in polar regions.Correspondingly, the non-dimensional wave attenuation rate, k i , is obtained from This representation allows for small-scale lab results to be examined together with dynamically similar field data on the same curve.Ultimately, this capability can facilitate the extraction of generalized attenuation information from lab and/or model simulations, which can then be used for the configuration of field-scale forecasts in preparation for future field expeditions.
Normalized attenuation results for the eight model configurations of Table 2 are presented graphically (in log-log format) in Figure 8, along with additional results from both lab and field experiments.The model results line up reasonably well with the linear best fit to the other lab and field data that was determined by [26], with modeled attenuation rates slightly lower than the trend line in all eight cases.which can then be used for the configuration of field-scale forecasts in preparation for future field expeditions.
Normalized attenuation results for the eight model configurations of Table 2 are presented graphically (in log-log format) in Figure 8, along with additional results from both lab and field experiments.The model results line up reasonably well with the linear best fit to the other lab and field data that was determined by [26], with modeled attenuation rates slightly lower than the trend line in all eight cases.

Wavenumber in Ice
The wavenumber in ice, r k , is also expected to vary with different wave and ice con- ditions.Experimental results from [6] indicate that the value of r k for a given open-water wavenumber generally increases with increasing ice thickness.These findings are supported by theoretical predictions of wave dispersion in ice based on a mass loading ice representation [34].However, an opposite trend (i.e., decreasing r k with increasing ice thickness) is predicted when an elastic plate ice representation is employed [35,36].
Using the surface elevation data at the two ends of the test range shown in Figure 1, the wavenumber in ice is estimated by following the procedure outlined in Parra et al. [35].For each trial, the time series from the two locations are compared, and the time  [6,21,32] and field data of [2,33], along with model output from L1-L6 (asterisks) and F1-F2 (squares).Best fit line was computed by [26] based on selected lab and field data.(Figure adapted from [6]).

Wavenumber in Ice
The wavenumber in ice, k r , is also expected to vary with different wave and ice conditions.Experimental results from [6] indicate that the value of k r for a given openwater wavenumber generally increases with increasing ice thickness.These findings are supported by theoretical predictions of wave dispersion in ice based on a mass loading ice representation [34].However, an opposite trend (i.e., decreasing k r with increasing ice thickness) is predicted when an elastic plate ice representation is employed [35,36].
Using the surface elevation data at the two ends of the test range shown in Figure 1, the wavenumber in ice is estimated by following the procedure outlined in Parra et al. [35].For each trial, the time series from the two locations are compared, and the time required (∆t) for a selected wave crest to travel the length of the test range (∆X test ) is determined.The wave phase speed is then computed as c = ∆X test /∆t.The wavenumber is then determined by making use of the relationship c = 2π f k r (12) in which f = 1/T is the wave frequency (1 Hz or 0.75 Hz, depending on the specific trial).
Values of k r are summarized for the eight model test cases in Table 3, along with their corresponding open-water values, k 0 = 2π/λ, where λ is the open-water wavelength.For each case, the value of ∆X test is fixed at 4λ, as indicated in Figure 1.The values measured for ∆t increase consistently with increasing ice drag C d .In cases L1-L3 as well as L4-L6, this results in a phase speed that decreases slightly as C d increases.This decreasing relative phase speed is reflected in the increasing values of k r in Table 3.For cases F1 and F2, both ∆X test and ∆t are significantly greater for the 10 s waves, and the drag C d remains constant, producing the expected smaller value of k r for the larger F2 waves.The wavenumber magnitudes for the cases where drag is weak (L1, L4, F1, and F2) are all close to their open-water values, whereas those for ice with greater drag (L2, L3, L5, L6) are somewhat larger.When normalized by ice thickness in the same manner as the attenuation coefficients, the wavenumber values for all eight cases are reasonably close to those predicted by theory for waves of their frequency passing through a viscous surface layer of comparable thickness, and they align well with several sets of lab measurements (Figure 9).Note, however, that for this lower range of normalized frequency, the wavenumber curve based on viscous-layer theory is essentially indistinguishable from an open-water wavenumber curve (see, for example, Figure 2a in Yu [26]).

Discussion
In Sections 3.2 and 3.3, the profiles of RMS velocity and Reynolds stress obtained from the adapted NHWAVE model are qualitatively comparable to theoretical results based on a viscous layer ice representation and to laboratory measurements.However, there are several areas in which the modeled profiles differ notably from the lab and theory.For normalized RMS velocities, the boundary layer velocity magnitudes of model cases L1-L6 (Figure 5) are consistently larger than the values obtained from the lab and theory (Figure 4).This difference is likely just an issue of model ice configuration and can easily be resolved through further adjustments of the ice drag coefficient, D C , and/or other model parameters.It does, however, highlight the need for a more comprehensive calibration of this model.As field data are not generally available, such a calibration would have to rely on both theory and measured data to provide "ground truth" values, as we have done above.Although it is not certain that these two data sources will be sufficient, we plan to at least attempt such a calibration in the near future and describe the

Discussion
In Sections 3.2 and 3.3, the profiles of RMS velocity and Reynolds stress obtained from the adapted NHWAVE model are qualitatively comparable to theoretical results based on a viscous layer ice representation and to laboratory measurements.However, there are several areas in which the modeled profiles differ notably from the lab and theory.For normalized RMS velocities, the boundary layer velocity magnitudes of model cases L1-L6 (Figure 5) are consistently larger than the values obtained from the lab and theory (Figure 4).This difference is likely just an issue of model ice configuration and can easily be resolved through further adjustments of the ice drag coefficient, C D , and/or other model parameters.It does, however, highlight the need for a more comprehensive calibration of this model.As field data are not generally available, such a calibration would have to rely on both theory and measured data to provide "ground truth" values, as we have done above.Although it is not certain that these two data sources will be sufficient, we plan to at least attempt such a calibration in the near future and describe the results in a companion article.
The model's k-ε turbulence formulation, in combination with the porous ice representation, produces micro-scale turbulent energy within the modeled ice layer (generated by fluid flow past individual "stems") that is not considered in most theoretical analyses utilizing viscous or viscoelastic ice.Chen et al. [37] conducted a similar study of wave attenuation in a porous viscoelastic ice layer in 2019, finding that the wave attenuation rate became non-monotonic when frictional turbulence was activated in their model.This type of turbulence generation is roughly comparable to that which might be seen under broken ice and frazil in the lab or field, although it is much more evenly distributed in the model ice layer than it would be in a real marginal ice zone.By including both porosity and friction, this hybrid ice representation may eventually contribute to explaining unresolved aspects of wave behavior, such as the roll-over of amplitudes that has been measured in the field and also noted in the analyses of fragmented ice layers by [13] and of porous ice layers by [37].The extent to which the modeled ice turbulence affects the fluid flow immediately beneath the ice has not yet been fully investigated and quantified, but it seems likely that it is contributing to the growth of the modeled boundary layer.
Though the larger 6 cm/1.34 s waves of model cases L4-L6 yield Reynolds stress profiles of the same order and shape as those from the theory and lab (compare Figures 6 and 7), the smaller 2 cm/1 s waves of cases L1-L3 produce rather stunted profiles (Figure 7, left panel), which incorrectly portray some details of the under-ice boundary layer.The rather choppy results obtained for the smaller-scale cases emphasize the need for a sufficiently high vertical resolution in the adapted wave-ice model in order to accurately represent velocity fields for very small waves.Velocity and stress profiles were not generated in this analysis for the significantly larger field-scale cases, F1-F2, as their vertical resolution (0.5-1.0 m) is much too low to capture the details of the boundary layer.
The extreme resolution requirements of this type of investigation could potentially make the model too computationally expensive and time-consuming for some users.Owing to the specific combinations of length and time scales associated with the wave-ice boundary layer, simulations with the adapted model will consistently require a fairly large amount of computation time.The eight tests described above each needed 1-3 days of parallel computation by 56 processors in order to complete roughly 60 s of simulation.Because the model is currently limited to working on a uniformly spaced grid, the high vertical resolution of lab-scale cases L1-L3 (dz = 2.5 mm), in a tank of sufficient length (52 m) and depth (1.5 m) to effectively measure attenuation in deep-water wave conditions, require a 2D grid consisting of nearly 1.58 million cells (2632 × 600).This requirement would increase considerably in a 3D simulation.It could be reduced somewhat if NHWAVE were properly configured to use exponential vertical spacing in its grids, but this option does not presently seem to be functional in the model.Although an exponential vertical grid is advertised as being a model option [28], the capability does not appear to have been fully validated.Exponential spacing has not generated physically realistic results when utilized for these boundary layer simulations.
Although somewhat costly, this model approach has a few advantages that may make its use appealing, particularly in highly resolved studies of the wave-ice boundary layer for which the ice behavior and properties must be as realistic as possible.As described in Section 2, the vegetation representation developed by [29] for NHWAVE is highly configurable and includes several relevant parameters, each of which may also be adjusted in this adapted version of the model to modify the effects of surface ice.In the analysis presented here, we have only allowed the drag coefficient to be varied, keeping the remaining quantities constant, including the vegetation density, stem width, virtual mass coefficient, and several others.In contrast to the unmeasurable pseudo-parameters of some other wave-ice models (e.g., viscoelasticity), the physically based ice parameters of this model offer an opportunity to configure characteristics of the surface ice layer that more directly match measured lab and field values.
By varying the vegetation density (λ) and stem width (b w ), for example, the model surface ice can be tuned to represent the effects of a broad range of ice fields, from tightly packed ice floes to scattered pancakes.A low setting for λ, paired with a value of around 0.5 m for b w in a 10 cm thick, low-drag ice layer will simulate a field of scattered pancakes.Larger values of both parameters and a thicker ice layer will create a region of tightly packed ice floes.
Adjustment of the virtual mass coefficient (C M ) alters the inertial force and effective form drag of the ice in response to the oscillating wave motion.A larger value of C M should simulate a denser ice cover, extracting a greater amount of energy from the waves, and a smaller value would allow waves to propagate more freely through the ice with less attenuation.
As noted earlier, although it is configured to follow the motion of the water surface in the vertical direction, the ice in the adapted model is horizontally fixed and does not presently move in response to horizontal wave fluxes.Yu [26] predicted that a mobile surface ice layer should experience a wave-induced Eulerian steady streaming that augments the Stokes drift effect to produce a net Lagrangian transport of the ice in the primary wave direction.Implementation of horizontally mobile ice in this adapted model would facilitate a more precise investigation and quantification of these two transport contributors.Available options for implementing such a capability in this adapted model are currently being examined.

Conclusions
A new formulation for simulating the effects of surface ice on propagating ocean waves has been implemented within the non-hydrostatic model NHWAVE.The formulation is based on a modified version of the model's existing representation of the effects of vegetation at the seabed.To limit the scope of the present analysis, only a single parameter, the ice drag coefficient, has been varied in these initial tests of the new ice formulation.
The results presented here demonstrate that the adapted model does a reasonably good job of qualitatively reproducing the profiles of RMS velocity and Reynolds stress that are predicted by viscous layer theory and measured in laboratory experiments for these types of waves in ice.They also show that the model can reproduce the attenuating effects of surface ice on wave amplitudes for a range of wave frequencies and ice configurations.
A number of important questions must still be addressed to enable more accurate and reliable prediction of wave attenuation in the MIZ.In particular, it is not clear to what degree other physical parameters (e.g., thickness, viscosity and elasticity of the ice layer, sub-ice wave orbital velocities, relative ice and water density, etc.) affect the attenuating waves, and if so, how their effects can be quantified.Although computationally expensive, this adapted wave model may ultimately contribute to resolving some of these continuing issues.This analysis will continue by testing the effects of varying other model parameters and calibrating the system with a range of lab and field datasets, with results to be presented in a follow-up companion paper.

Figure 1 .
Figure 1.General NHWAVE configuration for simulations of wave interaction with frazil ice region.The distance between wavemaker and the ice edge, wi X  , is 6 m for lab-scale waves and 600 m for field-scale waves (i.e., roughly 2.5 to 4 wavelengths, depending on the case).The length of the test range is four times the open-water wavelength,2 2 gT  = , in each test case.The total length of

Figure 1 .
Figure 1.General NHWAVE configuration for simulations of wave interaction with frazil ice region.The distance between wavemaker and the ice edge, ∆X wi , is 6 m for lab-scale waves and 600 m for field-scale waves (i.e., roughly 2.5 to 4 wavelengths, depending on the case).The length of the test range is four times the open-water wavelength, λ = gT 2 /2π, in each test case.The total length of the ice is at least twelve wavelengths.

Figure 2 .
Figure 2. Sample flow fields with vorticity for a surface ice section of cases L1 (top) and L6 (bottom) at simulation time t = 35 s.Longest velocity vector shown corresponds to 0.05 m/s in top panel and 0.17 m/s in bottom panel.

Figure 2 .
Figure 2. Sample flow fields with vorticity for a surface ice section of cases L1 (top) and L6 (bottom) at simulation time t = 35 s.Longest velocity vector shown corresponds to 0.05 m/s in top panel and 0.17 m/s in bottom panel.

Figure 3 .
Figure 3. Vertical profiles of RMS along-tank velocity computed using (8) for model test cases L1-L3 (left panel) and L4-L6 (right panel).Profiles are computed with time series data from location a X in Figure 1.Dashed line in each panel shows the expected velocity profile that would be ex- pected for those wave conditions in open water (based on linear theory).Cyan shading in each panel indicates region where ice was at least intermittently present.

Figure 3 .
Figure 3. Vertical profiles of RMS along-tank velocity computed using (8) for model test cases L1-L3 (left panel) and L4-L6 (right panel).Profiles are computed with time series data from location X a in Figure 1.Dashed line in each panel shows the expected velocity profile that would be expected for those wave conditions in open water (based on linear theory).Cyan shading in each panel indicates region where ice was at least intermittently present.

Figure 4 .
Figure 4. Vertical profiles of normalized RMS velocity.Left panel: Based on viscous layer theory [26], for different ratios of ice to water viscosity,

Figure 4 .
Figure 4. Vertical profiles of normalized RMS velocity.Left panel: Based on viscous layer theory [26], for different ratios of ice to water viscosity, ν i /ν w .Depths z are normalized by ice layer thickness, h.Third case (*) has 3× greater relative ice thickness.Right panel: From lab measurements by [6], for 4 cm, 1 s waves in surface ice of thickness 3.8 cm.Cyan shading in each panel indicates ice layer region (Figure adapted from [6]).

Figure 5 .
Figure 5. Normalized RMS velocity profiles for cases L1-L3 (left panel; T = 1.0 s, H = 2 cm) and L4-L6 (right panel; T = 1.34 s, H = 6 cm), obtained at same along-tank location as in Figure 3. Velocities are normalized by open-water velocity at interface, and depths are normalized by effective ice thickness h' (wave amplitude plus half ice thickness).Cyan shading indicates area partially or entirely occupied by ice layer.

Figure 5 .
Figure 5. Normalized RMS velocity profiles for cases L1-L3 (left panel; T = 1.0 s, H = 2 cm) and L4-L6 (right panel; T = 1.34 s, H = 6 cm), obtained at same along-tank location as in Figure 3. Velocities are normalized by open-water velocity at interface, and depths are normalized by effective ice thickness h' (wave amplitude plus half ice thickness).Cyan shading indicates area partially or entirely occupied by ice layer.
below to again establish a baseline for what might be expected from the adapted wave-ice model.As in the preceding section, these profiles are displayed in a non-dimensional form, with the water depth z normalized by the ice thickness h, and the Reynolds stress computed with u and v velocities that have been normalized by the openwater velocity amplitude at the level of the water-ice interface.The lab-measured stress profile is computed for the same test case as was used in Section 3.2, with monochromatic waves of period T = 1 s and height H = 4 cm and a surface ice thickness of 3.8 cm.J. Mar.Sci.Eng.2023, 11, x FOR PEER REVIEW 12 of 20returning to zero at greater depths.The vertical extent of the boundary layer and negative peak region is somewhat greater in the model results than in Reynolds stress profiles obtained from viscous layer theory, and in this respect, the modeled profiles are once again more like the experimental profile.In magnitude, however, the range of the modeled stress profiles this time much closer to that of the theoretical profiles, with consistently positive values in the shallowest region close to the ice interface and only a slightly negative minimum stress value 2 .

Figure 6 .
Figure 6.Left panel: Vertical profiles of normalized Reynolds stresses computed based on viscous layer theory, for which ratio of ice viscosity to water viscosity iw is 10, 100, or 100, with an ice layer that is approximately 3× thicker (*) than the other cases.Right panel (different axes): Normalized Reynolds stress profile computed from experimental data [6].Depths z are normalized by ice thickness h.For stresses ,

Figure 6 .
Figure 6.Left panel: Vertical profiles of normalized Reynolds stresses computed based on viscous layer theory, for which ratio of ice viscosity to water viscosity ν i /ν w is 10, 100, or 100, with an ice layer that is approximately 3× thicker (*) than the other cases.Right panel (different axes): Normalized Reynolds stress profile computed from experimental data [6].Depths z are normalized by ice thickness h.For stresses τ uv , velocities u and v are normalized by open-water velocity amplitude at interface.Cyan shading again indicates ice region (Figure adapted from [6]).

Figure 6 .
Figure 6.Left panel: Vertical profiles of normalized Reynolds stresses computed based on viscous layer theory, for which ratio of ice viscosity to water viscosity iw is 10, 100, or 100, with an ice layer that is approximately 3× thicker (*) than the other cases.Right panel (different axes): Normalized Reynolds stress profile computed from experimental data [6].Depths z are normalized by ice thickness h.For stresses ,

Figure 7 .
Figure 7. Vertical profiles of modeled Reynolds stresses for cases L1-L3 (left panel, T = 1 s, H = 2 cm) and L4-L6 (right panel; T = 1.34 s, H = 6 cm).Stress data are computed over a time range equal to 10 wave periods and normalized in same manner as in Figure 6.Depths are normalized by h in the same manner as in Figure 5.Note different horizontal scales.Cyan shading indicates ice layer region.

Figure 8 .
Figure 8. Normalized attenuation coefficient i i ice k k h = plotted against normalized frequency

Figure 8 .
Figure 8. Normalized attenuation coefficient k i = k i h ice plotted against normalized frequency ω = ω h ice /g on log-log axes.Values are shown for lab experiments of[6,21,32] and field data of[2,33], along with model output from L1-L6 (asterisks) and F1-F2 (squares).Best fit line was computed by[26] based on selected lab and field data.(Figure adapted from[6]).

Figure 9 .
Figure 9. Normalized real wavenumber, k r = k r h ice , plotted against normalized frequency, ω = ω h ice /g.Values are shown for model output from L1-L6 and F1-F2, along with a curve predicted by viscous layer theory [26] and three sets of laboratory cases (Lab 1-Lab 3) with error bars [6] (Figure adapted from [6]).

Table 3 .
Wavenumber in ice and open water.