  It is well-known that foreground subtraction in 21cm surveys removes large
scale power. We investigate associated systematic biases. We show that removing
line-of-sight fluctuations on large scales aliases into suppression of the 3D
power spectrum across a broad range of scales. This bias can be eliminated by
marginalizing over small k in the 1D power spectrum; however, the unbiased
estimator will have unavoidably larger variance. We also show that Gaussian
realizations of the power spectrum permit accurate and extremely rapid
Monte-Carlo simulations for error analysis; repeated realizations of the fully
non-Gaussian field are unnecessary. We perform Monte-Carlo maximum-likelihood
simulations of foreground removal which yield unbiased, minimum variance
estimates of the power spectrum in agreement with Fisher matrix estimates.
Foreground removal also distorts the 21cm PDF, reducing the contrast between
neutral and ionized regions. We show that it is the subtraction of large-scales
modes which is responsible for this distortion, and that it is less severe in
the earlier stages of reionization. It can be reduced by using larger
bandwidths for foreground removal. In the late stages of reionization, the
largest ionized regions (which consist of foreground emission only) provides
calibration points which potentially allow recovery of large-scale modes.
Finally, we also show that: (i) the broad frequency response of synchrotron and
free-free emission will smear out any features in the electron momentum
distribution and ensure spectrally smooth foregrounds; (ii) extragalactic radio
recombination lines should be negligible foregrounds.
