A Lagrangian theory of geostrophic adjustment for zonally invariant flows on a rotating spherical earth

We examine the late-time evolution of an inviscid zonally symmetric shallow-water flow on the surface of a rotating spherical earth. An arbitrary initial condition radiates inertia–gravity waves that disperse across the spherical surface. The simpler problem of a uniformly rotating (f-plane) shallow-water flow on the plane radiates these waves to infinity, leaving behind a nontrivial steady flow in geostrophic balance (in which the Coriolis acceleration balances the horizontal hydrostatic pressure gradient). This is called “geostrophic adjustment.” On a sphere, the waves cannot propagate to infinity, and the flow can never become steady due to energy conservation (at least in the absence of shocks). Nonetheless, when energy is conserved a form of adjustment still takes place, in a time-averaged sense, and this flow satisfies an extended form of geostrophic balance dependent only on the conserved mass and angular momentum distributions of fluid particles, just as in the planar case. This study employs a conservative numerical scheme based on a Lagrangian form of the rotating shallow-water equations to substantiate the applicability of these general considerations on an idealized aqua-planet for an initial “dam” along the equator in a motionless ocean. VC 2021 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http:// creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0054535


I. INTRODUCTION
Geostrophic adjustment theory describes the intricate manner in which an arbitrary initial rotating stratified flow transforms to a final, steady, and geostrophically balanced state. The paramount importance of this theory is clear in view of the dominance of geostrophic balance on scales larger than mesoscale (in the ocean) and synoptic scale (in the atmosphere). 1,2 The theory was first developed in the 1930s on the f-plane 3,4 and since then, several follow-up studies [5][6][7][8] have extended the original results to various initial conditions, including the classical fluid dynamical problem of a dam break in a fluid initially at rest. 7 Accounting for latitudinal variations in the Coriolis frequency (or background planetary vorticity), even in planar geometry through the b effect, greatly complicates the analysis and is the reason for the slow progress in advancing the geostrophic adjustment theory beyond the f-plane. In this study, we employ a Lagrangian formulation of the dynamics to study the adjustment process on a rotating spherical planet. We consider the simplest relevant model, namely, the one-dimensional, zonally symmetric shallow-water model that ignores longitudinal variations of both velocity and height. This approximation limits the applicability of the model to zonal flows.
As in all geostrophic adjustment theories, the analysis involves a derivation of the final, presumably steady, geostrophic or "balanced" state together with a description of the transient waves that transform the initial state to the final state. Since the surface of a sphere is a bounded domain, the transient waves cannot propagate to infinity and hence, in contrast to the classical studies mentioned above, the integral (over the entire sphere) of the initial state's energy must be conserved. The conserved energy must, therefore, be divided between the steady balanced state and the transient waves. As the energy in the balanced state is fixed, then so must be the energy in the waves. Thus, because the domain is bounded, the waves cannot escape (as in an unbounded plane), and in the absence of dissipation they may mask the underlying steady balanced flow for all time.
The Lagrangian approach adopted here is necessary to isolate the underlying steady balanced state having the same invariants as the initial state. That is, each fluid "particle" (here a latitude circle) in the balanced state must have the same mass and angular momentum as that in the initial state. This novel Lagrangian formulation allows one to derive the appropriate nonlinear balance equation (in a Lagrangian coordinate) determining the balanced state corresponding to any prescribed initial state. The Lagrangian formulation is also useful numerically for studying the transient flow evolution, as conservation of mass, angular momentum, and potential vorticity on each fluid particle can be exactly satisfied at all times. The Lagrangian approach is unique in that no high-order viscous terms are required to stabilize the numerical solution. This is of primary importance when studying the nonlinear rotating shallow-water equations in spherical coordinates that struggle with the coordinate singularity at the poles. 9,10 In more complex two-dimensional and layerwise-twodimensional flows, an (at least partially) Lagrangian formulation has proved highly beneficial in the study of fundamental aspects of geophysical fluid dynamics. The suite of "contour advection" methods [11][12][13][14][15] enforces material conservation of potential vorticity by tracking material contours in a Lagrangian manner, but otherwise use Eulerian (fixed-grid) methods for efficiency. These methods originate in the fully Lagrangian "contour dynamics" method, [16][17][18] first applied to two-dimensional flows in an unbounded domain. Contour dynamics requires an invertable linear relation between (potential) vorticity and the flow field (e.g., through a Green function). Contour advection relaxes this requirement by making use of Eulerian inversion methods, allowing access to a much wider range of fluid models. The Lagrangian formulation has been extensively employed in studies of inertial (particle) motion on a rotating sphere and the rotating spherical earth. 19,20 A comment is in order regarding our use of the term "spherical earth" in the title instead of "sphere." On a sphere that rotates at the Earth's rotation frequency, the centrifugal acceleration exceeds the Coriolis acceleration even for speeds of 100 ms À1 . However, the equatorial bulge that results from the slight eccentricity of the shape of the Earth's surface ($0:003) adds a poleward directed component of the gravitational acceleration that balances the equatorward directed centrifugal acceleration. Since the Earth's eccentricity is small, and since its effect on wave dynamics is second order, 21 the fractional error expected from the approximation of Earth's ellipsoidal surface by a sphere, while ignoring the centrifugal acceleration, is only Oð10 À5 Þ. Here, we consider this negligible.
Another important point is that we are not directly considering the problem of "balance," a hypothetical state (for a general twodimensional flow) in which inertia-gravity waves are absent. [22][23][24][25][26] Such a state is often constructed by "PV inversion," using the definition of potential vorticity alongside a pair of "balance relations" that effectively filter the inertia-gravity waves. In this study, the zonal symmetry means that there is an exact balanced state, a completely stationary state, entirely determined by the Lagrangian mass and angular momentum distributions. These distributions determine the potential vorticity distribution for zonal flows. The focus of this study is instead on the time-dependent problem, specifically the role played by inertia-gravity waves generated from the initial conditions while the slowly propagating Rossby waves are eliminated by the assumed zonal symmetry. Notably, all unsteadiness in zonal flows may be attributed to inertia-gravity waves, the imbalanced motions. Our objective is to examine how such flows adjust to the underlying steady balanced state, at least in a time-averaged sense. This paper is organized as follows. In Sec. II, the Lagrangian forms of the zonally symmetric spherical shallow-water equations are derived. From these, a single second-order nonlinear ordinary differential equation (ODE) is obtained for the balanced flow consistent with the initial state's distributions of mass and angular momentumthis appears to be new. The eigenvalue problem for the unsteady linear waves is also discussed, together with the numerical methods employed. The unsteady flow dynamics is described in detail in Sec. III for a specific but representative case. Here, it is shown that the time-averaged flow is nearly identical to the balanced flow, obtained independently. Balance is then the focus of Sec. IV, which explores the wider parameter space and discusses various flow properties. Our conclusions are offered in Sec. V.

II. FORMULATION OF THE LAGRANGIAN MODEL AND ITS NUMERICAL METHOD OF SOLUTION
The shallow-water equations possess a number of fundamental conservation properties due to underlying symmetries. Besides mass, angular momentum, and total energy (kinetic plus potential), potential vorticity (circulation) is locally conserved; that is, it remains constant on fluid particles (in the absence of forcing and dissipation). In a sense, mass is also locally conserved, but following infinitesimally small deformable fluid elements. For the zonally symmetric flows considered here, angular momentum also becomes locally conserved, that is, it does not change following each moving latitude.
It would therefore appear useful to cast the equations of motion in Lagrangian terms, to the extent possible (see below). In this way, conservation is built in directly (apart from energy, which is an integral invariant). This results in a pair of equations for each "particle" (really a latitude circle), one for the rate of change of its latitude and another for the rate of change of its meridional velocity. This is analogous to a particle in a force field, except here there are an infinite number of particles, and they interact locally through hydrostatic pressure variations.

A. Equations of motion
We next derive the Lagrangian forms of the zonally symmetric shallow-water equations from their standard Eulerian forms. These Lagrangian forms appear to be new and are especially useful for simulating the evolution of weakly to moderately nonlinear flows in a very nearly conservative manner. Their steady versions also enable one to link the underlying balanced state to any initial state.
To set notation, let u and v be the zonal and meridional velocity components and h be the local fluid height (here for simplicity above a flat bottom) divided by the mean fluid depth H. Let / stand for latitude and t for time. Let X be the background rotation rate, and consider the fluid motion in a frame of reference rotating at this rate. Let g be the acceleration due to gravity and c ¼ ffiffiffiffiffiffi gH p the short-scale gravity-wave speed. Without loss of generality, we take the radius of the planet to be 1.
The Eulerian forms of the equations read Physics of Fluids ARTICLE scitation.org/journal/phf where f ¼ 2X sin / is the Coriolis frequency (or planetary vorticity). It is convenient in what follows to define r ¼ cos / as the radial distance from the z-axis (connecting the two poles), and z ¼ sin /. We also use _ n ¼ @n=@t þ v@n=@/ to denote the material derivative or Lagrangian rate of change of any quantity n.
Equation (1) expresses conservation of angular momentum. Multiplying it by r and using _ / ¼ v, we find which shows that angular momentum U ¼ rðu þ XrÞ is materially conserved: _ U ¼ 0. This motivates introducing a Lagrangian label a for each fluid particle, and denoting the latitude of the particle by /ða; tÞ. The label itself could be the initial latitude, in which case /ða; 0Þ ¼ a. Material conservation of angular momentum is then equivalent to the statement that U depends only on a.
Equation (3) expresses conservation of mass: the mass dM between any two material latitudes /ða; tÞ and /ða þ da; tÞ is fixed for all time. As differential area on a spherical surface is proportional to rd/, it follows that dM hrd/ depends only on a. That is, Similarly, the angular momentum distribution U(a) determines the zonal velocity from uða; tÞ ¼ UðaÞ rða; tÞ À Xrða; tÞ: Potential vorticity conservation, Q ¼ QðaÞ, is a by-product of the conservation of M and U. From the definition, we have where f is the relative vorticity, by (7). It follows that by (5). This proves that Q is a function of a only, i.e., Q is materially conserved. The remaining Eq. (2) of the shallow-water set for the meridional velocity v, together with v ¼ _ /, provide the dynamical evolution equations in the Lagrangian formulation. Using (7) to replace u in (2), these equations are (The reader is reminded that r ¼ cos / and z ¼ sin /.) Only the hydrostatic acceleration Àc 2 @h=@/ is left in the Eulerian form. It can be converted using (6) for h and the chain rule @h=@/ ¼ ð1=/ 0 Þ@h=@a. This is a nonlocal term since it depends on / 0 and / 00 . It is through this term that all particles interact with one another. The boundary conditions are that / ¼ 6p=2 at the north and south poles, and v ¼ 0 there for consistency. Typically, the angular momentum U $ r 2 near each pole, so v remains zero at each pole only if @h=@/ ¼ 0 there.
Equations (11) and (12) conserve total energy (the Hamiltonian, an integral invariant), the sum of kinetic and (available) potential energy Àp=2 uða; tÞ 2 þ vða; tÞ 2 þ c 2 ðhða; tÞ À 1Þ À Á M 0 ðaÞ da: (13) This is actually the total energy divided by the mean depth H due to the scaling adopted for h. Numerically, Eqs. (11) and (12) are simple to solve as the only coupling term is the hydrostatic acceleration. We discretize the system into n þ 1 particles (latitudes / j ), indexed j ¼ 0; 1; …; n, equally spaced in the label a (reserving j ¼ 0 and n for the fixed polar values, / 0 ¼ Àp=2 and / n ¼ p=2). From a prescribed initial height distribution hð/; 0Þ, we calculate the fixed mass m j ¼ DM j between adjacent latitudes / j by integrating (5) at t ¼ 0 This is also required to hold for all subsequent times t > 0. Between latitudes / j ðtÞ and / jþ1 ðtÞ for j ¼ 0; 1; …; n À 1, the height h is represented by the quadratic spline with the coefficients h j , a j , and b j determined by: (a) enforcing continuity of h and @h=@/ at each internal latitude / j (j ¼ 1; 2; …; n À 1), (b) requiring (14) to hold in each interval (j ¼ 1; 2; …; n), and (c) requiring @h=@/ ¼ 0 at each pole. This gives 3n conditions to determine the 3n coefficients h j , a j , and b j (j ¼ 0; 1; …; n À 1). Note h n ¼ h nÀ1 þ a nÀ1 þ b nÀ1 . The resulting system can be put into a simple tri-diagonal form (see the Appendix) and is readily solved in O(n) operations. This provides @h=@/ at each / j , needed to calculate _ v j in (12). Besides the mass distribution, the initial meridional velocity vða; 0Þ and angular momentum distribution U(a) must be specified. The latter is calculated from the initial zonal velocity field uð/; 0Þ, taking /ða; 0Þ ¼ a for simplicity. Then UðaÞ ¼ rða; 0Þðuða; 0Þ þ Xrða; 0ÞÞ with rða; 0Þ ¼ cos ð/ða; 0ÞÞ ¼ cos a. The time integration is carried out using a fourth-order Runge-Kutta method, with a time step Dt taken sufficiently small to resolve the fastest gravity waves (ensuring numerical stability and excellent energy conservation).

B. Steady flow: Balance
The theory of geostrophic adjustment argues that any initial state will tend to a steady balanced flow having the same invariants M(a)

ARTICLE
scitation.org/journal/phf and U(a) as t ! 1. While this scenario is strictly prohibited on a bounded surface like a sphere, due to energy conservation, we will see below that the steady balanced flow is often (if not always) found upon time averaging over a sufficiently long period T. Indeed, if we time average (11), the LHS is bounded by p=T, and letting T ! 1 shows that the time-averaged meridional velocity v must be zero for all particles a in this limit. If we time average (12), however, we must assume that v(a, t) remains bounded pointwise to conclude that the LHS vanishes as T ! 1. We cannot guarantee this using energy conservation, but it is plausible. It is therefore of interest to compute the steady solutions of (11) and (12), for specified angular momentum and mass distributions U(a) and M(a). Then v ¼ 0, and using (6) to substitute for h(a) in (12), we obtain the nonlinear balance equation where r ¼ cos ð/ðaÞÞ and z ¼ sin ð/ðaÞÞ. This is a second-order nonlinear boundary value problem for /ðaÞ, with Dirichlet boundary conditions /ð6p=2Þ ¼ 6p=2. Notably, all steady zonal flows on a spherical earth satisfy (16). Once a solution /ðaÞ is obtained, the balanced height and zonal velocity distributions are found from (6) and (7), i.e., One can also calculate the relative vorticity f from (9). For M ¼ sin a and U ¼ X cos 2 a, (16) admits the solution / ¼ a. This corresponds to a flow at rest u ¼ 0 having a uniform height h ¼ 1. Otherwise, (16) is a highly nonlinear equation, and we must resort to a numerical solution. Here, (16) is solved on an equally spaced grid in a. Starting with an initial guess /ðaÞ ¼ a, subsequent guesses are obtained by linearizing (16) about the previous guess, approximating / 0 and / 00 by centered differences, and solving the resulting tri-diagonal problem for the perturbation D/. Underrelaxation is used by adding only half of D/ to the previous guess / to form the new guess (this is essential for convergence). When maxjD/j < 10 À12 , we accept /ðaÞ as the solution. Notably, this method even works under extreme conditions in which the fluid depth h nearly vanishes (incipient out-cropping).

C. Linear waves
Besides steady motions, Eqs. (11) and (12) admit propagating inertia-gravity waves, which cause particle oscillations (latitudinal displacements) and which fundamentally depend on both rotation and stratification (gravitational restoration of the free surface). A general initial condition can be thought of as a collection of these waves and an underlying steady balanced flow. Here we focus on the linear problem and consider only a basic state at rest. We shall see that the results are nonetheless relevant to the nonlinear problem.
A basic state at rest (u ¼ v ¼ 0, h ¼ 1) corresponds to the mass and angular momentum distributions M ¼ sin a and U ¼ X cos 2 a. Both the displacement dða; tÞ ¼ /ða; tÞ À a and the meridional velocity v(a, t) are assumed to be sufficiently small to be able to linearize (11) and (12). Taking dða; tÞ ¼dðaÞe Àixt and vða; tÞ ¼vðaÞe Àixt , where x is the frequency, (11) and (12) are readily combined into a single linear boundary value problem fordðaÞ with Dirichlet boundary conditions,dð6p=2Þ ¼ 0. Above, r ¼ cos a and f ¼ 2X sin a. This is a Sturm-Liouville problem with k ¼ x 2 =c 2 serving as the eigenvalue. In general, there are an infinite number of solutions, with 0 < k 1 < k 2 < … all real. Moreover, the k k depend only on the dimensionless parameter L d ¼ c=ð2XÞ, equal to the "Rossby deformation length" divided by the radius of the planet (here unity). Note the scaled frequency x=ð2XÞ ¼ ðx=cÞL d . Figure 1 shows how x=ð2XÞ depends on L d for the first 10 modes (k ¼ 1-10). For L d ) 1; x=ð2XÞ ! ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi kðk þ 1Þ p L d and the corresponding eigenmodesdðaÞ are associated Legendre polynomials P 1 k ðsin aÞ. In fact, this is a good approximation until kL d < 1. For small L d , the eigenmodes approach parabolic cylinder functions which are localized within a distance of OðL d Þ of the equator, while x=ð2XÞ ! ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð2k À 1ÞL d p as L d ! 0. These dispersion curves accurately reproduce those shown in Fig. 1 of Longuet-Higgins, 27 which were derived in an Eulerian formulation using standard methods. The excellent match between the results obtained by the two vastly different methods of solution attests to the applicability and accuracy of the Lagrangian formulation of the zonally symmetric rotating shallowwater equations. These results also agree with the explicit expression derived for these waves in the limit of large L d in Ref. 28.

III. DYNAMICS: THE DAM BREAK PROBLEM
We next examine the evolution of a flow initially at rest but with a difference in height h between the northern and southern hemisphere. Specifically, we consider the initial condition /ða; 0Þ ¼ a; vða; 0Þ ¼ 0 and hða; 0Þ ¼ 1 À Atanhðz=wÞ; where z ¼ sin a, w is the width of the step, and A is an amplitude parameter. Without loss of generality, it is sufficient to consider A > 0. Furthermore A < 1 is necessary or else the fluid would outcrop at some latitude. The corresponding mass distribution MðaÞ ¼ z ÀAw log ðcoshðz=wÞÞ. The angular momentum distribution is In what follows, we take X ¼ 2p without loss of generality. Then t ¼ 1 corresponds to one "day." Figure 2 depicts the flow evolution for h À 1, u and v (left to right) over the first day, for a dam break with A ¼ 0.05, w ¼ 0.1, and L d ¼ 0:2. [Here n ¼ 500 intervals in a were used, but the results are closely reproducible at half and double resolution, at least up to early times (t ¼ 10); thereafter solutions gradually diverge.] On the left, the initial step in hÀ1 is seen to split into two steps of nearly equal size (but opposite signs) propagating toward the poles. The one propagating toward the south pole is faster as the gravity-wave speed at the leading edge is approximately c ffiffi ffi h p and is larger where h is larger in the southern hemisphere. This leads to an asymmetry that grows in time. It results in more nonlinear behavior in the northern hemisphere, causing small-scale ripples in h À 1 and v. In the zonal velocity u, jets develop on either side of the equator but then collapse (only to re-emerge repeatedly at later times). The meridional velocity is initially positive everywhere until it reverses around t ¼ 0.5 then becomes Physics of Fluids ARTICLE scitation.org/journal/phf progressively more complex. At these early times, there is little sign of an underlying steady balanced flow. The flow remains unsteady with no significant change in v rms (in a running average) until the end of the simulation at t ¼ 250. Over this time, the total energy H defined in (13) remains conserved to within 0.8%, while the kinetic and potential components vary by 92.7% of H (in such a way as to nearly conserve total energy). Some of this energy loss is due to the formation of under-resolved small-scale gravity waves by nonlinearity. No explicit numerical damping or de-aliasing is included, but evidently the interpolation used for h results in weak damping under these circumstances (a smaller time step makes no difference).

ARTICLE
scitation.org/journal/phf Figure 3 shows the results of time-averaging the latitudinal displacement dða; tÞ /ða; tÞ À a, the height anomaly hða; tÞ À 1, and the zonal velocity u(a, t) vs a, in the same simulation depicted in Fig. 2 at early times. The solid curves show the time-averaged profiles and the shaded regions show the variation (61 standard deviation) about the average. First of all, this variation is large, indicating that the inertia-gravity waves play a major if not dominant role in the dynamics. The flow is never close to balanced. Notably, the time-averaged height distribution hhi is flat at the equator, where initially it had its steepest gradients. The corresponding zonal velocity distribution hui has a broad eastward jet in the northern hemisphere, and a broad westward jet in the southern hemisphere. The jets appear antisymmetric, but there is a slight asymmetry (see below).
We next compare the time-averaged profiles with the steady balanced flow obtained by solving the nonlinear boundary value problem (16) for /ðaÞ, using the same distributions U(a) and M(a) as in the simulation above. The balanced height and zonal velocity distributions are diagnosed from /ðaÞ using (17). Figure 4 shows the results. The agreement is excellent, with only a slight mismatch of the height distribution near the equator (this discrepancy is likely due to the finite resolution used and the finite period of time averaging). Remarkably, the energy H of the balanced flow (0.006 001) is just 42.2% of the total (initial) energy (0.014 212) in the dam-break simulation, confirming that the inertia-gravity waves dominate energetically, as already suggested by the large variability in Fig. 3.
To verify that the unsteady motions are purely inertia-gravity waves, the frequency power spectrum of v(a, t) was computed for each a over the period 0 < t 250, then averaged over a (after weighting the spectrum by cos a) to form P v ðxÞ, where x is the frequency. Figure 5 shows P v over an intermediate range of frequencies after log -log scaling. The prominent peaks starting around log 10 x ¼ 0:764 (or x ¼ 5:813) coincide almost exactly with the frequencies x k obtained from the linear analysis (18) for L d ¼ 0:2. This is despite the fact that the linear analysis assumes a basic state of rest-it would be more accurate to use the underlying balanced flow. The magenta dashed lines mark the frequencies of the odd modes x 1 , x 3 , etc., while the green lines mark those for the even modes. For the odd modes, the height anomalyĥ of the eigenmode is also odd in a, i.e., antisymmetric in latitude, while for the even modesĥ is even, see Fig. 6. Moreover, the kth mode exhibits k zeros in the eigenmode h À 1 as a function of a.
The most prominent frequency in Fig. 5 corresponds to x 1 and a mode with one zero in h À 1, just like the initial condition hða; 0Þ À 1 in (19). The next most prominent frequency is x 3 , also an odd mode. In fact all odd modes dominate in P v ; the leading even mode k ¼ 2 has less power than any of the odd modes shown. The presence of

ARTICLE
scitation.org/journal/phf even modes is partly due to the fact that the underlying balanced flow (the proper basic state) has a height anomaly that is not precisely antisymmetric in a, but it is also partly due to nonlinear wave-wave interactions. Those nonlinear interactions also give rise to additional peaks (sums and differences of two or more frequencies) seen between the marked frequencies. A weaker dam break [reducing A by 10 in (19)] removes nearly all of the secondary peaks and weakens those corresponding to the even modes (not shown).
In summary, a dam break evolves into a set of weakly interacting inertia-gravity waves riding on top of a steady balanced flow. The wave-wave interactions excite new waves at smaller scales, and this is expected to continue indefinitely, leading to a weak forward cascade of wave energy at small scales. This process eventually leads to numerical dissipation of energy, but viscosity would ultimately halt the energy cascade. This evolution has been seen across parameter space (varying A, w, and L d ), so long as A is not too large. When A exceeds 0.2 approximately, energy becomes poorly conserved and a strong forward wave cascade is seen, with many grid-scale variations. This may be an indication of shock formation, and indeed the local Froude number jvj= c ffiffi ffi h p À Á can exceed unity. Shocks however violate the longwave hydrostatic assumption underpinning the shallow-water model. A more accurate model capturing wave dispersion at small scales, such as the nonhydrostatic Green-Naghdi model, is required to go further. 29,30 The lack of dispersion at small scales in the shallow-water model is responsible for shock formation. 31 We hope to report on the dispersive model extension in future work.

IV. BALANCE: EXPLORATION OF PARAMETER SPACE
We next explore the diverse forms of the balanced flows having the same distributions of mass M(a) and angular momentum U(a) as the dam-break initial condition (19). Here, we take the width of the step w ¼ 0.01, much smaller than in the example presented in the previous section, to better approximate a step discontinuity in hða; 0Þ. Recall MðaÞ ¼ z À Aw log ðcoshðz=wÞÞ and UðaÞ ¼ Xr 2 , where r ¼ cos a and z ¼ sin a. For fixed w, the dimensionless balanced flow h and u=X depends only on the amplitude A of the step and the Rossby deformation length L d (the latter is already dimensionless since we have taken the radius of the planet to be 1). Figure 7 shows balanced flow profiles of a number of quantities for L d ¼ 0:2, and for a range of A. The top row shows displacement dðaÞ ¼ /ðaÞ À a, height h(a) and scaled zonal velocity uðaÞ=X plotted against the Lagrangian coordinate a, while the bottom row shows the scaled relative vorticity ffiffiffiffiffiffiffiffiffiffiffi ffi 1 À A p f=ð2XÞ, h(a) and uðaÞ=X plotted against actual latitude /ðaÞ. The sharp change in slope of the  v(a, t) averaged over a (with area weight cos a) for the simulation examined at early times in Fig. 2. Here, the spectrum was accumulated over 0 < t 250 using 10 6 time samples. The vertical dashed lines indicate the frequencies x k of linear waves on a basic state of rest, with the odd modes in magenta and the even modes in green. Note, P v continues to decay toward lower frequencies x, and also generally diminishes (with weaker peaks) toward higher frequencies.
FIG. 6. First three eigenmodes, displayed for the displacementd, height anomalyĥ ¼ Àr À1 ðrdÞ 0 , and scaled zonal velocityû=X ¼ 2zd, for a state of rest and for L d ¼ 0:2. The Lagrangian label, a, is plotted on the ordinate.  displacement dðaÞ ¼ /ðaÞ À a is associated with the small step width w. This causes kinks in u and in f, but not in h(a) when plotted as a function of actual latitude /ðaÞ. All flow profiles become increasingly asymmetric as A increases. For A ¼ 0.999, we see the most extreme situation in which h appears to vanish at some intermediate latitude in the northern hemisphere. In fact, h is very small but nonzero above this latitude. Only when A ¼ 1 exactly does the flow "outcrop": then there is no fluid north of the outcrop line (the limit A ! 1 is difficult to analyze; no results are presented here). As shown below, the proportion of energy in the balanced flow (as a fraction of that in the initial state) depends only weakly on A across the entire range of A. The relative vorticity f appears to diverge in the limit A ! 1 while u abruptly falls to zero at the outcrop line (consistent with divergent f). While this outcropping case is mathematically interesting, the local Froude number max juj= c ffiffi ffi h p À Á ¼ 5:2083 exceeds unity, and it is therefore likely that any disturbance to this steady flow will rapidly form shocks, violating the long-wave assumption underpinning the shallow-water model used here. 30 We next consider both smaller and larger Rossby deformation lengths L d . Recall that L d controls the "elasticity" of the free surface: when L d is large displacements in h tend to become small and broad scale; they also propagate fast (c ¼ 2XL d ). The opposite is true for small L d . Figure 8 shows balanced flow profiles for L d ¼ 1, five times larger than considered in Fig. 7. As expected, the profiles extend over a broader latitude range. They are also less asymmetric for smaller amplitudes A. The near outcropping latitude moves poleward (at larger L d outcropping does not occur), zonal speeds increase (exceeding the maximum rotational speed of the planet X for A >0.75 approximately), and the relative vorticity becomes larger. Moreover, the height anomaly h À 1 shows a strongly nonlinear dependence on A compared to the case with L d ¼ 0:2 in Fig. 7. All profiles migrate north when plotted against actual latitude /ðaÞ; hemispheric asymmetry grows with L d .

Physics of Fluids
Considering now smaller L d , Fig. 9 shows balanced flow profiles for L d ¼ 0:04, five times smaller than considered in Fig. 7. Now the near outcropping latitude moves equatorward and the flow variations are much more confined near the equator. The height anomaly h À 1 stays closely anti-symmetric until A is near its limiting value of 1. Zonal velocities are much weaker, but a bias toward a stronger (though narrower) northern hemisphere jet remains. This asymmetry is due to a shorter effective length scale c ffiffi ffi in the northern hemisphere where h < 1. It is not evident in u(a) vs a, but only in u(a) vs /ðaÞ as changes in / correspond to actual distance, whereas changes in a have no such meaning.
We conclude this section on balanced flows with a general summary of how flow properties like energy depend on L d and A, keeping w ¼ 0.01 fixed. For each balanced flow, we calculate the kinetic and potential energies from Àp=2 ðhða; tÞ À 1ÞÞM 0 ðaÞ da; (20b) whose sum gives the total energy H, see (13). Moreover, we diagnose the maximum (polar) Rossby number Ro max jfj=ð2XÞ and the maximum Froude number Fr max juj= c ffiffi ffi h p À Á .

Physics of Fluids
ARTICLE scitation.org/journal/phf The results are summarized in Fig. 10. In (a), H is seen to increase with A as expected, but also with L d until L d ¼ Oð1Þ. Thereafter, there is little dependence on L d . This is because the kinetic energy K dominates H for L d > 1 [see panel (b)], while the potential energy P diminishes as the free surface flattens and the flow becomes increasingly "barotropic." The maximum Rossby number Ro also increases with A and L d , and again shows little dependence on L d for L d > 1. The kinks in the contours in the middle of this panel are real: they are caused by max jfj occurring for f > 0 in one region and for f < 0 in another (see lower left panels in Figs. 7-9). The maximum Froude number Fr in panel (d) mainly increases with A but is also larger for small L d . Above the bold line indicating Fr ¼ 1 in this panel, the flow is subject to (spurious) shock formation if disturbed. This only occurs for A > 0.831 for small L d , and not at all for L d > 7:3 when the flow is strongly barotropic.
To better understand the importance of the transient waves in the dam break problem, we consider next the fraction of energy in the balanced flow as a function of L d and A. First of all, the initial dam break has only potential energy, and the total energy is readily computed from (19) and (13). Elementary integration gives In Fig. 11, we show the balanced state energy H, together with the potential P and kinetic K components, divided by H 0 . Only three selected values of A are displayed because there is almost no variation with A up to A ¼ 0.5 and remarkably little variation overall with A. At small L d , the balanced flow contains most of the energy of the dam break H 0 , and moreover this is mainly potential energy P. Thus, the waves contain only a small proportion of the total energy for small L d , and this proportion appears to reduce to zero as L d ! 0 (at least when w < L d ). For larger L d , even O(1) values, the situation is reversed: now the balanced flow contains little of the initial energy H 0 , and therefore the waves dominate, obscuring any underlying balance.

ARTICLE
scitation.org/journal/phf Another way of seeing this is to compare the balanced height h in Fig. 8 with the initial condition hða; 0Þ in (19), which is essentially a step function in a at the small value of w considered. For example, the A ¼ 0.2 curve (red) is much flatter than hða; 0Þ; the polar values of balanced h are just 1.0295 and 0.9595, compared to 1.2 and 0.8 for hða; 0Þ. This means that a large proportion of the initial (entirely potential) energy must be released into the waves, and this proportion only increases with L d . On the other hand, when L d ( 1 (see, e.g., Fig. 9 for L d ¼ 0:04), the balanced height h remains close to the initial condition except in a small region of OðL d Þ width centered on the equator. There is very little difference in potential energy between the balanced state and the initial state, and this difference goes to zero as L d ! 0 (at least when w < L d ). This is why so little adjustment takes place when L d ( 1: the balanced flow is already a close approximation to the initial one.

V. CONCLUSIONS
We have revisited Rossby's famous geostrophic adjustment problem 3,4 first solved in planar geometry under zonal symmetry. Geostrophic adjustment means that any localized initial condition will radiate inertia-gravity waves, all of which eventually escape to infinity, leaving behind a steady, zonal flow in geostrophic balance (a balance between the Coriolis acceleration and the horizontal hydrostatic pressure gradient). Moreover, this zonal flow is uniquely determined by the mass and momentum distributions of the initial fluid particles. These distributions are preserved as the fluid particles displace from the initial to the final, steady adjusted state.
In spherical geometry, the focus of this paper, there is an important difference: the domain is finite. Waves cannot escape, and cannot disappear without some form of damping (which is excluded from the geostrophic adjustment theory and would violate conservation of momentum). As a result, there is no possibility of a final steady flowindeed conservation of energy forbids it. Any underlying steady flow would have less energy than the initial flow, and the difference, which is entirely due to the presence of inertia-gravity waves, must be conserved. The initial flow already contains these waves; in time they spread through the domain and interact, but never change their energy content.
What then is "adjustment" on a spherical surface, or for that matter, on any finite surface? We argue that any confined flow, at least having zonal symmetry, can be decomposed into a steady "balanced" flow, and a residual unsteady flow consisting entirely of inertia-gravity waves. The balanced flow is the solution to the time-independent equations having the same mass and angular momentum distributions on fluid particles as the initial conditions. In general, this is determined from a highly nonlinear second-order ODE. Results of our timedependent simulations show that this balanced flow almost exactly corresponds to a long time average of the unsteady flow. We have gone further by comprehensively exploring the wide range of balanced flows permitted, these being determined by the solution of a novel highly nonlinear boundary value problem in a Lagrangian coordinate.
Though nonlinear phenomena such as shocks may develop in the rotating shallow-water equations 32 here the simulated fields after hundreds of days are accurately described by a combination of the balanced steady state and linear waves (the latter cannot leave the bounded domain). A caveat is that our simulations can only access the weakly to moderately nonlinear flow regime in which shocks do not form. The numerical method developed, which exactly enforces conservation of mass and angular momentum on fluid particles, is not appropriate when shocks form. On the other hand, shocks violate the long-wave assumption underpinning the shallow-water model. Arguably, nonhydrostatic dispersive effects should be taken into account at horizontal scales comparable to the fluid depth. [29][30][31] This is left to future work.
The idealization to zonally symmetric flows in this paper is of course not realistic. Two-dimensional flows also go through a form of adjustment, but the underlying balanced state is neither exact nor steady in general. 25 In such flows, balance is still an important concept, and is often used to explain general features of the large-scale atmospheric and oceanic circulation. 33 In effect, inertia-gravity waves are constantly produced in unsteady flows by spontaneous-adjustmentemission; 34 yet this is countered by adjustment toward an (unsteady) balanced flow. In contrast to the one-dimensional problem studied here in which Rossby waves are filtered out, these waves prevail in a two-dimensional problem. This is left for future work.

ACKNOWLEDGMENTS
The authors report no conflicts of interest.

APPENDIX: HEIGHT INTERPOLATION
Here we provide details of the method to interpolate the height field hð/; tÞ between discrete latitudes / j t ð Þ; j ¼ 0; 1; :::; n. Note / 0 ¼ Àp=2, while / n ¼ p=2 are the fixed polar values. In what follows, we suppress the dependence on t.
The numerical method used in Sec. III starts with the fixed mass m j between latitudes, given by the integral FIG. 11. Fraction of energy in the balanced flow relative to that in the initial dam break (19), as a function of the Rossby deformation length L d , and for three selected values of the amplitude A. Results are shown for the total, potential and kinetic energy divided by the total energy H 0 of the dam break initial condition (see the text).