The efficient simulation of multijet final states presents a serious computational task for analyses of LHC data and will be even more so at the HL-LHC. We here discuss means to accelerate the generation of unweighted events based on a two-stage rejection-sampling algorithm that employs neural-network surrogates for unweighting the hard-process matrix elements. To this end, we generalise the previously proposed algorithm based on factorisation-aware neural networks to the case of multijet merging at tree-level accuracy. We thereby account for several non-trivial aspects of realistic event-simulation setups, including biased phase-space sampling, partial unweighting, and the mapping of partonic subprocesses. We apply our methods to the production of Z+jets final states at the HL-LHC using the Sherpa event generator, including matrix elements with up to six final-state partons. When using neural-network surrogates for the dominant Z+5 jets and Z+6 jets partonic processes, we find a reduction in the total event-generation time by more than a factor of 10 compared to baseline Sherpa.