The present study proposes a methodology to generate stochastic slip field with desired probability structure and spectral characteristics. The characterization utilizes 100 tsunamigenic and 130 non-tsunamigenic slip models available in the SRCMOD database. First, the cumulative distribution function (CDF) of marginal slip distributions is obtained. Then, it compares the estimated CDFs with nine standard distribution functions based on model fitting criteria like Akaike Information Criterion, mean squared error, correlation coefficient, and QQ plots. The analysis showed that the majority of slip distributions follow generalized Pareto or truncated exponential distribution functions. For the parameters of CDF estimated from slip fields, an empirical relation is proposed as a function of magnitude (Mw) to predict the corresponding values. As the next step, to understand the correlation structure of the slip field, semi-variograms are estimated along strike and along downdip for all rupture models considered in the study. A total of nine theoretical variogram models is compared with the semi-variogram of each slip field. From the analysis based on model fitting criteria, most of the slip field is observed to follow a stable variogram model. An empirical equation is developed to predict the sill and range of stable variogram function for a given Mw. Then, an iterative method based on translational field theory is proposed to obtain an ensemble of non-Gaussian slip fields for a given Mw and type of source. The generated sample fields from the proposed method match with the prescribed non-Gaussian distribution and variogram structure. In the sequel, the method estimates the variogram and the corresponding spectral density of the underlying Gaussian field according to translation field theory. The scaling relations and the proposed iterative methodology can be used to obtain an ensemble of possible earthquake slip fields.