Abstract: The fission nuclear reactor is a typical complicated multi-physics coupling system because of the nonlinear coupled terms among different physical fields. The Newton-Krylov method is an effective method for solving the multi-physics coupling problem, featuring strong stability and a high-order convergence rate. Recently, the Newton-Krylov method with an explicit Jacobian matrix has become popular. Compared with the JFNK (Jacobian-free Newton-Krylov) method which doesn’t form the Jacobian matrix explicitly, it has a better preconditioner matrix (the Jacobian matrix itself) and can achieve a more stable and fast convergence. How to calculate the Jacobian matrix efficiently and automatically is a major challenge. The finite difference method is an effective way to compute the Jacobian matrix automatically and can avoid the derivation of matrix element expressions. Besides, the serial graph coloring algorithm has been utilized to achieve an efficient computation of Jacobian. By exploiting the sparsity of the Jacobian matrix, all structurally orthogonal columns can be computed simultaneously through one function evaluation. Thus, the Jacobian computational cost can be reduced by one order of magnitude. However, when solving the larger scale problem under a distributed-memory parallel environment, the parallel graph coloring algorithm is required because the Jacobian is distributed among all processors. In this study, a Distance-2 graph coloring algorithm, arising from the field of graph theory, is applied to color the Jacobian matrix in parallel. This algorithm is performed iteratively, with each iteration consisting of four stages. Local coloring: Each processor tentatively colors its diagonal submatrix using greedy coloring algorithms; Color transfer: The colors of the diagonal submatrix are transferred to other processors to update the off-diagonal submatrix. Conflict detection: Each processor concurrently checks whether the color of the off-diagonal submatrix conflicts with the diagonal submatrix. The conflicts are gathered back to the corresponding diagonal submatrix. Color reset: Reset columns with color conflicts, waiting for recoloring in the next iteration. The entire iteration process continues until each column is correctly colored. This parallel coloring algorithm was verified by solving the steady and transient neutronic/thermal-hydraulics coupling problem arising from a simplified pressurized water reactor model. The Jacobian coloring and Jacobian computing based on this method were tested under a different number of processors. The results indicate that the Jacobian matrix computed by this parallel coloring algorithm is completely correct. The parallel coloring number only slightly increases with the number of processors. The parallel Newton-Krylov method with explicit Jacobian is more efficient and stable than the parallel JFNK method. In a word, this parallel coloring algorithm is successfully applied to the parallel Newton-Krylov method for solving the large-scale neutronic/thermal-hydraulics coupling problem.