A great reference for your question is https://arxiv.org/abs/2009.14334, in particular the review in section 2. Another great reference is the asymptotic symmetries review by Compère https://arxiv.org/abs/1801.07064. The basic point is that a phase space is a pair $(\Gamma,{\pmb \Omega})$ where $\Gamma$ is a smooth manifold of even dimension and ${\pmb \Omega}$ is a two-form in $\Gamma$ which is closed, $\mathbf{d}{\pmb \Omega}=0$, and non-degenerate, i.e., ${\pmb \Omega}(\mathbf{X},\mathbf{Y})=0$ for all $\mathbf{Y}$ must imply $\mathbf{X}=0$ where $\mathbf{X}$ and $\mathbf{Y}$ are vector fields in $\Gamma$.
Now, there exists a prescription to construct, in a covariant manner, starting from some Lagrangian form $\mathbf{L}$ in spacetime, a phase space for a field theory. This prescription is the subject of the so called Covariant Phase Space (CPS) Formalism. As a set, $\Gamma$ is taken to be a subset of the space of solutions to the equations of motion. Now, since solutions to the equations of motion are in one-to-one correspondence with their initial data on Cauchy slices $\Sigma\subset M$, it turns out that one may use this data as coordinates on $\Gamma$. Vector fields in this space of fields are just standard field variations.
The (pre)-sympletic form ${\pmb \Omega}$ can be constructed from $\mathbf{L}$ in a precise way which is very well explained in the references I provided. It turns out that it is one integral over such Cauchy slice $\Sigma$, which of course does not depend on which slice you choose. The reason for calling it pre-sympletic form instead of sympletic form is at the heart of your question and I'll explain in a moment. Before that, let me just state what I said more precisely. There exists one object, called pre-sympletic density ${\pmb \omega}$ which is a $(d-1)$-form in spacetime and a two-form in $\Gamma$, from which we obtain ${\pmb \Omega}$ as follows:
$${\pmb \Omega}[\Phi,\mathbf{X}_\Phi,\mathbf{Y}_\Phi]=\int_{\Sigma}{\pmb\omega}[\Phi,\mathbf{X}_\Phi,\mathbf{Y}_\Phi],\tag{1}$$
where ${\pmb \omega}[\Phi,\mathbf{X}_\Phi,\mathbf{Y}_\Phi]$ means the two-form in $\Gamma$ at the point $\Phi\in\Gamma$ evaluated at the vectors $\mathbf{X}_\Phi,\mathbf{Y}_\Phi$ at $T_\Phi\Gamma$. The result is a $(d-1)$-form in spacetime which can be integrated over a $(d-1)$-dimensional slice which we take to be $\Sigma$. This gives the two-form ${\pmb \Omega}$ in $\Gamma$.
Now, in phase space, symmetries are canonical transformations. They are generated by vector fields in $\Gamma$ which are Hamiltonian. These obey $\mathscr{L}_{\mathbf{X}}{\pmb\Omega}=0$ which can be shown to be equivalent to the existence of a charge $Q_\mathbf{X}\in C^\infty(\Gamma)$ such that ${\pmb\Omega}(\mathbf{X},\mathbf{Y})=-\mathbf{Y}(Q_\mathbf{X})$. This is nothing but the differential geometry way of stating that the charge generates the symmetry through the Poisson bracket! In other words, this is just a fancy way of saying that $\{Q,f\}=\delta f$ where $\delta$ is the appropriate symmetry action.
In particular, a gauge transformation with parameter $\varepsilon$ corresponds to a vector field $\mathbf{X}_\varepsilon$ on $\Gamma$ and one may study ${\pmb \Omega}(\mathbf{X}_\varepsilon,\mathbf{Y})$ to study such symmetry and understand its charge. The key thing is that when one writes ${\pmb \Omega}$ as in (1) choosing a slice $\Sigma$, it turns out one has ${\pmb \Omega}(\mathbf{X}_\varepsilon,\mathbf{Y})$ as one integral over $\partial \Sigma$. This is a consequence of Noether's second theorem, which applies to local symmetries.
Now comes the important point. Two things might happen:
For some set of $\varepsilon$, ${\pmb \Omega}(\mathbf{X}_\varepsilon,\mathbf{Y})$ vanishes identically for all variations $\mathbf{Y}$. This is the statement that $\mathbf{X}_\varepsilon$ is a degenerate vector. Wait, but this means that ${\pmb \Omega}$ is a degenerate two-form! Therefore it does not qualify right now as a sympletic form. For that reason we called it pre-sympletic before. One must now get rid of these degenerate vectors to have an actual phase space. The way to do so is by choosing a gauge condition which eliminates all such gauge transformations. This is the rationale behind gauge-fixing in this language. These are the trivial gauge transformations.
So, the bottom line is that if your gauge-fixing does not eliminate all these gauge transformations it is not a good gauge-fixing and it must be supplemented by additional ones until all of these transformations are ruled out.
I urge you to observe that what classifies the gauge transformation as trivial is tied to its behavior at $\partial \Sigma$. In particular if the gauge parameter vanishes at $\partial\Sigma$, i.e., if the gauge transformation acts as the identity there, one has a trivial gauge transformation. Observe that since $\Sigma$ is a Cauchy slice, $\partial\Sigma$ lies at the boundary of spacetime. If you consider for example $\Sigma={\cal I}^+\cup i^+$ or $\Sigma={\cal I}^-\cup i^-$ in asymptotically flat spacetimes, then $\partial\Sigma = \cal I^+_-$ or $\partial\Sigma={\cal I}^-_+$. In either case, it is the behavior of the gauge transformation at the boundary which is important.
In summary this is the reason why residual gauge transformations must affect the boundary of spacetime. If there is some residual gauge transformation which is trivial at the boundary, you just have a bad gauge-fixing because the sympletic form is still degenerate.
For some other set of ${\varepsilon}$ we have ${\pmb \Omega}(\mathbf{X}_\varepsilon,\mathbf{Y})=-\mathbf{Y}(Q_\varepsilon)$ for some $Q_\varepsilon\in C^\infty(\Gamma)$ which is actually a codimension two integral over $\partial\Sigma$. These are not degenerate vectors, they are actual canonical transformations with well-defined Hamiltonians. These are the so-called asymptotic symmetries or large gauge transformations.
Finally, just to mention something that exists, but which is beyond the scope of your question. It might be the case that ${\pmb \Omega}(\mathbf{X}_\varepsilon,\mathbf{Y})$ is non-zero, is still a codimension two integral over $\partial \Sigma$, but is not of the form ${\pmb \Omega}(\mathbf{X}_\varepsilon,\mathbf{Y})=-\mathbf{Y}(Q_\varepsilon)$ because of an additional term which cannot be written as $\mathbf{Y}$ applied to something. In that case one has a non-integrable charge, something which is discussed for example in https://arxiv.org/abs/1801.07064.