Numerical simulation of multiphase flow in porous media is essential for many geoscience applications. Data-driven machine learning methods provide faster alternatives to traditional simulators by training neural network models with numerical simulation data mappings. Here we present U-FNO, a novel neural network architecture for solving multiphase flow problems with superior speed, accuracy, and data efficiency. U-FNO is designed based on the newly proposed Fourier neural operator (FNO) that learns an infinite-dimensional integral kernel in the Fourier space, which has shown excellent performance for single-phase flows. Here we extend the FNO-based architecture to a CO2-water multiphase problem, and proposes the U-FNO architecture to enhance the prediction accuracy in multiphase flow systems. Through a systematic comparison among a CNN benchmark and three types of FNO variations, we show that the U-FNO architecture has the advantages of both the traditional CNN and original FNO, providing significantly more accurate and efficient performance than previous architectures. The trained U-FNO predicts gas saturation and pressure buildup with a 6*10e4 times speed-up compared to traditional numerical simulators while maintaining similar accuracy. The trained models can act as a general-purpose simulator alternative for 2D-radial CO2 injection problems with wide ranges of permeability and porosity heterogeneity, anisotropy, reservoir conditions, injection configurations, flow rates, and multiphase flow properties.