Discrete Graphical Models (GMs) are widely used in Artificial Intelligence to describe complex systems through a joint function of interest. Probabilistic GMs such as Markov Random Fields (MRFs) define a joint non-normalized probability distribution while deterministic GMs such as Cost Function Networks (CFNs) define a joint cost function. A typical query on GMs consists in finding the joint state that optimizes this joint function, a problem denoted as the \textit{Maximum a Posteriori} or Weighted Constraint Satisfaction Problem respectively. In practice, more than one function of interest may need to be optimized at the same time. In this paper, we develop a two-phase scalarization method for solving bi-objective discrete graphical model optimization, with the aim of computing a set of non-dominated solutions --- the Pareto frontier --- representing different compromises between two GM-defined objectives. For this purpose, we introduce a dedicated higher-order constraint, which bounds the value of one GM-defined objective while minimizing another GM on the same variables. Discrete GM optimization is NP-hard, and its bi-objective variants are even harder. We show how existing GM global lower and upper bounds can be exploited to provide anytime bounds of the exact Pareto frontier. We benchmark it on various instances from Operations Research and Probabilistic Graphical Models.