In the context of the provision of clean renewable energy, the utilization of the subsurface of urban areas has been gaining more attention. Especially major urban centers like the city of Berlin, Germany, have moved into focus because of the major share in global CO2-Emissions. One of the most promising mitigation factors for at least part of these emissions, is to rely on geothermal energy resources. The subsurface of Berlin is characterized by a complex geological setting beneath which different, coupled heat transport processes interact and are also overprinted by anthropogenic activities to a certain degree. A major challenge herein is the specific present-day utilization scenario, featuring the production of groundwater from the shallow subsurface (drinking water supply) and extraction of shallow geothermal heat from the underground. Therefore, to predict and reconstruct observed temperature, pressure and mass distributions a precise knowledge of the geological configuration and physical parameters of the subsurface is key. In this thesis I attempt to provide answers to the aforementioned issues by relying on numerical modeling studies which are based on a detailed 3D reconstruction of the underground and on physical principles of heat-, fluid- and mass transport. I present models of increasing complexity and detail concerning all physical processes and boundary conditions at work, enabling detailed representations of the present state, and a reconstruction of the natural state of the system, that is, before any human intervention. Based on the understanding gained from these studies, I derive some conclusive asserts on sustainable utilization scenarios in the underground of Berlin, and based on similarities, of other large urban centers. An initial set of models investigated, for a first time, the effect of major surface water bodies (lakes and rivers) on the thermal and hydraulic configuration of the subsurface. The results of these models show that major surface water bodies significantly modify the shallow to intermediate geothermal and hydrogeological setting especially where connected to areas showing a certain degree of overprinting by human activities. Due to groundwater production and associated lower hydraulic heads near the surface water bodies, forced infiltration from the latter is predicted, which fits observations closely. The thermal modifications are on the order of interest for shallow geothermal installations and also connect to a shift in areas most promising for exploitation. The second set of models focused on reconstructing the natural state, that is, the state before pumping activities commenced. This step was deemed necessary in order to be able to quantify the modifications to the subsurface hydrothermal configuration as related to human activities. The reconstructed "natural" state shows a complete replenishment of the depression cones, resulting in a shift of recharge and discharge areas, where rivers and lakes display gaining conditions only. In contrast, models that integrate pumping activities, illustrate that the effects of subsurface production is larger, both in magnitudes and areal as well as depth extents, than previously captured by models that relied on fixed hydraulic surface boundary conditions only. The presence of active wells provides a more realistic representation of flow rates, the net results of which is a sensitive modification of predicted fluid pathways, in agreement with the monitored hydraulics in Berlin (as exemplary demonstrated for the site of Karolinenhöhe). In a latter stage, I carried out a quantitative study on the energetic potential of the underground, by conducting a systematic analysis of the geothermal potential stored at different levels beneath the city of Berlin. Two approaches were chosen to identify and quantify promising areas for geothermal utilization. The first approach makes direct use of modeled temperature distribution at different depth levels. These highly nonuniform distributions reflect a heterogeneous potential of heat in place in the underground. In the shallow subsurface, utilization for shallow geothermal heat production and storage could be opted for. In the deep subsurface heat and electricity production are deemed possible. Building on this, the second approach investigated two possible deep geothermal target horizons by simulating a hydrothermal power-plant utilizing the predicted temperatures as well as reservoir geometry and physical properties. The connected geothermal potential of the deeper subsurface then shows nominally promising results, depicting up to ~10 MWth for a single virtual power plant at the most promising locality, while locally 0 MWth are encountered as well. All of these model results show most importantly that the subsurface thermal, hydraulic and mass distribution is highly sensitive to the parameters under study, which highlights the amount of caution that should be given to any planned change in the utilization of this space. This relates to any planned geothermal utilization of the different groundwater compartments which should be studied in depth using the models of this thesis as starting conditions.