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)$?
 A: Yes.
To see why, begin with the definition of the variogram for a second-order stationary random field $Z$ as
$$\gamma(h;Z) = \frac{1}{2}E\left[(Z(u+h)-Z(u))^2\right]$$
for any $u$ where both $u$ and $u+h$ are in the support of $Z.$  (Stationarity implies these values do not depend on $u.$)
Suppose $W$ is another second-order stationary random field and let both $u$ and $u+h$ be in the supports of both $W$ and $Z.$  To abbreviate the notation, write
$$z(h) = Z(u+h)-Z(u);\quad w(h) = W(u+h)-W(u)$$
and note that $z$ and $w$ are both second-order covariance stationary with constant zero expectations.
In this case, the definition gives
$$\begin{aligned}
\gamma(h;Z+W)  &= \frac{1}{2}E\left[((Z+W)(u+h)-(Z+W)(u))^2\right] \\
&= \frac{1}{2}E\left[(z(h) + w(h))^2\right] \\
&= \frac{1}{2}E\left[z(h)^2 + w(h)^2 + 2z(h)w(h)\right] \\
&= \frac{1}{2}E\left[z(h)^2\right] +  \frac{1}{2}E\left[w(h)^2\right] + 2 \frac{1}{2}E\left[z(h)w(h)\right] \\
&= \gamma(h;Z) + \gamma(h;W) + E\left[z(h)w(h)\right]
\end{aligned}$$
All steps are just simple algebraic manipulations.  That's as far as we can get in general.  But when $z(h)$ and $w(h)$ are uncorrelated -- which is equivalent to $E[z(h)w(h)] = E[z(h)]E[w(h)]$ -- then this simplifies.  Moreover, this lack of correlation is guaranteed by the stronger independence hypothesis of the question, because independence of $Z$ and $W$ implies independence of $z$ and $w$ and finiteness of the variances of $Z(u),$ $Z(u+h),$ $W(u),$ and $W(u+h)$ then implies lack of correlation.
Indeed, the same analysis goes through when we replace $Z$ by $\lambda Z$ and $W$ by $\mu W$ for fixed constants $\lambda$ and $\mu,$ resulting in the conclusions

$\lambda Z+\mu W$ is second-order stationary and $\gamma(h; \lambda Z + \mu W) = \lambda^2 \gamma(h; Z) + \mu^2 \gamma(h; W).$

Repeated application of this result to, say, $Z_1$ and $Z_2 + Z_3$ yields the formula $\gamma=\gamma_1+\gamma_2+\gamma_3$ stated in the question.
Remarks
Because we never had to assume isotropy in the derivation, isotropy is an unnecessary assumption.  We can further deduce, though, that when all the $Z_i$ are isotropic, then so is any linear combination of them.  (This follows directly from any definition of isotropy.)
Beware of conditional simulations!  If each of your simulations is conditional on data, then it won't do to add them up, because then you obtain an expectation of three times each data value in the support of the data.  Instead, simulate exactly one of the fields conditional on the data (it doesn't matter which one) and simulate all the others conditional on all the observations being zero.  The sum will have the intended correlation structure (as shown above) and now it will have the correct conditional expectations at the data points.
A: 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$
