Determining an optimal transport velocity in the marginal ice zone using operational ice-ocean prediction systems

Knowledge of transport in the marginal ice zone (MIZ) is critical for operations in the Arctic and associated emergency response applications, for example, the transport of pollutants, such as oil, as well as predicting drift associated with search and rescue operations. This paper proposes a general transport equation for the MIZ that can be used for operational purposes in the MIZ. This equation is designed to use a mean velocity of the ice and water velocity, which is weighted by the ice concentration. A key component is the introduction of a leeway coefficient for both the ocean and ice components. 5 These leeway coefficients are determined by minimizing the velocity error between the transport model and observed drifter velocity in the MIZ. These leeway values are found to be 3% of the wind for the water leeway and 2% and 30◦ to the right of the wind for the ice leeway, which are consistent with "rule of thumb" values for surface drifters and sea ice respectively. This general transport model is compared with other transport models and the error is reduced by a factor of 2 compared with traditional transport models for 48 hour lead times. The inclusion of a leeway coefficient in the ice is the key component to 10 reduce trajectory errors in the MIZ. Copyright statement. The works published in this journal are distributed under the Creative Commons Attribution 4.0 License. This licence does not affect the Crown copyright work, which is re-usable under the Open Government Licence (OGL). The Creative Commons Attribution 4.0 License and the OGL are interoperable and do not conflict with, reduce or limit each other. ©Crown copyright 2021 15


Introduction
Estimating the drift and transport of material in the marginal ice zone (MIZ) is a crucial aspect of polar operations. This can be related to search and rescue operations (Rabatel et al., 2018) as well as the transport of oil and other contaminants (French-McCay et al., 2017;Nordam et al., 2019). The requirement of accurate drift and transport estimates for these processes is exacerbated by the remoteness of these regions. Thus, accurate predictions up to at least 48 hours are required to provide 20 sufficient time to coordinate response efforts. Due to the numerous transient processes in the ocean, these short-term predictions can be particularly difficult to predict accurately (Christensen et al., 2018).
The first step in predicting transport in the MIZ is accurate predictions of ice, water and wind velocities, and in the arctic there are a few operational products which provide forecasts of at least 48 hours. A couple examples of such Arctic systems being the Regional Ice-Ocen Prediction System (RIOPS) in Canada (Dupont et al., 2015) and the TOPAZ system as part of 25 the Copernicus Marine Environment Monitoring Service (CMEMS) (Sakov et al., 2012). These operational prediction systems typically will have the ocean and ice components coupled, and typically do not include forcing from surface waves. Rarely they are coupled with the atmosphere, but one example of an ice-ocean-atmosphere coupled model is the Canadian Arctic Prediction System (CAPS) which is the RIOPS system two-way coupled with the GEM (Global Environmental Multi-scale) atmospheric model (Côté et al., 1998b, a;Girard et al., 2014). The ice models used in these prediction systems assume a 30 continuous ice cover and a viscous-plastic rheology, which is typically derived from principles introduced by Hibler III (1979).
For ice concentrations of less than 80%, the ice is assumed to be in "free drift" and the rheology is expected to have a minimal impact on the dynamics (Hibler III, 1979). While free drift models exist, and they perform quite well when compared with observations, they do require data in order to tune the drag coefficients (Schweiger and Zhang, 2015) which can be difficult to obtain in the arctic. 35 It is common for operational prediction of drifting objects to include a leeway term, which is a fraction of the wind that is added to the predicted water velocity. This has been common for a long time in search and rescue (Breivik et al., 2011) and oil spill trajectory modelling . These leeway values can represent direct wind forcing on the object (Kirwan Jr et al., 1975), and/or wind-dependent physics which are not included in the operational prediction system. As most operational marine prediction systems don't include surface waves it is common to include the Stokes drift  40 which can be approximated by a leeway coefficient (Breivik et al., 2011;Sutherland et al., 2020). In addition, material at the surface, such as oil but also ice, will strongly attenuate surface waves creating an additional force on the material (Weber, 2001). Oil and ice also impact the "roughness" of the surface; slicks in the case of oil which reduce the roughness (Wu, 1983), and the presence of form drag in the MIZ which can increase roughness (Lüpkes et al., 2012;Tsamados et al., 2014). For oil, it's long been established that the competing wind and wave effects combine for a net leeway of about 3% of the velocity (Wu, 45 1983). For sea ice it is not clear whether a leeway term should be included in the drift, however, operational models such as CAPS and TOPAZ do not explicitly include processes associated with surface waves or form drag due to the ice roughness. Therefore, a non-zero leeway term in the MIZ may be necessary.
There are only a few examples of using a leeway coefficient with the drift of sea ice, but these are typically restricted to cases where there are no ocean current measurements and the ice speed is modelled as solely a function of the wind. Wilkinson and 50 Wadhams (2003) observed an ice leeway which was dependent on ice concentration. The leeway values ranged from 3.9% for ice concentrations less than 25% to 2.2% for ice concentrations greater than 75%. More recently, Lund et al. (2018) observed that the wind and sea ice motion is not always well correlated and that ice leeway coefficients span a large range of values from 3.5 to 5% across a wide range of ice concentrations in the arctic MIZ. In these papers, it is possible that some of the ice 2 https://doi.org/10.5194/tc-2021-289 Preprint. Discussion started: 8 November 2021 c Author(s) 2021. CC BY 4.0 License. motion was due to tidal and/or inertial motions, which will not directly scale with the instantaneous wind and would lead to 55 uncertainties in the above estimates.
While typically leeway coefficients are determined for particular objects using detailed field measurements (Breivik et al., 2011), recently it was shown by Sutherland et al. (2020) that leeway coefficients can also be predicted using operational prediction systems. The mean leeway coefficients using the operational prediction data were found to be consistent with detailed observations and independent of the choice of operational prediction system. However, these principles have only 60 been applied in open ocean conditions and have not been tested in the presence of sea ice. It would seem that there are certain analogues with the open ocean that would warrant the use of a leeway coefficient in the MIZ. First, most operational ice-ocean models don't include surface waves, which are not uncommon to the MIZ and will impact the drift (Weber, 1987;Squire, 2020). In addition, in the MIZ where ice concentrations are typically less than 80%, the ice is said to be in free drift and the ice motion will depend strongly on a accurate parameterization of the drag coefficient (Cole et al., 2017) and not on an accurate rheology of the internal ice stress (Hibler III, 1979). Form drag due to sea ice is also another physical process which will effect the ice motion and has a strong effect in the MIZ (Lüpkes et al., 2012).
In this paper, we compare the observed drift of four ice drifters in the MIZ and compare these with empirical leeway models.
These drifters span ice concentrations from 10 to 90%. Inputs for the ocean and ice velocities are obtained from two operational ice-ocean prediction systems focused on the arctic. The outline of the paper is as follows. Section 2 a leeway model is developed 70 for the MIZ. Section 3 describes the data, the ice-ocean prediction systems and the methodology used in calculating the leeway coefficients. Results are presented in section 4 followed by the conclusions.
2 Transport equations for the MIZ Ice and ocean surface velocities in the MIZ are strongly coupled. In the dynamical models of operational prediction systems, this is through the relative friction between the ocean and ice components. Most operational ice-ocean prediction systems will 75 have the ice and ocean components two-way coupled. Differences between ice and ocean surface velocities can arise due to ice inertia, which is generally negligible in the MIZ, and internal stresses within the ice, which are also negligible for ice concentrations less than 80% (Hibler III, 1979). So the ice and ocean surface velocities in the MIZ are predominantly in a steady-state free drift mode, with the magnitude, as predicted by operational ice-ocean prediction systems, being dependent on the drag coefficients and the ice concentration.

80
Given the uncertainties associated with modelling velocities in the MIZ, it makes sense to use both ice and ocean velocities to derive a general transport velocity. This is precisely the approach used for oil spill modelling in ice-covered waters (Nordam et al., 2019), where a mean transport velocity is calculated and the ice and surface ocean components are weighted by a function of ice concentration. The model used by the oil spill community uses some empirical estimates for the weighting function as well as for the leeway. We will use this idea of a mean transport model as a template, but will make it more general for a wider 85 use. First, we will go through the oil transport equations and then develop a more general model for transport in the MIZ.

Oil transport equation in the MIZ
The basic equation for the advection velocity in ice-covered water used by the oil spill community (Nordam et al., 2019) is where u o is the oil velocity, u w is the water velocity, u i is the ice velocity, U 10 is the wind velocity at 10 m, α w is a leeway 90 coefficient for oil in water which is typically about 3% Nordam et al., 2019) and k i is the ice transfer coefficient and a function of ice concentration A. Often k i is presented as a piece-wise linear function of ice concentration (Nordam et al., 2019), commonly referred to as the "80/30" rule, and defined as The arguments for the "80/30" rule originated from observations by Venkatesh et al. (1990) and are qualitative in nature. 95 Venkatesh et al. (1990) observed that for sea ice concentrations greater than 80% the oil appeared to drift with the surrounding ice. For ice concentrations less than 30%, the oil drifted with the water. In between these two limits a linear weighting is assumed. The functional form of (2) has not been investigated in detail (French-McCay et al., 2017;Nordam et al., 2019), but any form for k i will inevitably be a monotonically increasing function of ice concentration with the limits k i = 0 when A = 0 (no ice) and k i = 1 when A = 1 (all ice).

General transport equation in the MIZ
For the transport of material in the MIZ we propose a more general equation than (1) that allows for a non-zero leeway coefficient in the ice, where α i is the leeway coefficient in the ice. As mentioned previously, k i in (3) can be any monotonically increasing function 105 of ice concentration A that has the limits k i = 0 at A = 0 and k i = 1 at A = 1. The simplest parameterization of k i that satisfies these conditions is the linear relation k i = A and this will be used for this study.
The inclusion of a leeway coefficient in the sea ice has certain analogues with the open ocean. First, most operational ice-ocean models don't include surface waves, which are not uncommon to the MIZ and will impact the drift (Weber, 1987;Squire, 2020). There are also advantages to using both the ice and ocean components for transport. One potential advantage 110 is that the wind stress in coupled ice-ocean prediction systems (Sakov et al., 2012;Dupont et al., 2015), as well as for wave models (Rogers et al., 2016), is partitioned as a linear function of ice concentration and, therefore, biases in the ice concentration could lead to biases in the ice and ocean velocities. Using a weighted mean ice-ocean velocity ensures that all of the windgenerated currents will be present and biases can be adjusted via the leeway coefficient. Essentially, the assumption is that we will formulate a transport model which assumes that a weighted average of the ice and ocean velocities in the MIZ, along with a corresponding leeway coefficient to be determined, will provide a more accurate estimate for transport velocities than using either the ocean or the ice velocities separately. In some ways this is similar to how wave models calculate their solutions, that is they take a weighted mean of the forcing functions, weighted by ice concentration, and calculate one wave spectrum for both the ice and water components combined, even though on the sub-grid scale the waves are most likely different under ice than under no ice.

Ice-ocean prediction systems
The Canadian Arctic Prediction System (CAPS) is a coupled atmosphere-ice-ocean prediction system, with separate analyses for the atmosphere, ice and ocean components, which is run operationally at Environment and Climate Change Canada (ECCC)  (Côté et al., 1998b, a;Girard et al., 2014).
Snapshots of the ice concentration field (Figure 1) as well as the along-track time series of ice concentration ( Figure 2) are calculated for both CAPS and TOPAZ. A consistent feature between the two ice concentrations is that after approximately 25 September 2018 both CAPS and TOPAZ predict that there is no ice at the drifter locations. This will allow us to test the model of (3) for an object that travels between ice and open water. There are also clearly some differences between the predicted ice 155 concentrations, which are most likely due to the different assimilation cycles of each model. While the analysis of the ice field in CAPS is updated for each forecast (every 12 hours) the TOPAZ fields have a weekly analysis cycle, which could lead to the differences, especially for the two drifters which begin in lower ice concentrations (drifters 14435 and 14438).
A comparison of the model and observed velocities for each drifter is shown in Figure 3. Time series for ice, water and wind (2%) velocities, interpolated in time and space to available drifter observations, are shown. It is apparent from Figure 3 that 160 the observed drifter velocities are nearly always greater than either of the ice and water velocities from CAPS or TOPAZ. Ice and ocean velocities in CAPS vary in magnitude and direction more than those from TOPAZ. This is most likely due to ocean tides being included in CAPS, while these are not in TOPAZ, as this region is known to have large diurnal tidal currents due to a near-resonant, forced topographic wave propagating around the Yermak Plateau (Hunkins, 1986;Padman et al., 1992).

165
As the drifters were deployed on ice floes, it is assumed that the drift of these ice floes in the MIZ can be given by (3) with the drift of ice floes. However, even for ice floes over a wide range of ice concentrations, the weighting of (3) states that the ice floe will drift more like the water in low ice concentrations and more like the ice in high ice concentrations. The primary assumption for using the ice floes as a proxy for general transport is that the ice floe is small enough that inertial forces are 170 negligible. To investigate this assumption, the results from using (3) will be compared with results using (1) as well as an ocean-leeway model (k i = 0, α w = 0.03) and an ice-only model (k i = 1, α i = 0).
The optimal leeway coefficients to be used in (3)

Ice and water leeway coefficients
The MAE between the observed drifter velocity and the model velocity from (3) is minimized for each drifter and choice of ice-ocean prediction system. The optimal leeway coefficients, as well as the MAE in km/day, are located in Table 1 for using the CAPS forcing and Table 2 for the TOPAZ forcing. To demonstrate the relative sensitivity to the choice of leeway coefficients, 2-D slices of the MAE contours are made through a fixed leeway coefficient. That is, the downwind and crosswind contours 185 for the ice leeway are presented for a fixed water leeway and vice versa. Optimal leeway values are found to be approximately α w = 0.03 and α i = 0.02e −iπ/6 and these are used for the constant leeway coefficients used for the MAE slices. No leeway angle is used for α w , but for α i an angle of -30 • is found where the negative implies clockwise rotation relative to the wind as we are using complex notation for our vectors. The MAE for these fixed leeway coefficients are also shown in Tables 1 and 2, denoted MAE * in each    Contours of MAE for α i for the downwind and crosswind components at a fixed value of α w are shown in Figure 4. The MAE contours for each drifter are very similar when using the CAPS data or the TOPAZ data. There is a clear sensitivity related to the ice concentration as the drifter furthest in the ice (14432) has the sharpest contour gradients in MAE and the 195 drifter at the ice edge (14438) shows very little change in MAE for different values of α i . Not including drifter 14438, the optimal value for α i is approximately 0.02 and 30 • to the right of the wind using the CAPS data and 0.03 and 30 • to the right of the wind using the TOPAZ data. The difference between 0.02 and 0.03 are within the 1 km/day MAE contour, which is less than 10% of the total MAE.
Contours of MAE for α w for the downwind and crosswind components at a fixed value of α i are shown in Figure 5. Not 200 surprisingly, the MAE contours in Figure 5 have the opposite sensitivity as those in Figure 4, with sharper contour gradients in MAE for the lowest ice concentration and MAE showing minimal sensitivity to the choice of α w in high ice concentration.
Also, similar to Figure 4, the MAE values in Figure 5 are similar in magnitude when using either CAPS or TOPAZ data and range between 10 and 14 km/day.  and α w = 0.03 and α i = 0.02e −iπ/6 . For comparison, we test an ice-only transport (k i = 1, α i = 0) and ocean-only transport (k i = 0, α w = 0.03). The use of a leeway in the ocean is for consistency with 1).

Short-term trajectory predictions
A few examples of the observed and modelled trajectories for each drifter can be found in Figure 6. A more detailed presentation of the separation distance, d, after 48 hours between the observed and modelled trajectories as a function of 215 trajectory start time is located in Figure 7. A summary of the mean and standard deviation of the separation distance d can be found in Table 3).
For the drifter furthest in the ice, 14432, and the drifter second furthest in the ice, 14437, the separation distance after 48 hours is the smallest for forecasts beginning before 21 September and this does not vary much with choice of transport model or ice-ocean forcing (Figure 7a,b). This changes rapidly after 21 September where the predicted trajectories using the 80/30  and ice-only transport models perform poorly (large separation distance), while the linear model and ocean-leeway models do not produce such a dramatic change (Figure 7a,b). The time when the transport models diverge corresponds with the increase in wind speed ( Figure 3a). As the ice-only and 80/30 transport models do not have a leeway term (these two transport models are identical for A > 0.8), that the leeway term can improve 48 hour trajectories when using operational ice-ocean prediction systems. For the drifter in the lowest ice concentration, 14438, the 80/30, linear, and ocean-leeway models all have essentially the same trajectory error. As the ice concentration, as provided by the model, are less than 0.3 (with the exception of TOPAZ before 21 September) it is expected that these three transport models should reproduce the same trajectory.
In general, the linear model using the best-fit leeway coefficients produces the smallest separation distance as well as the 235 least amount of variation as a function of simulation start time (Table 3). The ice-only transport model had the largest trajectory errors even though we are using drifters on ice floes as our proxy for transport in the MIZ. The ocean-leeway had errors comparable to the linear model and the 80/30 model did not perform well in high ice concentration but did perform well in low ice concentration. These results were generally independent of the choice of CAPS or TOPAZ for the ice-ocean forcing.

240
Presented here is a general leeway model, which assumes a linear weighting of the ice and ocean velocities based on ice concentration and allows for a non-zero leeway in both the ice and ocean components and this leeway can vary in magnitude and direction. Optimal leeway coefficients are calculated by minimizing the error between observed drifter velocities in the waters (Nordam et al., 2019), but we allow for an ice leeway to possibly not be zero to account for missing physics and uncertainties in the ice model and assume a linear weighting for all ice concentrations and not only a subset.
By minimizing the error between the observed and modelled trajectories, optimal values for the water and ice leeway coefficients, α w and α i respectively, were determined. Optimal values for α w were found in the range 0.02 < α w < 0.04 and for α i they were 0.01 < α i < 0.03 with a mean direction of 30 • to the right of the wind. Slightly larger values for α i using the 250 TOPAZ data compared to the CAPS data. For α w , this range of values is consistent with 3% required for the prediction of surface drifters (Sutherland et al., 2020) and oil spill modelling (Nordam et al., 2019). The α i values are curiously consistent with canonical values used for icebergs of about 2% of the wind and 30 • to the right (Leppäranta, 2011).
To assess the quality of the general leeway model we compared the trajectory difference between several transport models for a series of short-term (48 hour) forecasts using the two different ice-ocean prediction systems (CAPS and TOPAZ). The 255 general leeway model with the linear k i = A is used with the optimal leeway values of α w = 0.03 and α i = 0.02 and 30 • to the right of the wind direction. This is compared with the 80/30 transport equation used by the oil spill community, as well as an ice-only transport (k i = 1, α i = 0) as well as an ocean-leeway model (k i = 0, α w = 0.03). Results were generally independent of the choice of ice-ocean forcing. The linear general leeway model consistently had the smallest trajectory error for all four of the drifters. Somewhat surprisingly the ocean-leeway model consistently had the second smallest errors, even for the drifter 260 in the highest ice concentration (14432). This is most likely due to the ice-ocean prediction systems being coupled so that the respective velocities will be strongly correlated in the MIZ when the internal ice stresses are small. The inclusion of the leeway is also key as it will compensate for missing physics, notably from surface waves which are not included in the ice-ocean prediction systems, and can as well compensate for any biases in the respective drag coefficients in the atmosphere-ice-ocean system. The 80/30 model had mixed results and generally had smaller prediction errors for smaller ice concentrations than for 265 large. This is most likely due to the lack of a leeway coefficient in the ice. The ice-only prediction had the largest errors, further emphasizing the importance of a leeway coefficient for accurate short-term prediction of drift in the MIZ.
It is not obvious why the errors associated with calculating the leeway coefficients are similar between the CAPS and TOPAZ forcing (Figures 4 and 5, while the separation errors are much smaller for CAPS than TOPAZ (Figure 7). One likely explanation is due to the drifter observations being available every 3 hours, hence velocities calculated using forward differences are also 270 every 3 hours, while CAPS and TOPAZ provide currents every hour. While averaging the CAPS and TOPAZ forcing would most likely reduce the magnitude of the MAE, it is not expected that this averaging would affect the relative sensitivity, i.e. the optimal leeway coefficients should not change. As CAPS includes tides, and therefore has more variability at these frequencies than TOPAZ (Figure 3), impacts on the magnitude of the MAE will probably be greater for CAPS than TOPAZ. However, simulation of the trajectories don't require the observed drifter velocities and should therefore be more accurate, i.e. the results 275 in Figure 7 and summarized in Table 3.
It is clear from the available data that the inclusion of an ice leeway improves short-term predictions in the MIZ. However, there is not enough information to ascertain whether this is due to missing physics in the prediction systems or correcting for biases that are correlated with the wind. For the ice leeway, α i , it is not clear the physical origin. It is not uncommon to observe "radiation stress" (Weber, 1987), which is the force due to the conservation of wave momentum that arises from the amplitude attenuation due to sea ice. This force can also affect push the ice edge and impact the thickness (Sutherland and Dumont, 2018), which will in turn affect how waves propagate into the MIZ (Sutherland et al., 2019). There are also large uncertainties associated with the drag coefficients (Heorton et al., 2019) that could also be part of the leeway coefficient. It should also be emphasized that while the use of an ice leeway improved predictions in the MIZ, there is little reason to expect that this would 285 be the case in the pack ice where internal stresses significantly impact the drift. In such a case, a more sophisticated expression for α i , which depends on ice concentration such that α i → 0 as A → 1 would most likely be required that also preserves some of the improvements shown in the MIZ. Author contributions. This article was written by GS with input from all co-authors. GS conceived of the idea, with feedback from all coauthors, and wrote the software for data analysis and production of figures. VA, LH and MD contributed to preliminary data analysis. JR 295 designed the drifters and was responsible for data collection and quality control. JR, LH and ØB contributed to drifter deployment.