[proofplan]
We apply the $\operatorname{EVI}_\lambda$ inequality for each flow against the other flow, using a standard Steklov averaging argument to justify the fact that the comparison point is moving in time. The two EVI inequalities have opposite energy terms, so after adding them the energy contributions cancel. This yields a distributional differential inequality for $W_2^2(\rho_t,\eta_t)$, and multiplication by the integrating factor $e^{2\lambda t}$ gives monotonicity. Taking square roots gives the desired contractivity estimate.
[/proofplan]
[step:Introduce the squared distance between the two flows]
Define the function
\begin{align*}
F:[0,\infty)\to[0,\infty)
\end{align*}
by
\begin{align*}
F(t):=W_2^2(\rho_t,\eta_t).
\end{align*}
Since $\rho$ and $\eta$ are locally absolutely continuous curves in the metric space $(\mathcal P_2(\mathbb R^n),W_2)$, the function $F$ is locally absolutely continuous on $[0,\infty)$. It is therefore enough to prove that
\begin{align*}
F'(t)+2\lambda F(t)\le 0
\end{align*}
for $\mathcal L^1$-a.e. $t>0$, because the resulting one-dimensional differential inequality implies
\begin{align*}
F(t)\le e^{-2\lambda t}F(0)
\end{align*}
for every $t\ge 0$.
[/step]
[step:Regularize the moving comparison by averaging in two time variables]
For $s,t\ge 0$, define
\begin{align*}
D:(0,\infty)\times(0,\infty)\to[0,\infty)
\end{align*}
by
\begin{align*}
D(t,s):=W_2^2(\rho_t,\eta_s).
\end{align*}
Fix $0<a<b<\infty$ and choose $\varepsilon\in(0,a)$. Define the Steklov-averaged function
\begin{align*}
F_\varepsilon:[a,b]\to[0,\infty)
\end{align*}
by
\begin{align*}
F_\varepsilon(t):=\frac{1}{\varepsilon^2}\int_0^\varepsilon\int_0^\varepsilon D(t+r,t+q)\,d\mathcal L^1(r)\,d\mathcal L^1(q).
\end{align*}
The local absolute continuity of $\rho$ and $\eta$ implies that $D$ is locally absolutely continuous in each variable, and hence $F_\varepsilon$ is absolutely continuous on $[a,b]$. Whenever the derivative in the first variable exists, define $\partial_1D(t,s)$ to be the derivative of the one-variable map $u\mapsto D(u,s)$ at $u=t$. Whenever the derivative in the second variable exists, define $\partial_2D(t,s)$ to be the derivative of the one-variable map $v\mapsto D(t,v)$ at $v=s$.
[guided]
The EVI inequality is formulated with a fixed comparison measure. We cannot directly insert $\eta_t$ as the comparison point in the EVI for $\rho_t$ without justification, because $\eta_t$ also depends on time. The averaging argument avoids this problem by first comparing $\rho_{t+r}$ with the fixed measure $\eta_{t+q}$ and $\eta_{t+q}$ with the fixed measure $\rho_{t+r}$, where $r$ and $q$ are frozen inside the integral.
For $s,t>0$, define
\begin{align*}
D:(0,\infty)\times(0,\infty)\to[0,\infty)
\end{align*}
by
\begin{align*}
D(t,s):=W_2^2(\rho_t,\eta_s).
\end{align*}
This function records the squared distance between the two curves at two independent times. The independence of the two time variables is the key device: the EVI for $\rho$ differentiates in the first variable while the second variable is fixed, and the EVI for $\eta$ differentiates in the second variable while the first variable is fixed.
Fix $0<a<b<\infty$ and $\varepsilon\in(0,a)$. Define
\begin{align*}
F_\varepsilon:[a,b]\to[0,\infty)
\end{align*}
by
\begin{align*}
F_\varepsilon(t):=\frac{1}{\varepsilon^2}\int_0^\varepsilon\int_0^\varepsilon D(t+r,t+q)\,d\mathcal L^1(r)\,d\mathcal L^1(q).
\end{align*}
Because both curves are locally absolutely continuous in $(\mathcal P_2(\mathbb R^n),W_2)$, the map $(t,s)\mapsto W_2^2(\rho_t,\eta_s)$ is locally absolutely continuous in each variable on compact subintervals of $(0,\infty)$. Therefore the averaged function $F_\varepsilon$ is absolutely continuous on $[a,b]$. Whenever the derivative in the first variable exists, $\partial_1D(t,s)$ denotes the derivative of $u\mapsto D(u,s)$ at $u=t$. Whenever the derivative in the second variable exists, $\partial_2D(t,s)$ denotes the derivative of $v\mapsto D(t,v)$ at $v=s$. These partial derivatives are locally integrable in their respective variables on compact rectangles because they are metric derivatives of locally absolutely continuous one-variable functions. The role of the averaging is to make the subsequent differentiation legitimate before passing to the diagonal limit $r,q\downarrow 0$.
[/guided]
[/step]
[step:Add the two EVI inequalities so the energy terms cancel]
Let $A_\rho:=A_\rho((a,b+\varepsilon))$ and $A_\eta:=A_\eta((a,b+\varepsilon))$ be the full-measure Borel sets from the strong EVI hypothesis. Since $\rho_{u},\eta_{v}\in\operatorname{Dom}(\mathcal E)$ for every $u,v>0$, the comparison measures $\eta_{t+q}$ and $\rho_{t+r}$ are admissible whenever $t+r,t+q\in(a,b+\varepsilon)$. By Fubini's theorem applied to the full-measure sets $\{(t,r,q):t+r\in A_\rho\}$ and $\{(t,r,q):t+q\in A_\eta\}$ in $[a,b]\times(0,\varepsilon)^2$, both EVI inequalities below hold for $\mathcal L^3$-a.e. triple $(t,r,q)$.
For such a triple, the $\operatorname{EVI}_\lambda$ inequality for $\rho$ with fixed comparison point $\eta_{t+q}$ gives
\begin{align*}
\partial_1 D(t+r,t+q)+\lambda D(t+r,t+q)\le 2\mathcal E[\eta_{t+q}]-2\mathcal E[\rho_{t+r}].
\end{align*}
The $\operatorname{EVI}_\lambda$ inequality for $\eta$ with fixed comparison point $\rho_{t+r}$ gives
\begin{align*}
\partial_2 D(t+r,t+q)+\lambda D(t+r,t+q)\le 2\mathcal E[\rho_{t+r}]-2\mathcal E[\eta_{t+q}].
\end{align*}
Adding these two inequalities cancels the two finite energy terms and yields
\begin{align*}
\partial_1D(t+r,t+q)+\partial_2D(t+r,t+q)+2\lambda D(t+r,t+q)\le 0.
\end{align*}
The local integrability of the partial derivatives and the nonnegativity of $D$ make the averaged terms integrable on $[a,b]\times(0,\varepsilon)^2$. Differentiating $F_\varepsilon$ under the integral sign for absolutely continuous functions and using Fubini's theorem gives, for $\mathcal L^1$-a.e. $t\in[a,b]$,
\begin{align*}
F_\varepsilon'(t)+2\lambda F_\varepsilon(t)\le 0.
\end{align*}
[guided]
The point of the strengthened EVI hypothesis is that the exceptional set of times does not depend on the comparison measure. Let $A_\rho:=A_\rho((a,b+\varepsilon))$ and $A_\eta:=A_\eta((a,b+\varepsilon))$. These are Borel subsets of $(a,b+\varepsilon)$ with full $\mathcal L^1$-measure. Since $\rho_u,\eta_v\in\operatorname{Dom}(\mathcal E)$ for every positive time, the measures $\eta_{t+q}$ and $\rho_{t+r}$ are valid comparison measures whenever $t+r,t+q\in(a,b+\varepsilon)$.
We now justify the phrase "for almost every triple." The set of triples $(t,r,q)\in[a,b]\times(0,\varepsilon)^2$ for which $t+r\notin A_\rho$ has zero $\mathcal L^3$-measure by Fubini's theorem, because for each fixed $(r,q)$ the bad set of $t$ is a translate of a null subset of $(a,b+\varepsilon)$. The same argument gives zero $\mathcal L^3$-measure for the triples with $t+q\notin A_\eta$. Hence, outside a null set of triples, both EVI inequalities may be applied simultaneously.
For such a triple, applying the EVI for $\rho$ at the time $t+r$ with fixed comparison point $\eta_{t+q}$ gives
\begin{align*}
\partial_1 D(t+r,t+q)+\lambda D(t+r,t+q)\le 2\mathcal E[\eta_{t+q}]-2\mathcal E[\rho_{t+r}].
\end{align*}
Applying the EVI for $\eta$ at the time $t+q$ with fixed comparison point $\rho_{t+r}$ gives
\begin{align*}
\partial_2 D(t+r,t+q)+\lambda D(t+r,t+q)\le 2\mathcal E[\rho_{t+r}]-2\mathcal E[\eta_{t+q}].
\end{align*}
The two right-hand sides are finite because both comparison measures lie in $\operatorname{Dom}(\mathcal E)$. Adding the two displayed inequalities cancels the energy terms and gives
\begin{align*}
\partial_1D(t+r,t+q)+\partial_2D(t+r,t+q)+2\lambda D(t+r,t+q)\le 0.
\end{align*}
Because $D$ is locally absolutely continuous in each variable, the partial derivatives appearing here are locally integrable on the compact averaged region. Therefore Fubini's theorem permits integration of the inequality over $(r,q)\in(0,\varepsilon)^2$, and the standard differentiation rule for Steklov averages of absolutely continuous functions gives
\begin{align*}
F_\varepsilon'(t)=\frac{1}{\varepsilon^2}\int_0^\varepsilon\int_0^\varepsilon \bigl(\partial_1D(t+r,t+q)+\partial_2D(t+r,t+q)\bigr)\,d\mathcal L^1(r)\,d\mathcal L^1(q)
\end{align*}
for $\mathcal L^1$-a.e. $t\in[a,b]$. Averaging the preceding pointwise inequality over $(r,q)$ gives
\begin{align*}
F_\varepsilon'(t)+2\lambda F_\varepsilon(t)\le 0.
\end{align*}
[/guided]
[/step]
[step:Pass from the averaged inequality to the diagonal inequality]
Since $\rho$ and $\eta$ are continuous as maps into $(\mathcal P_2(\mathbb R^n),W_2)$, for every $t\in[a,b]$ we have
\begin{align*}
F_\varepsilon(t)\to F(t)
\end{align*}
as $\varepsilon\downarrow0$. The convergence is uniform on $[a,b]$ because $\rho$ and $\eta$ are uniformly continuous on compact time intervals. Let $\varphi\in C_c^\infty((a,b);[0,\infty))$ be a nonnegative test function. The distributional form of
\begin{align*}
F_\varepsilon'(t)+2\lambda F_\varepsilon(t)\le 0
\end{align*}
is
\begin{align*}
-\int_a^b F_\varepsilon(t)\varphi'(t)\,d\mathcal L^1(t)+2\lambda\int_a^b F_\varepsilon(t)\varphi(t)\,d\mathcal L^1(t)\le 0.
\end{align*}
The uniform convergence $F_\varepsilon\to F$ on $[a,b]$ lets us pass to the limit in both Lebesgue integrals, giving
\begin{align*}
-\int_a^b F(t)\varphi'(t)\,d\mathcal L^1(t)+2\lambda\int_a^b F(t)\varphi(t)\,d\mathcal L^1(t)\le 0.
\end{align*}
Thus
\begin{align*}
F'(t)+2\lambda F(t)\le 0
\end{align*}
in the sense of distributions on $(a,b)$. Since $0<a<b<\infty$ were arbitrary and $F$ is locally absolutely continuous, the same inequality holds for $\mathcal L^1$-a.e. $t>0$.
[/step]
[step:Apply the integrating factor and take square roots]
Define
\begin{align*}
G:[0,\infty)\to[0,\infty)
\end{align*}
by
\begin{align*}
G(t):=e^{2\lambda t}F(t).
\end{align*}
For $\mathcal L^1$-a.e. $t>0$, the product rule for absolutely continuous functions gives
\begin{align*}
G'(t)=e^{2\lambda t}\bigl(F'(t)+2\lambda F(t)\bigr)\le 0.
\end{align*}
Hence $G$ is nonincreasing on $[0,\infty)$, and so for every $t\ge 0$,
\begin{align*}
e^{2\lambda t}F(t)\le F(0).
\end{align*}
Substituting the definition of $F$ gives
\begin{align*}
W_2^2(\rho_t,\eta_t)\le e^{-2\lambda t}W_2^2(\rho_0,\eta_0).
\end{align*}
Both sides are nonnegative, so taking square roots yields
\begin{align*}
W_2(\rho_t,\eta_t)\le e^{-\lambda t}W_2(\rho_0,\eta_0).
\end{align*}
This is the desired contractivity estimate.
[/step]