We developed a novel reservoir characterization workflow that addresses reservoir history matching by coupling a physics-informed neural operator (PINO) forward model with a mixture of experts' approach, termed cluster classify regress (CCR). The inverse modelling is achieved via an adaptive Regularized Ensemble Kalman inversion (aREKI) method, ideal for rapid inverse uncertainty quantification during history matching. We parametrize unknown permeability and porosity fields for non-Gaussian posterior measures using a variational convolution autoencoder and a denoising diffusion implicit model (DDIM) exotic priors. The CCR works as a supervised model with the PINO surrogate to replicate nonlinear Peaceman well equations. The CCR's flexibility allows any independent machine-learning algorithm for each stage. The PINO reservoir surrogate's loss function is derived from supervised data loss and losses from the initial conditions and residual of the governing black oil PDE. The PINO-CCR surrogate outputs pressure, water, and gas saturations, along with oil, water, and gas production rates. The methodology was compared to a standard numerical black oil simulator for a waterflooding case on the Norne field, showing similar outputs. This PINO-CCR surrogate was then used in the aREKI history matching workflow, successfully recovering the unknown permeability, porosity and fault multiplier, with simulations up to 6000 times faster than conventional methods. Training the PINO-CCR surrogate on an NVIDIA H100 with 80G memory takes about 5 hours for 100 samples of the Norne field. This workflow is suitable for ensemble-based approaches, where posterior density sampling, given an expensive likelihood evaluation, is desirable for uncertainty quantification.