Why are multivariate distributions often used for hierarchical model when errors are conditionally independent? Suppose you have data on $y_{it}$: outcome $y$ of unit $i$ at time $t$. You fit a varying intercept (hierarchical) model:
$$y_{it} \sim Normal(\mu_{i}, \sigma_y)\\
\mu_i \sim Normal(0, \psi_{\mu})\\
\sigma_y \sim Cauchy(0,1)\\
\psi_{\mu} \sim Cauchy(0,1)\\
$$
I often see specifications where the random effects $\mu_i$ follow some multivariate distribution. The stated rationale for this is typically that the random effects are correlated in some way. For example, https://cran.r-project.org/web/packages/brms/vignettes/brms_multivariate.html talks about the need to correlate a bunch of random intercepts and proposes a multivariate distribution for this.
This puzzles me since typically the justification for the varying intercepts is that the outcomes $y_{it}$ might be correlated between time periods $t$ but are uncorrelated within units $i$. Accordingly, including the random effects in the model now means that $y$ are independent of each other. If we thought the $u_i$ were correlated with each other, then we could add another set of intercepts to model those higher level correlations. However, it's typically not the case!
So, why don't people just use univariate normal likelihoods/priors instead of multivariate?

Your intuition works for univariate multilevel/hierarchical models, but the link you provide describes multivariate multilevel models.  So rather than

data on $_{}$: outcome $$ of unit $$ at time $$

the outcomes are a vector, say $[y_{it}, z_{it}]^{'}$ (variables tarsus and back in the posted example).  Thus, there are 2 random intercepts (for $y$ and $z$) rather than 1 (for $y$ only), say:

*

*$u_i$: the degree to which subject $i$'s average score on $y$ is higher/lower than the grand mean $\bar{y}$

*$v_i$: the degree to which subject $i$'s average score on $z$ is higher/lower than the grand mean $\bar{z}$
Conditioning on these random intercepts can indeed make the errors $[e^{(y)}_{it}, e^{(z)}_{it}]^{'}$ conditionally independent across occasions (e.g., $t$ with $t-1$) within each subject $i$.  However, that has nothing to do with correlations across variables at the within- or between-subject levels.

*

*Correlation of $u_i$ with $v_i$ indicates whether subjects above average on $y$ tend to be above/below average on $z$

*Correlation of $e^{(y)}_{it}$ with $e^{(z)}_{it}$ indicates whether subject $i$ being higher than usual on $y$ at Time $t$ also tends to be higher/lower than usual on $z$ at Time $t$
So, whereas a univariate MLM decomposes $y_{it}$ into subject-level ($\sigma^2_i$) and occasion-level ($\sigma^2_t$) variance components, a multivariate MLM decomposes the entire covariance matrix $\Sigma$ among multiple outcomes ($y$ and $z$) into a subject-level $\Sigma_i$ and an occasion-level $\Sigma_t$.  This is the basis for multilevel SEM, if you are familiar with such models, e.g.:
https://psu-psychology.github.io/psy-597-SEM/15_multilevel/multilevel_sem.html#multilevel-sem-msem-overview
So to clarify: across-variable correlations of random effects (at Level 2) or errors (at Level 1) are not a violation of conditional independence.  Conditional on the random effects, the error-vectors for rows of data are still (conditionally) independent, although there can still be correlation across columns of data (i.e., correlation of residuals of variable $y$ with residuals of variable $z$).
