How to write down a 3-level hierarchical model with a level-1 equation (reduced form)? Suppose longitudinal data on $y_{rct}$: the average temperature of region $r$ in country $c$ at time $t$. You $y_{rct}$ fit a 3-level hierarchical linear model with country-year and country deviations from an overall mean. Assuming $y_{rct}$ is centered at 0, a model could be (variance hyper-priors omitted):
$$y_{rct} \sim N(\mu_{ct}, \sigma_{rct}) \\
\mu_{ct} \sim (\mu_c,\sigma_{ct})\\
\mu_c \sim N(0,1)\\$$
This notation makes clear that the mean at level 1 is a function of the mean at level 2, which is a function of the mean at level 3. However, in popular software like the lme4 R package or Stata, the code to produce a 3-level model usually involves a single-level model with independent country and country-year offsets to the grand mean. For example, here or p.7 here. This seems to be equivalent to something like:
$$y_{rct} \sim N(\mu_{ct} + \mu_c, \sigma_{rct})$$
QUESTION
While Model 2 adds the same deviations to the mean of $y_{rct}$, are these models actually equivalent? My hunch is no: Model 1 shrinks the $\mu_{ct}$ toward its respective $\mu_c$ but Model 2 shrinks $\mu_{ct}$ toward the grand mean of 0. However, perhaps I am missing something about how these work that makes Model 2 the only way to specify nested intercepts in most software.
 A: In what follows I'm assuming that you follow the convention that $\mathcal{N}\!\left(\mu,\sigma\right)$ denotes a normal distribution with mean $\mu$ and standard deviation $\sigma$.

Using the law of total expectation we find
$$
\mathbb{E}\left(y_{rct}\right) =
\mathbb{E}\left(\mathbb{E}\left(y_{rct}\mid\mu_{ct}\right)\right) = \mathbb{E}\left(\mu_{ct}\right) = \mathbb{E}\left(\mathbb{E}\left(\mu_{ct}\mid\mu_c\right)\right) =
\mathbb{E}\left(\mu_{c}\right) = 0.
$$
Assuming fixed variance parameters $\sigma_{rct}^2,\sigma_{ct}^2$ (and $\sigma_{c}^2$) , the law of total variance delivers
$$
\begin{align}
\mathbb{V}\left(y_{rct}\right) &=
\mathbb{E}\left(\mathbb{V}\left(y_{rct}\mid\mu_{ct}\right)\right) + 
\mathbb{V}\left(\mathbb{E}\left(y_{rct}\mid\mu_{ct}\right)\right) \\&=
\mathbb{E}\left(\sigma_{rct}^2\right) + \mathbb{V}\left(\mu_{ct}\right) \\&= \sigma_{rct}^2 +
\mathbb{E}\left(\mathbb{V}\left(\mu_{ct}\mid\mu_{c}\right)\right) +
\mathbb{V}\left(\mathbb{E}\left(\mu_{ct}\mid\mu_{c}\right)\right) \\&=
\sigma_{rct}^2 + \mathbb{E}\left(\sigma_{ct}^2\right) + \mathbb{V}\left(\mu_{c}\right) \\&= 
\sigma_{rct}^2 + \sigma_{ct}^2 + \sigma_{c}^2,
\end{align}
$$
where $\sigma_{c}^2 \equiv 1$ in your model.
Similarly, the law of total covariance yields
$$
\begin{align}
&\mathrm{Cov}\left(y_{rct}, y_{\tilde{r}ct}\right) = \sigma_{c}^2 + \sigma_{ct}^2 &\text{for} \; \tilde{r} \neq r, \\
&\mathrm{Cov}\left(y_{rct}, y_{\tilde{r}c\tilde{t}}\right) = \sigma_{c}^2
&\text{for} \; \tilde{t} \neq t, \\
&\mathrm{Cov}\left(y_{rct}, y_{\tilde{r}\tilde{c}\tilde{t}}\right) = 0 &\text{for} \; \tilde{c} \neq c.
\end{align}
$$
The marginal model for the vector of average temperatures $\mathbb{y}$, i.e. the joint distribution of the $y_{rct}$s, is thus given by a multivariate normal distribution with mean $\mathbb{0}$ and covariance-matrix entries according to the (co)variances derived above. In a sense this gives you a "level-1 expression" for $\mathbb{y}$. In particular, we have
$$
y_{rct} \sim \mathcal{N}\!\left(0, \sqrt{\sigma_{c}^2 + \sigma_{ct}^2 + \sigma_{rct}^2}\right)
$$
for the marginal distribution of the individual $y_{rct}$s.
Such a model for $\mathbb{y}$ can be estimated in many mixed model packages.  To this end, it's useful to write the model, in analogy to Eoin's answer, as
$$
y_{rct} = \mu_c + \gamma_{ct} + \varepsilon_{rct},\\
\mu_c \overset{\mathrm{i.i.d.}}{\sim} \mathcal{N}\left(0, \sigma_c\right),\\
\gamma_{ct} \overset{\mathrm{i.i.d.}}{\sim} \mathcal{N}\left(0, \sigma_{ct}\right),\\
\varepsilon_{rct} \overset{\mathrm{i.i.d.}}{\sim} \mathcal{N}\left(0, \sigma_{rct}\right),\\
\mu_c \mathrel{\perp\!\!\!\perp} \gamma_{ct}, \mu_c \mathrel{\perp\!\!\!\perp} \varepsilon_{rct}, \gamma_{ct} \mathrel{\perp\!\!\!\perp} \varepsilon_{rct},
$$
where $\varepsilon_{rct}$ is the observation-level error term. As an example, the corresponding lme4::lmer syntax would then be
lmer(temp_centered ~ (1 | country) + (1 | country:time)).
Fixing $\sigma_{c}$ to $1$ would potentially (depending on the chosen software package) require more work, and considering the variance parameters as random variables in a fully Bayesian approach would likely indicate the usage of packages tailored to Bayesian (multilevel) regression analysis.
A: Your question, "While Model 2 adds the same deviations to the mean of $y_{rct}$, are these models actually equivalent?" isn't the same as the question title, "How to write down a 3-level hierarchical model with a level-1 equation (reduced form)?". @statmerkur has answered the question title, but I'm just going to answer the former question.

$$y_{rct} \sim N(\mu_{ct}, \sigma_{rct}) \\
\mu_{ct} \sim (\mu_c,\sigma_{ct})\\
\mu_c \sim N(0,1)\\$$
is not quite equivalent to
$$
y_{rct} \sim N(\mu_{ct} + \mu_c, \sigma_{rct})
$$
as the $\mu_{ct}$ terms do not mean the same thing in both cases.
However, your first equation is equivalent to
$$y_{rct} \sim N(\mu_c + \gamma_{ct}, \sigma_{rct}) \\
\gamma_{ct} \sim N(0, \sigma_{ct})\\
\mu_c \sim N(0, 1)\\$$
where $\gamma_{ct} = \mu_{ct} - \mu_c$ (e.g. how much country $c$ at time $t$ deviates from the overall mean for country $c$). It follows that shrinking $\gamma_{ct}$ towards $0$ is equivalent to shrinking $\mu_{ct}$ towards $\mu_c$.
