A variational method is proposed for calculating the percolation threshold, the real-space structure, and the ground-state energy of a disordered two-dimensional electron liquid. Its high accuracy is verified against exact asymptotics and prior numerical results. The inverse thermodynamical density of states is shown to have a strongly asymmetric minimum at a density that is approximately the triple of the percolation threshold. This implies that the experimentally observed metal-insulator transition takes place well before the percolation point is reached.The discovery that a two-dimensional (2D) electron liquid can be a metal at moderate electron density n e and an insulator at small n e was a major surprise that questioned our fundamental understanding of the role of disorder in such systems. Today, almost a decade later, it remains a subject of an intense debate. 1 One important reason why the conventional theory fails could be its flawed basic premise of the "good" metal, i.e., a uniform electron liquid slightly perturbed by impurities and defects. Indeed, modern nanoscale imaging techniques 2,3,4,5 unambiguously showed that low-density 2D electron systems are strongly inhomogeneous, "bad" metals, where effects of disorder are nonperturbatively strong. In particular, depletion regions (DR), i.e., regions where n(r) is effectively zero, exist. They appear when n e is too small to adequately compensate fluctuating charge density of randomly positioned impurities. As n e is reduced in the experiment, e.g., in order to approach the vicinity of the metal-insulator transition (MIT), the DRs are expected to grow in size and concentration and eventually merge below some percolation threshold n e = n p . An important and controversial issue is whether or not this percolation transition plays any role in the observed MIT. 1 To resolve it one needs to have a theory that is able to calculate n p and that can describe the inhomogeneous structure of the 2D metal at n e ∼ n p . A great progress in this direction has been achieved by Efros, Pikus, and Burnett, 6 whose paper is intellectually tied to earlier work on nonlinear screening by Efros and Shklovskii. 7 Still, analytical results remained scanty and numerical simulations 6,8 were the only known way to quantitatively study the n e ∼ n p regime. These simulations are very time consuming and redoing them in order to get any information beyond what is published 6,8 or to study novel experimental setups seems impracticable. Below I will show that a variational approach to the problem can be a viable alternative. Comparing it with the available numerical results for a typical model of the experimental geometry (Fig. 1), I establish that it correctly predicts the value of n p and accurately reproduces the energetics of the ground state, in particular, the density dependence of the electrochemical potential µ and of the inverse thermodynamical density of states (ITDOS), χ −1 = dµ/dn e , over a broad range of n e . The most striking feature of the resultant function χ −1 (n ...