The viscous evolution of two co-rotating vortices is analysed using direct two-dimensional numerical simulations of the Navier–Stokes equations. The article focuses on vortex interaction regimes before merging. Two parameters are varied: a steepness parameter n which measures the steepness of the initial vorticity profiles in a given family of profiles, and the Reynolds number Re (between 500 and 16 000). Two distinct relaxation processes are identified. The first one is non-viscous and corresponds to a rapid adaptation of each vortex to the external (strain) field generated by the other vortex. This adaptation process, which is profile dependent, is described and explained using the damped Kelvin modes of each vortex. The second relaxation process is a slow diffusion phenomenon. It is similar to the relaxation of any non-Gaussian axisymmetrical vortex towards the Gaussian. The quasi-stationary solution evolves on a viscous-time scale toward a single attractive solution which corresponds to the evolution from two initially Gaussian vortices. The attractive solution is analysed in detail up to the merging threshold a/b ≈ 0.22 where a and b are the vortex radius and the separation distance respectively. The vortex core deformations are quantified and compared to those induced by a single vortex in a rotating strain field. A good agreement with the asymptotic predictions is demonstrated for the eccentricity of vortex core streamlines. A weak anomalous Reynolds number dependence of the solution is also identified. This dependence is attributed to the advection–diffusion of vorticity towards the hyperbolic points of the system and across the separatrix connecting these points. A Re1/3 scaling for the vorticity at the central hyperbolic point is obtained. These findings are discussed in the context of a vortex merging criterion.