The recently developed multiphase geothermal flow simulator, Waiwera, is now able to handle unstructured grids. However, to reduce errors related to the grid, certain quality factors should be optimised. In applying the finite volume method, Waiwera assumes that for each pair of neighbouring tetrahedral cells the line connecting their centroids is perfectly orthogonal to their shared interface. Any deviation from orthogonality will introduce an error in the flow simulation. Our work aims to optimize a given grid so as to minimise departures from orthogonality as much as possible. Our algorithm is written in Python. It uses the MeshIO library to read and write various formats of mesh files, and the SciPy library for optimisation algorithms. Our approach is to minimise the sum of all deviations from orthogonality within the mesh. The inputs are the node coordinates and connectivity and any surface constraints (e.g., model boundaries, faults). Then, our algorithm moves the nodes around in space in an attempt to minimise deviations from orthogonality. Several challenges were addressed in developing the optimisation algorithm. First, we require that some nodes should be constrained to move only along boundaries (model edges, topography, faults) so that geometric features honoured by the original mesh are not altered. This reduces a node’s degrees of freedom from 3 – the number of coordinate directions - to (3 – n), where n is the number of boundaries the node lies on. As a result, the algorithm can handle plane and complex boundaries (including faults and topography), intersections between planes and between a plane and a complex surface. Second, we aimed to minimise the running time of the algorithm so that the approach is practical for full-sized meshes (up to one million blocks). This included implementing and passing analytical Jacobians to the optimiser, and the use of multiprocessing to subdivide and parallelise optimisation of mesh subsections. The mesh is divided into m sets of disjoint (non-overlapping, non- adjacent) node clusters, where 𝒎𝒎 is the number of parallel processors, and each cluster contains 𝒌𝒌 nodes. The coordinates of the 𝒌𝒌 nodes are optimised and updated before a new set of 𝒎𝒎 disjoint clusters are chosen. Calibrating tolerance parameters carefully, we are able to achieve convergence of an optimal mesh in less than 10 iterations (for a 10 4 node mesh we have achieved optimisation in just over 90 seconds). In this paper, we present the algorithm details, several test cases, and a challenge mesh: a one million block grid of the TVZ, with realistic topography, and approximate faults and basement contact