  The 2017 detection of a kilonova coincident with gravitational-wave emission
has identified neutron star mergers as the major source of the heaviest
elements, and dramatically constrained alternative theories of gravity.
Observing a population of such sources has the potential to transform
cosmology, nuclear physics, and astrophysics. However, with only one confident
multi-messenger detection currently available, modelling the diversity of
signals expected from such a population requires improved theoretical
understanding. In particular, models which are quick to evaluate, and are
calibrated with more detailed multi-physics simulations, are needed to design
observational strategies for kilonovae detection, and to obtain rapid-response
interpretations of new observations. We use grey-opacity models to construct
populations of kilonovae, spanning ejecta parameters predicted by numerical
simulations. Our modelling focuses on wavelengths relevant for upcoming optical
surveys, such as the Rubin Observatory Legacy Survey of Space and Time (LSST).
In these simulations, we implement heating rates that are based on nuclear
reaction network calculations. We create a Gaussian-process emulator for
kilonova grey opacities, calibrated with detailed radiative transfer
simulations. Using recent fits to numerical relativity simulations, we predict
how the ejecta parameters from BNS mergers shape the population of kilonovae,
accounting for the viewing-angle dependence. Our simulated population of binary
neutron star (BNS) mergers produce peak i-band absolute magnitudes $-20 \leq
M_i \leq -11$. A comparison with detailed radiative transfer calculations
indicates that further improvements are needed to accurately reproduce spectral
shapes over the full light curve evolution.
