We have developed a novel reservoir characterization workflow that solves the reservoir history matching problem coupling a physics Informed neural operator (PINO) forward model with a mixture of experts' approach termed cluster classify regress (CCR). The inverse modelling is computed via an adaptive Regularized Ensemble Kalman inversion (aREKI) method, and the method is optimal for rapid inverse uncertainty quantification tasks during reservoir history matching. We parametrize the unknown permeability and porosity fields for the case of sampling non-gaussian posterior measures by using, for instance, a variational convolution autoencoder and another instance using a denoising diffusion implicit model (DDIM). The CCR serves as exotic priors in the inverse modelling and as a coupled supervised model with the prior PINO surrogate for replicating the nonlinear Peaceman well equations in the forward model. The CCR is flexible, and any separate and independent machine learning algorithm can be used for each of the 3 stages. The loss function required for the PINO reservoir surrogate updating is obtained from the combined supervised data loss and losses arising from the initial conditions and residual of the governing black oil pde. The output from the coupled PINO-CCR surrogate is the pressure, water and gas saturations as well as the well outputs of oil, water, and gas production rates. The developed methodology is first compared to a standard numerical black oil reservoir simulator for a waterflooding case on the Norne field, and the similarity of the outputs with those from the PINO-CCR surrogate is very close. Next, we use this developed PINO-CCR surrogate in the developed aREKI history matching workflow. The overall workflow is successful in recovering the unknown permeability and porosity field and simulation is very fast with a speedup of up to 6000X to conventional numerical reservoir simulation approaches. Learning the PINO-CCR surrogate on an NVIDIA H100 with 80G memory takes approximately 5 hours for 100 training samples of the Norne field. The workflow is also suitable to be used in an ensemble-based workflow, where sampling the posterior density, given an expensive likelihood evaluation, to quantify uncertainty, as with most scientific computation, is desirable.