What do we mean when we say gravitational waves are non-linear and do not superpose like EM waves? I have read this question:

Now it's not actually true that general relativity obeys a law of superposition, but it is an extremely good approximation for a small-amplitude gravitational wave passing through the static curvature of an object like the earth.

Why does spacetime propagate gravitational waves?
And this one:

Gravity as described by general relativity is highly non-linear. Therefore it does not have any superposition principle. Gravitational waves do not have a superposition principle. However, at very large distances these waves can be approximated. And then this operator might be linear and you can reasonable speak of superpositions again.

Do all waves of any kind satisfy the principle of superposition?
My question is this, how does non-linearity cause gravitational waves not to superpose? The answer seems to be non-trivial, and somehow makes GWs special relative to EM waves.
Question:

*

*What do we mean when we say gravitational waves are non-linear and do not superpose like EM waves?

 A: The full Einstein equation is:
$$ G_{\alpha\beta} = 8\pi T_{\alpha\beta} $$
where $G_{\alpha\beta}$ is the Einstein tensor, which is a non-linear function of the metric. In general solving the equation directly is impossibly difficult.
However if the spacetime is approximately flat we can write the metric as:
$$ g_{\alpha\beta} = \eta_{\alpha\beta} + h_{\alpha\beta} $$
where $\eta_{\alpha\beta}$ is the flat spacetime metric and $h_{\alpha\beta}$ is a small perturbation. The importance of $h_{\alpha\beta}$ being small is that higher powers of $h_{\alpha\beta}$ are small compared to $h_{\alpha\beta}$ and can be neglected. If we do this we end up with the linearised Einstein equation, which can be written (in the Lorenz gauge) as:
$$ \Box \bar h^{\alpha\beta} = -16\pi T^{\alpha\beta}$$
where:
$$\bar{h}^{\alpha\beta}=h^{\alpha\beta}-\frac{1}{2}\eta^{\alpha \beta} h $$
Remember this is an approximation that works only when $h^{\alpha\beta}$ is small, but provided this is the case the approximation is very useful as it is much, much easier to solve than the full equation. In particular if we set $T_{\alpha\beta} = 0$ we get gravitational waves as solutions. Since the linearised Einstein equation is, well, linear, the solutions obey the superposition principle just like EM waves.
The problem is that this only works when higher powers of $h_{\alpha\beta}$ are so small they can be ignored. In practice this is always true for the sort of gravitational waves detected by LIGO, but very near the source of the waves the curvature would be so large that $h^2$, $h^3$ etc cannot be neglected and would need to be added to the linearised Einstein equation to keep the approximation good. When you do this the equation is no longer linear and the solutions no longer obey the superposition principle.
A: This simple definition is axiomatic:

The principle of superposition says:
When two or more waves cross at a point, the displacement at that point is equal to the sum of the displacements of the individual waves.
The individual wave displacements may be positive or negative. If the displacements are vectors, then the sum is calculated by vector addition.
'''
When the waves pass beyond a point of intersection, they separate out again and are unaffected.

Italics mine.
You ask:

What do we mean when we say gravitational waves are non-linear and do not superpose like EM waves?

That when crossing they do not obey "the displacement at that point is equal to the sum of the displacements of the individual waves." This is due to the mathematical solution of the gravitational wave in general relativity.
I like this definition of a wave!

What is a wave?
In the mathematical sense, a wave is
any function that moves.
To displace any function $f(x)$
to the right, just change its
argument from $x$ to $x-x_0$,
where $x_0$ is a positive number.


It includes all periodic phenomena discussed in the links you give.
For me, it is easy to understand that only when approximations to linearity work will one see simple phenomena of diffraction, etc. Now for gravitational waves, approximations can be used to fit observations, if interference phenomena are measured. The strict meaning of nonlinearity is that the sum of the individual wave disturbances at a point are not added vectorially.
A: The ordinary wave equation
$$\frac{\partial^2 u}{\partial t^2} = v_w^2 \frac{\partial^2 u}{\partial x^2}$$
where $u$ is the wave function and $v_w$ is the wave speed, is linear (in the sense of a linear operator - or "homogeneous and linear" in the sense of differential equations), in that if $u_1$ and $u_2$ - both functions of time as well as space - are solutions, then any weighted combination, i.e. $u = a_1 u_1 + a_2 u_2$, for constants $a_1$ and $a_2$, and that is to say, whose values are given by $u(t, x) = a_1 u_1(t, x) + a_2 u_2(t, x)$, is also a solution. This can be seen by noting that each side consists only of differentiations and no products, so that the expressions will "distribute" nicely over the sum.
In general relativity, if you have two metric waves $g_1$ and $g_2$, now tensor-valued instead of scalar-valued as above (so a "tensor wave", c.f. how EM waves are "vector waves"), while $g = a_1 g_1 + a_2 g_2$ may be another space-time metric tensor field, it is not a solution of the corresponding nonlinear wave equation, meaning that while "at any given spatial slice" it may be a valid instantaneous configuration of the geometry field, the propagation is all wrong.
That is because that the wave equations for gravitational waves are nonlinear wave equations. A simple example of such a wave equation would be
$$\frac{\partial^2 u}{\partial t^2} + v_w^2 \frac{\partial^2 u}{\partial x^2} + k u^2$$
for some $k > 0$. If you have $u_1$ and $u_2$ as solutions, then $u = u_1 + u_2$ is not a solution, for while at any given time it may be okay as a configuration of the wave system, the future evolution of that configuration is no longer merely the sum of the evolutions of the components - the squaring part introduces cross terms between the two, namely $2 u_1 u_2$, which cannot be created simply by taking the time derivative on the left, for it just distributes right across the sum, and introduces no new cross terms.
More specifically, the Einstein field equations are an equation in the Ricci tensor $R_{\mu \nu}$, and each component there takes (lots and lots of) products of partial derivatives of the metric tensor components. Hence they are all cross terms.
As to why this nonlinearity exists, you should note that even from elementary calculus, curvature of a simple curve given by a differentiable real valued function of one real variable is a nonlinear property in that simply adding the two curves won't add the curvatures. The Ricci tensor is a curvature for a general manifold of any dimensionality and not just such a curve, so it should be expected likewise to be similarly nonlinear in the "shape" thereof (here given by the metric field).
