Only the mean of means allows for a strict recovery of multiple simulations with limited storage capacities. However, Monte Carlo being an approximation method, approximations are always available.
The theoretical basis of Monte Carlo estimation is the Law of Large Numbers. This means that the empirical cdf
$$\hat{F}_i(x)=\dfrac{1}{n}\sum_{i=1}^n \mathbb{I}_{x_i\le x}$$
is an unbiased estimator of the true cdf $F$. Unfortunately, you cannot store this empirical cdf $\hat{F}_i$ without storing the entire simulation sample of the $n$ $x_i$'s. However, if you replace $\hat{F}_i$ with an approximation based on the empirical percentiles
$$\hat{\hat{F}}_i(x)=\dfrac{1}{99}\sum_{i=1}^{99} \mathbb{I}_{\hat{c}_i\le x}$$where $\hat{c}_i$ is the $i$-th empirical percentile, given by
$$\hat{F}_i(\hat{c}_i)=\dfrac{i}{100}$$
you only need to store $99$ values. (If the percentile precision is not sufficient, you can move to the permile precision or beyond.) You can then repeat simulations without undue pressure on storage by averaging the $\hat{\hat{F}}_i$'s. From this average, you can derive a convergent estimator of the quantiles of interest.
As you noticed, since $\text{cor}(X,Y)$ can be consistently estimated as $$\dfrac{\frac{1}{n}\sum_{i=1}^n x_iy_i -\frac{1}{n}\sum_{i=1}^n x_i\,\frac{1}{n}\sum_{i=1}^n y_i}{\left\{\frac{1}{n}\sum_{i=1}^n x_i^2-\frac{1}{n^2}(\sum_{i=1}^n x_i)^2\right\}^{1/2}\left\{\frac{1}{n}\sum_{i=1}^n y_i^2-\frac{1}{n^2}(\sum_{i=1}^n y_i)^2\right\}^{1/2}}$$you can easily update the five quantities involved in this expression
$$\sum_{i=1}^n x_iy_i,\ \sum_{i=1}^n x_i,\ \sum_{i=1}^n y_i,\ \sum_{i=1}^n x_i^2,\ \sum_{i=1}^n y_i^2$$
without increase storage pressure.
