Thermal motions in supercoiled DNA are studied by Brownian dynamics (BD) simulations with a focus on the site juxtaposition process. It had been shown in the last decade that the BD approach is capable of describing actual times of large-scale DNA motion. The bead model of DNA used here accounts for bending and torsional elasticity as well as the electrostatic repulsion among DNA segments. The hydrodynamic interaction among the beads of the model chain and the aqueous solution is incorporated through the Rotne-Prager tensor. All simulations were performed for the sodium ion concentration of 0.01 M. We first showed, to test our BD procedure, that the same distributions of equilibrium conformational properties are obtained as by Monte Carlo simulations for the corresponding DNA model. The BD simulations also predict with accuracy published experimental values of the diffusion coefficients of supercoiled DNA. To describe the rate of conformational changes, we also calculated the autocorrelation functions for the writhe and radius of gyration for the supercoiled molecules. The rate of site juxtaposition was then studied for DNA molecules up to 3000 bp in length. We find that site juxtaposition is a very slow process: although accelerated by a factor of more than 100 by DNA supercoiling, the times of juxtaposition are in the range of ms even for highly supercoiled DNA, about two orders of magnitude higher than the relaxation times of writhe and the radius of gyration for the same molecules. By inspecting successive simulated conformations of supercoiled DNA, we conclude that slithering of opposing segments of the interwound superhelix is not an efficient mechanism to accomplish site juxtaposition, at least for conditions of low salt concentration. Instead, transient distortions of the interwound superhelix, followed by continuous reshaping of the molecule, contribute more significantly to site juxtaposition kinetics.