I eventually found a derivation. The issue is finding the right relations to use for the metric tensor since they are so many. Instead of the notation used in the paper I used $q_{ij}$ to denote the 3-metric and $g_{\mu\nu}$ for the 4-dimensional one. Here it is in case someone else needs it:
From 3.9 a) we immediately obtain $g^{00}=-\frac{1}{N^2}$.
From the orthogonality relation $g^{\mu0}g_{\mu0}=1$, we obtain
\begin{equation}
g^{00}g_{00}+g^{i0}g_{i0}=1\implies -\frac{1}{N^2}g_{00} + g^{i0}N_i=1
\end{equation}
\begin{equation}
\implies g_{00}=-N^2(1-g^{i0}N_i).
\end{equation}
We now use $g_{\mu i}g^{\mu 0}=0$ to get
\begin{equation}
g_{0i}g^{00}+g_{ji}g^{j0}=0\implies N_i(-\frac{1}{N^2})+q_{ij}g^{j0}=0.
\end{equation}
Multiplying by the inverse 3-metric $q^{ik}$ and summing over $i$ gives
\begin{equation}
N^k(-\frac{1}{N^2})+\delta^k_jg^{j0}=0
\end{equation}
where $N^k=q^{ik}N_i$. Then
\begin{equation}
N^k(-\frac{1}{N^2})+g^{k0}=0
\end{equation}
\begin{equation}
\implies g^{k0}=\frac{N^k}{N^2}.
\end{equation}
This leads to
\begin{equation}
g_{00}=-(N^2-N^iN_i)
\end{equation}
Further, we have
\begin{equation}
g^{\mu j}g_{\mu i}=\delta^j_i \implies g^{0j}g_{0i}+g^{kj}g_{ki}=\delta^j_i \implies g^{0j}N_i + g^{kj}q_{ki}=\delta^j_i
\end{equation}
\begin{equation}
\implies g^{kj}q_{ki}=\delta^j_i - g^{0j}N_i
\end{equation}
Multiplying by $q^{il}$ and summing over $i$ gives
\begin{equation}
g^{kj}\delta^l_k=q^{jl}-g^{0j}N^l
\end{equation}
and we finally get, relabelling the indices $l \to i$ and using the symmetry property of the metric tensor,
\begin{equation}
g^{ij}=q^{ij}-\frac{N^iN^j}{N^2}.
\end{equation}