Abstract: We propose and construct a two-parameter perturbative expansion  around a Friedmann-Lemaître-Robertson-Walker geometry that can be used  to model high-order gravitational effects in the presence of nonlinear  structure. This framework reduces to the weak-field and slow-motion  post-Newtonian treatment of gravity in the appropriate limits, but also  includes the low-amplitude large-scale fluctuations that are important  for cosmological modeling. We derive a set of field equations that can  be applied to the late Universe, where nonlinear structure exists on  supercluster scales, and perform a detailed investigation of the  associated gauge problem. This allows us to identify a consistent set of  perturbed quantities in both the gravitational and matter sectors, and  to construct a set of gauge-invariant quantities that correspond to each  of them. The field equations, written in terms of these quantities,  take on a relatively simple form, and allow the effects of small-scale  structure on the large-scale properties of the Universe to be clearly  identified. We find that inhomogeneous structures source the global  expansion, that there exist new field equations at new orders, and that  there is vector gravitational potential that is a hundred times larger  than one might naively expect from cosmological perturbation theory.  Finally, we expect our formalism to be of use for calculating  relativistic effects in upcoming ultralarge-scale surveys, as the form  of the gravitational coupling between small and large scales depends on  the nonlinearity of Einstein’s equations, and occurs at what is normally  thought of as first order in cosmological perturbations.