  For the iterative decoupling of elliptic-parabolic problems such as
poroelasticity, we introduce time discretization schemes up to order $5$ based
on the backward differentiation formulae. Its analysis combines techniques
known from fixed-point iterations with the convergence analysis of the temporal
discretization. As the main result, we show that the convergence depends on the
interplay between the time step size and the parameters for the contraction of
the iterative scheme. Moreover, this connection is quantified explicitly, which
allows for balancing the single error components. Several numerical experiments
illustrate and validate the theoretical results, including a three-dimensional
example from biomechanics.
