The retrieval algorithms in remote sensing generally involve complex physical forward models that are nonlinear and computationally expensive to evaluate. Statistical emulation provides an alternative with cheap computation and can be used to quantify uncertainty, calibrate model parameters, and improve computational efficiency of the retrieval algorithms. Motivated by this, we introduce a framework of building statistical emulators by combining dimension reduction of input and output spaces and Gaussian process modeling. The functional principal component analysis (FPCA) via a conditional expectation method is chosen to reduce the dimension of the output space of the forward model. In addition, instead of making restrictive assumptions regarding the dependence structure of the high-dimensional input space, we apply the gradient-based kernel dimension reduction method to reduce the dimension of input space. A computationally efficient Gaussian process model is then constructed at the low-dimensional input and output spaces. Theoretical properties of the resulting statistical emulator are explored, and the performance of our approach is demonstrated with numerical studies with NASA’s Orbiting Carbon Observatory-2 (OCO2) data.