$\newcommand{\R}{\mathbb R}\newcommand{\tr}{\operatorname{tr}}$The answer is yes.
Indeed, it is easy to see (cf. e.g. Proposition 7 or the beginning of its proof) that the Wasserstein distance between $N(a,A)$ and $N(b,B)$ is
\begin{equation*}
W_2(N(a,A),N(b,B))=\sqrt{\|a-b\|^2+W_2(N(0,A),N(0,B))^2}, \tag{10}\label{10}
\end{equation*}
where $\|\cdot\|$ is the Euclidean norm. So,
\begin{equation*}
\|a-b\|\le W_2(N(a,A),N(b,B)). \tag{20}\label{20}
\end{equation*}
Let now $X\sim N(0,A)$ and $Y\sim N(0,B)$. Then for any unit vector $u\in\R^n=\R^{n\times1}$
\begin{equation*}
E\|X-Y\|^2\ge E(u^\top X-u^\top Y)^2\ge(\sqrt{u^\top A u}-\sqrt{u^\top B u})^2,
\end{equation*}
where the last inequality holds (say, by mentioned Proposition 7) because $u^\top X$ and $u^\top Y$ are real-valued Gaussian zero-mean random variables with variances $u^\top A u$ and $u^\top B u$.
So, by \eqref{10} and the definition of the $W_2$-distance, for any unit vector $u\in\R^n=\R^{n\times1}$,
\begin{equation*}
\begin{aligned}
W_2(N(a,A),N(b,B))&\ge W_2(N(0,A),N(0,B)) \\
&\ge|\sqrt{u^\top A u}-\sqrt{u^\top B u}| \\
&=\frac{|u^\top A u-u^\top B u|}{\sqrt{u^\top A u}+\sqrt{u^\top B u}} \\
&\ge\frac{|u^\top(A-B)u|}{\sqrt{\|A\|}+\sqrt{\|B\|}} \\
&=\frac{\|A-B\|}{\sqrt{\|A\|}+\sqrt{\|B\|}}
\end{aligned}
\end{equation*}
for some unit vector $u\in\R^n=\R^{n\times1}$, where $\|M\|$ is the spectral norm of a matrix $M$.
So,
\begin{equation*}
\|A-B\|\le(\sqrt{\|A\|}+\sqrt{\|B\|\,})\, W_2(N(a,A),N(b,B)). \tag{30}\label{30}
\end{equation*}
Thus, by \eqref{20} and \eqref{30}, we have the desired local Lipschitz property. $\quad\Box$
Remark: Inequality \eqref{20} is "exact in the limit", when $A$ and $B$ are each (close to) the zero matrix.
Inequality \eqref{30} turns into the equality when $n=1$ to, more generally, "exact in the limit" when $A$ and $B$ are (close to) commuting matrices of rank $1$ each. $\quad\Box$
User ABIM wrote in a comment: "@Justin_other_PhD OP claimed that $f$
is bi-Lipschitz but why is the upper-bound true?"
Let me answer this question as well. First of all, it is not true that $f$ is Lipschitz (and the OP did not claim that). Indeed, even when $n=1$, we have $W_2(N(0,A),N(0,B))=|\sqrt A-\sqrt B|$, so that there is no real $L>0$ such that $W_2(N(0,A),N(0,B))\le L|A-B|$ for all real $A,B>0$.
What the OP said, and what is true, is that $f$ is locally Lipschitz. Indeed,
let $X\sim N(0,A)$ and $Y:=B^{1/2}A^{-1/2}X$. Then $Y\sim N(0,B)$ and hence
\begin{equation*}
\begin{aligned}
W_2(N(0,A),N(0,B))^2&\le E\|X-Y\|^2 \\
&=\tr(A+B-B^{1/2}A^{1/2}-A^{1/2}B^{1/2}) \\
&=\tr[(A^{1/2}-B^{1/2})^2]=\|A^{1/2}-B^{1/2}\|_F^2 \\
&\le n\,\|A^{1/2}-B^{1/2}\|^2,
\end{aligned}
\tag{40}\label{40}
\end{equation*}
where $\tr$ denotes the trace and $\|\cdot\|_F$ is the Frobenius norm.
Next, for real $a\ge0$,
\begin{equation*}
a^{1/2}=\frac1\pi\int_0^\infty dt\,t^{-1/2}\Big(1-\frac t{t+a}\Big)
\end{equation*}
and hence the value of the derivative of the function $M\mapsto M^{1/2}$ at $A$ at a symmetric matrix $H$ is
\begin{equation*}
\begin{aligned}
(A^{1/2})'(H)&=\frac1\pi\int_0^\infty dt\,t^{1/2}(tI+A)^{-1}H(tI+A)^{-1} \\
& \le\frac1\pi\int_0^\infty dt\,t^{1/2}(t+c)^{-2}|H|
=\frac1{2c^{1/2}}\,|H|
\end{aligned}
\end{equation*}
provided that $A\ge cI$ for some $c\in(0,\infty)$ (where $I$ is the identity matrix), so that the operator norm of $(A^{1/2})'$ is $\le\dfrac1{2c^{1/2}}$.
So, by \eqref{10} and \eqref{40},
\begin{equation*}
W_2(N(a,A),N(b,B))\le\|a-b\|+W_2(N(0,A),N(0,B)) \\
\le \|a-b\|+\dfrac{n^{1/2}}{2c^{1/2}}\,\|A-B\|
\end{equation*}
provided that $A\ge cI$ and $B\ge cI$ for some $c\in(0,\infty)$. Thus, $f[=N]$ is locally Lipschitz. $\quad\Box$
