The iKID model treats each particle, or element, as a vertical column of ice that experiences drag from the ocean, atmosphere, sea ice, and seafloor (if grounded); normal forces, shear, and torques between elements ( 20 ); a force due to sea surface slope; a wave radiation force; and the Coriolis force. When the tensile stress on a bond exceeds the tensile ice strength, the bond breaks ( 21 ). While each element has a horizontal area and thickness, the iKID model is two-dimensional in the sense that there is only one vertical layer of elements and only horizontal forces are represented. Therefore, bending effects associated with grounding or changes in dynamic ocean topography over the length of the iceberg are neglected.
When the A68a experiment was run in serial on an Intel Xeon CPU ES-2697 v4, the wall clock time for the MTS scheme averaged about 10 s to simulate 1 day of iceberg evolution. Parallelized runs achieve similar wall clock times because all bonded particles comprising an iceberg conglomerate are transferred to all processing domains that the conglomerate overlaps (see Supplementary text) and computed redundantly. Therefore, the conglomerate with the greatest number of particles has a strong influence over how quickly a simulation will run. Nevertheless, iceberg A68a was the sixth largest iceberg on record ( 4 ), so we conclude that our bonded-particle model is computationally efficient enough for implementation within century-scale climate simulations. Further speed up may be possible through vectorization or by increasing the “short” MTS time step increments, which may be possible without sacrificing stability if particle size is increased or the Young’s modulus is decreased (see the “Tuning” section).
To increase computational efficiency, we developed a multiple time step (MTS) velocity Verlet scheme to integrate the equations of motion. In the MTS scheme, all forces are evaluated on a “long” time step increment (here, 30 min), except for the grounding drag and interactive forces between elements belonging to the same “conglomerate” of bonded elements, which are evaluated more frequently over a series of shorter substeps small enough to guarantee stability (here, 20 s). This scheme reduces how often each force must be evaluated, decreases the number of interpolations of gridded data to particles, and minimizes memory transfers between processing domains during parallel runs.
We forced the A68a simulation with OSCAR (Ocean Surface Current Analysis Real-time) near-surface ocean current velocities ( 22 ), SSALTO/DUACS (Segment Sol multimissions d’ALTimétrie, d’Orbitographie et de localisation précise/Data Unification Altimeter Combination System) sea surface heights ( 23 ), and NCEP/NCAR (National Centers for Environmental Prediction/National Center for Atmospheric Research) Reanalysis 1 10m vector winds ( 24 ). These fields were interpolated to the particles from a 1/8° background grid at the start of each half-hour time step. Guided by NASA Aqua MODIS (Moderate Resolution Imaging Spectroradiometer) imagery, we initialized the iceberg position within the OSCAR current at its observed longitude on 9 December 2020, and we assigned the iceberg an initial eastward velocity of 0.22 m/s. We arranged the bonded particles on a regular Cartesian lattice, i.e., square packing, where each particle has a maximum of four bonds, a constant radius of 1.5 km, and an estimated ice thickness of 200 m. Figure 2A shows this initial December 9 iceberg configuration, where the orange rectangle marks the grounding zone responsible for the first breakup event and is the only area where we activate the iceberg grounding drag. We manually delineated this zone to be slightly southwest of the observed grounding zone because the OSCAR ocean currents do not appear to flow close enough to the observed grounding zone.
Tuning
The primary tuning parameters that affect model behavior are the Young’s modulus (E), the horizontal (c o, h ) and vertical (c v, h ) ocean drag coefficients, the grounding drag coefficient (c g ), and the tensile bond strength (σ c ). A full description of these parameters is provided in the Supplementary text, and the values of all model parameters used for the A68a experiment are given in table S1. We caution that these values may not be applicable for all icebergs. Therefore, additional icebergs should be modeled in future studies to better constrain these values and estimate how they may vary between icebergs. We describe our tuning process here as a guide for future studies.
E. Pure, undamaged ice has a Young’s modulus between 1 and 10 GPa. However, there are both numerical and physical reasons to decrease E when modeling icebergs. Numerically, smaller values of E decrease the velocity of seismic waves, which increases computational efficiency by allowing longer time steps when evaluating bonded-particle interactive forces. Physically, E should decrease under the following conditions that may be applicable to icebergs: (i) when ice temperature increases and under strain rate effects over long loading times (E to 5 MPa, which is a large enough value to guarantee that iceberg behavior is visually stiff while allowing the bonded-particle scheme to be computationally efficient enough to use within climate models. We began the tuning process by determining an appropriate Young’s modulus,. Pure, undamaged ice has a Young’s modulus between 1 and 10 GPa. However, there are both numerical and physical reasons to decreasewhen modeling icebergs. Numerically, smaller values ofdecrease the velocity of seismic waves, which increases computational efficiency by allowing longer time steps when evaluating bonded-particle interactive forces. Physically,should decrease under the following conditions that may be applicable to icebergs: (i) when ice temperature increases and under strain rate effects over long loading times ( 25 ); (ii) when seawater or surface meltwater infiltrates into firn ( 26 27 ); and (iii) to account for crevassing, which decreases the ice thickness along which stresses are transmitted ( 28 ). Iceberg A68a exhibited substantial crevassing that was present when it was part of the Larsen C ice shelf ( 29 30 ). We setto 5 MPa, which is a large enough value to guarantee that iceberg behavior is visually stiff while allowing the bonded-particle scheme to be computationally efficient enough to use within climate models.
c o, h = 0.02136 and c o, v = 16.02. These coefficients differ from those typically used for unbonded, point-particle models of small icebergs (c o, h :c o, v that is used in point-particle iceberg models ( Next, we determined the ocean drag coefficients that yielded a modeled iceberg drift path that best matched observations:= 0.02136 and= 16.02. These coefficients differ from those typically used for unbonded, point-particle models of small icebergs ( 6 31 ), because here, each particle only constitutes a portion of a large bonded-particle conglomerate. Furthermore, the tuning may make up for error in the ocean current data or the fact that we force our model with ocean surface currents alone rather than currents that vary over the depth of the iceberg. When determining the optimal coefficients, we assumed for simplicity that the coefficients retained the 1:750 ratio forthat is used in point-particle iceberg models ( 6 31 ). Other values may yield a similar model response. Note that we do not pursue a similar tuning exercise for the wind drag, because wind contribution to the motion of large icebergs is small ( 32 ). Instead, we simply use the wind drag coefficients from a point-particle iceberg model ( 6 ). We do not tune the sea-ice drag coefficients because sea ice is absent in the vicinity of the iceberg in December 2020. We are able to attribute each breakup event to contact with either the seafloor or ocean currents because these are the primary processes that determine iceberg stresses and drift.
c g . Small values of c g will only slow the iceberg drift, while larger values resemble an impact event, which is more appropriate for A68a. We set c g = 104 kg m−2 s−1. Last, we determined the tensile ice strength, or tensile stress threshold for breaking bonds (σ c ), that gives the best match between modeled and observed iceberg breakup. We set σ c to 18 kPa, which is similar to the stress scale of order 10 kPa estimated in a previous study for rift calving of large tabular icebergs ( After finalizing the Young’s modulus and ocean drag coefficients, we tuned the grounding coefficient,. Small values ofwill only slow the iceberg drift, while larger values resemble an impact event, which is more appropriate for A68a. We set= 10kg m. Last, we determined the tensile ice strength, or tensile stress threshold for breaking bonds (σ), that gives the best match between modeled and observed iceberg breakup. We set σto 18 kPa, which is similar to the stress scale of order 10 kPa estimated in a previous study for rift calving of large tabular icebergs ( 33 ). This previous estimate neglected local stress amplification near rift tips, an assumption that is also likely applicable to the current study given the coarse resolution of the model, where each bond is ∼3 km long.
The above tuning exercise may not have arrived at the only combination of parameters, e.g., tuning of Young’s modulus and ocean drag coefficients, sufficient to model A68a, nor did we document the sensitivity to parameters, which could be nonlinear in response and thus make particular tuning nonunique or not robust. For example, it has yet to be shown that the tuning for an iceberg would be optimal as ice thickness changes over time. Note that tuning is not a trivial process; small changes to a single parameter can sometimes divert the drift trajectory of an iceberg drastically due to strong spatial variations in ocean and climate forcings. Despite tuning sensitivity, the existence of at least one set of parameters that appear to explain the A68a breakup suggests that we may be able to find parameters for other icebergs and events, and ultimately, one can imagine a general model for these parameters.