Simulation of continuous variables is an important aspect in the construction of numerical models of Earth Sciences variables. In many cases, there is a need for high resolution models either to model large volumes or surfaces, or to characterize on a very dense grid some spatial phenomena. Traditional simulation algorithms such as turning bands and sequential Gaussian simulation must be adapted to handle such large models and to perform the computations in a reasonable time. Parallelization requires a review of the algorithm design. We present the changes and performance improvements of parallelizing the turning bands algorithm using graphical processing units (GPU). Our implementation focuses on parallelization of the unconditional simulation step. One example is presented comparing the performance of a serial simulation implementation versus our parallelized implementation of turning bands. Speed ups of up to 60 are reached and an analysis of possible additional improvements is presented.