Dynamic traffic assignment models rely on a network performance module known as dynamic network loading(DNL), which expresses the dynamics of flow propagation, flow conservation, and travel delay at a network level. The DNL defines the so-called network delay operator,which maps a set of path departure rates to a set of path travel times. It is widely known that the delay operator is not available in closed form, and has undesirable properties that severely complicate DTA analysis and computation, such as discontinuity, non-differentiability, non-monotonicity, and computational inefficiency. This paper proposes a fresh take on this important and difficult problem, by providing a class of surrogate DNL models based on a statistical learning method known as Kriging.We present a metamodeling framework that systematically approximates DNL models and is flexible in the sense of allowing the modeler to make trade-offs among model granularity, complexity, and accuracy. It is shown that such surrogate DNL models yield highly accurate approximations (with errors below 8%) and superior computational efficiency (9 to 455 times faster than conventional DNL procedures). Moreover, these approximate DNL models admit closed-form and analytical delay operators, which are Lipschitz continuous and infinitely differentiable, while possessing closed-form Jacobians. The implications of these desirable properties for DTA research and model applications are discussed in depth.