Raster distortion from parallel kriging in R

Question

I am relatively new to parallel computing in R. I am using a high performance cluster to perform kriging for soil data points. I found this previous post which I used as a template:

How to achieve parallel Kriging in R to speed up the process?

This is the code I have been using, non-parallel, it works fine a produces kriged rasters. For clarity, I am using Random Forests to initially create the spatial predictions from the point data, then run kirging on the residuals.

qs_s5: spatial points data frame of the point data

p2: spatial pixels data frame of the covariate rasters

p <- as(p, Class= "data.frame")
p1 <- p[complete.cases(p),]
p2 <- na.omit(p1)
colnames(p2)[colnames(p2) == "s1"] <- "lat"
colnames(p2)[colnames(p2) == "s2"] <- "long"
v.rf <- variogram(qs_s5$resid.rf~1, data = qs_s5, width = 1000,cutoff=10000)
vm <- vgm(psill = 10, model = "Sph", range =4000, nugget = 20)
vmf.rf <- fit.variogram(v.rf, model = vm)
rf.s_ppm <- krige(resid.rf ~ 1, locations = qs_s5, newdata = p2, model = vmf.rf, debug.level = -1)
gridded(p2)<- ~lat+long
proj4string(p2) <- CRS(utm14)
rf.s_ppm$trend <- p.rf


rf.s_ppm$trend <- p.rf

rf.s_ppm$pred <- rf.s_ppm$trend + rf.s_ppm$var1.pred

rf.s_ppm$lower <- rf.s_ppm$pred - 1.645 * sqrt(rf.s_ppm$var1.var)
rf.s_ppm$upper <- rf.s_ppm$pred + 1.645 * sqrt(rf.s_ppm$var1.var)

writeGDAL(rf.s_ppm["pred"], fname = "./out10/s_ppm_parallel.tif", drivername = "GTiff", type = "Float32")

Here is the code with parallel integrated:

v.rf <- variogram(qs_s5$resid.rf~1, data = qs_s5, width = 1000,cutoff=10000)
vm <- vgm(psill = 10, model = "Sph", range =4000, nugget = 20)
vmf.rf <- fit.variogram(v.rf, model = vm)

gridded(p2)<- ~lat+long
proj4string(p2) <- CRS(utm14)

numCores <- detectCores() -1
numCores 

cl <- makeCluster(numCores)

parts <- split(x = 1:length(p2), f = 1:numCores)

clusterExport(cl = cl, varlist = c("qs_s5", "p2", "vmf.rf", "parts"), envir = .GlobalEnv)
clusterEvalQ(cl = cl, expr = c(library('sp'), library('gstat')))

parallelX <- parLapply(cl = cl, X = 1:numCores, fun = function(x) krige(resid.rf ~ 1, locations = qs_s5, newdata = p2[parts[[x]],], model = vmf.rf, debug.level = -1))

stopCluster(cl)

mergeParallelX <- maptools::spRbind(parallelX[[1]], parallelX[[2]])
for (j in 3:length(parallelX)) {
mergeParallelX <- maptools::spRbind(mergeParallelX, parallelX[[j]])
}

mergeParallelX <- SpatialPixelsDataFrame(points = mergeParallelX, data = mergeParallelX@data)

rf.s_ppm <- mergeParallelX

rf.s_ppm$trend <- p.rf

rf.s_ppm$pred <- rf.s_ppm$trend + rf.s_ppm$var1.pred

rf.s_ppm$lower <- rf.s_ppm$pred - 1.645 * sqrt(rf.s_ppm$var1.var)
rf.s_ppm$upper <- rf.s_ppm$pred + 1.645 * sqrt(rf.s_ppm$var1.var)

writeGDAL(rf.s_ppm["pred"], fname = "./out10/s_ppm_parallel.tif", drivername = "GTiff", type = "Float32")

Sdata<-raster("./out10/s_ppm_parallel.tif")
plot(Sdata)

When the raster is generated, the image is distorted to swirls. I understand that this must be do to the way the analysis is cut, distributed, and merged; but I dont know what in the code is the cause nor the solution.