The western margin of central South America is one of the best studied regions of subduction zones. The oceanic Nazca plate descends beneath the continental South American plate at an oceanward concave trench. This form of the trench is unique on Earth and lead to the formation of the Andes. Here, I focus on a three-dimensional thermomechanical finite element model along the north Chilean margin. At first a reference model with simplified geometry reaching down to 400 km depth is created that makes use of an elasto-viscoplastic rheology including shear heating along the plate interface. The subduction process is simulated for geologically short timescales to reflect the recent state of subduction and the model results are in first-order agreement with observations. The model is analyzed within the framework of a parametric study revealing the most important parameters to be convergence and slab pull velocity, type and number of elements of the finite element mesh, type of rheology, cohesion, friction coefficient and angle of internal friction. Because of its lower thickness the oceanic crust is generally more affected by changes in these parameters than the continental crust. The final model includes a realistic geometry and density distribution that is based on a well constrained density model, not considering topography and bathymetry. It simulates 100,000 years of subduction and accounts for elastic deformations. The geometry has a major impact on the model results causing significant differences between the northern and southern parts of the model. This, as well as the final distributions of shear stress, strain rate and heat flux is in good agreement with observations from global models except for the magmatic arc because inelastic deformations are neglected. The usage of a realistic geometry and density distribution for finite element simulations of subduction zones can thus help to significantly improve future geodynamic models.