The variational free energy profile (vFEP) method is extended to two dimensions and tested with molecular simulation applications. The proposed 2D-vFEP approach effectively addresses the two major obstacles to constructing free energy profiles from simulation data using traditional methods: the need for overlap in the reweighting procedure and the problem of data representation. This is especially evident as these problems are shown to be more severe in two dimensions. The vFEP method is demonstrated to be highly robust and able to provide stable analytic free energy profiles with only a paucity of sampled data. The analytic profiles can be analyzed with conventional search methods to easily identify stationary points (e.g., minima and first-order saddle points) as well as the pathways that connect these points. These “roadmaps” through the free energy surface are useful not only as a post-processing tool to characterize mechanisms but can also serve as a basis from which to direct more focused “on-the-fly” sampling or adaptive force biasing. Test cases demonstrate that 2D-vFEP outperforms other methods in terms of the amount and sparsity of the data needed to construct stable, converged, analytic free energy profiles. In a classic test case, the two-dimensional free energy profile of the backbone torsion angles of alanine dipeptide, 2D-vFEP, needs less than 1% of the original data set to reach a sampling accuracy of 0.5 kcal/mol in free energy shifts between windows. A new software tool for performing one- and two-dimensional vFEP calculations is herein described and made publicly available.