Can I simulate a random field with nested variogram by a summation of independently simulated random field for each componnet of the nested variogram? I want to simulate a random field $Z(u)$ that has a nested variogram, say $\gamma(h)=\gamma_1(h) + \gamma_2(h) + \gamma_3(h)$, assuming the variogram is isotropic. Whether can I simuate independently three random fields: $Z_1(u)$ with correlation structure $\gamma_1(h)$, $Z_2(u)$ with correlation structure $\gamma_2(h)$, $Z_3(u)$ with correlation structure $\gamma_3(h)$, and then sum them up $Z(u) = Z_1(u) + Z_2(u) + Z_3(u)$?

Yes.
This question is introduced in the book "Goovaerts, P., 1997. Geostatistics for natural resources evaluation. Oxford University Press", which is at page 95 with a title of "the linear model of regionalization".
The linear model of regionalization builds the random function $Z(u)$ as a linear combination of $L+1$ independent random functions $Y^l(u)$, each with zero mean and basic covariance function $c_l(h)$:
$Z(u)=\sum_{l=0}^{L}a^lY^l(u)+m \; \;(1)$
with

*

*$E\{Z(u)\}=m$

*$E\{Y^l(u)\}=0, \forall l$

*$Cov\{Y^l(u), Y^{k}(u+h)\}=c_l(h) \; if \; l=k, otherwise \; 0$
The decomposition shown in eq. (1) entails the covariance function of the RF $Z(u)$ is a linear combination of the $(L+1)$ basic covariance functions $c_l(h)$:
$\begin{align}
C_z(h)&=Cov\{Z(u),Z(u+h)\} \\ 
&=\sum_{l=0}^{L}\sum_{k=0}^{L}a^la^kCov\{Y^l(u),Y^k(u+h)\} \\
&=\sum_{l=0}^{L}a^la^lc_l(h)
\end{align}$
The covariance model $C_z(h)$ can be written as: $C_z(h)=\sum_{l=0}^{L}(a^l)^2c_l(h)$. The coviariance model should also be permissible, for which two requirements are:
1.the basic function $c_l(h)$ should be permissible
2.the sill $(a^l)^2$ of each basic covariance model $c_l(h)$ is positive
In addition, the decomposition equation also holds for semivariograms:
$\gamma (h)=\sum_{l=0}^{L}(a^l)^2g_l(h)$
Also, the a priori variance of random field $Z(u)$ can be derived from the coefficients:
$Var\{Z(u)\}=C(0)=\sum_{l=0}^{L}(a^l)^2$
