History matching and optimization problems often involve several, potentially conflicting, objectives. For example, we might seek to minimize a misfit function involving differences in reservoir pressure, multiphase production history and 4D time-lapse seismic data and these differences do not always change in tandem. It is a common practice to treat these differences as a single objective optimization problem by aggregating all objectives into a scalar function (weighted-sum), resulting in incomplete exploration of the solution space. The problem is particularly severe if the objectives are conflicting. In this paper we propose to use a Pareto-based multi-objective evolutionary algorithm (MOEA) focusing on finding a set of optimal solutions called Pareto optima. The MOEA makes direct use of the dominance relation for fitness assignment instead of a fitness score in one-dimensional objective space. The dominance concept can define levels of optimality without reduction of objective dimensionality to sort populations accordingly, and the given populations constitute typically several ranks (fronts) of classification for individuals. Because it uses a population of solutions in the search process and optimizes such that the ranks are minimized, the Pareto optima can provide a measure of uncertainty in prediction. We show how the MOEA identifies optimal solutions by examining the trade-off between multiple objectives from a set of plausible solutions. Specifically, we demonstrate that it outperforms the commonly used weighted-sum approach. For practical applications, we provide a novel history matching workflow with a grid connectivity-based transformation (GCT) basis coefficients as parameters for calibration using the gradient-free evolutionary optimization algorithms. The basis functions are obtained from a spectral decomposition of the grid connectivity Laplacian and avoid ad hoc redefinitions of regions while preserving the geologic heterogeneity. We demonstrate the power and utility of the proposed workflow using multiple examples. These include 2D synthetic examples for validation and a 3D field application for matching production and seismic data with uncertainty and conflicting information.