Abstract. In this paper, we propose a multiscale data-driven stochastic method (MsDSM) to study stochastic partial differential equations (SPDEs) in the multiquery setting. This method combines the advantages of the recently developed multiscale model reduction method [M. L. Ci, T. Y. Hou, and Z. Shi, ESAIM Math. Model. Numer. Anal., 48 (2014), pp. 449-474] and the datadriven stochastic method (DSM) [M. L. Cheng et al., SIAM/ASA J. Uncertain. Quantif., 1 (2013), pp. 452-493]. Our method consists of offline and online stages. In the offline stage, we decompose the harmonic coordinate into a smooth part and a highly oscillatory part so that the smooth part is invertible and the highly oscillatory part is small. Based on the Karhunen-Loève (KL) expansion of the smooth parts and oscillatory parts of the harmonic coordinates, we can derive an effective stochastic equation that can be well-resolved on a coarse grid. We then apply the DSM to the effective stochastic equation to construct a data-driven stochastic basis under which the stochastic solutions enjoy a compact representation for a broad range of forcing functions. In the online stage, we expand the SPDE solution using the data-driven stochastic basis and solve a small number of coupled deterministic partial differential equations (PDEs) to obtain the expansion coefficients. The MsDSM reduces both the stochastic and the physical dimensions of the solution. We have performed complexity analysis which shows that the MsDSM offers considerable savings over not only traditional methods but also DSM in solving multiscale SPDEs. Numerical results are presented to demonstrate the accuracy and efficiency of the proposed method for several multiscale stochastic problems without scale separation. 1. Introduction. In recent years, there has been an increased interest in the simulation of systems with uncertainties. Many physical and engineering applications involving uncertainty quantification can be described by stochastic partial differential equations (SPDEs). Several numerical methods have been developed in the literature to solve SPDEs, such as the stochastic finite element method [24], Wiener chaos expansion or generalized polynomial chaos (gPC) method [31,25,40,41,42,39,33,32], and stochastic collocation method [43,5,30,6]. These methods can effectively solve the SPDEs when the dimension of stochastic input variables is low. However, their performance deteriorates dramatically when the dimension of stochastic input variables is high. This so-called curse of dimensionality is one of the essential challenges in uncertainty quantification. Recently, some progress has been made to alleviate this difficulty by exploring the sparse structure of the solutions and constructing a problem-dependent stochastic basis to solve these SPDEs; see, e.g., the data-driven stochastic method [11,44] and dynamically biorthogonal method [9,10].