Parameterizing Tabular‐Iceberg Decay in an Ocean Model

Large tabular icebergs account for the majority of ice mass calved from Antarctic ice shelves, but are omitted from climate models. Specifically, these models do not account for iceberg breakup and as a result, modeled large icebergs could drift to low latitudes. Here, we develop a physically based parameterization of iceberg breakup based on the “footloose mechanism” suitable for climate models. This mechanism describes breakup of ice pieces from the iceberg edges triggered by buoyancy forces associated with a submerged ice foot fringing the iceberg. This foot develops as a result of ocean‐induced melt and erosion of the iceberg freeboard explicitly parameterized in the model. We then use an elastic beam model to determine when the foot is large enough to trigger calving, as well as the size of each child iceberg, which is controlled with the ice stiffness parameter. We test the breakup parameterization with a realistic large iceberg calving‐size distribution in the Geophysical Fluid Dynamics Laboratory OM4 ocean/sea‐ice model and obtain simulated iceberg trajectories and areas that closely match observations. Thus, the footloose mechanism appears to play a major role in iceberg decay that was previously unaccounted for in iceberg models. We also find that varying the size of the broken ice bits can influence the iceberg meltwater distribution more than physically realistic variations to the footloose decay rate.


Introduction
Icebergs serve as a major pathway for how freshwater from ice sheets is distributed into the ocean. About half of the current mass loss from both the Antarctic and Greenland ice sheets is attributed to calving icebergs (Depoorter et al., 2013;van den Broeke et al., 2016). After calving, icebergs can drift for thousands of kilometers before disappearing as a result of melting and breaking into smaller pieces. The freshwater released by icebergs as they melt can modify ocean circulation, enhance sea-ice formation, and affect biological primary productivity (e.g., Arrigo et al., 2002;Stern et al., 2015).
Icebergs could potentially trigger large shifts in climate. Heinrich events, which describe abrupt and massive discharges of icebergs from the Laurentide Ice Sheet during the last glacial period (Heinrich, 1988), weakened the Atlantic Meridional Overturning Circulation (AMOC), leading to reduced heat transport to the pole and regional cooling (e.g., Levine & Bigg, 2008). It is possible that a similar slowdown of the AMOC will occur over the next century, due to increases in meltwater from Greenland glaciers and icebergs (Golledge et al., 2019). In the Southern Hemisphere, Antarctic icebergs drifted farther north at the start of all glacial periods over the past 1.6 million years before decaying (Starr et al., 2021). The resulting redistribution of iceberg meltwater was associated with reorganizations in deep-water mass, which increased sequestration of CO 2 from the atmosphere and may have facilitated the development of ice age conditions. In recent years, many ocean and climate models have included iceberg representation such as tracking the iceberg drift, melt, and freshwater flux to the ocean (e.g., Jongma et al., 2009;Marsh et al., 2015;Martin & Adcroft, 2010). While these models have confirmed the role that iceberg meltwater plays in climate, a major shortcoming is that they severely underestimate the initial size of calved icebergs from the Antarctic ice sheet. Calving size influences climate as a major control on iceberg drift patterns, and therefore where iceberg meltwater is deposited (Stern et al., 2016). The current iceberg modules of climate models typically assign the calving size of new icebergs according to a prescribed size distribution, where the maximum possible horizontal iceberg area is ∼3.5 km 2 . However, the observed Antarctic calving size distribution follows a −3/2 power law (i.e., the probability that an iceberg has an area A is proportional to − 3 2 ), where giant icebergs with areas larger than 100 km 2 make up 89% of the total volume of all icebergs, despite being statistically rare .
Three primary issues have prevented the inclusion of large tabular icebergs (defined here as having a horizontal surface area greater than about 3.5 km 2 ) within climate models: (a) large icebergs may be subject to differing ocean forcings along their edges and can disrupt ocean circulation, so that the common approach of modeling icebergs as Lagrangian point particles that do not exert pressure on the ocean may be inappropriate; (b) the frequency of giant iceberg calving events and spatial variations in calving-size are poorly constrained; and (c) large icebergs modeled as Lagrangian point particles travel to unrealistically low latitudes, likely due to lack of a physically justified parameterization for iceberg fracture and breakup (Stern et al., 2016). Here, we address the third issue by developing a physically based parameterization, appropriate for climate models, of the "footloose mechanism"-a process by which large icebergs can rapidly deteriorate by calving smaller icebergs from their edges (Scambos et al., 2005;Wagner et al., 2014). We build on a previous approach to parameterize this mechanism, where icebergs are treated as elastic beams to calculate the size of each "child" iceberg-an iceberg that breaks off from a parent iceberg . However, unlike this previous approach, we also develop a physically based parameterization of the footloose calving rate, which depends on ocean-induced melt and erosion. Furthermore, we develop a technique to increase the computational efficiency of the parameterization by tracking these footloose ice bits (that we term "child icebergs") in clusters rather than as individual icebergs, thereby facilitating efficient implementation into coupled general circulation or Earth system models. We test our methods within the Geophysical Fluid Dynamics Laboratory global ocean-sea ice model (GFDL OM4), which serves as the ocean/sea ice/iceberg component for the GFDL climate and Earth system models (Adcroft et al., 2019), where the iceberg module (Martin & Adcroft, 2010) is largely based on . However, our methods may also be easily applied within most iceberg components of other ocean and climate models which follow the same  formulation (e.g., Marsh et al., 2015). This paper is outlined as follows: in Section 2, we overview the existing iceberg module of OM4; in Section 3, we describe the footloose parameterization; in Section 4, we carry out experiments to demonstrate the impact of the footloose parameterization on iceberg trajectories, decay, and meltwater flux; in Section 5, we discuss the results; and in Section 6 we offer concluding remarks.

The GFDL Iceberg Module
The GFDL iceberg module (Martin & Adcroft, 2010) is a sub-component of the Sea Ice Simulator version 2 (SIS2) code, which is coupled with the Modular Ocean Model version 6 (MOM6) in OM4 (Adcroft et al., 2019). The atmospheric state and runoff are prescribed using the JRA-55 surface atmospheric data set (Tsujino et al., 2018). Calving of icebergs from ice sheets is modeled by accumulating frozen runoff in coastal grid cells, where it is subdivided into a fixed number of specified iceberg size classes according to a statistical distribution. Whenever the accumulated ice mass reaches the given mass for a size class, it is released as an iceberg. Each iceberg is modeled as a cuboid Lagrangian particle, which is assigned an initial length-to-width ratio of 1.5 and a thickness provided by the respective iceberg size class. There are generally many more small icebergs than large icebergs. To reduce computational cost, small icebergs within the same size class are released in clusters. Each cluster is represented by a single particle, which takes the dimensions of a single iceberg in the cluster. A "mass scaling factor" defines the number of icebergs in the cluster, and is used to scale the meltwater flux to represent the entire cluster.
The model for iceberg evolution (Martin & Adcroft, 2010), described in the remainder of this section, is based on the equations given in  and , with updates for capsizing given by Wagner, 10.1029/2021MS002869 3 of 19 Stern, et al. (2017). The icebergs drift in response to oceanic and atmospheric forces, and their dimensions decrease over time according to a decay parameterization that includes melt and erosion terms. Iceberg meltwater is distributed to the ocean surface. The orientation of the icebergs is not tracked because they are treated as point particles during any interpolations to and from the background grid.

Momentum-Balance Formulation
The momentum balance for each iceberg is given by where is the total (Lagrangian) derivative, M is the iceberg mass, is the iceberg velocity, ⃗ = − × ⃗ is the Coriolis force (where f denotes the Coriolis parameter), ⃗ is the pressure gradient force due to the sea surface slope, ⃗ is the force from air drag, ⃗ is the water drag, ⃗ is the sea-ice drag, and ⃗ is the wave radiation force. The pressure gradient force is given as where η is the sea surface height. The drag forces are given by Subscripts a, o, and si are associated with air, ocean surface, and sea ice, respectively. Parameter ρ x (with x = a, o, si) denotes density, are velocities, and c x,h and c x,v are horizontal and vertical drag coefficients, respectively. We set c a,h = 0.0055, c a,v = 1.3, c o,h = 0.0012, and c o,v = c si,v = 0.9 . Sea-ice drag is assigned the same coefficient as water, following . The sea-ice force parameterization is likely incorrect in the limit that the iceberg is locked in strong sea ice (Morison & Goldberg, 2012). Parameters L, W, T, F, and D are iceberg quantities corresponding to length, width, thickness, freeboard, and draft. These quantities are defined so that L ≥ W, T = F + D, and = , where ρ is the iceberg density. Lastly, the wave radiation force is defined as where g is gravitational acceleration, = 0.010125|⃗ − ⃗ | 2 is the empirical wave amplitude, and c r is the wave drag coefficient Here, L c = 0.125L w is a lower cutoff length, and L t = 0.25L w is a upper length limit, where = 0.32|⃗ − ⃗ | 2 is an empirical wavelength.
We solve the momentum equations with a velocity Verlet time-stepping scheme used previously for a bonded-particle iceberg model . However, unlike the bonded-particle model, the point-particle model presented here does not account for interaction or collision between particles, and particles do not apply pressure to the ocean surface, which would require too small of time steps for climate modeling. Grounding of icebergs is also neglected.

Melting and Erosion Parameterization
Melt and erosion rates are parameterized following Martin and Adcroft (2010) and . These rates are used to calculate the latent heat and freshwater fluxes to the ocean surface from iceberg decay, and to evolve the dimensions of each iceberg. At the iceberg base, turbulence caused by the relative motion of the water and the iceberg generates a heat flux to the iceberg, and the resulting basal melt rate (m/day) is estimated by where ̃= −4 • C is the effective iceberg temperature and ̃ is the sea surface temperature. A heat flux also occurs at the iceberg side walls, which is related to buoyant convection caused by the temperature difference between the iceberg and ocean. The associated melt rate (m/day) is approximated by Further decay of the sides is caused by wave erosion and melting due to the contact with warm surface waters, which carve out a notch at the waterline of the iceberg known as a "wavecut." Once the wavecut grows large enough, the overhanging ice slab becomes unstable and breaks off. The rate of ice loss (m/day) from this erosion process is parameterized by where A i is the fractional sea-ice area and S s is the sea state, which is estimated by a fit to the Beaufort scale A percentage of the mass flux from wave erosion forms small child icebergs known as "bergy bits" (Martin & Adcroft, 2010). This percentage is poorly constrained, and we arbitrarily set it to 100% here. Since bergy bits can melt rapidly due to their small size, any wave erosion that should lead directly to melt, but that we convert into bits, is melted in due course. The bergy bits are assigned a cuboid geometry with side lengths equal to the shortest dimension of the parent iceberg or 40 m, whichever is smaller. Because they are so small, bergy bits melt according to Equations 8 and 9 alone, and we assume that they travel alongside their parent iceberg rather than representing them with a new particle. After iceberg dimensions are updated = − and = = − − , capsizing (rolling) of icebergs is enforced by swapping W and H whenever the ratio < √ 6 (1 − ) , where = (Wagner, Stern, et al., 2017).

The Footloose Mechanism
After an overhanging ice slab breaks off on an iceberg during the wave erosion process parameterized in Equation 10, a protruding underwater ice foot is left behind that is still attached to the parent iceberg (Figures S1a and S1b in Supporting Information S1). Buoyancy forces associated with this submerged foot induce a bending torque on the iceberg, producing internal stresses proportional to the length of the foot ( Figure S1c in Supporting Information S1). Once the foot grows long enough, a critical stress is reached that triggers calving of a smaller child iceberg from the parent iceberg (Figure S1d in Supporting Information S1). This process of calving is known as the "footloose" mechanism (Wagner et al., 2014) or "edge-wasting" (Scambos et al., 2005).

One-Dimensional Beam Approximation
Assuming that footloose calving is an elastic process, the critical length of the foot needed to trigger footloose calving can be estimated in one dimension using beam theory (Sergienko, 2013;Timoshenko & Goodier, 1970;Wagner et al., 2014) where σ y is the iceberg yield stress and is the buoyancy length, which describes the balance between the stiffness of the beam and loading from hydrostatic pressure. In Equation 13, B = ET 3 /12(1 − ν 2 ) is the iceberg stiffness, where E is Young's modulus and ν = 0.3 is Poisson's ratio. The length of the footloose child iceberg that is calved (Wagner et al., 2014) is given by For a typical tabular iceberg with a thickness of 200-300 m, is typically on the order of a few 100 m (Section 4.1).
A schematic of and is provided in Figure S2 in Supporting Information S1.
Equations 12-14 are formulated under the assumption that ≪ . This "small foot" assumption is most accurate for thick icebergs, such as giant tabular icebergs, which are generally over 200 m thick. Note that with all else being equal, icebergs with smaller l f will calve child icebergs more frequently that icebergs with larger , and are therefore associated with an greater rate of footloose mass loss. Similarly, with all else being equal, icebergs with larger l c will have a greater footloose mass loss rate than those with smaller , because each child iceberg calving event is larger.

Parameterization for the Iceberg Modules of Climate Models
We use Equations 9-14 to parameterize the footloose mechanism for the iceberg modules of climate models. We develop our parameterization under the assumption that a submerged foot grows along all sides of the parent berg, where child bergs are calved roughly regularly over time and randomly along the edges. This assumption is consistent with observations, where edge-wasting generally decreases the parent iceberg area without changing its shape (Scambos et al., 2005).
We track the area of the submerged foot, A f , for an iceberg as where A e is the horizontal area of ice lost from the side walls due to the accumulated ice eroded from the iceberg freeboard (calculated using Equation 10) and A v is the ice lost to melting along the side of the iceberg below sea-level (calculated using Equation 9). Generally, A v ≪ A e (Martin & Adcroft, 2010). A schematic of these areas is provided in Figure S3 in Supporting Information S1. Following , we assume that a newly calved child iceberg has horizontal dimensions of l c × 3l c and the same thickness as their parent. Consistently, we calve a child iceberg whenever A f is greater than or equal to 3l c l f . Upon calving, we subtract the area of the calved foot, 3l c l f , from A f , and we initialize the position of the child iceberg to a random point along the perimeter of the parent berg. Then, we adjust the shape of the parent iceberg to account for the mass lost to the calving event by subtracting the same distance, ds, from both the length (L) and width (W) of the parent berg: This adjustment to the parent area preserves the shape of the iceberg over time, and therefore is more realistic than the previous footloose scheme from , where only the length of the parent iceberg was adjusted to account for the mass lost to footloose calving.
The footloose mechanism does not apply to many small icebergs, even if L > l c , due to interference between the buoyancy forces from opposite sides of the iceberg. Following , we assume this interference inhibits footloose calving when the parent icebergs dimensions no longer satisfy L > 3l c . In addition, we also introduce the condition: if W ≤ 3l c and L > 3l c , we allow footloose calving only from the length-axis of iceberg, and not the width-axis. In this case, we only track the area of the submerged foot along the shortest edges, that is, perpendicular to the length-axis. This case is an exception to the idea that the iceberg length and width should 10.1029/2021MS002869 6 of 19 decrease evenly over time using Equation 16. Instead, we must reduce the length of the parent iceberg by 3 2 ∕ without modifying the width.

Tuning Process
The parameterization has two poorly constrained parameters: the Young's modulus (E), and yield stress (σ y ) for ice. Inferences from remote sensing observations suggest that the laboratory values of E are significantly higher than those required to produce the observed ice-shelf flexure (Sergienko, 2006). The laboratory experiments on the strength of ice are inconclusive about characteristic values of σ y (Murdza et al., 2021). As indicated by Equations 12 and 13, l f is proportional to σ y and inversely proportional to E 1/4 . In addition, from Equations 13 and 14, l c is proportional to E 1/4 . Given their influence on l f and l c , E and σ y together control the rate of parent iceberg mass loss due to footloose calving. Furthermore, with l c proportional to E 1/4 , then E also sets the parent iceberg length threshold that must be exceeded to allow footloose calving (L > 3l c ).

Increasing Computational Efficiency
Accounting for the footloose mechanism greatly increases the number of icebergs that must be tracked during the simulation. Tracking the drift and decay of each individual footloose child berg, as done in , is too computationally expensive for climate models (see Section 4). Therefore, we propose a scheme to reduce the number of particles in the simulation, which is similar to the small iceberg clustering scheme used when calving from ice sheets. Each parent iceberg is associated with a single "footloose bin," where the overall mass of child icebergs accumulates until the binned mass surpasses a given threshold, M th . When this threshold is exceeded, the cluster of binned icebergs are released from the parent as a single particle. Periodically converting the binned child icebergs into a new particle, according to M th , helps account for the spread of child icebergs from the parent bergs while still reducing the number of particles in the simulation.
The drift and decay of clustered child icebergs are modeled using the properties of a characteristic iceberg to represent the entire group. When a cluster is binned alongside the parent iceberg, the geometry of the characteristic child iceberg is always approximated as equal to that of a freshly calved child iceberg, according to the thickness of the parent iceberg. This geometry is needed for calculating the child melt and erosion (Equations 8-10), where the erosion is converted into bergy bits. The geometry must be approximated because new child icebergs are added to the bin at different times. Iceberg sizes within the bin will therefore vary depending on the calving size and decay history of each binned iceberg. Lower values of M th will release binned icebergs as new particles more often, thereby minimizing any size differences between binned icebergs, but at the cost of increasing the computational expense. The meltwater flux from the characteristic iceberg is scaled to represent the entire group according to a mass scaling factor, defined as the total mass of the binned icebergs over the mass of the characteristic iceberg. When a cluster is released from the bin as a new particle, it initially retains the dimensions of the characteristic iceberg. Then, it drifts and melts completely separately from its parent, like the clustered icebergs that calve from ice sheets (Martin & Adcroft, 2010).

Experiments Setup and Results
In this section, we compare the performance of the iceberg model within OM4, with and without the footloose parameterization. We test several calving-size distributions, and a range of realistic values for the Young's modulus and yield stress used in the footloose parameterization. We compare the modeled iceberg trajectories and sizes to observations, and analyze the differences in meltwater fluxes between the experiments.

Experimental Setup
For all experiments, the MOM6 ocean model and SIS2 sea-ice model are run on the same horizontal C-grid with 35 isopycnal layers in the vertical. The horizontal resolution is 0.5°, where some additional refinement, mostly in the meridional direction, is applied at the poles and equator following the "OM4p5" configuration detailed in Adcroft et al. (2019). Each experiment is split into two 60-year runs, where 60 years  is the length of the JRA-55 data set used to prescribe runoff and surface atmospheric forcing. The first of the two runs serves as a spin-up, where the iceberg mass on the ocean increases from zero at the start of the run to an equilibrium state after 60 years, which is the same time to equilibrium obtained in a previous study (Martin & Adcroft, 2010). The final states of the iceberg, sea-ice, and ocean components at the end of the first run then serve as the respective initial states for the second run.
All experiments use a calving-size distribution in the Northern Hemisphere based on   (Table  S1 in Supporting Information S1). In the Southern Hemisphere, we test two distributions mentioned in the introduction: the "Gladstone distribution" (Table S2 in Supporting Information S1;  and the "power law distribution" (Table S3 in Supporting Information S1; . The Gladstone distribution does not include icebergs with horizontal areas over 3.5 km 2 , thereby severely underestimating the size of Antarctic tabular icebergs, which can have areas several orders of magnitude larger. However, this distribution has widely been used in ocean models (Marsh et al., 2015;Martin & Adcroft, 2010). We therefore test it, without the footloose parameterization, for comparison with the experiments that use the more realistic power law distribution and the new footloose parameterization. The power law distribution is based on the observed −3/2 power law , and includes iceberg areas ranging from 0.3 to 1,000 km 2 . Following , icebergs with areas greater than 1,000 km 2 are neglected under the assumption that these rare icebergs break up into pieces smaller than 1,000 km 2 before they enter the open ocean, or soon afterward.
The experiments are summarized in Table 1, which gives the experiment name, the Southern Hemisphere distribution used, and five additional parameters for experiments that use the footloose parameterization: the Young's modulus (E), the yield stress (σ y ), the foot length to trigger calving (l f ), the length of new child icebergs (l c ), and dimensionless measure of the rate of mass loss from the footloose mechanism (l c /l f ). These last three parameters are normalized to the highest values of all experiments, which eliminates any dependence on ice thickness. However, for a typical 250 m thick tabular iceberg, the range of values for l f and l c vary between about 33-93 m and 216-384 m, respectively, between all experiments. These values are calculated using Equations 12-14, where we use constant densities for the icebergs and ocean of ρ = 850 kg/m 3 and ρ o = 1,025 kg/m 3 , respectively. All footloose experiments are run with the child iceberg binning scheme using M th = 1 × 10 12 kg. Experiment E100_S.5 is additionally run without the binning scheme to compare computational expense. We label this "no binning" run as "E100_S.5_nb." As indicated in Table 1, we test a range of values for E (10 MPa ≤ E ≤ 100 MPa) that are much smaller than typically assumed for non-damaged ice (1 GPa ≤ E ≤ 10 GPa). Modeling studies have found that such low values for E are necessary to simulate observed rampart-moat profiles (Mosbeux et al., 2020;Scambos et al., 2005;Wagner et al., 2014). In part, a lower E is expected because E decreases as ice temperature increases, and under strainrate effects over long loading times (Sinha, 1978), such as those associated with the formation of submerged feet. However, the primary justification for lowering E is likely related to crevassing, which effectively decreases the thickness of ice that is able to support bending stresses. In both the one-dimensional footloose beam approximation and two-dimensional full-Stokes models, decreasing E has a similar effect as decreasing the ice thickness (Mosbeux et al., 2020). Crevassing also reduces σ y because the ice is damaged and will therefore fracture more easily. While σ y is typically set on the order of 1 MPa for undamaged ice, we set σ y as low as 0.1 MPa to account for the presence of already-damaged ice, following van der Veen (1998).
All experiments use the  iceberg size distribution in the Northern Hemisphere. b Defined as l c /l f , and reported relative to experiment E10_S.1.

Table 1
Descriptions of the OM4 Iceberg Experiments

Results
The results described below use the iceberg areas and melt fluxes analyzed for the final 60 model years, that is, after the spin-up period.

Southern Hemisphere
We first compare modeled iceberg areas to observations collected by the National Ice Center (NIC). The NIC detects large icebergs in satellite imagery of the Southern Hemisphere and reports their weekly positions, lengths, and widths. The majority of observations are from the early 1990s to the present, and were compiled into a single database by Budge and Long (2018), to which we compare our results. Figure 1 shows the average area of icebergs in the Southern Hemisphere that drift within 100 km of each grid point from 1990 to 2018 for each experiment (Figures 1a-1e) and the NIC (Figure 1f). Here, average area is calculated for icebergs with areas greater than 200 km 2 (the minimum iceberg area consistently detected by the NIC) and less than 1,000 km 2 (the upper limit on iceberg area in the simulations). Therefore, the icebergs described by the Gladstone distribution are excluded because it does not include icebergs with areas greater than ∼3.5 km 2 . For all experiments, the pattern of simulated iceberg drift generally matches observations, where large icebergs travel westward (counter-clockwise) in the Antarctic Coastal Current (ACoC) before escaping into the eastward-flowing (clockwise) Antarctic Circumpolar Current (ACC), primarily through the Weddell Sea "iceberg alley." The NIC data set (Figure 1f) also shows that many large icebergs exit the ACoC near the Kerguelen Plateau and the Ross Sea. However, this behavior is less apparent in the simulated iceberg drift, likely due to difficulty in resolving coastal currents and the simplified representation of large icebergs as non-interactive levitating point particles (see Section 5). Nevertheless, it is evident that the footloose mechanism has a strong impact on the spatial spread of large icebergs. The large icebergs from the "No footloose" power law experiment (Figure 1e) reach far lower latitudes than the large icebergs in the footloose experiments (Figures 1a-1d), which have trajectories that more closely resemble the NIC observations (Figure 1f). The drift distance of icebergs in the footloose experiments is inversely correlated with the footloose decay rates given in Table 1. Simulations with slower footloose decay rates, such as E10_S.25 (Figure 1c), show icebergs that survive further into lower latitudes than runs with faster footloose decay rates, such as E10_S.1 (Figure 1d). Note that the two experiments used to test the binning scheme, E100_S.5 (with child binning) and E100_S.5_nb (no child binning), generally show similar spatial extents of large icebergs, as well as similar average areas in well-trafficked areas such as along the Antarctic coast and the Weddell Sea. The few differences between these runs can be explained by the random differences in behavior of a few single icebergs. Most notably, one large iceberg escapes the ACoC near Kerguelen Plateau in E100_S.5, but not in E100_S.5_nb. This iceberg drifted north until being deflected to the east by the ACC near 60°S.
Overall, experiments E10_S.25, E100_S.5, and E100_S.5_nb best replicate the observed spatial spread of large icebergs from NIC observations. The average iceberg areas of these experiments generally match observations as well, especially in the ACoC and the Weddell Sea region where average iceberg area consistently ranges between 400 and 600 km 2 . These regions are well-trafficked by icebergs, which facilitates comparisons of average iceberg areas because the influence of outliers is minimized. Unfortunately, a more detailed comparison in area between the experiments and NIC observations is challenging because the size, frequency, and locations of tabular calving events are inconsistent between modeled and observed tabular calving events from the Antarctic ice sheet. Furthermore, comparing areal changes in individual icebergs between the simulated and NIC data sets is not possible, because NIC areas are calculated by hand from satellite data by estimating the shape of observed icebergs as rectangular, and are not accurate enough to capture gradual decay via the footloose mechanism. Nevertheless, all power law experiments yield a significantly improved match to observations of large iceberg spatial extent and areas when the physically based footloose parameterization is activated, suggesting that the footloose mechanism is an important process in iceberg decay that was previously neglected in iceberg drift models.
The average melt flux (mm/day) over the final 60 years for all experiments is given in Figure 2 on a logarithmic scale. First, note that nearly identical meltwater fluxes are obtained when using the footloose binning scheme (E100_S.5; Figure 2a) versus evolving each child iceberg individually (E100_S.5_nb; Figure 2b), with the exception of the meltwater input from the single iceberg that escaped the ACoC near the Kerguelen Plateau only in the E100_S.5 experiment. However, as reflected in year for E100_S.5 is about an order of magnitude faster than E100_S.5_nb and two times faster than when using the Gladstone distribution.
When comparing all Power Law experiments, the melt flux distribution (Figure 2) does not appear to be strongly correlated with the spatial spread of large icebergs shown in Figure 1, with the exception of the "No footloose" case. In the "No footloose," both the large iceberg drift (Figure 1e) and melt flux (Figure 2e) are concentrated  (Figure 2b). This increased melting further from the coast is especially apparent in the three most common regions for icebergs to exit the ACoC: the Kerguelen Plateau and Ross and Weddell Seas. The similarity between the E10_S.25 and E10_S.1 melt patterns occurs despite the fact that E10_S.25 shows the greatest spatial spread of large icebergs from the coast of all footloose runs (i.e., the slowest footloose mass loss rate; Figure 1c) while E10_S.1 shows the least spatial spread (i.e., the fastest footloose mass loss rate; Figure 1d). We conclude that the footloose mass-loss rate and drift patterns of parent icebergs are not the primary mechanisms that drive the differences between the meltwater flux results, despite their influence on where child icebergs are calved. Instead, these differences in melt pattern appear to be driven by the size of the footloose child bergs upon calving, l c . According to Table 1, the l c associated with experiments E10_S.25 and E10_S.1 is nearly half the l c for experiments E100_S.5 and E100_S.5_nb. A previous study that used the same iceberg module (Stern et al., 2016) showed that smaller icebergs are more easily driven away from the coast by katabatic winds, while larger icebergs tend to stay entrained in coastal currents and only rarely are driven north, by ocean current meanders caused by topographic features. This notion that the drift of differently sized icebergs are sensitive to different forcings is well-known, and can be inferred from the iceberg momentum balance (Equation 1). First, note that under approximate geostrophy in the ocean, the pressure gradient force can be rewritten as = × ⃗ . Furthermore, note that mass increases with iceberg size more than surface area does. Consequently, for small icebergs, strong wind drag forces, which are proportional to surface area, are often larger than the "Coriolis-related" forces that are proportional to mass, F c and F p . Conversely, the "Coriolis-related" forces are dominant for large icebergs. This effect is discussed in detail in several studies Stern et al., 2016;.
In Figure (Figures 1a and 1b). In other words, larger child icebergs are more likely to follow the drift of their parents. Additionally, Figure 3 further reflects that the large child iceberg experiments show decreased melt in the ACC downstream of the Kerguelen Plateau compared to the Gladstone experiment, given that there are fewer small icebergs that can drift out of the ACoC into this region.
Figures 4 and 5 give the child and parent melt fluxes, respectively, for the footloose experiments. The percentage of the total melt flux attributed to child or parent icebergs alone is given in each of the respective figures, for each experiment. In all cases, the child iceberg melt dominates over the parent melt flux proportionately to the footloose decay rates (Table 1). For the small child iceberg experiments, melt is increased within the ACC downstream of the Kerguelen Plateau (north of 60°S and between 90° and 135°E) for not only the child icebergs (Figures 4c and 4d), but the parents as well (Figures 5c and 5d). This increased parent melt in the ACC occurs because footloose calving is deactivated when parent L ≤ 3l c . Therefore, when l c is small, parent icebergs are allowed to lose more mass from the footloose mechanism alone as compared to when l c is large. The resulting smaller parent icebergs are then able to escape the ACoC more easily, as their smaller size means their drift is more sensitive to Katabatic winds.

Northern Hemisphere
The average area of icebergs in the Northern Hemisphere that drift within 100 km of each grid point from 1990 to 2018 is given for each experiment in Figure 6. Here, average area is calculated using icebergs with areas greater than 0.5 km 2 (note that all experiments use the Bigg calving size distribution in the Northern Hemisphere). The Gladstone result is not reported here because it is essentially identical to the "No footloose" experiment; both experiments use the same calving-size distribution without footloose in the Northern Hemisphere. Furthermore, unlike in the Southern Hemisphere, an observational data set is not available.
Like the Southern Hemisphere results, the spatial spread of icebergs in the Northern Hemisphere footloose experiments (Figures 6a-6d) is decreased compared to the case without footloose (Figure 6e). However, the differences in spatial spread of large icebergs between footloose experiments is not as strongly dependent on the footloose calving rate as it is the Southern Hemisphere ( Figure 1). Instead, large iceberg spread is mostly dependent on child iceberg calving size l c because the iceberg size distribution is smaller, so that the parent length threshold on   E100_S.5,E100_S.5_nb,and Gladstone deactivating the footloose mechanism is the dominant control on how far icebergs drift. In experiments where l c is low (E10_S.25; Figure 6c and E10_S.1; Figure 6d), more footloose calving is allowed before footloose is deactivated. Consequently, the parent bergs decay quickly, and both the small parent and child icebergs melt rapidly once they encounter the high sea surface temperatures south of Greenland. In other words, these small l c icebergs are unable to survive to as low of latitudes as the icebergs in the larger l c experiments (E100_S.5; Figure 6a and E100_S.5_nb; Figure 6b). Thus, in the Northern Hemisphere, the smaller l c icebergs deposit more of their meltwater at higher latitudes than the larger l c icebergs. This melt pattern is shown in Figure 7, which gives the  (Figures 7b-7d). While all footloose experiments show increased/decreased meltwater at higher/lower latitudes as compared to the "No footloose" run, these differences are more exaggerated for the small l c experiments (Figures 7d and 7e) compared to the large small l c experiments (Figures 7b and 7c).
This strong dependence on l c in the Northern Hemisphere is further supported when considering the child and parent melt fluxes separately (Figures 8 and 9, respectively). The percentage of the total melt flux attributed to child or parent icebergs alone is given in each of the respective figures, for each experiment. Unlike in the Southern Hemisphere, the parent melt flux is dominant in the Northern Hemisphere for all experiments because the smaller Bigg size distribution causes parent icebergs to reach the length threshold to deactivate footloose calving more quickly, and less footloose calving occurs overall. However, experiments characterized by smaller values for l c can calve more footloose child icebergs, and therefore show a relatively high percentage of child melt (Figures 8c and 8d) as compared the experiments with larger l c (Figures 8a and 8b).

Discussion
Of all experiments, the E10_S.25 experiment (Figure 1c) appears to yield the closest match to observed (Figure 1f) drift and areas of large icebergs. Consequently, it may be reasonable to choose the values for the footloose parameters used in this experiment as the recommended tuning: E = 10 MPa and σ y = 0.25 MPa. However, note that it is possible to assign other values for these parameters that give a similar footloose decay rate, and therefore similar drift patterns and areas, but with a different l c that changes the melt flux distribution. Therefore, we consider the suite of footloose experiments presented in this paper as a sensitivity test for the footloose parameters rather than a rigorous tuning exercise. Further constraining the values for these parameters could involve using high resolution observations to determine a range of target sizes for l c . However, such an approach will be complicated in the likely case that E and σ y vary significantly between different icebergs. Nevertheless, combining the modified power law calving-size distribution with the footloose parameterization in this study, regardless of the exact tuning chosen, has significantly improved upon the representation of large icebergs as compared to the previously used Gladstone distribution without breakup. Therefore, the methods developed here are ready to be incorporated within the iceberg modules of climate models.
Although implementation of a physically based parameterization of the footloose mechanism has improved the representation of large icebergs, several processes remain unaccounted for in the model. For example, we do not properly account for icebergs larger than 1,000 km 2 despite the fact that icebergs as large as 11,000 km 2 , while rare, have been observed. Our assumption, following , that these giant icebergs break into fragments no greater than 1,000 km 2 shortly after calving, is a limitation of the current study. Occasionally, giant icebergs with areas over 1,000 km 2 (e.g., iceberg A68a) do enter the open ocean. Properly accounting for giant icebergs with areas over 1,000 km 2 would likely require an ice shelf calving parameterization that can predict when and where these infrequent, but massive, icebergs are created. Furthermore, a Lagrangian point particle model is not suitable for representing the dynamics of giant icebergs, which should exert pressure on the ocean and sea ice, and endure differing ocean forcings along their edges (see Section 1). While bonded-particle methods have been developed that account for these ocean/iceberg interactions , they require time step sizes of fractions of a second and are therefore not yet feasible in climate models, which typically use time step sizes on the order of an hour. One known improvement to our representation of large icebergs in the current model would be to replace the ocean surface forcing in the iceberg momentum equations with an ocean forcing integrated over the depth of the iceberg, which has been shown to displace thick icebergs further from the coast in both the North Atlantic (Marson et al., 2018) and Southern Ocean (Merino et al., 2016). Implementing a depth-integrated ocean forcing may therefore help large icebergs exit the ACoC more often in the Ross Sea and near the Kerguelen Peninsula, a process which the iceberg model currently under-represents (Section 4.2.1). Accounting for grounding and increasing model resolution to better resolve the coastal currents would also likely help to realistically divert large icebergs out of the ACoC. While these additions are likely to improve model fidelity, we think they would not change the sensitivity results presented here concerning the footloose parameterization.

Summary
We present a parameterization of footloose calving for the iceberg components of climate models. We adapt a 1-D elastic beam approximation to calculate the size of the submerged ice foot needed to trigger calving and the size of the calved child iceberg, according to specified values for the Young's modulus (E) and yield strength (σ y ) of ice. Growth of the submerged foot is tracked according to empirical models for the melt and erosion of the iceberg freeboard, thereby ensuring physically based and realistic footloose calving rates that can respond to changes in climate. We test the parameterization within the GFDL OM4 sea-ice-ocean model using realistic calving-size distributions for each hemisphere. Without the footloose parameterization, large tabular icebergs drift to unrealistically low latitudes in the Southern Hemisphere. However, activating the footloose parameterization yields modeled iceberg areas and trajectories that closely match observations, which suggests that the footloose mechanism plays a major role in iceberg decay that was previously unaccounted for in climate models.
Our experiments reveal several fundamental aspects of iceberg dynamics. The trajectories of large parent icebergs determines where child footloose icebergs are released. However, for the range of footloose parameters (E and σ y ) tested here, the overall iceberg meltwater distribution is more sensitive to the variations in child iceberg size (as controlled by E) than variations in large parent iceberg trajectories, because smaller child icebergs more easily escape boundary currents. Furthermore, smaller child icebergs also allow parent icebergs to become smaller via footloose calving alone, that is, an increase in total mass loss due to footloose calving, because footloose is only active along parent dimensions that exceed three times the length of the potential child iceberg. This effect is particularly influential in determining how rapidly the smaller parent icebergs in the Northern Hemisphere decay entirely.