This paper proposes a simulation method to model the program Vth distribution of 3-D vertical channel TLC/QLC charge-trapping NAND flash memory. The program Vth distribution can be calculated by considering ISPP noise, WL-WL interference, and the RTN effect of tunneling oxide and poly Si, which are the major physical factors affecting the width of program Vth distribution. Then, the program Vth distribution shapes with different ISPP incremental voltage steps are compared, and the results are found to be consistent with the experimental results. Code and layer-dependent coupling coefficients of WL-WL interference in 3-D vertical channel NAND flash memory are considered. The effect of RTN on the program Vth distribution is comprehensively studied. The program Vth distribution of a WL is calibrated with the measurement, and a good agreement is obtained, validating the array program Vth distribution simulation method. The simulation method can help in improving the reliability of 3-D TLC NAND flash memory and provides guidance for the design and optimization of 3-D QLC NAND flash memory technology.