Motivated by a large ground-level ozone dataset, we propose a new computationally efficient additive approximate Gaussian process. The proposed method incorporates a computational-complexity-reduction method and a separable covariance function, which can flexibly capture various spatio-temporal dependence structure. The first component is able to capture nonseparable spatio-temporal variability while the second component captures the separable variation. Based on a hierarchical formulation of the model, we are able to utilize the computational advantages of both components and perform efficient Bayesian inference. To demonstrate the inferential and computational benefits of the proposed method, we carry out extensive simulation studies assuming various scenarios of underlying spatio-temporal covariance structure. The proposed method is also applied to analyze large spatio-temporal measurements of ground-level ozone in the Eastern United States.