This paper investigates the three-dimensional stability of a Lamb–Chaplygin columnar vertical vortex pair as a function of the vertical wavenumber kz, horizontal Froude number Fh, Reynolds number Re and Schmidt number Sc. The horizontal Froude number Fh (Fh = U/NR, where U is the dipole travelling velocity, R the dipole radius and N the Brunt–Väisälä frequency) is varied in the range [0.033, ∞[ and three set of Reynolds-Schmidt numbers are investigated: {Re = 10 000, Sc = 1}, Re = 1000, Sc = 1}, {Re = 200, Sc = 637}. In the whole range of Fh and Re, the dominant mode is always antisymmetric with respect to the middle plane between the vortices but its physical nature and properties change when Fh is varied. An elliptic instability prevails for Fh > 0.25, independently of the Reynolds number. It manifests itself by the bending of each vortex core in the opposite direction to the vortex periphery. The growth rate of the elliptic instability is reduced by stratification effects but its spatial structure is almost unaffected. In the range 0.2 < Fh < 0.25, a continuous transition occurs from the elliptic instability to a different instability called zigzag instability. The transitional range Fhc = 0.2–0.25 is in good agreement with the value Fh = 0.22 at which the elliptic instability of an infinite uniform vortex is suppressed by the stratification. The zigzag instability dominates for Fh [les ] 0.2 and corresponds to a vertically modulated bending and twisting of the whole vortex pair. The experimental evidence for this zigzag instability in a strongly stratified fluid reported in the first part of this study (Billant & Chomaz 2000a) are therefore confirmed and extended. The numerically calculated wavelength and growth rate for low Reynolds number compare well with experimental measurements.

The present numerical stability analysis fully agrees with the inviscid asymptotic analysis carried out in the second part of this investigation (Billant & Chomaz 2000b) for small Froude number Fh and long wavelength. This confirms that the zigzag instability is related to the breaking of translational and rotational invariances. As predicted, the growth rate of the zigzag instability is observed to be self-similar with respect to the variable Fhkz, implying that the maximum growth rate is independent of Fh while the most amplified dimensional wavenumber varies with N/U. The numerically computed eigenmode and dispersion relation are in striking agreement with the analytical results.