Physical interpretation of the nil ADM mass of gravitational waves The ADM mass is defined for any asymptotic flat spacetime.  Using cartesian coordinates:
\begin{equation}\tag{1}
  E_{\text{ADM}} = -\: \frac{1}{16 \pi G} \, \lim_{r \, \rightarrow \, \infty}\oint_{\mathcal{S}_r} (\, \partial_j \, h_{ij} - \partial_i \, h_{kk}) \, dS_i.
\end{equation}
I consider linearized theory: $g_{\mu \nu} = \eta_{\mu\nu} + h_{\mu \nu}$, and gravitational waves emitted from an arbitrary localized source.  In the $\textsf{TT}$ (transverse-traceless) gauge, we have the following cartesian components:
\begin{align}
  h_{00}^{\textsf{TT}} = h_{kk}^{\textsf{TT}} &= 0, \tag{2} \\[1em]
  \partial_i \, h_{ij}^{\textsf{TT}} &= 0, \tag{3}
\end{align}
so (1) gives trivially $E_{\text{ADM}} = 0$. (take note that all waves quantities must be evaluated at some retarded time)
Now, the physical interpretation.  I know that the ADM mass doesn't give the wave's energy moved to infinity.  The ADM mass is the energy of the source (which should be computed with other linear terms from the metric, since (1) is linear).
But then how can I properly interpret the result $E_{\text{ADM}} = 0$ for the gravitational waves?  Am I right in saying that it's actually a physical inconsistency of the first order approximation, and that the weak waves don't have a "mass" contribution to the source since the waves themselves don't generate other waves (linearized theory, waves too weak to have a feedback on the sources), while the full non-linear theory should imply a mass contribution from the strong waves?
This question is related to this one, for which the answers aren't satisfying:
The ADM energy of gravitational waves?
 A: Let us start with electromagnetic flat space analogy.
Consider the following quantity:
$$
\mathcal{Q}=-\frac{1}{4\pi}\oint_S\left(\frac1c \frac{\partial\mathbf{A}}{\partial t} +\nabla\Phi\right)\cdot \mathbf{n}\,dS,
$$
where the integration is done over the “large” sphere containing all charges within the Minkowski space at specific moment $t=t_0$.
Now, let us impose the Coulomb gauge condition on the EM potential: $\mathbf{\nabla} \cdot \mathbf{A}=0$ and write the potentials $(\Phi_\text{rad},\mathbf{A}_\text{rad})$ of an EM radiation from an isolated source.
By analogy with OP we could ask ourselves:

*

*What would be the quantity $\mathcal{Q}$ calculated for the potentials $(\Phi_\text{rad},\mathbf{A}_\text{rad})$?


*If it is zero, what is the “deepest” reason for that?


*If instead of electromagnetism we consider some non-Abelian gauge theory, would corresponding generalization of the radiation fields $(\Phi_\text{rad},\mathbf{A}_\text{rad})$ produce nonzero $\mathcal{Q}$?
The quantity $\mathcal{Q}$ is just the charge of the system (in Gaussian units), so obviously it would be zero for purely radiative electromagnetic fields  $(\Phi_\text{rad},\mathbf{A}_\text{rad})$. One obvious reason for that is the fact that EM field itself does not carry charge. But we must also point to some vagueness in OP's expression for ADM energy (and also in our definition of $\mathcal{Q}$): rather than integrating over a single unspecified “large sphere” we must consider a limiting procedure over a sequence of spheres $S_R$ with increasing radii $R$:
$$
\text{instead of} \, \oint_S \,\text{ we must use } \lim_{R\to \infty} \oint_{S_R}.
$$
But the time is held fixed in this limiting procedure, and since EM radiation propagates with the speed of light, radiation fields (and their derivatives) for any source radiating finite total energy would be zero in this limit:
$$
\lim_{R\to \infty} (\Phi_\text{rad},\mathbf{A}_\text{rad})\Big|_{|\mathbf{r}|=R,\,t=t_0}=0
$$
Another point to remember is that our ability to write the total charge as an integral over a distant sphere instead of volume integral is because the electromagnetic theory due to the gauge invariance has constraints equations (that is equations without time derivatives). For field strengths that would be equations $\nabla\cdot \mathbf{E}=4\pi \rho$ and $\nabla\cdot \mathbf{B}=0$, while for the potentials one equation would be the Coulomb gauge condition and the other is  equation for the scalar potential:
$$
\Delta \Phi = - 4\pi \rho.
$$
One way to interpret this equation is to think of $\Phi$ as “propagating” with infinite speed. And we also can write the expression for charge $\mathcal{Q}$ in terms of $\Phi$ only:
$$
\mathcal{Q}=-\frac{1}{4\pi}\lim_{R\to\infty}\oint_{S_R}\left(\nabla\Phi\right)\cdot \mathbf{n}\,dS.
$$
So the “moral” reason for why radiative fields do not contribute to $\mathcal{Q}$ could be stated as follows: Charge $\mathcal{Q}$ could be written as a surface integral in the limit $R\to \infty$, $t=t_0$, and there could be no radiative fields in this limit. The only field at this limit contributing to $\mathcal{Q}$ is the scalar potential that propagates with infinite speed.
Such reasoning will also hold for nonabelian gauge theories where radiation could carry charges. Radiative fields in such theories still propagate with finite speed, so they could not contribute at the limit of $R\to\infty$, $t=t_0$.
Now, back to ADM energy
The argument presented above generalizes also to general relativity and asymptotically flat spacetimes. Let us formulate it for the full (nonlinear) GR using language  independent of specific coordinate systems.
We start with the spacelike hypersurface $\Sigma$ of our spacetime. On that hypersurface we have the induced metric $h_{ij}$ and sectional curvature $K_{ij}$ that must satisfy the Einstein constraint equations. Those equations could be seen as generalizations of Poisson's equation for scalar Newtonian and vector gravitomagnetic potentials.  ADM mass is defined as a limit of integrals over “spheres” of increasing radii (defined via the asymptotically Euclidean coordinates).  This limiting procedure corresponds to spacelike infinity $\mathcal{i}^0$. But there are no dynamic radiative degrees of freedom “living” at spacelike infinity: all gravitational radiation from an isolated source resides  at (future) null infinity $\mathcal{I}^{+}$. Therefore even in full nonlinear general relativity gravitational radiation would not be contributing to ADM mass, this result being not an artifact of linear approximation, but a direct consequence of finite speed of propagation of gravitational radiation and gauge freedom of general relativity leading to the constraint equations.
Note also the link between this argument and the “no-hair” theorems of the black hole physics: mass, angular momentum and charge are precisely the quantities that could be observed from “behind” the event horizon of a black hole. But these are the quantities that could be found from the fields that “propagate” instantaneously in a suitable gauge, by virtue of the corresponding constraint equation: scalar and vector Einstein constraints and Gauss law. While all the other “information”, that requires dynamical degrees of freedom to transmit ends up hidden behind the event horizon.
References:

*

*Wald, R.  General relativity. U. of Chicago press, 1984, ch. 10, 11.


Update. If one is interested in the definition of energy (–momentum) of spacetime that would reflect the loss of energy and momentum carried away by gravitational radiation, the necessary quantity would be the Bondi energy (and also momentum and mass). These quantities are defined using the asymptotic behavior along the null hypersurfaces reaching to the null infinity $\mathcal{I}^{+}$, which with a suitable choice of coordinate would be the surfaces of constant retarded time $u$.
