Stochastic methods have been typically used for the design and operations of hydraulic infrastructure. They allow decision makers to evaluate existing or new infrastructure under different possible scenarios, giving them the flexibility and tools needed in decision making. In this paper, we present a novel stochastic streamflow simulation approach able to replicate both temporal and spatial dependencies from the original data in a multi-site basin context. The proposed model is a multi-site extension of the modified Fractional Gaussian Noise (mFGN) model which is well-known to be efficient to maintain periodic correlation for several time lags, but presents shortcomings in preserving the spatial correlation. Our method, called Weighted-mFGN (WmFGN), incorporates spatial dependency into streamflows simulated with mFGN by relying on the Cholesky decomposition of the spatial correlation matrix of the historical streamflow records. As the order in which the decomposition steps are performed (temporal then spatial, or vice-versa) affects the performance in terms of preserving the temporal and spatial correlation, our method searches for an optimal convex combination of the resulting correlation matrices. The result is a Pareto-curve that indicates the optimal weights of the convex combination depending on the importance given by the user to spatial and temporal correlations. The model is applied to a number of river basins in Chile, where the results show that the WmFGN approach maintains the qualities of the single-site mFGN, while significantly improving spatial correlation.