Definitions
Denoting the Haar measure of some function $f\left(x\right)$ over $d$-dimensional unitaries as $\int_{\mathrm U\left(d\right)}f\left(x\right)d\mu\left(x\right)$, twirling some arbitrary channel $\varepsilon$ can be defined as the operation $\varepsilon \mapsto\int_{\mathrm U\left(d\right)}U^\dagger\varepsilon U dU$, which, when $\varepsilon$ is acting on some density matrix $\rho$ gives $$\rho\mapsto\int_{\mathrm{U}(d)} U^{\dagger} \varepsilon\left(U \rho U^{\dagger}\right) U d U.$$
While the above applies to intergrating over the $d$-dimensional unitary group, the same idea can be applied to (weighted) gatesets $\mathrm D_w$, consisting of elements of unitaries $U_i$ along with associated probabilities $p_i$ such that the weighted gate-set channel is $\mathcal U_{\mathrm D_w}:\rho\mapsto\sum_{U_i\in\mathrm D}p_iU_i\rho U_i^\dagger$. A special case of this are the unweighted gatesets $\mathrm D$ where $p_i = \frac{1}{\left|\mathrm D\right|}\,\forall\, i$. This allows us to define the Twirling operation over a gateset as the operation $\varepsilon\mapsto\sum_{U_i\in\mathrm D}p_iU^\dagger_i\varepsilon U_i$.
The Pauli twirling operation is then defined as $\varepsilon\mapsto\sum_{P_i\in\mathrm P_n}p_iP_i\varepsilon P_i$ where $\mathrm P_n$ is the Pauli group over $n$ qubits, while the Clifford twirl can either be defined as $\varepsilon\mapsto\sum_{U_i\in\mathrm C_n}p_iU^\dagger_i\varepsilon U_i$ or $\varepsilon\mapsto\sum_{U_i\in\mathrm C_1^{\otimes n}}p_iU^\dagger_i\varepsilon U_i$.1
Unitary Designs
There is a difference between what these two twirls do and it lies in the theory of unitary t-designs2: a set $\mathrm D_w$ is a $t$-design if and only if $$\sum_{U\in \mathrm D}\left(U^\dagger\right)^{\otimes t}\otimes \left(U^T\right)^{\otimes t} = \int_{\mathrm U\left(d\right)}\left(V^\dagger\right)^{\otimes t}\otimes \left(V^T\right)^{\otimes t}\,dV$$ and if $\mathrm D_w$ is a $t$-design, it's also a $\left(t-1\right)$-design for any positive integer $t$. This is equivalent to the statement that $$\sum_{U \in D}(U)^{\otimes t} \rho\left(U^{\dagger}\right)^{\otimes t}=\int_{U(d)}(V)^{\otimes t} \rho\left(V^{\dagger}\right)^{\otimes t} d V,$$ which gives the $t^{th}$ statistical moment of $\rho$.
Using Schur's Lemma
For a pure state (density matrix) in 2 dimensions $\rho$, $\int_{\mathrm U\left(d\right)} V\rho V^\dagger\,dV = \frac{1}{2}I$ and so, this fully depolarises (i.e. maximally mixes) the state and as the single qubit Pauli twirl is also the depolarisation channel, the Pauli group is a $1$-design. However, by Schur's lemma, it's not a $2$-design. In other words, as explained in section IIA of Gross, Audenaert and Eisert, twirling essentially writes the state in terms of a sum of irreducible subspaces, which, for a 2D 2-design, are the symmetric $\frac{1}{2}\left(I+SWAP\right)$ and antisymmetric $\frac{1}{2}\left(I-SWAP\right)$ subspaces and this is not what the Pauli twirl does. However, it is what the Clifford twirl does. Furthermore, as shown both here and here, the Clifford group forms a 3-design.
Randomised Benchmarking
This is useful and even relatively important from a quantum computing perspective when it comes to benchmarking - it turns out that local random circuits are approximate designs and can be created in polynomial time
Vectorisation
This all relates to the papers linked in the question by the process of vectorisation and the property that $|A B C\rangle \rangle=A \otimes C^T|B\rangle \rangle$, which when $\varepsilon$ is written as a sum of Kraus operators gives $$\int_{\mathrm{U}(d)} U^{\dagger} \varepsilon\left(U \rho U^{\dagger}\right) U d U \mapsto\sum_{K} \int_{U(d)}\left(U^{\dagger} \otimes U^T\right) \cdot\left(K \otimes K^{*}\right) \cdot\left(U \otimes U^{*}\right)|\rho\rangle \rangle d U.$$
If we then go a step further and start with a random state, we can effectively integrate over $\rho$ and repeat the vectorisation process to give $$\int_{\mathrm{U}(d)} U^{\dagger} \varepsilon\left(U \rho U^{\dagger}\right) U d U \mapsto\int_{U(d)} U^{\dagger} \otimes U^T \otimes U^T \otimes U^{\dagger} d U \sum_{K}\left|K \otimes K^* \right\rangle \rangle$$ and we now have something that looks more like a 2-design.
Overall, we've gone from twirling over some channel to an operation that looks something similar to 2-design acting on the vectorisation of the Kraus operator form of that channel and finally...
The difference between Pauli and Clifford Twirls
Rewriting $K = \sum_j|u_j^{\left(K\right)}\rangle\langle v_j^{\left(K\right)}|$ gives $K\otimes K^* = \sum_{j, l}|u_j^{\left(K\right)}\rangle|u_l^{*\left(K\right)}\rangle\langle v_j^{\left(K\right)}|\langle v_l^{*\left(K\right)}|$ and so $|K\otimes K^*\rangle\rangle = \sum_{j, l}| v_j^{\left(K\right)}\rangle| v_l^{*\left(K\right)}\rangle|u_j^{\left(K\right)}\rangle|u_l^{*\left(K\right)}\rangle$
From what's above, the Clifford twirl then permutes this to get $\sum_{j, l}|u_l^{*\left(K\right)}\rangle| v_j^{\left(K\right)}\rangle| v_l^{*\left(K\right)}\rangle|u_j^{\left(K\right)}\rangle$, sums the projections onto the symmetric and anti-symmetric subspaces and unpermutes again. Overall, information about the second statistical moment of the noise channel can be learned from repeated measurements.
By comparison, the Pauli twirl generally doesn't do this. It does, in the case of a single qubit, at least give the diagonal elements of the Process matrix, as per the paper they refer to, which happens to be the parameters they're looking for.
1 I'm not sure off the top of my head if there's a diffrerence between $\mathrm C_n$ and $\mathrm C_1^{\otimes n}$, although it seems clear that $\mathrm C_1^{\otimes n}\subset\mathrm C_n$
2 Also known as unitary k-designs
