I wrote in Matlab a code which generates 1D unconditional random fields following section 2.1.2 of “Application of conditional random fields and sparse polynomial chaos expansion to geotechnical problems” and it seems to work well. For example, I use it to model spatial variability of cohesion dependent on the depth in the soil layer. Attached a picture with 4 realizations of a depth dependent soil cohesion created with my code. The question is, how to generate random fields for 2D spatial variability? What I would like to have is for example soil cohesion which is dependent on the position (x,y) and not only on the depth y. More precisely, I don’t understand which size the covariance matrix is going to have. In 1D spatial variability, if my domain is discretized by 10 points, then my covariance matrix is going to be 10X10, but what if my domain is discretized by 10 point in x and y direction?