A computational method based on the non-linear Gaussian process (GP), known as deep Gaussian processes (deep GPs) for uncertainty quantification & propagation in modelling of flow through heterogeneous porous media is presented. The method is also used for reducing dimensionality of model output and consequently emulating highly complex relationship between hydrogeological properties and reduced order fluid velocity field in a tractable manner. Deep GPs are multi-layer hierarchical generalisations of GPs with multiple, infinitely wide hidden layers that are very efficient models for deep learning and modelling of high-dimensional complex systems by tackling the complexity through several hidden layers connected with non-linear mappings. According to this approach, the hydrogeological data is modelled as the output of a multivariate GP whose inputs are governed by another GP such that each single layer is either a standard GP or the Gaussian process latent variable model. A variational approximation framework is used so that the posterior distribution of the model outputs associated to given inputs can be analytically approximated. In contrast to the other dimensionality reduction, methods that do not provide any information about the dimensionality of each hidden layer, the proposed method automatically selects the dimensionality of each hidden layer and it can be used to propagate uncertainty obtained in each layer across the hierarchy. Using this, dimensionality of the full input space consists of both geometrical parameters of modelling domain and stochastic hydrogeological parameters can be simultaneously reduced without the need for any simplifications generally being assumed for stochastic modelling of subsurface flow problems. It allows estimation of the flow statistics with greatly reduced computational efforts compared to other stochastic approaches such as Monte Carlo method.