The paper introduces a non-linear version of the process convolutionformalism for building covariance functions for multi-output Gaussianprocesses. The non-linearity is introduced via Volterra series, one series pereach output. We provide closed-form expressions for the mean function and thecovariance function of the approximated Gaussian process at the output of theVolterra series. The mean function and covariance function for the jointGaussian process are derived using formulae for the product moments of Gaussianvariables. We compare the performance of the non-linear model against theclassical process convolution approach in one synthetic dataset and two realdatasets.