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.
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.
