In three-dimensional (3D) integrated circuits, the effect of process variations on clock skew differs from 2D circuits. The combined effect of inter-die and intra-die process variations on the design of 3D clock distribution networks is considered in this article. A statistical clock skew model incorporating both the systematic and random components of process variations is employed to describe this effect. Two regular 3D clock tree topologies are investigated and compared in terms of clock skew variation. The statistical skew model used to describe clock skew variations is verified through Monte-Carlo simulations. The clock skew is shown to change in different ways with the number of planes forming the 3D IC and the clock network architecture. Simulations based on a 45-nm CMOS technology show that the maximum standard deviation of clock skew can vary from 15 ps to 77 ps. Results indicate that simply increasing the number of planes of a 3D IC does not necessarily lead to lower skew variation and higher operating frequencies. A multigroup 3D clock tree topology is proposed to effectively mitigate the variability of clock skew. Tradeoffs between the investigated 3D clock distribution networks and the number of planes comprising a 3D circuit are discussed and related design guidelines are offered. The skew variation in 3D clock trees is also compared with the skew variation of clock grids.