Approximating the variance of the integral of a white noise Gaussian process

The disparity arises from the fact that your discretization of the continuous process does not assign the appropriate variance to the $X_i$. Here's the key (for heuristics, see here, Section 3.2):

If $\{X(t)\}_{t\in\mathbb{R}}$ is a continuous Gaussian process such that $$\begin{align}&E[X(t)]=0\\ &E[X(t)X(t')]=\sigma^2 \delta(t-t')\end{align} $$ then a discrete sampling of the process, with a uniform sampling interval $\Delta$, viz., $$X(t),X(t+1\Delta),X(t+2\Delta),...$$ is simulated as an i.i.d sequence of Gaussian$(\text{mean}=0,\text{variance}=\frac{\sigma^2}{\Delta})$ random variables, the quality of the simulation increasing as $\Delta\to 0$.

(The power spectral density of the i.i.d. sequence approaches that of the continuous process that it simulates -- flat and equal to the same constant value $\sigma^2$ -- except that for the i.i.d. sequence it is "band limited", i.e. vanishing outside of a finite-width interval.)

Thus, in the present case, to simulate $$I = \int_0^L X(t) \, dt \approx \sum_{i=1}^N X(i\Delta)\,\Delta $$ where $\Delta=\frac{L}{N}$, one would use $$\hat{I} = \Delta\sum_{i=1}^N X_i $$ with $X_1,X_2,\ldots$ i.i.d. $\text{Gaussian}(\text{mean}=\mu,\text{variance}=\frac{\sigma^2}{\Delta})$. Then $$\begin{align} \text{Var}[\hat{I}] = \Delta^2\,N\,\frac{\sigma^2}{\Delta} = (N\Delta)\sigma^2 = L\sigma^2. \end{align} $$

I think somewhere between your first StackExchange post reference in your CrossValidated post there's been some confusion about the process in question.

It seems to me the answer to your CrossValidated post regards a Wiener process/Brownian motion with independent increments in the process. Or perhaps the person providing that answer mistook the integral of the process you describe for such a process.

To me your numerical result seems correct.