I would like to generate data with "Model 1" and fit them with "Model 2". The underlying idea is to investigate robustness properties of "Model 2". I am particularly interested in the coverage rate of the 95% confidence interval (based on the normal approximation).
I often use the width of confidence intervals as a quick-and-dirty way to determine the number of iterations needed.
Let $p$ be the true coverage rate of the 95 % confidence interval when data from "Model 1" is fitted to "Model 2". If $X$ is the number of times that the confidence interval covers the true parameter value in $n$ iterations, then $X\sim {\rm Bin}(n,p)$.
The estimator $\hat{p}=X/n$ has mean $p$ and standard deviation $\sqrt{p(1-p)/n}$. For large $n$, $\hat{p}$ is approximately normal and $\hat{p}\pm 1.96\sqrt{\hat{p}(1-\hat{p})/n}$ gives you an approximately 95 % confidence interval for $p$. Since you know (would gess) that $p\approx 0.95$, it follows that the width of this interval is approximately $2\cdot 1.96\sqrt{0.95\cdot 0.05/n}$.
If you think that a confidence interval with width $0.1$ (say) is acceptable, you find the approximate number of iterations $n$ needed for this by solving the equation $$0.1=2\cdot 1.96\sqrt{0.95\cdot 0.05/n}.$$
In this way you can find a reasonable $n$ by choosing the accuracy that you are looking for.
Based on your follow up comment it sounds like you are trying to estimate the coverage probability of a confidence interval when you assume constant error variance when the true error variance is not constant.
The way I think about this is that, for each run, the confidence interval either covers the true value or it doesn't. Define an indicator variable:
$$ Y_i = \begin{cases} 1 & {\rm if \ the \ interval \ covers} \\ 0 & {\rm if \ it \ does \ not } \end{cases}$$
Then the coverage probability you're interested in is $E(Y_i) = p$ which you can estimate by the sample proportion which I think is what you're proposing.
How do I set the number of iteration runs?
We know that the variance of a Bernoulli trial is $p(1-p)$, and your simulations will generate IID bernoulli trials, therefore the variance of your simulation based estimate of $p$ is $p(1-p)/n$, where $n$ is the number of simulations. You can choose $n$ to shrink this variance as much as you want. It is a fact that $$p(1-p)/n \leq 1/4n$$
So, if you want the variance to be less than some pre-specified threshold, $\delta$, then you can ensure this by choosing $n \geq 1/4\delta$.
In a more general setting, if you're trying to investigate properties of the sampling distribution of an estimator by simulation (e.g. it's mean and variance) then you can choose your number of simulations based on how much precision you want to achieve in an analogous fashion to that described here.
Also note that, when the mean (or some other moment) of a variable is the object of interest, as it is here, you can construct a confidence interval for it based on the simulations using the normal approximation (i.e. the central limit theorem), as discussed in MansT's nice answer. This normal approximation is better as the number of samples grows, so, if you plan on constructing a confidence interval by appealing to the central limit theorem, you will want $n$ to be large enough for that to apply. For the binary case, as you have here, it appears this approximation is good even when $np$ and $n(1-p)$ are pretty moderate - say, $20$.
Is it true that larger than necessary replications may result in spurious biases? If so, how is that?
As I mentioned in a comment - this depends on what you mean by spurious. Larger numbers of simulations will not produce bias in the statistical sense, but it may reveal an unimportant bias that is only noticeable with an astronomically large sample size. For example, suppose the true coverage probability of the misspecified confidence interval were $94.9999\%$. Then, this isn't really a problem in a practical sense, but it you may only pick up this difference if you ran a ton of simulations.
If you are doing a simulation the minimum number of required runs depends on your objective (What are you trying to estimate and with what accuracy?). If you are trying to estimate the average response then the standard deviation of the sample average is the $\dfrac{\text{Population Standard Deviation}}{\sqrt{n}}$. So if $d$ is the required half-width for $95\%$ confidence interval for the mean you want $d= 1.96 \times \dfrac{\text{Pop.Std.Dev}}{\sqrt{n}}$ or $n=\dfrac{ (1.96 \times\text{Pop.Std.Dev})^2}{d^2}$.
Doing more simulations (assuming all samples arre generated by a random process) does nothing to hurt the estimation in terms of accuracy or bias.
The coverage of an approximate confidence interval will differ from the exact $95\%$ coveraged desired and the error in coverage should decrease with increasing $n$. As mentioned by Macro and MansT, you can bound the Monte Carlo estimate of coverage based on the variance of the binomial proportion being $\dfrac{p(1-p)}{n}$.
