  While the origin of neutrino masses remains unknown, several key neutrino
mass generation models result in a non-unitary three-neutrino mixing matrix. To
put such models to test, the deviations of the mixing matrix from unitarity can
be measured directly through neutrino oscillation experiments. In this study,
we perform a Bayesian analysis of the non-unitary mixing model using the recent
public data from atmospheric and reactor neutrino experiments - namely
IceCube-DeepCore, Daya Bay, and KamLAND. The novelty of our approach compared
to the preceding global fits for non-unitarity is in the detailed treatment of
the atmospheric neutrino data, which for the first time includes the relevant
flux and detector systematic uncertainties. From the Bayesian posteriors on the
individual mixing matrix elements, we derive the non-unitarity constraints in
the form of normalisations and closures of the mixing matrix rows and columns.
We find comparable constraints for electron and tau row normalisations as other
similar studies in literature, and additionally reveal strong correlations
between muon and tau row constraints induced by the atmospheric systematic
uncertainties. We find that the current data is well described by both unitary
and non-unitary mixing models, with a strong preference for the unitary mixing
indicated by the Bayes factor. With the upcoming IceCube-Upgrade and JUNO
detectors, both featuring superior energy resolution compared to the current
atmospheric and reactor neutrino experiments, our constraints on the row
normalisations are expected to improve by a factor of 2 (3) in the tau (muon)
sector and by nearly 25% in the electron sector. In the future, our approach
can be expanded to include solar and long-baseline neutrino experiments, with
the aim to provide more stringent constraints while keeping track of the
nuisance parameters that may be degenerate with non-unitarity.
