The objective of this work is the development of a novel finite element formulation describing the contact behavior of slender beams in complex 3D contact configurations involving arbitrary beam-to-beam orientations. It is shown by means of a mathematically concise investigation of well-known beam contact models based on point-wise contact forces that these formulations fail to describe a considerable range of contact configurations, which are, however, likely to occur in complex unstructured systems of thin fibers. In contrary, the formulation proposed here models mechanical contact interaction of slender continua by means of distributed line forces, a procedure that is shown to be applicable for any geometrical contact configuration. The proposed formulation is based on a Gauss-point-to-segment type contact discretization and a penalty regularization of the contact constraint. Additionally, theoretical considerations concerning alternative mortar type contact discretizations and constraint enforcement by means of Lagrange multipliers are made. However, based on detailed theoretical and numerical investigations of these different variants, the penalty-based Gauss-point-to-segment formulation is suggested as the most promising and suitable approach for beam-to-beam contact. This formulation is supplemented by a consistently linearized integration interval segmentation that avoids numerical integration across strong discontinuities. In combination with a smoothed contact force law and the employed C 1 -continuous beam element formulation, this procedure drastically reduces the numerical integration error, an essential prerequisite for optimal spatial convergence rates. The resulting line-to-line contact algorithm is supplemented by contact contributions of the beam endpoints, which represent boundary minima of the minimal distance problem underlying the contact formulation. Finally, a series of numerical test cases is analyzed in order to investigate the accuracy and consistency of the proposed formulation regarding integration error, spatial convergence behavior and resulting contact force distributions. For one of these test cases, an analytical solution based on the Kirchhoff theory of thin rods is derived, which can serve as valuable benchmark for the proposed model but also for future beam-to-beam contact formulations. In addition to these examples, two real-world applications are presented in order to verify the robustness of the proposed formulation when applied to practically relevant problems.