As you already pointed out in the comments, the question reduces to simulating a Brownian sheet. This can be done by generalizing simulation of Brownian motion in a straightforward way.
To simulating the Brownian motion, one can take an i.i.d. mean-0 variance-1 time series $W_i$, $i = 1, 2, \cdots$, and construct the normalized partial sum process
$$
X_n(t) = \frac{1}{\sqrt{n}} \sum_{i = 1}^{[nt]} W_i.
$$
As $n \rightarrow \infty$, $X_n$ convergence weakly (in the sense of Borel probability measures on a metric space) to the standard Brownian $B$ on the Skorohod space $D[0,1]$.
The i.i.d. with finite second moment case is the simplest way to simulate. The mathematical result (Functional Central Limit Theorem/Donsker's Theorem/Invariance Principle) holds in much greater generality.
Now to simulating the (say, two-dimensional) Brownian sheet, takes i.i.d. mean-0 variance-1 array $W_{ij}$, $i,j = 1, 2, \cdots$, and construct the normalized partial sum process
$$
X_n(t_1, t_2) = \frac{1}{ n } \sum_{1 \leq i \leq [nt_1] , 1 \leq j \leq [nt_2]} W_{ij}.
$$
As $n \rightarrow \infty$, $X_n$ convergence weakly to the standard Brownian sheet on the Skorohod space $D([0,1]^2)$ on the unit square.
(The proof is a standard weak convergence argument:
Convergence of finite dimensional distribution follows from the Levy-Lindeberg CLT.
Tightness on $D([0,1]^2)$ follows from a sufficient moment condition that holds trivially in the i.i.d. finite second moment case---see, e.g. Bickel and Wichura (1971).
)
Then, by the continuous mapping theorem
$$
X_n(t_1, t_2) - \prod_{j=1}^2 t_j X_n(t_1, t_2)
$$
converges weakly to the two-dimensional Brownian bridge.