While many studies point towards the existence of many-body localization (MBL) in one dimension, the fate of higher-dimensional strongly disordered systems is a topic of current debate. The latest experiments as well as several recent numerical studies indicate that such systems behave many-body localized - at least on practically relevant timescales. However, thus far, theoretical approaches have been unable to quantitatively reproduce experimentally measured MBL features - an important requirement to demonstrate their validity. In this Letter, we use fermionic quantum circuits as a variational method to approximate the full set of eigenstates of two-dimensional MBL systems realized in fermionic optical lattice experiments. Using entanglement-based features, we obtain a phase transition point in excellent agreement with the experimentally measured value. Moreover, we calculate the filling-fraction-dependent MBL phase diagram. We argue that our approach best captures the underlying charge-density-wave experiments and compute the mean localization lengths, which can be compared to future experiments.