Thermal energy storage technologies are of current interest in order to improve the integration of renewable energy sources as well as energy efficiency. Numerical simulations of thermochemical heat storage are especially challenging and time consuming due to the complexity of the mathematical description of the strongly coupled and highly nonlinear processes characteristics for such systems. These difficulties are exacerbated once practically relevant complex or large geometries are considered as they can occur around heat exchangers or due to internal heterogeneities of the reactive bed. To allow a computationally efficient simulation of such applications, an existing finite element implementation of a thermochemical heat storage model was parallelised using PETSc routines. Input/output, global assembly and the linear solver all work in a distributed fashion. The approach is implemented into the open source framework OpenGeoSys. The performance of the present parallelisation approach is tested by simulating the discharge of a heat store based on calcium oxide and water as a benchmark problem. The algorithm is tested on 2D as well as 3D meshes. The computational time required for the simulation could be reduced significantly. For example, a 3D model running almost 7 days on a single core could be solved in less than 1 h on 120 cores using the developed framework. The results strongly depend on linear solver and preconditioner settings.