Abstract

The authors study the splitting of a coherent vortex by a large-scale baroclinic background current. A criterion for the splitting of the vortex is defined, and the process is then studied numerically and analytically in a 2½-layer reduced-gravity model in which the vortex is represented by a potential vorticity (PV) patch in each layer. Three effects are important for the process: 1) the vortex “coherence,” which is a measure of the advective effect induced by the PV patches on each other; 2) the background current shear, which tears the vortex; and 3) the baroclinic β effect, associated with the background current PV gradient, which is shown to counteract the shear. When the baroclinic β effect is neglected, it is shown that PV patches oscillate around an equilibrium state, and they separate when the oscillation amplitude is larger than the splitting criterion. This model also shows that vortex core deformations play a (moderate) role when the vortex radius is larger than the first baroclinic radius of deformation. The baroclinic β effect substantially compensates the advective tearing and drastically reduces the oscillation amplitude. Thus, the vortex is able to resist much higher shear when the current PV gradient is taken into account. On the other hand, the baroclinic β effect also induces dispersion of the vortex, which is essential when the shear is strong enough. It is shown that, in fact, a vortex is generally scattered by Rossby waves before it is split.

1. Introduction

Mesoscale and submesoscale coherent vortices are observed in all ocean basins. Their typical scale is much smaller than the domain size: in the Atlantic Ocean, for instance, their radii range from about 150 km for Agulhas rings to 60 km for Gulf Stream rings and 30 km for Mediterranean water lenses (or even a few kilometers for eddies associated with convective plumes), sometimes with strong variations among each vortex type (e.g., McWilliams 1985; Olson 1980, 1991; Joyce and McDougall 1992; Olson and Evans 1986; van Ballegooyen et al. 1994; Pingree and Le Cann 1992, 1993a,b, 1994; Armi et al. 1989; Richardson et al. 1989; Chérubin et al. 1997; Paillet et al. 1999). They are able to trap fluid inside their cores and to propagate coherently for several thousand kilometers, thus providing a significant contribution to the global heat and salt transport. Their maximum vorticity is less than the local Coriolis frequency f0 (0.1–0.5f0 or so), and it decreases during their lifetime, which is on the order of one to several years.

The influence of the planetary vorticity gradient, referred to as β effect, on the dynamics of vortices has been extensively analyzed. Its role on the vortex propagation has been studied in McWilliams and Flierl (1979), Sutyrin (1987), Reznik and Dewar (1994), Sutyrin and Flierl (1994), Sutyrin and Morel (1997), Morel and McWilliams (1997), and Reznik et al. (2000). These studies have shown how the background potential vorticity (PV) field is distorted by the vortex circulation and how this modification of the background PV generates a secondary circulation called the “β gyre,” which in turn induces displacement and deformations of the vortex. A similar process is likely to take place whenever a vortex evolves in a background PV gradient.

In this study, we analyze the interaction of a vortex with “large-scale” background currents (currents with constant velocity field in the horizontal plane) associated with a PV gradient. When the background currents have no horizontal shear, barotropic currents (currents without vertical shear) have a simple effect: they merely advect the vortex at the current speed. Baroclinic background currents are vertically sheared and thus have a more complicated effect. In a stratified fluid, they are also associated with sloping isopycnals and corresponding horizontal PV gradients, because PV depends on the layer thickness variations. In the following, the PV gradient associated with background baroclinic currents will be called “baroclinic β effect” because of its similarity to the planetary β effect, and the β gyre will more generally refer to the secondary circulation that arises when the background PV gradient (in particular, associated with a background baroclinic current) is deformed by a vortex circulation. The influence of a current on the dynamics of a vortex has been studied by Meacham et al. (1990, 1994), Walsh (1995), Walsh and Pratt (1995), and, recently, Vandermeirsh et al. (2001, hereinafter VMS01). VMS01 have studied the influence of a large-scale baroclinic current on the vortex propagation. They showed that the advective effect of the current is always compensated by the baroclinic β effect, so that baroclinic large-scale currents only have a weak influence on the propagation of oceanic vortices.

A vertical shear also provides a mechanism for the destruction of oceanic vortices, which is the focus of this study. The tearing of a vortex by a vertically sheared current has been studied analytically by Hogg and Stommel (1990, hereinafter HS90) in a two-and-one-half-layer reduced-gravity model with point vortices and by Marshall and Parthasarathy (1993, hereinafter MP93) in a two-layer model with piecewise constant PV patches. In considering a vertically sheared mean current, they find regimes in which parts of the vortex core represented by point vortices (or PV patches) are not able to stay coherent and drift apart from each other when the current shear is strong enough in comparison with the vortex strength. However, the baroclinic β effect has been neglected in the HS90 and MP93 analytical calculations. MP93 considered also a few numerical examples that included the effect of nonuniform ambient potential vorticity. They found only a slight effect of the secondary flow on the vortex tearing process for small current shear but did not investigate this effect in detail.

The influence of the baroclinic β effect on the vortex tearing by baroclinic currents therefore remains unclear. In this paper we present an analytical theory and numerical simulations in a two-and-one-half-layer reduced-gravity configuration that demonstrate that the secondary flow associated with nonuniform ambient potential vorticity plays an essential role in the vortex tearing process: it increases the ability of the vortex to resist the current tearing effect, but it also provides a way to scatter the vortex.

In section 2, we define the equations and the model configuration. In section 3 we describe the numerical and analytical models used in our study, as well as the criterion for the splitting of vortices. In section 4 we consider the dynamics of vortex splitting without the baroclinic β effect, and in section 5 we discuss the influence of the baroclinic β effect. Application of the results to oceanic vortices, discussion of the model limits, and our general conclusions are presented in section 6.

2. The model configuration and the equations

As in HS90, we here consider a quasigeostrophic model with two active layers that overlie an infinitely deep and resting lower layer (see Fig. 1, in which Hk is the depth at rest of the kth layer and ρk is its density). We also take into account a background baroclinic current, with velocities Uk, that is zonal when the planetary β effect is taken into account. In each layer, PV consists of the background part βkY and vortex PV anomalies PVAk:

 
formula

where f0 is the Coriolis frequency; g1,2 = g(ρ2ρ1)/ρ1 and g2,3 = g(ρ3ρ2)/ρ1 are the reduced gravity of the interfaces between layers 1 and 2 and 2 and 3, respectively; βp is the gradient of the Coriolis parameter; Y is the meridional coordinate; Ψk and is the streamfunction associated with the vortex signature in the kth layer. Notice that the mean current PV gradient βUk varies between layers and represents the baroclinic β effect in addition to the planetary β effect (see VMS01).

Fig. 1.

Configuration used in this study. We consider a 2-layer system above an infinitely deep bottom layer. The vortex is composed of two PV patches (one in each layer) with radii R1 and R2 and strengths Q1 and Q2 in layers 1 and 2, respectively. Here, η represents the separation between the PV patch centers. A background current (U1, U2) is also taken into account

Fig. 1.

Configuration used in this study. We consider a 2-layer system above an infinitely deep bottom layer. The vortex is composed of two PV patches (one in each layer) with radii R1 and R2 and strengths Q1 and Q2 in layers 1 and 2, respectively. Here, η represents the separation between the PV patch centers. A background current (U1, U2) is also taken into account

The initial vortex structure is depicted in Fig. 1. The vortex core consists of two PV patches, one in each layer (with respective PV anomalies Q1 and Q2 and radii R1 and R2 in the upper and middle layers). The separation between the centers of each circular PV contour is denoted by η and represents the tearing of the vortex core.

To focus on the physics of the tearing process, we consider regimes suitable to quasigeostrophic dynamics; that is to say, we only consider vortices associated with small PV anomalies |Qk| ≪ f0. Further, we nondimensionalize all equations using the first baroclinic radius of deformation Rd1 (defined below) as the horizontal length scale and the inverse PV anomaly |Q1|−1 as the timescale. Thus, the equations of motion are (VMS01; Pedlosky 1987, chapter 6, section 16)

 
formula

Here, J(A, B) = ∂xAyB − ∂xByA is the Jacobian of A and B; t = T|Q1| is the nondimensional time; x = X/Rd1 and y = Y/Rd1 are the nondimensional east and poleward coordinates; β = βpRd1/|Q1| and βUk = βUkRd1/|Q1| are the nondimensional β coefficients, which are small for strong vortices; and Uk = Uk/Rd1|Q1| is the nondimensional large-scale current in layer k.

In this study, the background stratification is fixed corresponding to the following typical oceanic conditions: f0 = 10−4 s−1, H1 = H2 = 500 m, g1,2 = 10−2 m s−2, and g2,3 = 0.5 × 10−2 m s−2, yielding internal radii of deformation Rd1 ≃ 29 km and Rd2 ≃ 12 km. The corresponding nondimensional stratification is thus given by F+1 = 1.7, F2 = 1.7, and F+2 = 3.4.

In our basic configuration, the background current will be nonzero only in the upper layer (U2 = 0). In this case, the interface between the lower (infinitely deep) and middle layers is horizontal, and the interface between the upper and middle layers slopes and the PV gradients in the upper and middle layers then have opposite signs (βU1 = −βU2, exactly as in MP93). Without the planetary beta, baroclinic instabilities develops (Pedlosky 1987, chapter 7) and can perturb the interaction between the large-scale flow and the isolated coherent vortex in the long-time-duration run. To ensure stability of the background flow in our 2½-layer configuration, either an appropriate current in the middle layer or a strong enough planetary β has to be taken into account. This will be discussed in section 5c.

In this study, we consider a simple PV anomaly structure with Q1 = Q2 and R1 = R2. Without loss of generality, we consider anticyclonic vortices, so that the nondimensional PV anomalies are Δk = Qk/|Q1| = −1. In this configuration, the vortex characteristics strongly depend on the vortex radius: the maximum vortex velocity increases while the rotation rate at the center decreases with the vortex radius (see Fig. 2). Note also that, although PV anomalies are the same in both layers, the vortex circulation is weaker in the middle layer than in the upper layer because of the asymmetrical stratification. This is closer to reality than the configuration with F+2 = 0 considered in MP93.

Fig. 2.

(a) Maximum azimuthal velocity Vmax and (b) rotation rate Ω as a function of the vortex radius. The solid line is associated with the upper layer, and the dashed line is associated with the middle layer. Note how the rotation rate decreases with the vortex radius

Fig. 2.

(a) Maximum azimuthal velocity Vmax and (b) rotation rate Ω as a function of the vortex radius. The solid line is associated with the upper layer, and the dashed line is associated with the middle layer. Note how the rotation rate decreases with the vortex radius

In the following, the variable parameters of the study are the background current shear U = U1U2 and vortex radius r1 = R1/Rd1.

3. The approaches

a. The numerical model

In the following, numerical results are obtained with a pseudospectral code, described in Dewar and Flierl (1987). The domain is square and biperiodic with a width of about 12 vortex radii to avoid interaction of the vortex with the periodic continuation of the velocity field. The horizontal resolution is 128 × 128, giving a grid size Δx ≃ 0.1r1. A weak biharmonic vorticity diffusion term is used with a nondimensional coefficient (ν ≃ 4 × 10−8) to suppress small-scale numerical noise.

b. The analytical model

The analytical model used in the study is a generalization of MP93 in which the baroclinic β effect can be taken into account but the PV contours remain circular (as in HS90 and MP93). The different analytical steps are given in the appendix, and the basic idea behind its derivation is that three processes intervene in the evolution of the PV center separation: mutual PV patch interactions, advection by the background current, and β-gyre development. Each process is evaluated independently, and the evolution of the PV center separation (written in complex form η = xc + iyc, where i is −1⁠) is a simple superimposition of these three effects and can be written (see the appendix)

 
/dt = Ω(|η|) + U + βterm(t),
(3)

where the following apply:

  1. The first term represents the mutual advection effect of the PV patches. It is simply associated with the velocity induced by a patch at the center of the other one and is calculated under the hypothesis that each patch remains circular.

  2. Parameter U = U2U1 is the background shear effect.

  3. Parameter βterm(t) is the β-gyre effect, which is associated with the background PV deformation. It is calculated by assuming that this deformation is simply due to advection by the vortex initial symmetric circulation. Once the deformation of the background PV field has been evaluated, it can be inverted to derive the associated β gyre (see VMS01).

Notice that the Rossby wave dynamics are not represented in Eq. (3), because βterm is associated with a single azimuthal mode (mode 1; see the appendix and VMS01), which is the most important for the vortex propagation. When βterm is neglected, Eq. (3) reduces to the MP93 model.

c. Methods

The numerical and analytical models presented above will be used to calculate the critical shear necessary for the splitting of a vortex and to understand the physics of this process. Indeed, a comparison between the results of different models permits assessment of the importance of each effect.

  1. As in Eq. (3), all processes are taken into account except the deformation of the PV contours and the general Rossby wave dynamics. The differences between analytical model and numerical solutions can be attributed to vortex deformation and/or Rossby wave effects.

  2. In a similar way, neglecting βterm in Eq. (3) and then comparing with the results of the full equation permits an assessment of the influence of β-gyre development. Also, a comparison with numerical solutions in which βUk has been artificially set to 0 allows examination of the influence of vortex deformation alone (without the influence of the Rossby waves, which are absent in that case).

d. A criterion for splitting

To analyze the sensitivity of vortex splitting to different parameters, we also have to define a precise criterion that indicates when a vortex core consisting of two PV patches will be considered as split.

The results presented below show that splitting is, in general, very rapid, and we found that an appropriate criterion can be approximated as

 
ηr1→2Vmax + r2→1Vmax for t ∈ [0, 100],
(4)

where η is the separation between PV centers, t is time1 nondimensionalized such that Δ1 = −1 as in VMS01 (thus t = 100 time units represents about 10 vortex rotation periods), r1→2Vmax is the radius at which the velocity induced by the PV patch in layer 1 on layer 2 reaches its maximum (V1→2max; see Fig. 3), and r2→1Vmax is the radius at which the velocity induced by the PV patch in layer 2 on layer 1 is maximum. The physical ground behind the criterion in Eq. (4) is associated with the decrease of the velocity field induced by a PV patch in the adjacent layer when rrVmax (see Fig. 3). Indeed, if a particle located beyond rVmax is displaced outward by an external process (such as a background current), the influence of the vortex on the particle becomes weaker and the separation is then likely to increase further.

Fig. 3.

Schematic diagram of the velocity induced by a PV patch in the adjacent layer. The maximum velocity Vkk±1max is reached at a radius rkk±1Vmax that can sensibly differ from rk, the radius of the PV patch

Fig. 3.

Schematic diagram of the velocity induced by a PV patch in the adjacent layer. The maximum velocity Vkk±1max is reached at a radius rkk±1Vmax that can sensibly differ from rk, the radius of the PV patch

It can be shown that, for small vortices (r1 ≤ 0.5 or so), r1→2Vmax is much higher than r1. This means that, for small vortices, splitting occurs when the PV patches are fairly far apart. This is illustrated in Fig. 4 for a vortex radius r1 = r2 = 0.15. For this radius, rVmax ≃ 1.7r1, much higher than the PV patch radii. For such structure, we found that the critical shear for splitting is Ucrit = 0.0052. Figure 4 is associated with a shear U = 0.0051 slightly below the latter, and the vortex stays coherent even if the separation distance reaches about 4 vortex radii (which roughly represents the previous criterion). However, when r1 ≥ 0.5, which is the case of most oceanic coherent vortices, r1→2Vmax is roughly equal to r1, and the vortex can be considered as split when the PV patches no longer overlap. This result shows that strong tilting of vertical axis of a structure can be achieved before it loses its coherence. In the ocean, significant tilting has been observed for vortices that did remain coherent (Richardson et al. 1989; Walsh et al. 1996).

Fig. 4.

Upper (solid line) and lower (dashed line) PV contour evolution at time t = 20, 300, 600, and 800 (nondimensional units) for a vortex structure Δ1 = Δ2 = −1 and r1 = r2 = 0.15 and for a background current (U1, U2) = (0.0051, 0). Note how the structure stays coherent even if the separation reaches almost 4 vortex radii

Fig. 4.

Upper (solid line) and lower (dashed line) PV contour evolution at time t = 20, 300, 600, and 800 (nondimensional units) for a vortex structure Δ1 = Δ2 = −1 and r1 = r2 = 0.15 and for a background current (U1, U2) = (0.0051, 0). Note how the structure stays coherent even if the separation reaches almost 4 vortex radii

4. Vortex splitting without the baroclinic β effect

As a preliminary consideration, we revisit MP93 results for our 2½-layer configuration and artificially set βUk = βterm ≡ 0. The realism for the this set up will be discussed in section 6. The numerical and analytical models provide a way to explore the sensitivity of the vortex splitting to different parameters and to understand the governing physics.

a. A sufficient condition

When the baroclinic β effect is zero, one can expect the splitting criteria to be reduced to a simple kinematic condition. Indeed, when βterm = 0, a sufficient condition for splitting can be derived from Eq. (3):

 
U > max[|η|Ω(|η|)].
(5)

This condition states that the maximum velocity induced by the PV patches on each other has to be weaker than the background velocity shear. Also notice that Eq. (5) is associated with the condition for the existence of a stationary state (HS90). Indeed, for a background shear weaker than max[|η|Ω(|η|)], there exists a steady configuration (∂tη = 0 for η = ηs) given by

 
sΩ(|ηs|) = −U.
(6)

When all β effects are neglected, we can thus expect oscillations around the latter steady state when it exists, and splitting otherwise. We will see below that the dynamics are not so simple.

b. Results

Figure 5 represents the critical shear U necessary to split the vortex as a function of the core radius r1 and with βUk = 0. Crosses represent results from numerical experiments. For a given radius r1, each cross represents the critical shear beyond which Eq. (4) is met and the vortex is split. Below this value, the vortex remains coherent. The solid line corresponds to the analytical model: Eq. (3) predicts that Eq. (4) is met and the vortex is split2 for background shears below the solid line, and it remains coherent for background shears above it. The dashed line represents the sufficient condition in Eq. (5).

Fig. 5.

Critical shear U necessary to split the vortex, as a function of its radius r1. The baroclinic β effect is neglected. The solid line is the result for the analytical model, crosses represent numerical results, and the dashed line represents the sufficient condition for splitting

Fig. 5.

Critical shear U necessary to split the vortex, as a function of its radius r1. The baroclinic β effect is neglected. The solid line is the result for the analytical model, crosses represent numerical results, and the dashed line represents the sufficient condition for splitting

Figure 5 shows that all models are qualitatively consistent: increasing the vortex radius (or equivalently decreasing the stratification or increasing the layer depths to decrease Rd1) increases the critical shear necessary to split the structure. This is expected because the vortex coherence increases with the vortex size.

c. Interpretation

Results from the analytical model in Eq. (3) are close to numerical results when r1 ≤ 1, but some deviation appears for larger vortex radius. Because the PV patch deformations are neglected in Eq. (3), the difference between the solid line and crosses can be attributed to the deformation of the vortex. Indeed, numerical results reveal that for large radii, the deformation is very strong and sometimes leads to filamentation of one of the vortices, which thus becomes weaker. As a result, the vortex coherence decreases, and a weaker shear is necessary to split the structure.

As seen in Fig. 5, the sufficient condition in Eq. (5) overestimates the critical shear necessary for splitting: the condition in Eq. (5) is sufficient but far from necessary, which also means that the existence of a stationary state does not warrant coherence of the structure. In fact, this emphasizes the influence of the initial condition (the PV patches are vertically aligned) and subsequent oscillation induced by the background shear. Indeed, because the PV patches are initially aligned (η = 0 at t = 0), the separation does not tend to the stationary state ηs but oscillates around it (see also MP93). The oscillation amplitude corresponds to the maximum vortex separation and can be larger than ηs. It can thus reach values higher than r1→2Vmax + r2→1Vmax, and, in such a case, splitting is expected.

We thus conclude that PV patch deformation favors splitting but only plays a (moderate) role for large vortices. In addition, the oscillation of the structure is essential for the splitting process, and with the aligned initial state, splitting is possible even if the background current shear is smaller than the PV patch mutual advection.

5. Vortex resistance with βterm ≠ 0

In the following, we study the influence of the baroclinic β effect by comparing the previous results with the solutions obtained when the β-gyre development is taken into account in the numerical and the analytical models.

a. Results

Figure 6 represents the same results as Fig. 5 does but with the baroclinic β terms taken into account. The solid line represents the results obtained with the analytical model, and the crosses indicate numerical results. The previous results without baroclinic β effect have been superimposed (dashed line) for a better comparison. It is obvious that the critical shear necessary to split the structure is now much higher. When vortex radius r1 ≥ 1 or so, the critical shear is almost 2 times as large.

Fig. 6.

Same as Fig. 5 except that the baroclinic β effect is taken into account. The solid line is the result from the analytical model, crosses represent the validation by numerical tests, and the dashed line is for the previous results without β effect. Notice how the critical shear is modified and almost doubled. This is a feature of the compensation of the current advective properties by the baroclinic β effect

Fig. 6.

Same as Fig. 5 except that the baroclinic β effect is taken into account. The solid line is the result from the analytical model, crosses represent the validation by numerical tests, and the dashed line is for the previous results without β effect. Notice how the critical shear is modified and almost doubled. This is a feature of the compensation of the current advective properties by the baroclinic β effect

Notice that the results from numerical experiments are only shown for r1 ≤ 1. It was, in fact, difficult to interpret the numerical solutions in terms of splitting for r1 > 1 or so, because in this case Rossby wave dynamics become very important. This problem is discussed in detail below.

b. Interpretation

The differences between the solutions without and with the baroclinic β effect for r1 ≤ 1 shows that the latter reduces the effectiveness of tearing by the background shear. This effect is general and was expected. VMS01 have indeed shown that baroclinic large-scale currents have weak influence on the propagation of coherent vortices because the baroclinic β effect compensates the advective effect of the current [VMS01 have indeed shown that βterm(t) → −U when t → ∞], so that the effectiveness of tearing by the background current decreases with time. The reason for this is associated with the distortion of the background PV field by the differential vortex advection, as long as the vortex is stronger than the background current and can be considered to be coherent. Indeed, as discussed in VMS01, there is a direct correspondence between background PV and the background velocity field tearing the vortex. As shown in Sutyrin and Flierl (1994) and VMS01, the circulation induced by the vortex strongly distorts the background PV field and, after a few vortex turnovers, leads to the development of small-scale PV structures, rolled up into a spiral, where opposite sign PV filaments alternate along a vortex radius. As the streamfunction or current field averages out the small-scale PV structure, the background current becomes weaker and weaker in all layers where the vortex signature is strong, so that its tearing properties weaken (again, we refer to VMS01 for a detailed discussion).

As a result, the shear necessary to split the vortex is so strong that the vortex can actually no longer be considered to be coherent: the nonlinear term in Eq. (2a) becomes smaller than the β terms and can be neglected. The dynamics of the vortex then boil down to the evolution of a Rossby wave packet with wavelengths on the order of the vortex radius and smaller. Because Rossby waves are very dispersive, the vortex is scattered before it can be split. Because the baroclinic β effect varies vertically, the vertical structure of each wave is more complicated than on the planetary β plane (for which it boils down to the usual baroclinic modes associated with the stratification): it depends on the wavelength, the vortex, and the background current vertical structure (Dewar 1998).

Earlier studies have shown that a vortex is dispersed into Rossby waves if the change of background PV over the vortex core is nonnegligible in comparison with the vortex PV anomaly (Flierl 1977; McWilliams and Flierl 1979; Thierry and Morel 1999). In our study, the former is associated with the background current and is measured by βU1r1, which gives a measure of dispersion. When βU1r1 ≥ 1, the vortex represents only a perturbation of the background PV. There are no closed PV contours, and this perturbation is scattered by Rossby waves. When βU1r1 is small, the vortex has closed PV contours and is strong, and fast fluid rotation prevents the vortex from being dispersed.

To evaluate whether the vortex is scattered by Rossby waves before it is split by the current shear, we calculated βU1r1 using the critical shear for splitting given by the solid line in Fig. 6. The results are shown in Fig. 7, which represents ε = βUc1r1 (where βUc1 = F+1Ucrit) as a function of the vortex radius. Therefore, ε measures the vortex dispersion when it is close to splitting. The value of ε above which the vortex is dispersed by Rossby waves is not well defined, but notice that ε is no longer small when r1 ≥ 1 or so. Figure 8 shows the total potential vorticity evolution in layer 1 (Fig. 8a) and layer 2 (Fig. 8b) of a vortex with radius r1 = 1.5 and with a background current U = 0.75Ucrit ≃ 0.15. The dispersion coefficient in this case is ε ≃ 0.25. Notice how the vortex is distorted and dispersed by Rossby waves even if it stays aligned. Also notice that the propagation is weak, which illustrates the compensation of the background current advection by the baroclinic β effect.

Fig. 7.

Measurement of the background PV field dispersion effect ε as a function of r1 and for a background current that corresponds to the critical shear given by the solid line in Fig. 6. Notice that when r1 ≥ 1 the dispersion effect can no longer be considered weak so that, for such shear, the vortex is scattered

Fig. 7.

Measurement of the background PV field dispersion effect ε as a function of r1 and for a background current that corresponds to the critical shear given by the solid line in Fig. 6. Notice that when r1 ≥ 1 the dispersion effect can no longer be considered weak so that, for such shear, the vortex is scattered

Fig. 8.

Evolution of the total PV (vortex + background current) in (a) layer 1 and (b) layer 2. Times shown are t = 0, 20, 40, and 60. The vortex radius is r1 = 1.5, and the background current is (U1, U2) = (0.75Ucrit, 0) = (0.15, 0). Note how the vortex is dispersed by Rossby waves

Fig. 8.

Evolution of the total PV (vortex + background current) in (a) layer 1 and (b) layer 2. Times shown are t = 0, 20, 40, and 60. The vortex radius is r1 = 1.5, and the background current is (U1, U2) = (0.75Ucrit, 0) = (0.15, 0). Note how the vortex is dispersed by Rossby waves

The baroclinic β effect therefore drastically changes the physics of the interaction between a vortex and a background current. It strongly increases the critical shear for splitting, and, for vortices with radii larger than the internal radius of deformation (r1 ≥ 1), splitting becomes impossible. This result is summarized in Fig. 9, which represents the behavior of the vortex as a function of its radius and background shear. Three regions exist in which the vortex either remains coherent, is split, or is scattered. The solid line is the same as in Fig. 6 and represents the critical shear beyond which the vortex is split. It has been calculated from the analytical model in Eq. (3). We did not define precisely when scattering seriously alters the vortex, but from numerical solutions we can estimate that scattering becomes important when ε ≥ 0.2–0.3. The hatched zone corresponds to ε ∈ [0.2, 0.3] and is associated with this fuzzy frontier between scattering and coherence. Note that, for a given vortex strength, large-scale vortices are more easily destroyed than are vortices with moderate radius (r1 ≃ 1), because they are more sensitive to dispersion. A similar result was obtained by Thierry and Morel (1999) for the scattering of vortices by topographic Rossby waves.

Fig. 9.

Vortex evolution diagram as a function of r1 and the background shear U. When the baroclinic β effect is taken into account, the dispersion effect of the current has to be taken into account. As a result, vortices with radii above the first internal radius of deformation are dispersed before they can be split. The frontier between coherence and dispersion roughly corresponds to the line ε = 0.2–0.3 but is not well defined. The hatched region accounts for this fuzziness, and numerical experiments show that below the hatched region the PV patches stay coherent (at least for the time period considered here: 100 nondimensional time units), whereas above they are dispersed

Fig. 9.

Vortex evolution diagram as a function of r1 and the background shear U. When the baroclinic β effect is taken into account, the dispersion effect of the current has to be taken into account. As a result, vortices with radii above the first internal radius of deformation are dispersed before they can be split. The frontier between coherence and dispersion roughly corresponds to the line ε = 0.2–0.3 but is not well defined. The hatched region accounts for this fuzziness, and numerical experiments show that below the hatched region the PV patches stay coherent (at least for the time period considered here: 100 nondimensional time units), whereas above they are dispersed

c. Influence of the middle-layer background PV gradient

As already mentioned, in our configuration the background PV gradients in the upper and middle layers are opposite, and the configuration is unstable to baroclinic perturbations. In this section, we consider the additional effect of the planetary β or/and of a current in the middle layer, and we study how it affects the vortex splitting. These additional effects can be easily taken into account in the analytical model: as shown in the appendix [see Eq. (A7)] the equation for the general βterm is simply the sum of the planetary and baroclinic β effects.

First consider the influence of the planetary β. For eastward (positive) flows in the upper layer, βU is positive (negative) in the upper (middle) layer so that this additional effect strengthens (weakens) the background PV gradient [see Eqs. (1a)–(1b)]. Our investigations show that the planetary β favors destruction of the vortex both by splitting and dispersion (Vandermeirsch 1999). As one can expect, the addition of a planetary β effect increases the dispersive properties of the background current: a weaker current is now necessary to scatter the vortex. For instance, when βp = 0.1 and r1 = 1, we found that the critical shear needed to split the vortex is reduced from 0.14 to 0.11. The vortex also becomes more sensitive to splitting when r1 ≤ 1. This can be explained by the fact that the rotational advections induced by the vortex are different in each layer, and the corresponding planetary β gyres are different, too. This induces different propagation speed for PV patches, resulting in increasing the vortex separation η.

For the configuration with an eastward flow in the upper layer, stability of the background flow is achieved when the PV gradient in the middle layer is zero or becomes positive as in the upper layer, β2 = βp + βU2 ≥ 0, which limits the background current strength. For instance, when βp = 2 × 10−11 m−1 s−1 and with the stratification considered in this paper, stable background currents should have shear equal to or less than 2 cm s−1. Adding an eastward flow in the middle layer, however, permits the consideration of stronger shears. Indeed, as shown by Eqs. (1a)–(1b) and (2d), even without the planetary β, βU1 and βU2 have the same sign if 1 ≤ U1/U2 ≤ 1 + F+2/F2. In this case, the background current is baroclinically stable. For the stratification considered in this paper, the flow is thus stable if the velocity in the middle layer is one-third of the upper-layer velocity (U2U1/3). The middle-layer thickness is then constant, and there is no baroclinic β effect in this layer. In this case, the results described above—in particular, the compensation of the current-tearing properties by the baroclinic β effect in the upper layer—remain generally valid (Vandermeirsch 1999).

6. Discussion and conclusions

In this paper, we studied the resistance of a coherent vortex to a vertically sheared current. Numerical and analytical models are compared to gain insight into the mechanisms influencing vortex splitting.

We first derived a splitting criterion and validated it in a model in which the baroclinic β effect was artificially neglected. We showed that large tilting can be achieved before the vortex actually splits. We also showed that the vortex deformations play a role in the splitting process when the vortex radius is larger than the Rossby radius of deformation but that the effect remains modest.

We then took the influence of baroclinic β effect into account to revisit MP93 preliminary investigations. In their numerical solutions with nonuniform ambient potential vorticity, MP93 did notice some influence of the baroclinic β effect on the behavior of vortices, but this process was only briefly discussed (see the last paragraph of their section 6), and they advocated for a more detailed study on this subject. Our results show that the baroclinic β effect plays a major role in the dynamics of the interaction between a vortex and a vertically sheared current, whenever it is not zero within the vortex core. The effect is indeed associated with the development of a secondary flow generated by the background current PV gradient, which drastically reduces the effectiveness of tearing by the current. Thus, relatively, strong shears are necessary to split the vortex, and with such shears the vortex is seriously affected by the Rossby wave dispersion. As a result, the vortex is generally dispersed before it can be split by the current, at least when its radius is above the first internal radius of deformation, which is the case of many coherent vortices in the ocean. The planetary β effect is shown to favor splitting and dispersion for an eastward current in the upper layer.

In practice, in the ocean, most large-scale currents have a moderate vertical shear that is below the critical shear for the splitting or dispersion of newly formed coherent vortices. For instance, in the case of Gulf Stream rings, the maximum velocity of the ring is V ≃ 0.5 m s−1 and is reached at a radius Rm ≃ 50 km (Olson 1991). The radius of deformation in this region is Rd ≃ 25 km so that we can estimate r1 ≃ 2. The dimensional potential vorticity anomaly associated with such vortices is Q1 = V/RdVmax (r1 = 2), where Vmax = 0.5 as estimated using Fig. 2a. This yields Q1 ≃ 4 × 10−5 s−1. Figure 9 shows that such vortices can only be destroyed by scattering effects of Rossby waves when UU2U1 ≃ 0.1Rd1Q1. This corresponds to a current shear of U ≃ 10 cm s−1. For a smaller vortex with a similar PV but a smaller radius Rm ≃ Rd, the critical background shear would be even larger: U ≥ 0.14Rd1Q1 ≃ 14 cm s−1. Such intense shears can only be found in intense baroclinic jets.

Vortices with very small radii, such as eddies that are associated with convective plumes and have dimensional radii R1 ≃ 1–2 km (Gascard 1978, 1991), have a different behavior. Their nondimensional radii roughly correspond to 0.1–0.2 Rd1 (if we consider Rd1 ≃ 10 km as in the Mediterranean Sea), and their potential vorticity anomaly is mainly associated with their relative vorticity, which is close to |Q1| ≃ f0 = 10−4 s−1 (Legg and Marshall 1993). Because of their small size, the baroclinic β effect does not play a significant role (notice that the plain and dashed lines in Fig. 6 are close when r1 is small), so that the results for βUk ≡ 0 can be applied. Our results predict a minimum dimensional shear of U = 0.0052Rd1|Q1| ≃ 0.5 cm s−1 (with Rd1 ≃ 10 km), which can be easily achieved. In fact, even though their PV anomaly is fairly strong, the coherence of these vortices is very low (their maximum speed is proportional to their radius and is therefore small), and we therefore do not expect these small-scale vortices to resist even modest vertical shears.

This study thus gives some insights on the life of mesoscale oceanic vortices. The way these structures are generated and propagate has been extensively studied, but not much is known in regard to their erosion. This is, however, important, because coherent vortices play an important role in the thermohaline circulation, and they release their heat and salt content through dissipative processes. We have shown that, during the early stages of their life, they are insensitive to the tearing by large-scale baroclinic currents (VMS01 also showed that the latter have weak influence on their propagation). Unless they interact with intense jets (such as the Gulf Stream, the Azores current, etc.), or large bottom topography (see, e.g., Richardson et al. 1989), their erosion is thus initially slow and is likely to be associated with dispersion by Rossby waves. Their radius and strength, however, decrease with time (see, e.g., Armi et al. 1989) and, after a few years, when they have become weak enough, our analysis shows that their remaining heat and salt anomalies are rapidly destroyed by surrounding currents. The distance over which they can transport tracers is, however, very large, in particular owing to the compensation of the background current tearing effect by the baroclinic β effect, which increases their lifetime.

Note that in this paper we focused on a particular 2½-layer configuration with ambient PV gradients of opposite sign within the vortex core and weaker vortex rotation in the middle layer than in the upper layer. There do exist configurations in which the background baroclinic current may have no PV gradients within some part of the vortex core (e.g., βU2 = 0 as discussed in section 5c) or even within the whole vortex core (see Fig. 10). In these configurations, the reduction of the tearing effectiveness within the vortex core due to the baroclinic β effect depends on the stratification and vortex structure. For instance, if the vortex circulation is weak and can be considered to be inactive in the layers above and below the vortex core, the β gyres in those layers will also be weak and the compensation of the background current tearing properties will not be effective. To our knowledge, these special configurations can only be obtained if the layers above and below the vortex core are very deep or if very strong background PV exists in these layers, scattering the motion induced by the vortex. Chassignet and Cushman-Roisin (1991) have shown that increasing the lower-layer depth to obtain negligible vortex circulation leads to unrealistic situations. On the other hand, strong PV gradients in some layers are possible. Thierry and Morel (1999) have indeed shown that a steep bottom slope (or small-scale topography with steep slopes) can rapidly scatter the vortex signature in the lower layer so that it can be considered to be at rest (see Fig. 10b).

Fig. 10.

Configurations in which the current baroclinic β effect is negligible when interacting with a vortex. The background current shear must be constant within the vortex core (if the stratification is constant), and the layers above and below must be associated with strong PV gradients. This is possible for (a) strong current shear or (b) strong current shear in the upper layer and steep bottom topography in the lower layer

Fig. 10.

Configurations in which the current baroclinic β effect is negligible when interacting with a vortex. The background current shear must be constant within the vortex core (if the stratification is constant), and the layers above and below must be associated with strong PV gradients. This is possible for (a) strong current shear or (b) strong current shear in the upper layer and steep bottom topography in the lower layer

Last, the analytical and numerical models we have used are based on the quasigeostrophic equations, and comparison of the results with oceanic vortices is thus questionable because the latter often have large PV anomalies and Rossby numbers (e.g., Sutyrin 1989). For strong vortices, quantitative differences can be expected, but we believe the physics and the general results summarized below are still valid when considering more general primitive equation models.

Acknowledgments

We thank Drs. William Dewar and Glenn Flierl, who generously provided their quasigeostrophic pseudospectral model. Frédéric Vandermeirsch and Yves Morel are supported by the Service Hydrographique et Océanographique de la Marine (SHOM: French Navy Hydrographic and Oceanographic Service). Georgi Sutyrin is supported by the U.S. Office of Naval Research (Grant N000149710139) and the National Science Foundation (Grant ATM-9905209). This paper is a contribution to the SEMAPHORE and SEMANE research programs conducted by SHOM (Programme d'Etudes Amonts 982401).

REFERENCES

REFERENCES
Abramowitz
,
M.
, and
I. A.
Stegun
,
Eds.,
.
1970
:
Handbook of Mathematical Functions, with Formulas, Graphs, and Mathematical Tables.
Applied Mathematics Series, Vol. 55, National Bureau of Standards, 1046 pp
.
Armi
,
L.
,
D.
Hebert
,
N.
Oakey
,
J. F.
Price
,
P. L.
Richardson
,
H. T.
Rossby
, and
B.
Ruddick
,
1989
:
Two years in the life of a Mediterranean salt lens.
J. Phys. Oceanogr.
,
19
,
354
370
.
Chassignet
,
E.
, and
B.
Cushman-Roisin
,
1991
:
On the influence of a lower layer on the propagation of nonlinear oceanic eddies.
J. Phys. Oceanogr.
,
21
,
939
957
.
Chérubin
,
L.
, and
Coauthors
.
1997
:
Descriptive analysis of the hydrology and currents on the Iberian shelf from Gibraltar to Cape Finisterre: Preliminary results from the SEMANE and INTERAFOS experiments.
Ann. Hydrogr.
,
21
,
5
69
.
Dewar
,
W. K.
,
1998
:
On “too fast” baroclinic planetary waves in the general circulation.
J. Phys. Oceanogr.
,
28
,
1739
1758
.
Dewar
,
W. K.
, and
G.
Flierl
,
1987
:
Some effects of the wind on rings.
J. Phys. Oceanogr.
,
17
,
1653
1667
.
Flierl
,
G.
,
1977
:
The application of linear quasi-geostrophic dynamics to Gulf Stream rings.
J. Phys. Oceanogr.
,
7
,
365
379
.
Gascard
,
J-C.
,
1978
:
Mediterranean deep water formation, baroclinic instability and oceanic eddies.
Oceanol. Acta
,
1
,
315
330
.
Gascard
,
J-C.
,
1991
:
Open ocean convection and deep water formation revisited in the Mediterranean, Labrador, Greenland, and Weddell seas.
Proceedings of the Workshop on Deep Convection and Deep Water Formation in the Ocean, P. C. Chu and J. C. Gascard, Eds., Elsevier Science, 157–182
.
Hogg
,
N.
, and
H.
Stommel
,
1990
:
How currents in the upper thermocline could advect meddies deeper down.
Deep-Sea Res.
,
37
,
613
623
.
Joyce
,
T. M.
, and
T. J.
McDougall
,
1992
:
Physical structure and temporal evolution of Gulf Stream warm-core ring 82B.
Deep-Sea Res.
,
39
,
(Suppl.),
.
19
44
.
Legg
,
S.
, and
J.
Marshall
,
1993
:
A heton model of the spreading phase of open-ocean deep convection.
J. Phys. Oceanogr.
,
23
,
1040
1056
.
Marshall
,
J. S.
, and
B.
Parthasarathy
,
1993
:
Tearing of an aligned vortex by a current difference in two-layer quasi-geostrophic flow.
J. Fluid Mech.
,
255
,
157
182
.
McWilliams
,
J.
,
1985
:
Submesoscale, coherent vortices in the ocean.
Rev. Geophys.
,
23
,
165
182
.
McWilliams
,
J.
, and
G.
Flierl
,
1979
:
On the evolution of isolated, nonlinear vortices.
J. Phys. Oceanogr.
,
9
,
1155
1182
.
Meacham
,
S. P.
,
G. R.
Flierl
, and
U.
Send
,
1990
:
Vortices in shear.
Dyn. Atmos. Oceans
,
14
,
333
386
.
Meacham
,
S. P.
,
K. K.
Pankratov
,
A. S.
Shchepetkin
, and
V. V.
Zhmur
,
1994
:
The interaction of ellipsoidal vortices in background shear flows in a stratified fluid.
Dyn. Atmos. Oceans
,
21
,
167
212
.
Morel
,
Y.
, and
J. C.
McWilliams
,
1997
:
Evolution of isolated interior vortices in the ocean.
J. Phys. Oceanogr.
,
27
,
727
748
.
Olson
,
D. B.
,
1980
:
The physical oceanography of two rings observed by the cyclonic ring experiment. Part II: Dynamics.
J. Phys. Oceanogr.
,
10
,
514
528
.
Olson
,
D. B.
,
1991
:
Rings in the ocean.
Annu. Rev. Earth Planet. Sci.
,
19
,
283
311
.
Olson
,
D. B.
, and
R. H.
Evans
,
1986
:
Rings of the Agulhas.
Deep-Sea Res.
,
33
,
27
42
.
Paillet
,
J.
,
B.
Le Cann
,
A.
Serpette
,
Y.
Morel
, and
X.
Carton
,
1999
:
Real-time tracking of a Galician meddy.
Geophys. Res. Lett.
,
26
,
1877
1880
.
Pedlosky
,
J.
,
1987
:
Geophysical Fluid Dynamics.
Springer-Verlag, 710 pp
.
Pingree
,
R. D.
, and
B.
Le Cann
,
1992
:
Anticyclonic eddy X91 in the southern Bay of Biscay, May 1991 to February 1992.
J. Geophys. Res.
,
97
,
14353
14367
.
Pingree
,
R. D.
, and
B.
Le Cann
,
1993a
:
A shallow meddy (a smeddy) from the secondary Mediterranean salinity maximum.
J. Geophys. Res.
,
98
,
20169
20185
.
Pingree
,
R. D.
, and
B.
Le Cann
,
1993b
:
Structure of a meddy (Bobby 92) southeast of the Azores.
Deep-Sea Res.
,
40
,
2077
2103
.
Pingree
,
R. D.
, and
B.
Le Cann
,
1994
:
Winter warming in the southern Bay of Biscay and Lagrangian eddy kinematics from a deep-drogued Argos buoy.
J. Mar. Biol.
,
74
,
107
128
.
Reznik
,
G. M.
, and
W. K.
Dewar
,
1994
:
An analytic theory of distributed axisymmetric barotropic vortices on the β-plane.
J. Fluid Mech.
,
269
,
301
321
.
Reznik
,
G. M.
,
R.
Grimshaw
, and
E. S.
Benilov
,
2000
:
On the long-term evolution of an intense localized divergent vortex on the beta-plane.
J. Fluid Mech.
,
422
,
249
280
.
Richardson
,
P. L.
,
D.
Walsh
,
L.
Armi
,
M.
Schröder
, and
J. F.
Price
,
1989
:
Tracking three meddies with SOFAR floats.
J. Phys. Oceanogr.
,
19
,
371
383
.
Sutyrin
,
G. G.
,
1987
:
The beta-effect and the evolution of a localized vortex.
Sov. Phys. Dokl.
,
32
,
791
793
.
Sutyrin
,
G. G.
,
1989
:
The structure of a baroclinic eddy.
Oceanology
,
29
,
139
144
.
Sutyrin
,
G. G.
, and
G. R.
Flierl
,
1994
:
Intense vortex motion on the beta-plane: Development of the beta gyres.
J. Atmos. Sci.
,
51
,
773
790
.
Sutyrin
,
G. G.
, and
Y. G.
Morel
,
1997
:
Intense vortex motion in a stratified fluid on the beta-plane. An analytical model and its validation.
J. Fluid Mech.
,
336
,
203
220
.
Thierry
,
V.
, and
Y.
Morel
,
1999
:
Influence of a strong bottom slope on the evolution of a surface intensified vortex.
J. Phys. Oceanogr.
,
29
,
911
924
.
van Ballegooyen
,
R. C.
,
M. L.
Grüdlingh
, and
J. R. E.
Lutjeharms
,
1994
:
Eddy fluxes of heat and salt from the southwest Indian Ocean into the southeast Atlantic Ocean: A case study.
J. Geophys. Res.
,
99
,
14053
14070
.
Vandermeirsch
,
F.
,
1999
:
Interaction entre un jet et un tourbillon océanique.
Application au courant des Açores et aux meddies de la campagne SEMAPHORE. Ph.D. thesis, Université de Bretagne Occidentale, 210 pp. [Available from F. Vandermeirsch, IFREMER, B.P. 70, 29280 Plouzane, France.]
.
Vandermeirsch
,
F.
,
Y.
Morel
, and
G.
Sutyrin
,
2001
:
The net advective effect of a vertically sheared current on a coherent vortex.
J. Phys. Oceanogr.
,
31
,
2210
2225
.
Walsh
,
D.
,
1995
:
A model of a mesoscale lens in large scale shear. Part I: Linear calculations.
J. Phys. Oceanogr.
,
25
,
735
746
.
Walsh
,
D.
, and
L. J.
Pratt
,
1995
:
The interaction of a pair of point vortices in uniform shear.
Dyn. Atmos. Oceans
,
22
,
135
160
.
Walsh
,
D.
,
P. L.
Richardson
, and
J.
Lynch
,
1996
:
Observations of tilting meddies.
J. Phys. Oceanogr.
,
26
,
1023
1038
.

APPENDIX

Analytical Models

Following HS90 or MP93, we assume that both PV patches remain circular and we first neglect the baroclinic β effect. In that case the PV patch centers (x1, y1) and (x2, y2) verify

 
formula

where (x1, y1) and (x2, y2) are the position of the centers in layer 1 and 2, respectively; r = [(x2x1)2 + (y2y1)2]1/2 is the distance between the centers; and V2→1(V1→2) is the velocity induced by the second (first) layer PV patch in the first (second) layer.

This yields for the vortex separation η = (x, y) = (x2x1, y2y1):

 
formula

or in complex form, with η = x + iy,

 
formula

where

 
formula

[see Eq. (A9) below for an algebraic expression of Ω].

As described in VMS01, the β-gyre development is important, too, and can drastically modify the vortex evolution. An additional term thus has to be considered in Eq. (A3), which becomes

 
/dt = Ω(|η|) + U + βterm,
(A5)

where βterm accounts for the baroclinic β effect.

VMS01 developed an analytical model in which the latter is taken into account. The PV patch separation evolution is calculated, too [see Eqs. (A16)–(A18) in their appendix A] and is similar to Eq. (A5). It can indeed be written

 
/dt = Ωo + U + βterm(t),
(A6)

where Ωo = Ω(r1) is now constant and βterm(t) is a time function independent of η and has been calculated in VMS01:

 
formula

where Ωk represents the rotation rate of the initially axisymmetric vortex in layer k and is written

 
formula

Here, Ω(r) represents the sum of the rotation rate induced by a PV patch in the adjacent layer and is thus given by

 
formula

In these formulas, P(n) = [P(n)1, … , P(n)k, … , P(n)N] is the nth vertical eigenmode associated with the stretching matrix 𝗙𝗿; −γ2n is its corresponding eigenvalue. The matrix α with coefficients α(n)l is the inverse of the matrix 𝗣 whose columns are the vectors P(n). Here, G(n)1 is the Green function associated with the Helmholtz operator [rr(rr) − 1/r2γ2n] and is expressed in terms of modified Bessel functions K1 and I1 (Abramowitz and Stegun 1970, 231–233). When γn ≠ 0,

 
formula

If there exists a barotropic mode (e.g., if F+2 = 0) with eigenvalue γ0 = 0, G(0)1(r|r′) becomes

 
formula

Equations (A6)–(A7) are linear, and η(t) can be calculated explicitly. This model gives good results over a few tens of vortex turnover time; however, strictly speaking, the model is limited because it is based on a linearization of the evolution equations. Its limits are discussed in VMS01, but let us underline that it is only precise when η is small in comparison with the vortex radius. Because this hypothesis is not valid when the structure splits, in this study we expect nonlinearity of Ω to be important. To achieve quantitatively good results, the evolution equation we consider is thus Eq. (A5), where βterm is given by Eq. (A7).

To clarify the calculation of all coefficients appearing in the previous formulas, we give their explicit expressions in the case of the 2½-layer system considered in this paper. The matrices we have defined above, and their corresponding elements, are therefore given by

 
formula

Footnotes

Corresponding author address: Yves Morel, EPSHOM/CMO, B.P. 426, Brest, Cedex 29275, France. Email: morel@shom.fr

Notice, however, that for vortices with small radii (rk ≃ 0.1 or so) the time period has to be increased because the background shear necessary to split the vortex is low and the time period necessary to get a significant separation is thus large. However, for vortices with “reasonable” sizes (rk ≥ 0.5 or so), t ∈ [0, 100] yields good results.

Notice that the validity of Eq. (4) has been verified using both of these numerical and analytical solutions.