**1. Introduction**

Many aspects of geophysical fluid dynamics are most easily modeled in a frame of reference that moves along with the fluid. For example, we show below that in such a framework three relatively simple ordinary differential equations predict large-scale fluid motions and local changes to fluid tracers (Section 2.3). While *partially* Lagrangian models, in which fluid particles are tracked in flow fields simulated by Eulerian methods, are becoming increasingly popular for both the atmosphere and the oceans [1–3], few *fully* Lagrangian general circulation models have been used for either the oceans or the atmosphere. This paper discusses the development of a fully Lagrangian model for simulating global ocean circulations.

One particular application the author has in mind for this model, is to couple it to his fully Lagrangian atmospheric model (LAM) [4]. The resulting Lagrangian coupled model is expected to be useful for weather and climate prediction as well as climate dynamics experiments. The LAM has already had success in simulating the Madden Julian Oscillation (MJO) [4–6], a planetary scale tropical weather disturbance [7] that is poorly represented in many atmospheric models [8–10], and which has global impacts on weather and climate [11–13]. The MJO also represents a component of the global atmospheric circulation that is predictable at much longer time scales than mid-latitude baroclinic waves [14]. The LAM has a unique spherical geometry [4], which as we note below is shared by the new global Lagrangian ocean model (GLOM), which makes the GLOM well suited for coupling to the LAM. However, there are many other potential advantages and applications for a fully Lagrangian global ocean model.

One such advantage is that the GLOM can circulate water without producing numerical tracer diffusion. A tracer value for a given water parcel does not change unless a tracer source is included, or a parameterization of tracer mixing is applied [15]. This contrasts with the behavior of z-coordinate ocean models, which generate an uncertain amount of spurious numerical mixing associated with the advection of tracers [16]. While isopycnal models [17] remove spurious numerical mixing of tracers *across* isopycnal surfaces—essentially by using a Lagrangian vertical coordinate—they also generate numerical mixing of tracers in modeling flow *along* isopycnals.

Another advantage of Lagrangian modeling relevant to this study relates to the parameterization of convection. In nature, when a fluid is unstable, convective plumes form that transport dense fluid downward and less dense fluid upward. The horizontal scale of the convective circulations is typically much smaller than that of the grid spacing for a large-scale model. What geophysical convective plumes ultimately accomplish is transporting fluid to a different layer. For example, in the atmosphere, tropical deep convection moves air from the boundary layer to the upper troposphere [18], and mid-level air to the boundary layer [19]. In the ocean, deep convective plumes move dense water from near the surface to the ocean bottom. As noted in Section 2.2, the Lagrangian convective parameterization produces the same vertical transport—by changing the vertical positioning of parcels—without attempting to model the details of the small-scale circulation in the convective plumes. The author believes that this unique convective parameterization has helped the LAM to simulate MJOs with realistic vertical structures and life cycles [4–6].

A third advantage of the Lagrangian approach to weather and climate modeling is that it provides fluid trajectory information for every mass element in the oceans and the atmosphere. Each component of fluid mass has an identification number that does not change during the course of a simulation [4,15]. The modeler can look up the position of a given mass element at all previous times for which model data is saved with no additional computations. It is hard to overemphasize the utility of this information. For example, it is well known that in the ocean, water mass characteristics are intimately connected to the locations in which they form, which is why water masses are given names like "North Atlantic Deep Water", "Antarctic Bottom Water", and "Antarctic Intermediate Water". The Lagrangian ocean modeler can easily determine where each water parcel last had contact with the surface (e.g., see Section 3.5). Similarly, in the atmosphere, the air temperature and moisture of a given air parcel are strongly dependent on previous locations of the parcel, with terms like "Artic Air" or "Gulf Moisture" (indicated moisture coming from the Gulf of Mexico) frequently used by meteorologists, and transports of moisture largely determining where moist convective systems form and move [4].

Finally, we note that perhaps the most important benefit of developing fully Lagrangian weather and climate modeling components, is that they will increase the genetic diversity of models. By perturbating the planet's radiative forcing, humans are essentially using the earth as a laboratory for a giant climate experiment with unknown consequences. We need a diversity of tools to map out the possible impacts of this experiment. Much in the same way that different kinds of engines (e.g., standard combustion, diesel, electric motors) have different niches for which they are best suited, so it is for different kinds of models of the atmosphere and oceans. The models and parameterizations that are the most fundamentally different from the others make the largest contribution to diversity. The Lagrangian tools presented here fit into that category. Much in the same way that the unique animals living in Australia are a grea<sup>t</sup> treasure to the biologist, Lagrangian modeling components are unique, and they have the potential to provide a fresh and fundamentally different perspective on ocean and atmosphere modeling and dynamics.

*Climate* **2019**, *7*, 41

The work presented in this paper builds on a number of previous studies involving Lagrangian models of lakes and oceans. Haertel and Randall [20] developed a Lagrangian numerical method for the oceans in which a body of water was represented as a collection of conforming water parcels referred to as "slippery sacks". This method was closely related to two Lagrangian numerical methods used at that time: partical in cell (PIC; [21]) and smoothed particle hydrodynamics (SPH; [22]). However, in geophysical applications, PIC had only been used for one or two layers of fluid (e.g., [23]), and the slippery sacks method differed from SPH in that two conforming water parcels could not occupy the same physical space whereas two SPH particles could. Consequently, as slippery sacks moved around, they maintained a fairly uniform distribution over time, whereas regions highly concentrated with particles or nearly void of particles could develop under SPH. Haertel et al. [24] adapted horizontal and vertical mixing schemes to the slippery sacks method and created a three-dimensional model of a large lake. Haertel et al. [25] developed an idealized Lagrangian model of the North Atlantic Ocean. Haertel and Fedorov [15] added an Antarctic circumpolar to the idealized Atlantic Ocean Model. Applications of these models have included large lake upwelling [24], simulation of Stommel and Munk solutions [25], modeling of the circulation and stratification in the North Atlantic Ocean [15,25,26], and simulations of idealized equatorial oceans and tropical instability waves [27].

This study advances our previous Lagrangian ocean modeling efforts in several ways. In particular, for the first time global spherical geometry and realistic bottom topography are included in an ocean model. This means that for the first time a fully Lagrangian ocean general circulation model is presented that can be used in global weather and climate modeling and climate dynamics experiments. This model (the GLOM) was developed by combining components of the Lagrangian basin-scale ocean model used by Haertel and Fedorov [15] and the Lagrangian atmospheric model developed by Haertel et al. [4] This paper presents the first simulations conducted with the GLOM, which not only provide evidence of its usefulness as a climate modeling tool, but also help to illustrate the role of interior tracer mixing in maintaining global ocean circulation and stratification. This paper is organized as follows. Section 2 describes the components of the GLOM. Section 3 presents our first fully Lagrangian simulations of global oceans. Section 4 discusses the results presented in light of other studies.

#### **2. Materials and Methods**

The primary data source for this study is the Global Lagrangian Ocean Model (GLOM). It was developed by combining components of several previous Lagrangian ocean and atmosphere models. Each of these model components is briefly described in this section with the goal of providing the reader with an intuitive understanding of how they work. More complete technical details are available from the original references that discuss the development of particular model components.

## *2.1. Fluid Parcels*

The GLOM represents a body of water as a collection of flexible water parcels [20]. Each parcel has a horizontal mass distribution function that remains fixed in the parcel's frame of reference (Figure 1).

**Figure 1.** Parcel vertical thickness function (i.e., the outline of a parcel on a level surface).

The mass distribution function defines the amount of mass per unit of horizontal area between the top and bottom surface of the parcel. Since density variations in the ocean are small, this function also provides a good approximation of the shape of the vertical thickness of a parcel, and it is also referred to as the "vertical thickness" function. Here, we use a third-order polynomial in each horizontal dimension (latitude and longitude) to define the mass distribution function following [24]. Constructing the mass distribution function in this manner is computationally efficient, and it helps with conserving energy when calculating the pressure force, which also guarantees numerical stability [20]. Although this means that the horizontal projection of a parcel is a square, it turns out that most of the contours of vertical thickness are circular, as is illustrated in Figure 2b of [27].

Each column of mass within a given parcel can move up or down independently of neighboring columns as the parcel slides over variable bottom topography or other parcels. So, although all parcels have the same horizontal mass distribution or vertical thickness function (Figure 1), parcels have a variety of shapes owing to vertical shearing (e.g., Figure 2a). As parcels slide over one another, they conform, so there are no gaps between parcels. Since each parcel has a fixed amount of mass associated with it, parcel centers maintain an approximately uniform distribution. When two parcels meet, the parcel with a greater density slides under the parcel with a lesser density, so the fluid maintains a neutral or stable stratification (Figure 2a; note that darker shading denotes denser parcels).

**Figure 2.** A schematic illustration of the Lagrangian convective parameterization. (**a**) A collection of water parcels, each having the same vertical thickness function, with darker shading denoting denser parcels. (**b**) The density of the top middle parcel (P) increases (e.g., from a temperature decrease or a salinity increase) yielding an unstable stratification. (**c**) The convective parameterization changes the stacking order of the parcels, so parcel P sinks to the level at which its density is the same as neighboring parcels. The long dashed line shows the path of the parcel center. Surrounding parcels rise slightly with short dashed lines illustrating the paths of parcel centers.

#### *2.2. Lagrangian Convective Parameterization*

On occasion, cooling or an increase in salinity can cause a near surface parcel to become more dense than parcels beneath it (Figure 2b). The Lagrangian convective parameterization then changes the stacking order of the parcels, so a neutral or stable stratification is restored (Figure 2c). Note that as the dense parcel sinks, neighboring parcels rise around it (Figure 2c; dashes lines denote vertical displacements of individual parcel centers). In this study, there is no mixing between convective updrafts and downdrafts (i.e., rising and sinking parcels). In atmospheric applications, we have found that carefully setting convective mixing (i.e., entrainment) is a key to accurately simulating moist convective systems such as the Madden Julian Oscillation [4–6].

Note that Figure 2 is a schematic illustration of the Lagrangian convective scheme, and it is intended to provide the reader with a mental picture, or intuitive understanding of how the parameterization works. For simplicity, the parcels have been drawn in two dimensions in an extremely coarse-resolution box-shaped ocean. A depiction of parcel interfaces for the simulations in this paper would include many more lines, a varying free surface elevation, and irregular bottom topography. Moreover, the GLOM does not actually keep track of parcel interfaces, but only the positions of parcel centers. The horizontal position vector for each parcel is a prognostic model variable (see below), whereas the vertical position of a given parcel depends on the vertical thicknesses of parcels beneath it. Near the beginning of each model time step, the vertical positions are diagnosed one parcel at a time starting with the lowest (most dense) parcel. The way that the convective scheme is implemented in the model is parcels are sorted by density after each time step, so density is a non-increasing function of the parcel array index, which determines the stacking order.

There are several limitations to this convective scheme. First, it is primarily intended to represent the net transports by deep, buoyancy-driven convection. The depth to which a parcel sinks is determined by its density relative to that of the environment (Figure 2). It is not intended to represent near-surface wind driven convection due to turbulent mixing, which can determine the depth of the surface boundary layer. Moreover, this scheme does not predict convective vertical velocities or accelerations of parcels, but rather the integrated effects of transports by updrafts and downdrafts. While mixing between rising and sinking parcels (i.e., entrainment) has been implemented and tested in atmospheric applications, it has not ye<sup>t</sup> been tested in ocean simulations.

## *2.3. Model Equations*

The following equations describe the motions of a fluid parcel:

$$\frac{d\mathbf{x}}{dt} = \mathbf{u} \tag{1}$$

$$\frac{d\mathbf{u}}{dt} + f\mathbf{k} \times \mathbf{u} = \mathbf{A}\mathbf{p} + \mathbf{A}\mathbf{m} \tag{2}$$

where **x** is horizontal position, **u** is horizontal velocity, *f* is the Coriolis parameter, **Ap** is the acceleration due to pressure, and **Am** is the acceleration resulting from the mixing of momentum (i.e., viscosity). Evaluating the pressure acceleration involves approximating an integral of pressure over the surface of a parcel with a Riemann sum. Haertel et al. [24] developed a method to efficiently evaluate this sum for every parcel, and for technical details the reader is referred to that study. The pressure force is then divided by the parcel mass to yield the pressure acceleration vector. Vertical motion occurs when a parcel slides up or down variable topography, and/or when there is mass convergence or divergence beneath a parcel (i.e., as in conventional hydrostatic fluid models).

A fluid tracer *q* changes only in response to sources and mixing:

$$\frac{dq}{dt} = Sq + Mq\tag{3}$$

where *Sq* is the source of the tracer and *Mq* denotes mixing. In this study, the only tracer we consider is temperature, because we use the following simple equation of state:

$$
\rho = 1029 \text{ kg } \text{m}^{-3} \text{ } (1 - 0.0002 \text{ } T) \tag{4}
$$

where *ρ* is density, and *T* is temperature. The only temperature source for parcels is in the surface layer, where a restoring temperature is applied (see Section 2.7). The author chose not to include salinity's effects on density for several reasons. First, he judged the transition from basin-scale ocean simulations with idealized topography to global-scale simulations with realistic topography, along with the required rewrite of much of the model code, to be a task of sufficient complexity in and of itself without using a new equation of state. Second, a major goal of this paper is to extend the results of [15] to global scales and they did not use an equation of state with salinity. Third, the author's next planned application for the GLOM involves studying the impacts of air–sea coupling on the dynamics of the Madden Julian Oscillation, which does not require the use of an equation of state with salinity variations. In nature, of course, salinity variations make important contributions to density, so in interpreting the simulations presented here and comparing to observed ocean structure, it is probably best to think of temperature as a proxy for density.
