*2.4. Mixing*

For the purposes of computing horizontal and vertical mixing, the domain is divided into isopycnal layers and columns respectively [24,25]. Columns of points are treated like columns of points in a finite difference model following [24]. Within isopycnal layers, parcels are allowed to mix momentum with their nearest neighbors following [25]. The Gent Mcwilliams [28] parameterization of mixing by eddies is not explicitly included, but Haertel and Fedorov [15] found that random motions of parcels generate a similar kind of mixing in this kind of Lagrangian model. The vertical viscosity is set to 10 cm<sup>2</sup> s<sup>−</sup>1, and the horizontal viscosity is set to 3 × 10<sup>4</sup> m<sup>2</sup> s<sup>−</sup>1. For the simulation with interior mixing, the vertical tracer diffusivity is 1 cm<sup>2</sup> s<sup>−</sup><sup>1</sup> (it is zero in the simulation without interior mixing). The horizontal tracer diffusivity is set to zero for both simulations.

#### *2.5. Bottom Topography and Spherical Geometry*

The Lagrangian method requires a specification of the bottom surface elevation over a horizontal grid of points on which Riemann sums are evaluated to calculate the pressure force [24]. For the simulations presented in this paper, the surface elevation field is constructed in the following way. First, actual mean surface height of the earth's land surface and ocean bottom is averaged into 1-degree bins. Then, ocean depths are truncated at 4900 m, positive land surface values are set to 300 m, and the land mass of Central America is widened slightly. The resulting data are smoothed repeatedly using a 1-2-1 filter. Finally, the smoothed 1-degree data are averaged over the spherical grid of the ocean model. A comparison of the unsmoothed one-degree bathymetry (dashed) and the GLOM bottom topography (solid) is shown in Figure 3a for 30 W. Sharp topographic features are smoothed to be consistent with the large parcels used in this study, and the slope at the land ocean interface is reduced, which helps to improve numerical accuracy in low resolution simulations [25]. The land surface elevation is set to 300 m to prevent the smoothing of the bottom surface as well as the enhanced free surface height perturbations resulting from the use of gravity wave retardation [29] from overly distorting the locations of continental boundaries (e.g., from opening a water pathway through Central America).

The GLOM uses a spherical geometry in which grid boxes have a constant meridional width (in both degrees and meters), and a zonal width that is roughly the same in meters but increases in degrees longitude with distance from the equator. This geometry was first developed for a global Lagrangian atmospheric model by [4]. A mercator projection of the grid boxes is shown in Figure 3b for North America.

**Figure 3.** Bottom topography and spherical geometry. (**a**) Model bottom surface along 30 W (solid black line), the observed bathymetry (dashed line), and the undisturbed free surface height (blue). (**b**) Spherical geometry for North America. Each grid box used for calculating the pressure force is outlined in black. (**c**) Regions of land (gray) and ocean (blue) as determined by the zero height contour. The solid black line denotes the actual position of the land–sea interface.

Here, grid boxes are shaded gray for z > 0 m, and blue for z < 0 m. A similar shading of model bathymetry is shown for the entire world in Figure 3c, with the z = 0 contour for actual bathymetry delineated with a black line. Notice that small peninsulas and islands are largely smoothed out in the GLOM's bottom topography, but that the gross structure of continents and ocean basins is retained.

For the simulations discussed in this paper, grid boxes are 1 by 1 degrees wide at the equator, but have a larger width in degrees longitude at higher latitudes. Note that this grid is used to evaluate the pressure force and to create plots of layer thickness. In contrast, parcels span multiple grid boxes and move freely throughout the model domain. While it is difficult to precisely characterize the equivalent Eulerian resolution of a Lagrangian model, the author estimates it to be roughly 2 or 3 degrees at low latitudes, and somewhat coarser at high latitudes for the simulations presented in this paper.

#### *2.6. Merging and Dividing Parcels*

Following [25], the water in the oceans is divided into a collection of equally sized water mass elements (WMEs), which may be considered to be the building blocks of water parcels. The target size for parcels is one WME in the upper layer (0–700 m), two WMEs in the middle layer (700–2100 m), and four WMEs in the lower layer (depths greater than 2100 m). When a vertically thick parcel rises to

a higher layer, and it is larger than the target size, it is sliced in half and the resulting two parcels begin moving independently (Figure 4).

**Figure 4.** Dividing and merging parcels. (**a**) A single parcel comprising two water mass elements. (**b**,**<sup>c</sup>**) The parcel divides into two smaller parcels each comprising one water mass element, which then move independently. Viewing the panels in reverse order illustrates the merging process.

When a vertically thin parcel sinks to a lower layer, it maintains its size until it overlaps and is sufficiently close to another vertically thin parcel, at which time the two parcels are combined (Figure 4, view panels from the bottom up). Notice that each individual WME has its own identification number, and it can be tracked through the dividing and merging process. While merging generates a small amount of numerical mixing, Haertel and Fedorov [15] found little difference between simulations with dividing and merging when compared to simulations with vertically thin parcels everywhere. Presumably, this is because merging only occurs when parcels go through substantial vertical displacements (i.e., it is a rare event). While dividing and merging can be turned off to completely remove numerical mixing of tracers, the author selected to use it here because it dramatically speeds up the model (by a factor of 2 or more for the simulations presented in this paper). Owing to the dividing and merging process, the vast majority of parcels in the upper layer comprise a single water mass element and have a vertical thickness of approximately 78 m, with a factor of 2 (4) increase in vertical thickness in the middle (lower) layer. All parcels have a radius of 333 km in each horizontal dimension (i.e., 3 degrees latitude and 3 degrees longitude near the equator). There are approximately 154,000 water mass elements in each simulation, which allows the model to be run for millennial time scales on a single processor (the GLOM has not ye<sup>t</sup> been coded to run in parallel) .
