The Reynolds-averaged Navier-Stokes (RANS) simulations are commonly used in industrial applications due to their computational efficiency. However, the linear eddy viscosity model (LEVM) used in RANS often fails to accurately capture the anisotropy of Reynolds stress in complex flow conditions. To enhance RANS predictive accuracy, data-driven closure models, such as Tensor Basis Neural Network (TBNN) and Tensor Basis Random Forest (TBRF), have been proposed. However existing models, including TBNN and TBRF, have limitations in capturing the nonlocal patterns of turbulence models, resulting in irregular and unsmooth predictions. Convolutional neural networks (CNNs) are considered as an alternative approach, but their reliance on discretization poses challenges when dealing with arbitrarily designed meshes in RANS simulations. In this study, we propose a nonlinear convolutional neural operator as the RANS closure model. Our model satisfies Galilean invariance, can learn nonlocal physics, and recovers high-resolution physics even when trained on undersampled grids. The model outperforms existing TBNN and TBRF models, successfully predicting smooth fields of Reynolds stress in flows with adverse pressure gradients, separations, and streamline curvature, where existing models struggle or fail to provide accurate predictions.