Whether the autocorrelation structure of random field $Z(u) + Y(u)$ is equal to the autocorrelation structure of $Z(u)$ plus that of $Y(u)$?

Question

I want to simulate a Gaussian random field (RF) with correlation structure (represented by the geostatistic tool 'semivariogram' $\gamma (h) \: +\: pure \: nugget \: effect$). I want to know whether can I do the simulation by independently simulate a Gaussian RF with correlation structure $\gamma (h)$ and another rand field with pure nugget effect, and then sum them up. ($u$ in the title denotes spatial location $(x, y)$). I do a simulation experiment in Matlab:

clear;
clc;
rng(0);
mu = 0; 
sigma2 = 1; % sill
covar_type = 3; % Spherical variogram model
range = 20; 
knn = 12; % number of neighbors to consider
xsize = 100; % number of columns
ysize = 100; % number of rows

field1 = J_Simu2DSGS(xsize, ysize, [], mu, sigma2, covar_type, range, knn);

% field 2
nugget = 0.3;
purenugget = sqrt(nugget)*randn(xsize, ysize);

% field1 + field2
comRF = field1 + purenugget;
Sampling from simulated fields:
samplesz = 500;
[xcoord, ycoord] = meshgrid(1:xsize, 1:ysize);
xycc = [xcoord(:), ycoord(:), field1(:), comRF(:)];
nbcell = xsize*ysize;
idx = randperm(nbcell, samplesz);
obsxy = xycc(idx, 1:2);
obscc = xycc(idx, 3:end);

The function used for simulation 'J_Simu2DSGS' above is equivalent to the function 'MGSimulSGS.m' here: https://github.com/GAIA-UNIL/MATLAB_Geostat_Utilities/blob/master/MATLAB_UTILS/MGSimulSGS.m

I calculate the semivariograms of the simulated RFs. The result is:

We can see that the correlation structure of the composite RF (variogram at lower right corner) is similar to that of the RF1 (variogram at lower left corner), in spite of a general shift on the y-axis representing the nugget effect.

But I also want to validate it mathematically. Namely:

Assuming second-order stationarity holds, if Gaussian random field $Z(u)$ has correlation structure $C_1 (h)$, Gaussian random field $Y(u)$ has correlation structure $C_2 (h)$, then $Z(u) + Y(u)$ has correlation structure $C_1 (h) + C_2 (h)$.

My math is poor, could someone help demonstrate it?

Note: $\gamma (h)$ is only associated with the separated distance in space since we make a second-order stationarity assumption. More info can be found in wiki: https://en.wikipedia.org/wiki/Variogram