In　practical　engineering,　cross-correlated　random　fields　are　often　used　to　model　structural　materials　or　random　loads　containing　multiple　correlations.　Effective　and　accurate　simulation　of　these　cross-correlation　fields　is　an　important　prerequisite　for　subsequent　reliability　analysis　and　uncertainty　quantification　of　complex　systems.　Therefore,　this　paper　proposes　a　new　simulation　method　for　non-stationary　non-Gaussian　cross-correlated　random　fields　to　satisfy　realistic　engineering　requirements.　In　this　new　method,　L-moments-based　Hermite　polynomials　model　(LHPM)　is　extended　to　non-stationary　non-Gaussian　cross-correlated　fields.　Then,　a　transformation　model　from　non-stationary　non-Gaussian　correlation　matrix　function　(CMF)　to　the　underlying　Gaussian　CMF　is　explicitly　given.　Multivariate　K-L　expansion　is　further　used　to　approximate　the　underlying　Gaussian　cross-correlated　fields,　where　the　eigenvalue　and　eigenfunctions　are　estimated　via　the　Nystrom　method　with　uniform　weights.　The　approximation　permits　the　application　of　few　random　variables　to　characterize　the　entire　Gaussian　cross-correlated　fields.　Ultimately,　the　simulated　underlying　Gaussian　cross-correlated　fields　are　mapped　back　to　the　target　nonstationary　non-Gaussian　cross-correlated　fields　based　on　LHPM.　Three　typical　examples,　including　exponential　kernel,　Wiener　process　kernel　and　spatially　varying　non-Gaussian　and　nonstationary　seismic　ground　motions,　are　used　to　validate　the　effectiveness　of　the　proposed　method.　The　source　code　is　readily　available　at:　https://github.com/zhaozhao23/Simulation-ofnon-stationary-and-non-Gaussian-cross-correlated-fields.