  Numerical relativity is an essential tool in studying the coalescence of
binary black holes (BBHs). It is still computationally prohibitive to cover the
BBH parameter space exhaustively, making phenomenological fitting formulas for
BBH waveforms and final-state properties important for practical applications.
We describe a general hierarchical bottom-up fitting methodology to design and
calibrate fits to numerical relativity simulations for the three-dimensional
parameter space of quasicircular nonprecessing merging BBHs, spanned by mass
ratio and by the individual spin components orthogonal to the orbital plane.
Particular attention is paid to incorporating the extreme-mass-ratio limit and
to the subdominant unequal-spin effects. As an illustration of the method, we
provide two applications, to the final spin and final mass (or equivalently:
radiated energy) of the remnant black hole. Fitting to 427 numerical relativity
simulations, we obtain results broadly consistent with previously published
fits, but improving in overall accuracy and particularly in the approach to
extremal limits and for unequal-spin configurations. We also discuss the
importance of data quality studies when combining simulations from diverse
sources, how detailed error budgets will be necessary for further improvements
of these already highly accurate fits, and how this first detailed study of
unequal-spin effects helps in choosing the most informative parameters for
future numerical relativity runs.
