Let $(\Omega\text{, }\mathcal{F}\text{, }\mathbb{P})$ be a probability space and $\{B_t\}_{t\geq0}$ be a standard Brownian Motion.
Consider $\text{Var}\left(\displaystyle \int_{0}^t B_sds\right)$. We have: $$ \text{Var}\left(\displaystyle \int_{0}^t B_sds\right)=\mathbb{E}\left(\displaystyle\int_0^t B_sds \displaystyle\int_0^t B_udu\right)\tag{1} $$ At this point, I know that one can state that: $$ \mathbb{E}\left(\displaystyle\int_0^t B_sds \displaystyle\int_0^t B_udu\right)=\displaystyle\int_0^t ds\displaystyle\int_0^t du \hspace{0.2cm} \mathbb{E}\left(B_sB_u\right)\tag{2} $$
My question is: why does $(2)$ hold true? Is that an application of Tonelli-Fubini Theorem? If so, how is it specifically applied? How does the 'swapping integrals' process precisely occur?
Fubini applies if we can show that $$\int_{\Omega \times [0,t] \times [0,t]} |B_s(\omega)B_u(\omega)|(\mathbb{P} \otimes \lambda\otimes \lambda)(d\omega, ds,du) < \infty$$
However, now the integrand is positive and for positive functions we can use Fubini for positive functions (= Tonelli) freely:
$$\int_{\Omega \times [0,t] \times [0,t]} |B_s(\omega)B_u(\omega)|(\mathbb{P} \otimes \lambda\otimes \lambda)(d\omega, ds,du) $$ $$= \int_0^t \int_0^t \mathbb{E}[|B_s B_u|] ds du$$
The last integral is finite, which I leave to you as a good exercise on Brownian motion.