Einstein proposed his famous theory of gravity, general relativity, in 1915, but not until a century later have physicists been able to directly observe the long-predicted gravitational waves. To date, there have been about 90 detection events of gravitational waves sourced by the most energetic phenomena in the Universe, the mergers of two black holes. The success of gravitational-wave astronomy is inseparable from the advances in numerical relativity. Nowadays, numerical relativists are able to simulate the spacetimes of binary-black-hole mergers routinely and extract accurate waveforms from these simulations. Meanwhile, there are still numerous open questions in this field awaiting answers. In this thesis, we present our solutions to three important issues in simulating binary black holes. First, we set about improving efficiency of simulations. Because of the significant computational expense of high-spin simulations, such simulations are sparsely scattered in current public waveform catalogs. However, waveform modeling for detection analysis requires many more data points in this regime to accurately capture the behaviors of high-spin binary-black-hole systems. We develop new numerical techniques that can reduce the cost of high-spin simulations by a factor of 2 compared to those using traditional methods. These new techniques can certainly accelerate the expansion of the parameter space into the high-spin regime. Second, we concern ourselves with gravitational-wave analysis. For detection purposes, gravitational waves are usually analyzed in the frequency domain rather than in the time domain. Because a general numerical waveform has different values at both ends, it must be preprocessed in the time domain before being Fourier transformed. Otherwise, the value mismatch would lead to spurious contents that contaminate the frequency spectrum. The common scheme of preprocessing in literature has been very successful in producing frequency spectra of those waveforms without gravitational-wave memory, but it fails to render the expected low-frequency spectrum of a waveform with memory. To solve this issue, we propose a new preprocessing scheme that produces robust and accurate spectra of memory waveforms. As two applications of this new scheme, we carefully inspect the detailed structures of the spectrum of a memory waveform and briefly survey the detectability of the memory effects in both current- and next-generation detectors. Third, we shift our attention to the strongest-field region in a simulation, the black-hole horizons. There has been evidence in the literature that the multipole moments on the common horizon, a set of time-dependent values that quantifies the horizon shape, are described by Kerr perturbation theory. This is a surprising result, since the common horizon forms as the black holes merge. We would expect this to correspond to a highly non-linear phase of the merger, unlikely to be describable by linear perturbation theory. However, such evidence is based on either a head-on collision of two black holes or a definition of multipole moments that disregards the connection among quasilocal horizons on different time slices. In contrast, we construct a set of covariantly defined multipole moments on the common horizon formed by the merger of two orbiting black holes. We find that these multipole moments are indeed described by the fundamental quasinormal modes at late times, similar to the ringdown gravitational waves. By including overtones, we also find an excellent quasinormal description of the dominant multipole moment at all times after the merger. This signifies the perhaps remarkable capability of the Kerr perturbation theory.