$\newcommand{\M}{\mathcal{M}}$
$\newcommand{\ep}{\epsilon}$
I will sketch here a different approach than the one given by Robert's answer. (It is loosely based on Alex's answer). We want to prove that around every point $p \in \M$ there exist a local frame for $\bigwedge^k(T^*\M)$ whose elements are closed and co-closed.

For the Euclidean metric this is immediate: We have the standard (constant) frame $dx^I=dx^{i_1} \wedge \ldots dx^{i_k} $. Since every metric is locally close to being Euclidean on small neighbourhoods, the idea is to use an approximation argument:

Given a Riemannian metric $g$, we denote the space of $g$-harmonic forms of degree $k$ by $H^k_{g}$.

We view $H^k_{g}$, as a subspace of $\Omega^k(\M)$ which is "changing continuously" with the metric $g$. Suppose $g_{\ep} \to g_0$ in the $C^1$-norm where $g_0$ is the Euclidean metric; Then $H^k_{g_{\ep}} \to H^k_{g_0}$ in the following sense: there exist a family of bases of $H^k_{g_{\ep}}$, which converges to a basis of $H^k_{g_{0}} $ in $C^1$; this basis of $H^k_{g_{0}}$ forms a local frame for $\bigwedge^k(T^*\M)$. Since being a frame is an open condition, those bases for $H^k_{g_{\ep}}$ are local frames for sufficiently small $\ep$.

For the full details, see Appendix A in my paper here.

*Some more details:*

Even though the claim is local, and the approximation scheme is also inspired by a local phenomena, the implementation of the proof is based on a combination of local and global arguments. The reason is that on a **closed manifold**, being closed and co-closed is equivalent to being harmonic, and the dimension of the space of harmonic forms is a finite number which is a topological invariant of the manifold; it does not depend on the chosen metric.

Thus, given a family of metrics $g_{\ep} \to g_0$ on a closed manifold $\M$, we consider the behaviour of the finite-dimensional subspaces $H^k_{g_{\ep}}$ (all of the same dimension) as $\ep \to 0$.

That is, we look at the map $g \to H^k_{g}=\ker \Delta_g$. It turns out that this map is continuous in some appropriate sense; this relies on a certain **"stability property of kernels of linear operators"**. It turns out that a crucial factor in the existence of such a stability phenomenon is the assumption that all the kernels have the same finite dimension. The convergence of kernels does not always hold when the dimensions are not equal or infinite.