  We fit the next-to-leading order unitarity conditions to the
Two-Higgs-Doublet model with a softly broken $\mathbb{Z}_2$ symmetry. In doing
so, we alleviate the existing uncertainty on how to treat higher order
corrections to quartic couplings of its Higgs potential. A simplified approach
to implementing the next-to-leading order unitarity conditions is presented.
These new bounds are then combined with all other relevant constraints,
including the complete set of LHC Run I data. The upper $95\%$ bounds we find
are $4.2$ on the absolute values of the quartic couplings, and $235$ GeV ($100$
GeV) for the mass degeneracies between the heavy Higgs particles in the type I
(type II) scenario. In type II, we exclude an unbroken $\mathbb{Z}_2$ symmetry
with a probability of $95\%$. All fits are performed using the open-source code
HEPfit.
