Some thoughts:
A theoretical perspective
From a theoretical perspective, the depolarizing channel is the 'standard' (if there is such a thing) or by some means the most applicable.
Because the Paulis (together with the identity operator) form a basis for $SU(2)$, if a code can correct the $X, Y$ and $Z$ flips on a certain qubit (and it it is able to correctly identify no error (i.e. '$I$-flip') having happened), it can correct all errors on that qubit. A theoretical analysis of a code cares a little less about the relative probabilities of these three flips happening, so we just as well can set them all to $\frac{p}{3}$, thereby obtaining the depolarizing channel.
The depolarizing channel $\Lambda_{\mathrm{depo}}$ can be written as:
$$
\Lambda_{\mathrm{deph}}\left(\begin{bmatrix}a & b \\ b^{*} & d\end{bmatrix}\right) = \begin{bmatrix} (1-\frac{2p}{3})a + \frac{2p}{3}d & (1-\frac{4p}{3}) b \\ (1-\frac{4p}{3})b^{*} & (1-\frac{2p}{3})d + \frac{2p}{3}a\end{bmatrix}.
$$
Since $d = 1-a$, we can rewrite this to:
$$
\Lambda_{\mathrm{deph}}\left(\rho_{\mathrm{in}}\right) = (1-\frac{4p}{3})\rho_{\mathrm{in}} + \frac{4p}{3} \frac{I}{2},
$$
which is a convex combination of the input $\rho_{\mathrm{in}}$ and the maximally mixed state $\frac{I}{2}$.
Furthermore, if you would equate $p$ to the elapsed time, $p$ would asymptotically go to $\frac{3}{4}$, thereby obtaining the maximally mixed state; therefore the depolarizing channel is in some way the 'worst' noise channel: it destroys both all quantum- (i.e. coherent superpositions) and classical information (there's literally only noise left) in the qubit.
A physical perspective
If you want your simulation to be more true to the physical world, the depolarizing channel is not a very good model, as much as theorists might like it to be.
A good first model for noise in qubits is the combination of two channels, the dephasing channel $\Lambda_{\mathrm{deph}}$ and the amplitude damping channel $\Lambda_{\mathrm{amp}}$.
Loosely speaking, the dephasing channel map destroys the coherent phase between the $|0\rangle$ and $|1\rangle$ state:
$$
\Lambda_{\mathrm{deph}}\left(\begin{bmatrix}a & b \\ b^{*} & 1-a\end{bmatrix}\right) = \begin{bmatrix}a & e^{-\frac{t}{T_{2}}} b \\ e^{-\frac{t}{T_{2}}}b^{*} & 1-a\end{bmatrix},
$$
where $T_{2}$ is known as the characteristic qubit dephasing time.
The Kraus operators of $\Lambda_{\mathrm{deph}}$ are $A_{1} = \sqrt{1-p}I$ and $A_{2} = \sqrt{p}Z$, so it is still a Pauli channel, which can help in the analysis or simulation.
$\Lambda_{\mathrm{amp}}$ is a little bit trickier: it simulates the relaxation of the excited (by convention $|1\rangle$) state, and maps it to the $|0\rangle$ state:
$$
\Lambda_{\mathrm{amp}}\left(\begin{bmatrix}a & b \\ b^{*} & 1-a\end{bmatrix}\right) = \begin{bmatrix}a & e^{-\frac{t}{2T_{1}}} b \\ e^{-\frac{t}{2T_{1}}} b^{*} & e^{-\frac{t}{T_{1}}}(1-a)\end{bmatrix},
$$
where $T_{1}$ is known as the qubit relaxation time.
The Kraus operators are also a little trickier: $B_{1} = \begin{bmatrix}1 & 0 \\ 0 & \sqrt{1-p}\end{bmatrix}$ and $B_{2} = \begin{bmatrix}0 & \sqrt{p} \\ 0 & 0\end{bmatrix}$. This means that amplitude damping channel is not a Pauli channel (allthough its Kraus operators can of course be written as linear combinations of the Paulis).
A simulation of a system undergoing both an amplitude damping and dephasing is a good start. There is one big caveat, however: this does not take leakage errors into account whatsoever. Depending on the physical system you are trying to simulate, this might range from either not a very large issue (e.g. for quantum dots) to a gross oversimplification (e.g. for transmon qubits).