  The Two Higgs Doublet Model (2HDM) with spontaneously broken $Z_2$ symmetry
predicts a production of domain walls at the electroweak scale. We derive
cosmological constraints on model parameters for both Type-I and Type-II 2HDMs
from the requirement that domain walls do not dominate the Universe by the
present day. For Type-I 2HDMs, we deduce the lower bound on the key parameter
$\tan\beta > 10^5$ for a wide range of Higgs-boson masses $\sim$ 100 GeV or
greater close to the Standard Model alignment limit. In addition, we perform
numerical simulations of the 2HDM with an approximate as well as an exact $Z_2$
symmetry but biased initial conditions. In both cases, we find that domain wall
networks are unstable and, hence, do not survive at late times. The domain
walls experience an exponential suppression of scaling in these models which
can help ameliorate the stringent constraints found in the case of an exact
discrete symmetry. For a 2HDM with softly-broken $Z_2$ symmetry, we relate the
size of this exponential suppression to the soft-breaking bilinear parameter
$m_{12}$ allowing limits to be placed on this parameter of order $\mu$eV, such
that domain wall domination can be avoided. In particular, for Type-II 2HDMs,
we obtain a corresponding lower limit on the CP-odd phase $\theta$ generated by
QCD instantons, $\theta \ \stackrel{>}{{}_\sim}\ 10^{-11}/(\sin\beta
\cos\beta)$, which is in some tension with the upper limit of $\theta \
\stackrel{<}{{}_\sim}\ 10^{-11}$--$10^{-10}$, as derived from the
non-observation of a non-zero neutron electric dipole moment. For a
$Z_2$-symmetric 2HDM with biased initial conditions, we are able to relate the
size of the exponential suppression to a biasing parameter $\varepsilon$ so as
to avoid domain wall domination.
