The dynamics of a counter-rotating pair of columnar vortices aligned parallel to a stable density gradient are investigated. By means of numerical simulation, we extend the linear analyses and laboratory experiments of Billant & Chomaz (J. Fluid Mech. vol. 418, p. 167; vol. 419, pp. 29, 65 (2000)) to the fully nonlinear, large-Reynolds-number regime. A range of stratifications and vertical length scales is considered, withFrh< 0.2 and 0.1 <Frz< 10. HereFrh≡U/(NR) andFrz≡Ukz/Nare the horizontal and vertical Froude numbers,UandRare the horizontal velocity and length scales of the vortices,Nis the Brunt–Väisälä frequency, and 2π/kzis the vertical wavelength of a small initial perturbation. At early times withFrz< 1, linear predictions for the zigzag instability are reproduced. Short-wavelength perturbations withFrz> 1 are found to be unstable as well, with growth rates only slightly less than those of the zigzag instability but with very different structure. At later times, the large-Reynolds-number evolution diverges profoundly from the moderate-Reynolds-number laboratory experiments as the instabilities transition to turbulence. For the zigzag instability, this transition occurs when density perturbations generated by the vortex bending become gravitationally unstable. The resulting turbulence rapidly destroys the vortex pair. We derive the criterion η/R≈ 0.2/Frzfor the onset of gravitational instability, where η is the maximum horizontal displacement of the bent vortices, and refine it to account for a finite twisting disturbance. Our simulations agree for the fastest growing wavelengths 0.3 <Frz< 0.8. Short perturbations withFrz> 1 saturate at low amplitude, preserving the columnar structure of the vortices well after the generation of turbulence. Viscosity is shown to suppress the transition to turbulence for Reynolds numberRe≲ 80/Frh, yielding laminar dynamics and, under certain conditions, pancake vortices like those observed in the laboratory.