The vehicle routing problem is one of the most important areas of logistics management. This study considers two-echelon vehicle routing problem with simultaneous pickup and delivery (2E-VRPSPD) which is a variant of vehicle routing problem. In the 2E-VRPSPD, the pickup and delivery activities are performed simultaneously by the same vehicles through depot to satellites in the first echelon and from satellites to customers in the second echelon. To solve the problem, firstly, a node-based mathematical model is proposed and three valid inequalities from the literature are adapted to strengthen the model. Because of the NP-hardness of the 2E-VRPSPD, secondly, a hybrid heuristic algorithm based on variable neighborhood descent (VND) and local search (LS), called VND_LS, is developed to solve medium- and large-size instances of the 2E-VRPSPD. We conduct an experimental study to investigate the effectiveness of the valid inequalities on the mathematical model and also to evaluate the effectiveness and efficiency of the VND_LS. Computational results show that valid inequalities have significant effect to strengthen the mathematical formulation. Furthermore, the VND_LS finds good solutions for the problem efficiently. Finally, we apply the VND_LS to compare single- and two-echelon distribution systems for a supermarket chain located in Turkey. The results indicate that the VND_LS can easily be applied for real-world problems.