There is a subtle difference between the photon statistics from a single-mode thermal light source with mean photon number $\bar{n}$ and the photon statistics from $N_m$ independent thermal light sources with total mean photon number $\bar{n}$.
In a single-mode thermal light source, the probability distribution of the photon number $n$ is
$$P(n) = \frac{\bar{n}^n}{(\bar{n} + 1)^{n+1}}.$$
It has a variance of
$$\Delta^2 n = \langle n^2 \rangle - \langle n \rangle^2 = \bar{n} + \bar{n}^2,$$
as stated by equation $(1)$ in the question. In the limit of large intensity $\bar{n}^2 \gg \bar{n}$, we can neglect $\bar{n}$ and therefore have
$$\Delta^2n \approx \bar{n}^2 \sim \bar{I}^2.$$
If you average over several independent thermal light sources, the fluctuations will wash out. In a moment where one light source has a high intensity, another one might have a low intensity and vice versa. Let's say each independent light source has a mean photon number of $\bar{n_i} = \bar{n} / N_m$. Each has a variance of
$$\Delta^2 n_i = \bar{n_i} + \bar{n_i}^2.$$
Since the different modes are independent, the overall variance is simply the sum of the individual variances
$$\Delta^2 n = \sum_i \Delta^2 n_i = N_m \left( \bar{n_i} + \bar{n_i}^2 \right) = \bar{n} + \frac{\bar{n}^2}{N_m}.$$
This is equation $(2)$ in the question. If you now take the classical limit by using a large number of independent thermal light sources $N_m$, the result is
$$\Delta^2 n \approx \bar{n} \sim \bar{I},$$
like for a light source with Poissonian statistics! The reason is, as stated above, the independence of the different modes. Therefore, the photons from different modes arrive independently on the detector. DanielSank wrote a nice answer how this leads to a Poissonian distribution.