  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 time scales. However, thus far, theoretical
approaches have been unable to quantitatively reproduce experimentally measured
MBL features -- an important requirement to demonstrate their validity. In this
work, 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, an important feature which has not been addressed in previous
literature. 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.
