The problem of land subsidence due to groundwater pumping in highly compressible, randomly heterogenous aquitards with stress-dependent parameters is analyzed through one-dimensional Monte Carlo numerical simulations of nonlinear groundwater flow and consolidation. The natural logarithm of hydraulic conductivity (Y) and compression index (Cc) are considered spatially cross-correlated random variables. A stochastic simulation method based on a bivariate normal copula was applied to generate positive, negative, and uncorrelated realizations of the two random variables. The results show that when Y and Cc are correlated random parameters, the ensemble average total settlement is consistently larger (between 13% and 40% for the synthetic cases tested) than the total settlement obtained when the parameters are deterministic (perfectly known). In addition, the results of negative and no correlation between Y and Cc are similar in terms of the ensemble average and variance of total settlement, while positive correlation leads to larger values of total settlement and its variance. The stochastic approach employed here allows incorporating the uncertainty in the parameters in the simulation of the subsidence process, which is relevant for geotechnical engineering, design of structures and groundwater management.