  We study the double Higgs boson production processes $e^+e^- \to hh f\bar{f}$
($f\neq t$) with $h$ being the 125 GeV Higgs boson in the two-Higgs-doublet
model with a softly-broken $Z_2$ symmetry. The cross section can be
significantly enhanced, typically a few hundreds percent, as compared to the
standard model prediction due to resonant effects of heavy neutral Higgs
bosons, which becomes important in the case without the alignment limit. We
find a strong correlation between the enhancement factor of the cross section
and the scaling factor of the $hf\bar{f}$ couplings under constraints from
perturbative unitarity, vacuum stability and current experimental data at the
LHC as well as the electroweak precision data.
