Large-scale three-dimensional (3D) seismic modeling is considered the foundation of imaging and inversion. Parallelization strategies are crucial to solving such a computationally intensive problem. In this paper, we focus on two factors, decomposition direction and decomposition dimension, that significantly affect the computational performance. The decomposition direction determines the cache hit ratio during register addressing, and the decomposition dimension influences the communication size. We thoroughly analyze these two factors by performing time-space domain staggered-grid finite-difference (SGFD) modeling with a set of decomposition strategies. Four metrics, including computation time, speedup ratio, strong scaling property, and memory usage, are introduced to evaluate the computational performance of each trial. After theoretical analysis and metrics testing, we conclude that the optimized domain decomposition strategy is: decomposing models at two dimensions, the decomposition directions are perpendicular to the fastest and the second fastest dimensions, here we refer the dimension in which data are continuously saved as the fastest dimension. Three examples further verify the feasibility and efficiency of the optimized parallel scheme. Considering that domain decomposition-based 3D seismic parallel simulation packages are seldom available in the public domain, we provide a program template for the optimized domain decomposition strategies as an open-source package.