One class of examples in a compact, connected Lie group $G$ of rank $r$ are the conjugacy classes of codimension $r$. Taking an element $a\in G$ whose centralizer is a maximal torus (a generic condition), the conjugacy class of $a$ is a submanifold $C_a\subset G$ of codimension $r$ whose normal plane at $a$ is the tangent at $a$ to $Z_a$, the centralizer of $a$ in $G$, which is a flat, totally geodesic submanifold. By the $G$-homogeneity of $C_a = G/Z_a\subset G$, the fact that this holds at $a$ implies that it holds at all point of $C_a$.
This gives an $r$-parameter family of mutually noncogruent examples.
There are other examples: When $G=\mathrm{SU}(n)$, and $a = \mathrm{diag}(\lambda_1,\lambda_2,\ldots,\lambda_n)$ is a diagonal element for which the $\lambda_i^2$ are all distinct, the submanifold
$$
M_a = \{\ g_1 a g_2\ |\ g_1,g_2\in\mathrm{SO}(n)\ \}\subset\mathrm{SU}(n)
$$
is homogeneous under the isometry group of $\mathrm{SU}(n)$ (endowed with its biïnvariant metric), and its tangent plane at $a$ is orthogonal to the diagonal maximal torus $T\subset\mathrm{SU}(n)$, so, by homogeneity, its tangent plane at any point is orthogonal to a flat, totally geodesic submanifold of $\mathrm{SU}(n)$. Hence it has an abelian normal bundle.
This gives another $r=n{-}1$ parameter family of mutually noncongruent examples of codimension $r$, distinct from the conjugacy classes.
Using the methods of exterior differential systems, one can show that, when $n=3$, these two families account for all of the codimension $2$ submanifolds of $\mathrm{SU}(3)$ with abelian normal bundle, in the sense that any connected submanifold $M^6\subset\mathrm{SU}(3)$ with abelian normal bundle is, up to ambient isometry, an open subset of one of the examples listed above. The argument that I have written out is a calculation using exterior differential systems and the moving frame, but, when I have time, I can sketch the proof, if there is interest.
Addendum 1: I had a flight with some time to look at the other two compact simple rank 2 groups. It turns out that all of the connected codimension two submanifolds with abelian normal bundle in $\mathrm{SO}(5)$ and $\mathrm{G}_2$ are (open subsets of) homogeneous compact ones as well. In each case, there is one additional $2$-parameter family of examples beyond the principal conjugacy classes.
Addendum 2: I also checked the rank $r=3$ case $G = \mathrm{SU}(4)$, and found that every connected codimension $3$ submanifold of $G$ with abelian normal bundle is an open subset of one of the two types of homogeneous examples listed above. Since any codimension $2$ submanifold of $\mathrm{SU}(4)$ that is foliated by codimension $3$ submanifolds of $\mathrm{SU}(4)$ with abelian normal bundle will have abelian normal bundle, it follows that there are many non-homogeneous codimension $2$ submanifolds of $\mathrm{SU}(4)$ with abelian normal bundle.
In fact, it now seems likely that, for any compact simple group $G$ of rank $r>1$, there are two $r$-parameter families of homogeneous codimension $r$ submanifolds with abelian normal bundle and every connected codimension $r$ submanifold with abelian normal bundle is, up to ambient isometry, an open subset of a homogeneous one. The two families are as follows: The first family is the family of principal conjugacy classes in $G$, and the second family is the set of principal orbits of $K\times K$ acting by left and right multiplication in $G$ where $G/K$ is a symmetric space of rank $r$. For example, when $G=\mathrm{SU}(n)$ $(n\ge3)$, $K=\mathrm{SO}(n)$; when $G = \mathrm{SO}(n)$ $(n\ge 5)$, $K = \mathrm{SO}(p)\times\mathrm{SO}(n{-}p)$ where $p = \bigl[\tfrac12 n\bigr]$; when $G=\mathrm{Sp}(n)$ $(n\ge3)$, $K=\mathrm{U}(n)$; when $G=\mathrm{G}_2$, $K=\mathrm{SO}(4)$; when $G=\mathrm{F}_4$, $K=\mathrm{Sp}(3)\mathrm{SU}(2)$; when $G=\mathrm{E}_6$, $K=\mathrm{Sp}(4)$; when $G=\mathrm{E}_7$, $K=\mathrm{SU}(8)$, and when $G=\mathrm{E}_8$, $K=\mathrm{SO}'(16)$.
If you take the simplest case, in which $Q$ is a $3$-dimensional Riemannian manifold, then there are plenty of extrinsically flat surfaces $M\subset Q$. In fact, if one chooses a 'generic' curve $\gamma$ in $Q$ and a 'generic' normal vector field $\nu$ along $\gamma$, then there will be a unique 'extinsically flat' surface $M$ containing $\gamma$ that has $\nu$ as its normal along $\gamma$.
In Cartan's language, the extrinsically flat surfaces in $Q$ depend on two functions of one variable. This does not depend on any knowledge of the curvature of $Q$.
Addendum 1: (7/23/22) I thought a little more about this and realized that there is a natural conjecture about the existence of 'extrinsically flat' submanifolds. Here is the statement.
Conjecture: For each $r\ge2$, let $D_r = r^2(r^2{-}1)/12$ be the rank of the Riemann curvature tensor in dimension $r$. For a Riemannian manifold $Q$ of dimension $d\ge r + D_r$, the PDE for 'extrinsically flat' submanifolds of dimension $r$ is involutive. In particular, such local submanifolds are just as 'plentiful' in $Q$ as they are in flat $\mathbb{R}^d$.
The intuition for the Conjecture is this: 'Extrinsic flatness of $M\subset Q$ is the requirement that, at each point $x\in M$, the Riemann curvature tensor of $M$ at $x$ should be equal to the restriction of the Riemann curvature tensor of $Q$ at $x$ to $T_xM$. This is $D_r$ second order PDE on $M$ as a submanifold of $Q$. Since the $r$-dimensional submanifolds of $Q$ can be represented locally as graphs of $d{-}r$ functions of $r$ variables, this PDE system will be formally 'determined' if $d{-}r = D_r$ and 'overdetermined' if $d{-}r<D_r$. If the symbol of the PDE system is 'non-degenerate', in the appropriate sense, when $d{-}r = D_r$, then, at least in the real-analytic category, there will be 'plenty' of local solutions. In fact one would expect the local generality of the solutions in this case to be $2D_r$ functions of $r{-}1$ variables. Meanwhile, for a 'generic' $Q$ of dimension $d<r + D_r$, one would expect that it would not have any 'extrinsically flat' submanifolds of dimension $r$.
The above conjecture is easily established for $r=2$. (The case $(r,d)=(2,3)$ is just the statement in the first paragraph.) The case $(r,d)=(3,9)$ (note that $9 = 3 + D_3$) seems to work, but I haven't checked all of the details, as the algebra is a little tricky.
Also, note that, for the 'overdetermined' non-existence result mentioned above, it is essential that one assume that the metric on $Q$ be 'generic'. As Cartan showed, if $d=2r$ and $Q$ has constant sectional curvature $c$, then the PDE for $r$-dimensional submanifolds of constant sectional curvature $c$ is involutive, and $2r < r + D_r$ when $r>2$.
Best Answer
Try pages 71 and 72 of Cartan's "La géométrie des groupes de transformations" https://eudml.org/doc/235668