You are correct that the standard explanation of "filling up single-electron orbitals" is confusing. That's because it makes two key simplifying assumptions which are rarely stated explicitly:
First, it neglects the Coulomb interaction between the elections, so that the Hamiltonian can be decomposed as
$$H_\text{full} = \sum_{i=1}^n H^{(1)}_i,$$ where $H^{(1)}$ represents a single-electron Hamiltonian (e.g. the hydrogen atom Hamiltonian). In this very special case, it can be shown that the eigenfunctions can all be represented as Slater determinants of single-electron eigenfunctions, so that we really can meaningfully talk about the wave functions of the individual electrons without having to measure all of them at once. In the general entangled case where we fully incorporate the Coulomb interaction, we can't do this, and the "orbital" picture breaks down, and as you say, the physical consequences of the Pauli exclusion principle become very hard to intuit.
In practice, we very often use a hybrid approach called the "Hartree-Fock" approximation (which works surprisingly well and is ubiquitous in quantum chemistry). It's a variational approximation in which we try to minimize the energy of the exact interacting Hamiltonian, but only over the space of Slater determinants of single-particle wavefunctions. In this case it turns out that the best energies come from giving different electrons effective hydrogen-like orbitals, but with different effective nuclear charges that are less than the true nuclear charge $Ze$. Physically, this represents the fact that the interelectron repulsion is being approximately incorporated into a "screening" effect that the inner electrons have on the outer ones, by partially cancelling out the nuclear charge. (Moreover, the best effective nuclear charge turns out to depend on the angular momentum quantum number $l$ (although not on $m$). This breaks the energy degeneracy between orbitals with different values of $l$ that one finds in the hydrogen atom.) But it is inherently just an approximation; in the exact ground-state wave function, you can't talk about individual electron wave functions.
Within the HF approximation, we can assume that each electron has a well-defined orbital, but any arbitrary set of orbitals (that respects Pauli exclusion) is a valid eigenstate. Why do we always assume that they get filled up from lowest to highest energy? Because of the second implicit assumption, which is that the electrons are in thermal equilibrium at zero temperature, so that they are in the ground state of the full multi-electron Hamiltonian. This is an excellent approximation: except in exotic high-temperature systems like plasmas, the electrons are almost always found to be in their ground state. (This is lucky, because it turns out that for the exact exited states, the Hartree-Fock approximation works much less well than for the exact ground state.)