In this paper, Isomap and kernel Isomap are used to dramatically reduce the dimensionality of the output space to efficiently construct a Gaussian process emulator of parametrized partial differential equations. The output space consists of spatial or spatio-temporal fields that are functions of multiple input variables. For such problems, standard multi-output Gaussian process emulation strategies are computationally impractical and/or make restrictive assumptions regarding the correlation structure. The method we develop can be applied without modification to any problem involving vector-valued targets and vector-valued inputs. It also extends a method based on linear dimensionality reduction to response surfaces that cannot be described accurately by a linear subspace of the high dimensional output space. Comparisons to the linear method are made through examples that clearly demonstrate the advantages of nonlinear dimensionality reduction.