The Netherlands has set the ambitious goal to be CO2 neutral by 2050 and signed the Paris Treaty in 2015. The contribution of geothermal energy to reaching this goal are outlined in the Masterplan in 2018 which attempts to reduce CO2 emissions. It is imperative to enhance geothermal participation in renewable energy resources, so this thesis proposes a consideration of large scale geothermal field development in order to meet these requirements. For these large scale projects, large scale geological heterogeneities must be taken into account in order to propose a development strategy that honours subsurface variability in properties like porosity or permeability of an aquifer. Furthermore, the nature of the large scale operations, inherently requires the consideration and application of well patterns typically used in oil and gas developments. Operating in these extensive domains, carries a lot of uncertainty in the final economical output of the project, so modelling the process could indicate the optimal conditions that would deliver the best possible operational outcome. Conceptual 2D model approaches were adopted to demonstrate the main ideas behind large-scale geothermal well pattern optimisation. The main objective of this project is to model, evaluate and optimise the performance of large scale geothermal field development. The proposed strategy is based on the use of well patterns as are frequently used in the oil industry. The heterogeneity in geological properties that may be expected to be encountered at larger spatial scales is addressed by the concept of a flexible well density function. This density function allows the well patterns to be resized (or ultimately, reshaped) and adapt to spatial variations in geological characteristics. The flexible well pattern is fed to an objective function created with embedded simulator. The output of the function is the net present value (NPV) of the project. Four test cases are created, starting with a homogeneous static model and building up more heterogeneous aquifer models, aiming to test the performance of the flexible well density function. The aquifer property models are representative of the West Netherlands Basin and specifically the Delft Sandstone Member. Per each aquifer model, line drive and 5-spot development strategies are assessed. Each development scenario is modeled in an objective function and optimised. The optimisation algorithm chosen is the Simplicial Homology Global Optimisation, suitable for black-box functions that show multiple local optimum solution and among them, a global optimum pattern size is found. The NPV of each project realisation is calculated based on the energy recovered and an economic model under Dutch fiscal conditions. The results suggest that, the flexible well placement is successfully aligning with the different aquifer geological properties. Pattern size is inversely correlated to porosity which corresponds to higher volume of pore fluid from which heat can be recovered. The optimisation algorithm managed to identify the global optimum solution of pattern size that delivers the highest possible positive NPV. The most efficient in terms of profitable strategy, is suggested to be the 5-spot pattern. The optimal pattern size ranges between 500 - 2160m depending on the aquifer model. The sweep efficiency, in terms of energy recovery, is also assessed per aquifer model and development strategy. The most efficient is the 5-spot pattern. The performance of the adopted optimisation algorithm, on the fully homogeneous aquifer, is tested with an exhaustive response curve of NPV. It is confirmed that the algorithm manages to identify the global optimum. The character of the NPV as a function of well pattern size/density proved the complexity of the system with respect to the amount and well types introduced. Different sensitivity analyses in the context of the uncertainty of the aquifer thickness, thermal conductivity, thermal heat capacity and full economic model are conducted in order to show the impact on the optimal patter size. The performance of the optimisation algorithm is assessed as well, indicating that further investigation on the tuning meta-parameters could potentially lead to better global optimum solutions in more heterogeneous aquifer models.