I don't know of an exact reference, and in general this sort of result (transfering a "classical" result from the analysis of PDE in $\mathbb{R}^n$ to Riemannian manifolds) is often quite hard to track down. Instead, I can explain how you derive the inequalities you want from the standard Calderon-Zygmund estimates for domains in $\mathbb{R}^n$. What follows is quite long for a MathOverflow answer, but I've tried to be as complete as possible. It is a "standard" patching argument to turn estimates from $\mathbb{R}^n$ into estimates on a manifold. I don't remember where I learned this sort of argument from (maybe Wells' book on analysis on complex manifolds?) but it is ubiquitous in geometric analysis.
Take a deeep breath, and here we go.
There are two important hypotheses which mean the result holds. Firstly the system you are interested in is elliptic. Secondly there is no kernel, i.e. no solutions to the homogeneous equation. I will explain the result in this sort of generality (linear elliptic systems over a compact Riemannian manifold) since it's less notation and no more extra thought.
So, let $(M,g)$ be a compact Riemannian manifold with vector bundles $E, F \to M$ and let $L$ be a linear differential operator taking sections of $E$ to sections of $F$. In your case, $M=S^2$ is the shpere, $E=T^*S^2$ is the cotangent bundle and $F$ is the sum of two real line bundles: $F = \underline{\mathbb{R}} \oplus \Lambda^2T^*S^2$. Your operator is:
$$
L(\alpha) = (\mathop{div}(\alpha), \mathop{curl}(\alpha))
$$
We aassume that $L$ is elliptic. This means that the principal symbol in any given cotangent direction is an isomorphism $E \to F$. This is the case for your $L$ because the symbol in the direction $\beta$ is given by
$$
\sigma_L(\beta)(\alpha) = ( g(\alpha,\beta), \beta \wedge \alpha)
$$
(To see that $\sigma_L$ is an isomorphism it's enough to check it has no kernel.)
The first estimate I will explain is that for any elliptic $L$, any $k \in \mathbb{N}$ and any $p \in (1,\infty)$, there is a constant $C$ (which depends on $L,k,p$) such that for any section $s$ of $E$,
$$
| s |_{L^p_{k+r}} \leq C \left( | L(s)|_{L^p_k} + |s |_{L^p}\right)
$$
Here $r$ is the order of $L$ (so $r=1$ in your cases). Meanwhile, $L^p_k$ means the Sobolev norm given by taking the sum of the $L^p$ norms of derivatives up to order $k$, all defined using the Riemannian metric $g$. We also need to fix a choice of fibrewise innerproduct in both $E$ and $F$,as well as connections in the bundles. (Your question only asks about $k=1,2$.)
To prove this estimate, use a partition of unity to write $s = \sum \phi_i s$ where each $\phi_i$ is supported in a small coordinate chart. Suppose moreover that $E$ and $F$ are both trivial over each of these coordinate patches. Now on any single coordinate chart, $\Omega_i$, you can treat everything as if it were taking place on vector valued functions on $\mathbb{R}^2$. Apply the classical Calderon-Zygmund inequality to obtain:
$$
| \phi_i s |_{W^p_{k+r}(\Omega_i)} \leq C| L(\phi_i s)|_{W^p_k(\Omega_i)}
$$
where here $W^p_k(\Omega_i)$ means using the Euclidean $(p,k)$-Sobolev norm for vector valued functions on $\Omega_i \subset \mathbb{R}^n$. It is important here that $\phi_i s$ is compactly supported in $\Omega_i$.
To get a global inequality over the whole of $M$, we need to do two things: Firstly, we need to compare the $W^p_k(\Omega_i)$ norm and $L^p_k(\Omega_i)$ norm for sections supported in $\Omega_i$. It turns out that there is a constant $K>0$ such that for any section of $E$ or $F$ supported in $\Omega_i$,
$$
\frac{1}{K} | s |_{L^p_k(\Omega_i)} \leq | s|_{W^p_k(\Omega_i)} \leq K | s|_{L^p_k(\Omega_i)}.
$$
The constant $K$ will depend on $\Omega_i$, because we need to compare $g$ to the Euclidean metric on each chart $\Omega_i$ (and also our connections on $E$ and $F$ to the trivial connection). When $M$ is compact, there are only finitely many charts, so we can use the same $K$ for all of them. (When $M$ is not compact, one needs to assume some uniform geometric bounds here.)
From here we see that for any section $s$ at all,
$$
| s|_{L^p_{k+r}(M)}
\leq K \sum | \phi_i s|_{W^p_{k+r}(\Omega_i)} \leq CK^2 \sum | L(\phi_i s)|_{W^p_k(\Omega_i)}
$$
where the second inequality is Calderon-Zygmund in $\Omega_i$.
Next, we need to relate $L(\phi s)$ to $\phi_i L(s)$. We do this via the Leibniz rule. The difference $L(\phi_i s) - \phi_i L(s)$ is an expression in which we see terms with at most $r-1$ derivatives of $s$. This means that (with a little more work),
$$
| L(\phi_i s) |_{W^p_k(\Omega_i)} \leq B \left(| \phi_i L(s)|_{W^p_k(\Omega_i)} + | s |_{W^p_{r+k-1}(\Omega_i)}\right)
$$
Multiplication by a smooth function is a bounded map on Sobolev spaces and so
$$
| \phi_i L(s)|_{W^p_k(\Omega_i)} \leq A |L(s)|_{W^p_k(\Omega_i)}
$$
Putting this together we find that
$$
| s|_{L^p_{k+r}(M)} \leq D \sum \left(|L(s)|_{L^p_k(\Omega_i)} + |s|_{L^p_{k+r-1}(\Omega_i)}\right)
$$
where $D$ combines all our previous constants $C,K,A,B$.
We also choose our coordinate charts to be locally finite: any point of $M$ is in at most $m$ of the $\Omega_i$. This means that
$$
\sum |L(s)|_{L^p_k}(\Omega_i) \leq m |L(s)|_{L^p_k(M)}
$$
and similarly for the $L^p_{k+r-1}$ term. So we conclude that
$$
| s|_{L^p_{k+r}(M)} \leq C \left( |L(s)|_{L^p_k(M)} +|s|_{L^p_{k+r-1}}\right)
$$
for some new constant $C$. We now get the estimate we really want (with $L^p_{k+r-1}$ replaced by $L^p$ in the last term) by induction on $k$. (When $r=1$ the base of the induction is simple enough, but when $r>1$ it needs a little more thought, which I leave as an exercise - code for I didn't have time to think about this bit very much!)
This is almost the estimate you want in your original quesiton - which asked about the cases $k=1,2$ for $L = (\mathop{div}, \mathop{curl})$ - except there is an additional $L^p$ norm of $s$ on the right-hand side. To get rid of this we need an additional hypothesis: we assume that $L$ has no kernel. In your case this is true because the first Betti number of $S^2$ vanishes. (The estimate would be FALSE on the torus for example. You can see this by taking, in your notation, $\alpha$ to be a non-zero parallel 1-form. Then $\mathop{div}(\alpha)=0=\mathop{curl}(\alpha)$.)
With this additional assumption, let us prove that there exists $C$ for which
$$
|s|_{L^p(M)} \leq C|L(s)|_{L^p_k(M)}
$$
If not then, by contradiction, there is a sequence $s_j$ of sections of $E$ with $|s_j|_{L^p} =1$ and with $L(s_j) \to 0$ in $L^p_k(M)$. By the estimate we have just proved, $s_k$ is bounded in $L^{p}_{k+r}$. By Rellich's Lemma, a subsequence converges in $L^p_{k+r-1}$. The limit $s$ must be a solution of $L(s)=0$ and so must vanish, but this contradicts the fact that $|s|_{L^p}=1$. This means that, when $\mathrm{Ker}(L)=0$, we can absorb the $|s|_{L^p}$ term on the right-hand side in our previous estimate into the $|L(s)|_{L^p_k}$ term.
This completes the proof that when $L$ is elliptic and has no kernel, there is a constant $C$ (depending on $k,p$) such that
$$
|s|_{L^p_{k+r}(M)} \leq C |L(s)|_{L^p_{k}(M)}
$$
This is what you asked for, when you set $k=1,2$ and $L = (\mathop{div},\mathop{curl})$ acting on 1-forms over the 2-sphere.
