## Abstract

The dynamics of current retroflection and rings shedding are not yet fully understood. In this paper, the authors develop an analytical model of the Agulhas Current retroflection dynamics using three simple laws: conservation of volume, momentum balance, and Bernoulli’s principle. This study shows that, for a retroflecting current with a small Rossby number, this theoretical model is in good agreement with numerical simulations of a reduced-gravity isopycnal model. Otherwise, the retroflection position becomes unstable and quickly propagates upstream, leaving a chain of eddies in its path. On the basis of these findings, the authors hypothesize that the westward protrusion of the Agulhas retroflection and the local “zonalization” of the Agulhas Current after it passes the Agulhas Bank are stable only for small Rossby numbers. Otherwise, the retroflection shifts toward the eastern slope of the Agulhas Bank, where its position stabilizes due to the slanted configuration of the slope. This study shows that this scenario is in good agreement with several high-resolution numerical models.

## 1. Introduction

### a. General background

Recent investigations of the role of the Agulhas system in the variability of the Atlantic meridional overturning circulation (AMOC) and the global climate show that the Agulhas leakage variability can impact the strength of the AMOC on several time scales (Weijer et al. 2002; Knorr and Lohmann 2003; Biastoch et al. 2008a; Beal et al. 2011). In turn, the Agulhas leakage itself strongly depends on the position of retroflection (van Sebille et al. 2009).

According to Doglioli et al. (2006) and van Sebille et al. (2010), 35%–45% of the Agulhas leakage is carried within rings (the remainder is direct leakage and leakage carried by filaments and cyclonic eddies). Observations indicate that the leakage flux into the South Atlantic is about 10–15 Sverdrups (Sv; 1 Sv ≡ 10^{6} m^{3} s^{−1}) (see Table 1.1 from van Sebille 2009; see also Gordon 1986; Gordon et al. 1987, 1992; Ganachaud and Wunsch 2000; Garzoli and Goni 2000; Boebel et al. 2003; Richardson 2007). Therefore, the outflow via anticyclonic eddies is about 4–6 Sv.

Typically, Agulhas rings are shed at a frequency of 5–6 yr^{−1}. As a rule, the retroflection protrudes westward before shedding a recurrent eddy and abruptly shifts eastward after shedding (Lutjeharms and van Ballegooyen 1988a; Dencausse et al. 2010a). According to Dencausse et al. (2010a), the position of the Agulhas retroflection typically moves from 15° to 17°E during the eddy shedding to 20° to 22°E just after the shedding.

There were also periods of almost 6 months when no shedding event was observed (e.g., Gordon et al. 1987; Byrne et al. 1995; Schonten et al. 2000; Lutjeharms 2006; van Aken et al. 2003; Dencausse et al. 2010a,b). This increased length of the shedding period may be associated with a retroflection farther to the east (de Ruijter et al. 2004). Lutjeharms and van Ballegooyen (1988b) and Lutjeharms (2006) showed anomalous and more occasional eastward shifts of the Agulhas retroflection occurring two to three times per year with durations of 3–6 weeks. Also, very irregular, so-called early (upstream) retroflection events were observed in 1986 (Shannon et al. 1990) and 2000/01 (Quartly and Srokosz 2002; de Ruijter et al. 2004), when the Agulhas Current retroflected east of the Agulhas Plateau. However, these events are uncommon.

Often, the Agulhas Current protrudes westward from the Cape of Good Hope. Such a necklike protrusion is connected with the formation of a local “zonalization” of the Agulhas influx after it passes the Agulhas Bank between 17.5° and 20°E, as shown in Fig. 1 [adapted from Beal et al. (2011, their Fig. 1) as a fragment], where the zonalization and retroflection areas are indicated by dashed ellipses [see also a similar configuration in Lutjeharms (2006, their Fig. 1.2)]. This configuration is the most favorable for Agulhas ring shedding because the necklike protrusion is almost free of topographic effects, and the eddies are shed directly into the South Atlantic. Therefore, it is important to know under what conditions this protrusion is stable, and when the Agulhas Current is more likely to retroflect directly after passing the slanted eastern slope of the Agulhas Bank.

### b. Numerical background

The Agulhas retroflection dynamics have been studied in realistic numerical models. According to Biastoch et al. (2009), the retroflection occurs at 17°E during shedding and 23°E (i.e., east from the Agulhas Bank) afterward. Similar scenarios can be seen in the Backeberg et al. (2009) and Tsugawa and Hasumi (2010) simulations.

In the ° global Hybrid Coordinate Ocean Model (HYCOM) (www7320.nrlssc.navy.mil/GLBhycom1-12/agulha.html), the usual position of retroflection is nearly 20°–22°E during eddy shedding and 23°–25°E after shedding, which is 4° eastward compared to the Dencausse et al. (2010a) data. The aforementioned eastward shift is probably due to the HYCOM’s inability to reproduce the zonalization of the Agulhas influx. It is still not clear which of the simulated parameters are responsible for this inconsistency between observations and simulations.

We note, however, that the fine westward protrusion configuration was simulated recently by Loveday et al. (2014, see their Fig. 3).

### c. Theoretical background

Widely accepted theoretical models indicate that Agulhas rings are shed primarily because of inertial and momentum imbalances. Nevertheless, the exact mechanism of shedding is still under discussion. For example, according to Ou and de Ruijter (1986), the Agulhas Current front, moving slowly southwestward, forms a loop soon after separating from the coast, due to the coastline curvature. The loop occludes and forms a ring, whose shedding is accompanied by the instantaneous eastward retreat of the front. This mechanism was further discussed by Lutjeharms and van Ballegooyen (1988a) and Feron et al. (1992).

A purely inertial shedding mechanism was proposed by Nof and Pichevin (1996) and discussed by Pichevin et al. (1999). According to this mechanism, ring shedding from the retroflection area is necessary to circumvent the so-called retroflection paradox. Specifically, the “rocket force” caused by the westward-propagating eddies balances the nonzero, combined zonal momentum flux (or flow force) of incoming and outgoing currents. Nof and Pichevin (1996) did not address the formation of the ring in the retroflection area and focused instead on the detaching phenomenon.

An alternative idea is that the momentum fluxes of currents are compensated by the Coriolis force, which is caused by the “ballooning” of the basic eddy (a generic term to designate the retroflection area), a mechanism suggested by Nof and Pichevin (2001) for the outflows and elaborated further by Nof (2005) for the case of a plane.

However, Zharkov and Nof (2008a) pointed out a “vorticity paradox”; in the case of retroflecting currents, both the momentum and mass conservation equations can only be satisfied for small Rossby numbers. To circumvent this paradox, the authors considered the incoming current retroflecting from a coastline with a slant greater than a threshold value of ~15°. Zharkov and Nof (2008b) and Zharkov et al. (2010) elaborated on the effect of coastal geometry on ring shedding for 1.5-layer models with slanted and “kinked” coastlines and showed that, in the case of a rectilinear coast, there is a critical slant above which there is almost no shedding. These results are in agreement with the numerical runs in Pichevin et al. (2009).

Nevertheless, recent theoretical models still have shortcomings. The basic equations are (i) the mass conservation, (ii) the momentum balance, and (iii) the Bernoulli’s principle. It is assumed that these three equations can be satisfied together when (i) the basic eddy radius is much larger than the widths of upstream/downstream flows and significantly grows, and (ii) the potential vorticity (PV) of the basic eddy is constant. In fact, when the near-linear dynamics of approximately geostrophic incoming/outgoing currents are transformed into the dynamics of a radially symmetric basic eddy of constant PV, the Bernoulli integral cannot be conserved because an additional portion of energy is required to raise the basic eddy’s nonlinear advection [this forcing term is as strong as the Coriolis force in the basic eddy momentum equation; see Eq. (17) below]. As a consequence, there are a number of contradictions in the previous analytical models. Two of these contradictions are reported in appendix A.

### d. Present approach

To address these contradictions, we develop a new model of current retroflection and eddy shedding. Contrary to previous models of retroflection, in which the vorticity of the basic eddy, as well as the vorticity of the incoming and outgoing currents, were treated as constants, we allow the vorticity to vary. (Actually, in the text, we consider the Rossby numbers characterizing the relative vorticities, although, in principle, it does not matter which type of vorticity we mean because the Coriolis force is assumed almost constant at the scale of retroflection area.) By doing so, we aim to overcome the aforementioned paradoxes.

The main goal of this study is to use this model to understand why the Agulhas Current retroflection protrudes westward from the Agulhas Bank. We also want to clarify the relation of this westward protrusion to the ring-shedding event and establish a simple criterion for the stability of each configuration of the Agulhas Current.

The paper is organized as follows: In section 2, we introduce the governing equations for the basic eddy development and derive an analytical model of a zonal retroflecting current. We also discuss the reasonable intervals for the initial value of the current’s Rossby number. In section 3, we compare this theoretical model with a numerical model for several values of the Rossby number. On the basis of this investigation, we propose a criterion for the stability of the westward-protruding configuration of the Agulhas Current being dependent on its relative vorticity. In section 4, we confirm this hypothesis using data from eddy-resolving numerical models. Last, we summarize and discuss our results in section 5. For convenience, we define all the variables both in the text and in appendix E.

## 2. Theoretical model

In this section, we derive an analytical model of Agulhas Current retroflection. We first describe the incoming and outgoing currents, then the basic eddy. All these elements are then connected using three equations expressing Bernoulli’s principle, the momentum balance, and the conservation of volume. We allow the main variables of the model to vary in time and attempt to describe their temporal evolution. Finally, we consider the stability of this system and establish a criterion for eddy shedding.

### a. Description of the currents

Consider the situation depicted in Fig. 2; a boundary current flows along a zonal coast (in the Southern Hemisphere) and retroflects at some point. To describe this system, we consider the Cartesian coordinate system , where marks the limit between the incoming and outgoing currents, and is placed at the center of the basic eddy at (a stagnation point of the retroflection).

As in preceding papers by Nof et al., we use a reduced-gravity model (the 1.5-layer shallow-water approximation). The incoming and outgoing currents of density flow above an infinitely deep, stagnant layer, whose density is . To avoid the “thickness paradox” (see appendix A, section a), we assume that the thickness of the upper (active) layer vanishes outside the currents, that is, at , where denotes the widths of both currents. This assumption can probably be justified by the fact that the Agulhas Current carries much warmer (although saltier) water from the tropical Indian Ocean, implying that there occurs an additional pycnocline at the interface between Agulhas Current water and surrounding waters, and this pycnocline really outcrops at the borders of the current [the border is strongly pronounced, e.g., in Fig. 1 of van Sebille (2009)]. We impose a linear shear profile such that the velocity in the currents is given by

with as the Rossby number of the currents and as the linear approximation of the Coriolis parameter. The full Coriolis parameter is (in the usual -plane approximation). The variation of over a distance is evaluated using a parameter

For a current of width km and with s^{−1} and m^{−1} s^{−1}, we have . Henceforth, we use this small parameter when we perform a Taylor expansion of a variable.

Away from the retroflection point in the positive direction, the meridional velocities in the currents vanish, so the momentum equation in the meridional direction is

We integrate Eq. (4), retaining the first-order term in :

where is the upper-layer thickness at . We adjust such that vanishes at the borders of incoming and outgoing currents, that is, at the leading order in , . Therefore,

We can also introduce the deformation radius of the currents

so that Eq. (6) can be rewritten as

which results in . It is easy to verify that at the first order as well.

We then calculate the incoming and outgoing volume fluxes and according to

Thus, at the leading order of approximation, the incoming flux scales with . For the typical Agulhas Current parameters ( Sv, s^{−1} and m s^{−2}), we obtain km.

where

is a constant small parameter and, for the Agulhas Current conditions,

As we see from Eq. (11), the net volume transport going to the eddies in the retroflection area is nonzero only at the first order of [i.e., scales like ] and therefore is only due to . Using Eqs. (11)–(13), it is expressed as

Equation (14) implies that the leakage of a retroflecting current is small compared to the incoming flux. This affirmation is verified hereafter. For the Agulhas Current parameters,

Note that we used the -order approximation only for the calculation of . The use of the leading-order approximation for the other parameters is verified by comparing our model outputs with numerical calculations (see section 3). In particular, keeping as a fixed parameter in our further analysis, we can consider to be constant as well.

### b. Description of the basic eddy

In the retroflection area, the incoming and outgoing currents are connected through the basic eddy. At the leading order of approximation, we assume that its shape is circular, so that its dynamics can be described in a polar coordinate system bound with the movement of the basic eddy’s center. Following Nof and Pichevin (2001), we express the basic eddy’s orbital velocity as

with twice the basic eddy’s Rossby number (we note here that characterizes a zero PV eddy). The radial symmetry condition is

which corresponds to the momentum equation for the azimuthal velocity in the new coordinate system. We can compute the eddy thickness profile by integrating Eq. (17) using Eq. (16) and the boundary condition (since vanishes at the basic eddy’s rim):

The eddy volume is

### c. System of equations

At this stage, there are three unknowns that may vary in time: , , and . The first two, being twice the Rossby number of the basic eddy and Rossby number of the currents, respectively, express the variation of vorticity (see section 1d); and describes the development of the basic eddy before its detachment from the retroflection area. The variable can be obtained from Eq. (8), and we have Eq. (12) for . The other parameters are kept constant. Therefore, we need three more equations to close the system. These equations express the connection between the basic eddy and the incoming/outgoing currents and are derived from three basic principles we describe in the next subsections: Bernoulli’s principle, the momentum balance, and the conservation of volume.

#### 1) Bernoulli’s principle

This principle states that, for a steady state, the value of

is constant along a streamline. Thus, if we choose the free streamline (line of outcropping between the “gray” area of nonzero and “white” area in Fig. 2), we obtain

Using Eq. (8), we can rewrite this equation as

According to Eq. (23), there are two unknowns left: and . The time rate of the basic eddy evolution is

#### 2) Momentum balance equation

The derivation of an equation expressing the balance between the momentums of joint incoming/outgoing currents and the Coriolis force of the growing basic eddy is given in appendix B (see also Nof and Pichevin 2001). Finally, we have

where is the basic eddy center propagation speed in the meridional direction. Henceforth, we assume that . At the leading-order approximation, the calculation of the integral in Eq. (25) using Eqs. (1) and (5) and substitution of Eq. (19) for yield

The physical interpretation of Eq. (27) is that the momentums of incoming and outgoing fluxes are balanced by the Coriolis force, owing to the off-wall movement of the basic eddy center due to the growth of its radius.

#### 3) Conservation of volume

The last equation that we use to close the system is the conservation of volume:

The expression for the left-hand side, , is derived from Eq. (19):

After substituting Eq. (14) into the right-hand side of Eq. (28), and Eqs. (29), (23), and (26) into the left-hand side, we obtain

Interpreting Eq. (30) physically, we expect to be negative because the volume transport going to the basic eddy growth is small, so that the strong growth of should be, for the most part, compensated for by the decrease in and therefore in the basic eddy thickness [see Eq. (18)]. Hence, the “spreading” of the basic eddy (i.e., increasing its diameter but decreasing thickness) is the dominating process compared to its volume increase—this is a principal difference of our model from preceding ones. Also, we show in the next subsection that, despite the smallness of , its account in Eq. (30) is essential.

Now, Eqs. (27) and (30) [with substitution of Eq. (12) for ] form a system of two ordinary differential equations (ODEs) that can be solved numerically. To do so, we need to specify the initial value for the Rossby number of the incoming flux:

where stands for the “initial condition.” We also assume that , so the second initial condition is

The range of validity of is discussed in the next subsection.

### d. Stability of the system

We just derived a system of equations describing the evolution of a retroflecting current and the basic eddy. Let us recall the hypothesis made to obtain the final set of equations.

The transport is accurate if the value of is small, that is, the variation of the Coriolis parameter is small compared to its average value . This condition limits the width of the incoming current to km so that remains smaller than 0.1.

The model is also valid under the assumption of a small Rossby number for the incoming current. Here, we remind the reader [see our comments for Eqs. (27) and (30)] that the basic eddy radius should grow to ensure the balance between the Coriolis force and the momentum forces of the currents. The radius can grow on account of both mass flux from the currents and spreading of the basic eddy. However, according to Eq. (8), the adaptation of the system to the basic eddy formation requires the alteration of the currents’ vorticity, implying the change in the widths of the currents. To be more precise on the limit of validity of the model, we introduce the nondimensional function characterizing the relation between these two processes:

which corresponds to the ratio of the time rates of and change. For our convenience, we also introduce . The situation where corresponds to the case where the basic eddy expands slower than the width of the currents, a situation that is physically unacceptable. Henceforth, we refer to as the reasonable interval and as the unreasonable interval. To obtain an analytical value of , we expand the solutions , , and in the Taylor series around . After some algebra introduced in appendix C, we obtain

It is seen from Eq. (35) that, if (or ), is more than 3 for any value of , tending to infinity at . Physically, this means that the widths of incoming and outgoing flows increase faster than the radius of the basic eddy. Reasonable values of should be chosen such that . From Eq. (35), we conclude that should be less than , which is the root of

For small values of , we obtain the condition

This is an important condition that defines the limit of validity of the model, and we make use of it henceforth. Using Eqs. (37) and (33), we obtain an estimate for the upper bound of reasonable values of :

In section 3, we compare the solutions of the system of Eqs. (27) and (30) with numerical simulations. The goal is to investigate (i) whether is indeed a reasonable condition, (ii) what happens in simulations of our models and numerics when is inside and/or outside this interval, and (iii) whether or not we can extend this interval.

### e. Eddy-shedding criteria

We finish our theoretical analysis with some considerations on the conditions of the basic eddy detachment from the retroflecting current.

#### 1) First criterion

To be shed, the ring should be large enough to escape the front propagating behind it. Calculating the average front speed

we obtain

Moreover, according to Nof (1981) and Zharkov and Nof (2008b), the propagation speed for the open ocean lenses is

Using Eq. (23) for the value of and equating [Eq. (40)] with [Eq. (41)], we obtain the condition of the eddy detachment:

Remember that, initially, , and, according to our simulations (described in section 3), both and decrease with time. However, for Eq. (42) to be satisfied at some moment of time, should decrease faster. As we show in appendix C, this occurs in a reasonable interval of (when ).

#### 2) Second criterion

An alternative choice would be equating the currents’ frontal velocity with the eddy rear front velocity. As seen in Fig. 3, the zonal propagation speed of the basic eddy rear front changes from for to for , where is the angle in the polar coordinate system. The average rear front speed is

Equating [Eq. (43)] with [Eq. (40)] and using Eq. (41) for and Eq. (22) for [see section 2c(2)], we obtain

Since both and are less than one and we consider the case where decreases faster than , the right-hand side of Eq. (44) is expected to decrease quickly and become less than one, while the left-hand side grows and becomes positive and of order one. Therefore, we expect Eq. (44) to be satisfied at some time (the second criterion thereafter).

We end the development of the theoretical model here. In summary, the time evolution of the system of incoming and outgoing currents and the basic eddy in the retroflection area are characterized by Eq. (8) for , Eq. (11) for and in terms of , Eq. (12) for , Eq. (23) for , and a system of Eqs. (27) and (30) for and with initial conditions [Eqs. (30) and (31)]. The instant of the eddy detachment is defined by Eqs. (42) (first criterion) or (44) (second criterion). We also showed that the solution of the system has physical meaning when initially the parameter does not exceed the value defined by Eq. (38). In the following section, we proceed to time integrations and compare them to a 1.5-layer shallow-water model.

## 3. Numerical simulations

In this section, we compare the theoretical model developed in the previous section with several numerical simulations. We use a modified version of the Bleck and Boudra (1986) reduced-gravity isopycnal model with a passive lower layer and the Orlanski (1976) second-order radiation conditions for the open boundary.

### a. Model setup

The parameters are m s^{−2}, s^{−1}, and Sv. The other parameters are listed in Table 1. We either use a *regular* value for ( m^{−1} s^{−1}) or a *magnified* value for ( m^{−1} s^{−1}). The initial condition is taken between 0.04 and 0.25. The size of the domain is adjusted so that the boundary current remains away from the influence of the boundary of the domain.

The northeastern section of the domain is filled by land, and the position of the corner of the coastline is given in Table 1. We use the free-slip boundary condition at the “coastal” walls.

The spatial resolution is km; the time step is 30 s. We use a Laplacian viscosity coefficient m^{2} s^{−1} (the minimal value for the stability of long time simulations). Such a choice of parameters is adequate because the diffusion speed () is 0.1 m s^{−1}, which is small compared to the rings’ orbital speed (of the order 1 m s^{−1}). Our simulations are significantly nonlinear since, according to Dijkstra and de Ruijter (2001a,b), viscous effects become dominant when m^{2} s^{−1} (for their forced model of Agulhas retroflection). This conclusion is also in agreement with Boudra and Chassignet (1988).

All the simulations are run for more than 700 days and initialized with a fully retroflecting current at . We use an open-channel configuration with both the incoming current (adjacent to the zonal wall) and the outflow having streamlines parallel to the zonal wall. The initial cross-channel velocity profile is linear, and the thickness profile is parabolic. The inflow and outflow are connected by circular streamlines. The “center” of retroflection is initialized at km [ km for experiment (expt) 1].

### b. Snapshots of the numerical simulations

In Figs. 4–7, we plot the snapshots of the four experiments listed in Table 1. The time intervals between these snapshots are 150, 50, 50, and 10 days, respectively. We adjust these intervals according to the eddy-shedding periods in each experiment (several hundred days for , 170–180 days for , and 30–35 days for ).

The figures show the upper-layer contours through 100-m increments. Some contours are marked in meters. The scales on the coordinate axes are in kilometers.

In all these experiments, at least one eddy forms and detaches from the current. We note that the position and the propagation speed of detached eddies do not change significantly as we move the meridional wall relative to the retroflection, even in the case when the retroflection is strongly protruding into the “open ocean.”

### c. Comparison with the analytical model

Figures 8–11 show the comparison of the results of our theoretical model with the outputs of the corresponding numerical model for the four experiments (see Table 1).

The curves of and (solid lines in Figs. 8a, 9a, 10a, and 11a) are obtained by integrating the system of Eqs. (27), (30), and (12) using a fourth-order Runge–Kutta method. They can be directly compared to the curves of and obtained in the numerical model and fitted by polynomials of degree 5 (dashed lines in Figs. 8a, 9a, 10a, and 11a). In the latter case, these values are obtained by dividing the averaged values of the vorticity by the Coriolis parameter.

Similarly, we plot with solid lines (theoretical model) and dashed lines (numerical model) for and theoretical lines for in Figs. 8b, 9b, 10b, and 11b; , , and in Figs. 8c, 9c, 10c, and 11c; and and in Figs. 8d, 9d, 10d, and 11d. (In appendix D, we explain the plotting of numerics in detail and give an example of a plot incorporating individual data points for Fig. 10. We also plot a numeric curve for there; the reason we do not do it in all the Figs. 8–11 is discussed below in section 3d.) We restrict the comparison to a time interval defined by the instant where the first eddy detaches from the current. Indeed, after this event, the system theoretically returns to its initial condition (which is not the case for the system of ODEs that has to be manually reset). However, the initialization of our numerical simulations is not reproducible. After the eddy shedding, the new initial value of (for the development of the next eddy) is not quite clear in numerics since, owing to viscosity, the process of shedding is shadowed (i.e., the system does not return to the initial condition instantly after shedding the first eddy, and it is not clear how to define the exact instant for starting the formation of the next one). To define the eddy detachment periods, instead of solving Eqs. (42) and (44), we consider the intersection of and (Figs. 8c, 9c, 10c, and 11c) and and (Figs. 8d, 9d, 10d, and 11d). Note that we derived Eqs. (42) and (44) by equating with and , respectively.

In Fig. 8, we compare the two models for the configuration of expt 1 (). We find an excellent agreement in both and (Fig. 8a) and and (Fig. 8c). The theoretical plot of (Fig. 8b) looks like “time-averaged” numerical data, so the agreement is adequate as well; the value of the radius in the theoretical model is overestimated (10%) compared to the numerical model. In Fig. 8b, we see that grows from 300 km initially to 550 km after 1000 days, whereas is almost constant. This behavior coincides with (see section 2d), which is in the interval of validity of our model. The lower panels show that the period of detachment should be about 210–260 days, according to the first criterion (Fig. 8c), and 320–250 days, according to the second criterion (Fig. 8d). The numerical simulation suggests a higher value of 800 days (probably associated with the viscous “jamming” of the first eddy in the numerical simulation). (Also, the simulation probably does not perform the “complete” separation owing to the viscosity and the eddy spreading kind. Instead, a “numerical quasi separation” occurs, when the maximal depth of the upper layer in the “neck” between eddies is about 10 times less than the depth of the eddy itself.)

In Fig. 9, we compare the results of expt 2 . The variations of are well captured by the theoretical model, whereas is underestimated (Fig. 9a). From Fig. 9b, we see grows faster than (i.e., ), which means we might be outside the validity range of our model. Also, the agreement between the theoretical and experimental curves in Figs. 9c and 9d is worse than for expt 1 (Fig. 8).

To remain in the validity range, we used a magnified value of (expt 3). The results of this experiment are plotted in Fig. 10. In this case, we have a fairly good agreement for both and . In Fig. 10b, we see that the growth rate of and is comparable during the first 50 days ( is approximately one). In other words, exceeds but not significantly. After 50 days, keeps growing while remains stable (). According to this behavior, we define a “transient” interval where

Unfortunately, the values of (Fig. 10c) and (Fig. 10d) in the theoretical model are significantly larger than in the numerics, probably because of the viscous effects (here we note that, in our preceding investigations, the -simulated eddy propagation speed in the open ocean was about two times or more lower in numerics than that predicted by theory). For this case, the theoretical curves of and (Fig. 10c) and and (Fig. 10d) do intersect, showing the values of the detachment period according to our first and second criteria. The numerical curves rather touch each other, probably because of the fitting errors.

The results of expt 4 are plotted in Fig. 11. While the tendency for is correctly captured by the theoretical model, is strongly underestimated (Fig. 11a). The values of current width increase much faster than values of the basic eddy radius (see orange and blue solid lines in Fig. 11b), at least at the initial stage of the basic eddy development, corresponding to . The theoretical curves of and (Fig. 11c) and and (Fig. 11d) do not intersect, implying that the shedding criteria [Eqs. (42) and (44)] are not met. This is not in agreement with Fig. 7, where we see at least two eddies shedding in 60 days. The experimental curves do intersect and predict the detachment period of about 16 days by both criteria, but the agreement between theoretical and experimental curves is poor.

### d. Physical interpretation

On the basis of our comparisons, we suggest that, without loss of generality, the agreement of our theoretical model with numerics can be characterized by the values of instead of . This agreement is strong for (reasonable interval), moderate for (transient interval), and weak for (unreasonable interval). For convenience, we introduce as a value of at which . Figure 12 shows the plots of versus for vanishing (i.e., plane) and regular and magnified . The Roman numerals *I*, *II*, and *III* denote reasonable, transient, and unreasonable intervals, respectively. Based on Fig. 12, we conclude that (expt 1) is exactly at the border between reasonable and transient intervals for regular . The value of (expts 2 and 3) falls into the transient interval only for magnified and in the unreasonable interval for a regular value of . The higher value of (expt 4) is in the unreasonable interval.

We suggest that, when is in the unreasonable interval, large values of indicate that the Coriolis force momentum of a single shed eddy is still insufficient to balance the momentums of incoming and outgoing fluxes. Therefore, the shedding of several rings during the period defined by Eq. (42) (first criterion) or Eq. (44) (second criterion) is required to resolve the retroflection paradox. Indeed, as we can see from numerics for expts 2 and 4 (Figs. 5, 7), the retroflection does not stay at the same location after shedding the first eddy. Rather, it propagates upstream, forming a chain of eddies on its path. (This is the reason we avoided plotting the numerical values of ; it is not quite clear what is the width of the current that is gradually transforming into eddies, for which has no sense.) The larger the theoretical value of , the faster is this upstream propagation of the retroflection. So, the retroflection of a zonal current under these conditions looks unsteady. We assume that the retroflection of the jet flowing along a slanted coast is more stable because, in this case, smaller current momentums need to be balanced by the basic eddy Coriolis force momentum. In addition, eddies (with finite radii) can “fill” only a very limited area between incoming and outgoing currents. Conversely, when is in the reasonable interval, the shedding of each recurrent eddy occurs approximately at the same position, and the shed eddies propagate westward (see Fig. 4). So, the reasonable and unreasonable intervals of can be reclassified “stable” and “unstable,” respectively.

Our main hypothesis is that the unsteadiness of the zonal current retroflection with insufficiently small relative vorticity could cause the Agulhas retroflection to shift eastward from the open South Atlantic basin to some location east of the Agulhas Bank, where the incoming current flows along the slanted continental shelf. This hypothesis is based on the fact that, in our numerical simulations, the retroflection propagates upstream when the Rossby number falls into the unstable interval.

## 4. Comparison with more realistic numerical models

In this section, we apply the results obtained in the previous sections to the interpretation of the realistic numerical model simulations.

We first consider the study by van Sebille et al. (2009). The authors investigated the relation between the Agulhas retroflection position and the magnitude of the Agulhas leakage to the South Atlantic. The ° ocean numerical model (Biastoch et al. 2008a,b) was used. The retroflection location was defined based on the sea surface height isolines. The authors show that the Agulhas leakage is a function of the latitude of the retroflection point (see their Fig. 9). This figure is redrawn here as Fig. 13, upper panel. It shows the Agulhas leakage transport versus the position of the Agulhas retroflection. Van Sebille et al. (2009) highlighted a significant correlation between these two variables.

Assuming that approximately 40% of Agulhas leakage is carried by rings (Doglioli et al. 2006; van Sebille et al. 2010), we multiply the values of by 0.4 to obtain approximate values of and then calculate using Eq. (15). The results are shown in Fig. 13, lower panel. Here, as in Fig. 12, the Roman numerals *I*, *II*, and *III* denote stable, transient, and unstable intervals, respectively, for the regular value of . The borders between intervals are the dashed–dotted lines, where and . The blue diamonds correspond to the red circles in the upper panel (data from van Sebille et al. 2009), and the red line corresponds to the linear regression. If the regime is nearly stable ( is nearly ), the value of does not significantly change with time, as seen in Fig. 8. Therefore, we can compare the theoretical values of and directly with obtained values of . Most of the experiments with the retroflection location to the west of the Cape of Good Hope (~18°E) are concentrated in the stable interval (*I*); the regression line shows that characteristic values of the Rossby number are in the stable interval when the location is west of 16°E and in the unstable interval when the location is east of 19°E, which is close to the Agulhas Bank. This confirms our hypothesis; the retroflection is likely to shift toward a slanted coast when the Rossby number gets into the unstable interval.

Last, we consider the ° global HYCOM simulation (www7320.nrlssc.navy.mil/GLBhycom1-12/agulha.html). As could be inferred from the animation, the typical period of the Agulhas basic eddy rotation is about 4 days. Therefore, the characteristic value of is , so that is not less than 0.1, which is larger than values of and . Hence, according to our hypothesis, the zonal flow configuration for the retroflection is unstable, which is in agreement with the fact that, in simulations, the current retroflects directly from the slanted eastern slope of the Agulhas Bank.

## 5. Conclusions and discussion

### a. Summary

In this paper, we have presented a new model of a zonal current retroflection and eddy shedding. In contrast to the preceding models of retroflecting currents, this one satisfies, (i) the mass conservation equation, (ii) the momentum balance equation (i.e., to resolve the retroflection paradox), and (iii) the time-dependent, altered Bernoulli’s principle (i.e., the continuity of velocity along the free streamline bounding the currents and the retroflection area). Our main assumption is that the net mass flux going into the eddies is only due to the effect, circumventing the thickness paradox by avoiding the discontinuity of the upper-layer thickness between the rim of the retroflection area (basic eddy) and the incoming current at the zonal coastline (see appendix A). In this case, the increase in the basic eddy volume is very weak, so the momentum of the Coriolis force due to the basic eddy growth cannot compensate for the net momentums of incoming and retroflected currents. However, the vorticity variation spreads the basic eddy even without a significant volume increase, so that the center moves off the wall much faster than it does with a constant PV (i.e., owing only to the basic eddy volume growth), which is enough to balance the current momentums.

As shown, the system of aforementioned equations yields reasonable solutions for small values of the Rossby number (defined in the text as the stable interval). In this case, we obtained good agreement between our theoretical solutions and numerical simulations (Figs. 8, 10). Otherwise, the Rossby number is in the unstable interval. In the theoretical model, this unstable interval is characterized by an unphysical increase of the width of the incoming and outgoing currents (remember that, in our problem statement, we assumed the incoming mass flux and the current velocity to be constant). In this case, the agreement between our model and numerics significantly deteriorates (Figs. 9, 11). In numerics, we observe a strong eastward shift of the retroflection and the formation of a chain of eddies in its path (Figs. 5, 7). Based on these findings, we hypothesize that the retroflection of a zonal current becomes unstable for larger Rossby numbers and therefore the retroflection point is expected to shift toward the slanted coast (in the case of the Agulhas Current and the eastern slope of the Agulhas Bank).

This hypothesis is in agreement with the data of van Sebille et al. (2009). Using their data and Eq. (15), we show (Fig. 13) that the Rossby number most likely gets to the stable interval (defined in our theoretical model) when the Agulhas retroflection is located west of the Cape of Good Hope and to the unstable interval when it is located near and east of the Agulhas Bank. It is possible that some eddy-resolving models do not simulate the local “zonalization” of the Agulhas incoming flux after it passes the Agulhas Bank (which, in turn, can lead to underestimation of the leakage into the Atlantic) because the Rossby number is overestimated numerically, possibly because of an underestimation of mesoscale processes or a low viscosity value. The viscous friction is probably the main cause of the vorticity variation—we leave this mechanism as a subject for future investigations.

### b. Another physical interpretation

Our analysis of the Agulhas Current uses the Rossby number as the main variable. However, one can also describe the stability of the current using dimensional variables, namely, the width of the currents and the deformation radius . In fact, Eq. (35) can be rewritten with dimensional parameters in the form

A plot of as a function of and with s^{−1} and m^{−1} s^{−1} is shown in Fig. 14. The gray area corresponds to the unstable domain. We hypothesize that for these values of and , the Agulhas Current must retreat in a region where the coast is slanted. The white area is the stable region. For the corresponding values of and , the Agulhas Current can shed an eddy and remain in the zonal area.

The limit of stability (at leading order of approximation) can be reformulated with dimensional parameters:

or, in view of Eq. (11),

Using Fig. 14, we give a new interpretation of the retroflection position in realistic models; we anticipate that a wider Agulhas Current can help stabilize the current and give it the opportunity to elongate in the zonalization region. Since the viscosity effect is a possible mechanism of the vorticity alteration (decrease in ), this widening of the current can be achieved by increasing the viscosity coefficient. Also, according to Fig. 14, a lower value deformation radius would help stabilize the Agulhas Current, although in realistic assimilative runs like HYCOM simulations, the definition of is probably not as clear as in our reduced-gravity model.

A few comments should be made here. First, the widening of the current probably implies decreasing the current velocity. However, according to Eqs. (1) and (8), the velocity at the borders of the currents reaches , and, using Eqs. (11), (12), and (38), we conclude that the maximum velocity that can be reached when is the upper boundary of the stable interval should be

which is 1.1 m s^{−1} for the Agulhas Current parameters used in our model. Overall, this value is reasonable. We do not consider the entire current as a part of Indian Ocean Gyre, but focus only on the area not far upstream from the retroflection. According to Duncan and Schladow (1981), the current velocities not far from the Agulhas retroflection are in excess of 1.5 m s^{−1}. Gordon et al. (1987) gave an estimate of about 1.1 m s^{−1}. Therefore, the velocities are usually only a bit above the upper boundary of the stable interval, implying that the westward protrusion of the retroflection area is not as common as the retroflection east of Agulhas Bank but nevertheless is not rare.

Second, the viscosity in numerical models can affect the Agulhas leakage in different ways. Although, as we mentioned before, the viscosity increases the widths of the currents in our simulations, the main effect we consider is protrusion of the retroflection toward the South Atlantic, possibly resulting in the increase of Agulhas leakage, which is carried by the rings. We do not talk here about any other parts of leakage. On the other hand, Weijer et al. (2012) also show that, in their CCSM4 numerical model, the high viscosity simulation widens the current and indeed overestimates the Agulhas leakage. However, the reason is that a simulation of the Agulhas Current that is too viscous lacks the inertia necessary for retroflection and ring shedding. Instead, outflow to the South Atlantic is mainly due to direct leakage through a viscous sublayer. Such a scenario has no relation to the subject of our paper.

Overall, the current widening is due to an internal process that regulates the vorticity, which is not necessarily the effect of viscosity.

### c. Model applicability

This model can only be applied to the dynamics of Agulhas Current retroflection. Other retroflecting currents, such as the Brazil Current, North Brazil Current, and East Australian Current, have no zonalization because they retroflect directly from slanted coasts. Also, our present approach is not valid for meddies (Pichevin and Nof 1996) because, in that situation, the incoming and outgoing flows are separated by a wall, and the upstream retreat of the retroflection area is impossible.

The foremost shortcoming of the present model is that, for a reasonable Rossby number (getting into a stable interval), eddies are unrealistic (with radii about 300 km and more) and the period of shedding is extremely long (one or even several years). For a higher value of , we obtain an eddy detachment period of about 100 days, which is acceptable. The simulation with magnified rather than regular seems closer to reality because viscous damping compensates for the artificial acceleration. Concerning the large radii, observed radii are usually measured as distances between the center and point(s) of maximum radial velocity. The modeled linear distribution of radial velocity in eddies is idealized, and, naturally, the velocity profile is smoothed at the eddy’s periphery, which leads to underestimation of the eddy radius.

As for now, it is not quite clear how to describe the direct leakage and filaments in the momentum equations, so, in our theory, we assume that the net mass flux is completely converted into rings. Nevertheless, we test our hypothesis with observational data, assuming that 40% of the leakage is carried by rings. A possible explanation for this contradiction is that in some area close to the retroflection, the eddy shedding is probably the primary (and overwhelming) mechanism for westward transport of the water. At the same time, the direct leakage and filaments can be formed in a larger area containing both incoming and outgoing fluxes.

In our model, we focus only on the basic eddy detachment time scale and do not consider periodic solutions that would allow tracing the possible retreat at longer time scales. In the future, it would be interesting to develop the theoretical model aimed on the description of the upstream propagation of the retroflection.

Finally, our theory is developed only for lenses. Applying this theory to a finite-thickness upper layer is difficult in that, if the thickness is significant, increased eddy volume due to radial growth cannot be compensated for by the decrease in vorticity (and eddy thickness). This, in addition to the viscous effect, can cause the gradual collapse of the initialized 1.5-layer structure in numerical simulations.

## Acknowledgments

The study was supported by NASA Doctoral Fellowship Grant NNG05GP65H, LANL/IGPP Grant (1815), NSF (OCE-0752225, OCE-9911342, OCE-0545204, OCE-0241036, and OCE-1434780), BSF (2006296), and NASA (NNX07AL97G). Part of this work was done during Wilton Arruda’s visit to the Department of Earth, Ocean and Atmospheric Science at Florida State University supported by a post-doctoral fellowship from CNPq of the Ministry for Science and Technology of Brazil (Proc. 201627/2010-8). This research is also part of the activities of the Instituto Nacional de Ciência e Tecnologia do Mar Centro de Oceanografia Integrada e Usos Múltiplos da Plataforma Continental e Oceano Adjacente (INCT-Mar-COI, Proc. 565062/2010-7). We are grateful to Steve van Gorder for helping in the numerical simulations and to Erik van Sebille for presenting the data we used in Fig. 13. We also thank Donna Samaan for helping in preparation of the manuscript and Kathy Fearon and Meredith Field for assistance in improving the style.

### APPENDIX A

#### Contradictions in Previous Models of Retroflecting Currents

In this appendix, we describe two paradoxes (or inconsistencies) that were encountered in previous models of retroflections.

##### a. Thickness paradox on an plane

Equation (28) suggests that for a growing basic eddy. However, since we assume both incoming and outgoing fluxes are approximately in geostrophic balance (on an plane), we have

inside both. Therefore, by definition of incoming and outgoing fluxes,

and, to allow the eddy growth, we should have . According to Bernoulli’s principle, the thickness at the eddy rim is equal to , implying that, at the point where the eddy rim detaches from the wall, a jump in the thickness occurs (or, at least, the zonal gradients of both thickness and velocity are unnaturally large, which, in turn, leads to strong distortion of the basic eddy structure). This is a common case for all the preceding models of retroflection. We call this the “thickness paradox.” To resolve it, we suggest considering the retroflection on a plane instead of an plane and assuming that the net flux going to the eddies in the retroflection area is only due to .

##### b. Slowly varying problem statement

According to the mass conservation and momentum balance principles, the basic eddy develops by analogy with the ballooning outflow bulge. Its volume increases linearly with time, and its radius growth is [Nof and Pichevin 2001, Eq. (3.7)]. Under the assumption of constant vorticity, the basic eddy orbital velocity growth is too. Because the basic eddy kinetic energy is proportional to , its growth is , implying that the energy flux from the current into the basic eddy is not constant but essentially nonlinear in time. This conclusion contradicts the energetic sense of Bernoulli’s principle in terms of the slowly varying problem statement.

### APPENDIX B

#### Derivation of the Momentum Balance Equation

Taking into account the continuity equation, we write the momentum equation in the along-coast direction in the moving coordinate system as

Here, and are the velocities; (both negative) are the basic eddy center propagation speed projections on the axes of the fixed coordinate system. We do not consider the equation in the cross-coast direction because it involves an unknown force on the wall.

After integrating Eq. (B1) over the rectangular contour ABCDEA (Fig. 2) and using Stokes’ theorem, we get

where is the streamfunction, defined by , and is the area enclosed by the integration contour.

We use the slowly varying approximation, introducing the “fast” time scale (associated mostly with the basic eddy rotation) and the “slow” time scale . The ratio of the time scales is [see Eq. (13)]. Nof and Pichevin (1996, 2001) showed that, in similar problems, the terms with time derivatives in the momentum equations are negligible. Moreover, the second term in the left-hand side of Eq. (B2) is zero because at least one of the multipliers , , or vanishes at each segment of the contour . The main contribution to the third term in the left-hand side of Eq. (B2) comes from the segment *AB* (i.e., cross sections of both incoming and outgoing flows). Because both currents are approximately geostrophic (but not exactly geostrophic; see van Leeuwen and de Ruijter 2009; Nof et al. 2012), is approximately , which is [see Eq. (2)] and negligible compared to . Finally, in the last term of the left-hand side of Eq. (B2), we can neglect the change of the Coriolis parameter at the scale of the area , assuming . Also, with error , we assume the basic eddy to be circular. Under this assumption, we have .

### APPENDIX C

#### Reasonable Interval Criterion Derivation

We want to assess the limit of validity of the theoretical model. This is done via the characterization of the initial evolution of the system. We compute the expansion rate of the current and of the eddy near . The model is physically relevant when the eddy expands faster than the current ().

Near , we expand , , and in the Taylor series

where we used the initial conditions given in Eqs. (31)–(33). Even if is a function of , we keep the notation in Eq. (C3) for clarity. Using the two main equations of the model [Eqs. (27) and (30)], we obtain

and

and

Therefore,

Moreover, we noted in section 2e(1) that should decrease faster than . It follows from Eqs. (C1), (C2), and (C9) that

Therefore, the criterion is satisfied when , that is, again when .

### APPENDIX D

#### Plotting of Numerical Data in Figs. 8–11

Values of and are calculated from the average vorticity over the eddy and over the upstream currents far from the retroflection area, then fitted with polynomials;

is estimated by averaging the distance from the center to four points at the eddy boundary (north, south, east, and west);

is more complicated to directly measure from the model, so we calculate it from the maximum depth of the currents and use ; ;

is calculated through , where is the

*x*coordinate of the basic eddy center at time*t*;is the only parameter we differentiate from the polynomial for , since discrete differentiating of was noisy;

is also directly calculated from the data, using the coordinate of the easternmost point with maximum depth; and

is calculated according to Eq. (43) with using data for and .

As an example, Fig. D1 shows the same plot as Fig. 10 but with incorporated data points for all the parameters and the numerical data for (which is convenient for plotting because is not in the unreasonable interval). It is seen that the data for , , , , and are almost not dispersed, but there is a strong dispersion in , , and at the beginning of the basic eddy development process, which is strongly smoothed by fitting. Also, for this case, the theoretical and numerical data for are in good agreement.

### APPENDIX E

#### List of Abbreviations and Symbols

Coefficient with in expansion at

Segment of the integration contour

AMOC Atlantic meridional overturning circulation

Bernoulli’s integral

Coefficient with in expansion at

The basic eddy center propagation speed projections on the axes of the fixed coordinate system

Current front velocity

Eddy rear front velocity

Contour enclosing the integration area

Coriolis parameter

Absolute value of near the basic eddy center

Reduced gravity

Upper-layer thickness

Value of at the currents’ interface ()

HYCOM Hybrid Coordinate Ocean Model

Width(s) of current(s) at the zero-order approximation

Addition to at the first-order approximation (without factor)

Lower index meaning number of point number in numerical dataset (for , , , , etc.)

ODE Ordinary differential equation

PV Potential vorticity

Radius of the basic eddy

Deformation radius of currents

Polar radius in the coordinate system connected with the basic eddy center

Incoming current mass flux

Outgoing current mass flux

Ratio of the time rates of and change

Initial value of

Sv Sverdrup ( m^{3} s^{−1})

Time

Generation period

Mass transport going to the eddies

Particle velocity vector

Particle velocity projections on the horizontal axes

Volume of the basic eddy

Basic eddy orbital velocity

Cartesian coordinate system

Twice the time-dependent Rossby number of the basic eddy

Time-dependent Rossby number of the currents

, Initial values of and , respectively

Upper bound of reasonable values of

Value of at which

Meridional gradient of the Coriolis parameter

Difference in density between lower and upper layers

Horizontal grid

Small parameter ()

Initial value of

Coefficient with in expansion at

Small parameter ()

Viscosity coefficient (in numerical simulations)

Polar angle in the coordinate system connected with the basic eddy center

Upper-layer density

Streamfunction

Integration area around the basic eddy