Natural convection of water plays a key role in shaping the thermal structure of the subsurface. As a consequence of water density decreasing with temperature and temperature increasing with depth in the Earth’s crust, an unstable situation develops: cold, dense fluid overlies hot, less dense fluid. Given a sufficiently permeable and vertically extensive host, hot water will buoyantly rise towards the surface as cold water sinks to replace it, forming a convective circulation cell. Rock units with permeability far too low to facilitate convection may, however, contain high- permeability faults or fractures capable of sustaining convective fluid flow. We investigate the effects and emergent behavior of fracture convection by means of numerical simulation using the Complex Systems Modeling Platform (CSMP++). Despite being confined to a narrow, planar feature, natural convection in a fault or fracture significantly impacts the thermal structure of the subsurface via conduction heat exchange with the host rock and enhances fluid mixing. We compare convection in wide, filled fault zones and narrow, open slots and show that convection patterns within the fracture and the strength of this subsequent thermal perturbation are a function of transmissivity – permeability times thickness – rather than either permeability or thickness alone. Fluid circulation may lead to temperature changes of several tens of degrees in the host up to a distance roughly equal to half the height of the fracture, significantly changing the local geothermal gradient. Above the fracture, these thermal perturbations influence conductive heat flow and take on a diffuse pattern, appearing much wider than the originating fracture or fault. By varying host rock permeability, we show that fluid exchange between fracture and host varies with transmissivity, but has little effect on convection patterns. Conversely, Rayleigh convection within the host rock dominates and overprints convection within the fracture. We build on these results by investigating convection in five elliptic fractures in an en echelon array. Thermal perturbations arising from convection in each fracture eventually connect and align, thereby imposing mutual feedback of convection in neighboring fractures, despite a lack of hydraulic connectivity between them. This interaction leads to a self-organizing pattern of convection cells, with regions of upflow and downflow aligning with the upflows and downflows in their neighbors, inducing a regular pattern of alternating cooled and heated host rock regions between the fractures that is aligned (sub-) perpendicular to the en echelon fracture array. By varying fracture spacing, we determine that this “synchronization” effect may occur iii when fracture spacing is less than convection cell width. Even in fractures with heterogeneous permeability, convection cells continue to self-organize despite zones of low-transmissivity. Fractures which alone cannot facilitate convection may nevertheless host weaker fluid circulation cells in response to the thermal perturbation caused by convection in a neighboring fracture. Linear elasticity is then introduced to our fracture convection model. Using the Barton- Bandis relationship, we vary fracture transmissivity within the fracture as a function of effective normal stress and depth and show that convection strength and vertical extent depend on fracture stiffness, with convection becoming increasingly confined to shallow depths as stiffness decreases. We then include thermo-elastic strain, allowing temperature changes due to convection to induce thermal strains and alter effective normal stress and transmissivity, resulting in high-transmissivity channels which facilitate downward flow of cold water. Using a simplified fracture propagation model, we show that these channels may induce contractive thermal strains sufficient to induce tensile failure in the host rock below, propagating and deepening the fracture and allowing further infiltration of cold water. We then model injection and production of cold water into a single fracture, focusing on thermo-elastic strains and their associated stresses. Injection of cold water into a fracture induces contractive thermal strains behind the thermal front, generally resulting in decreased normal stress and enhanced transmissivity. We identify a zone of increased normal stress in a ring around this thermal front, decreasing fracture transmissivity. This “stress-ring” is shown to be responsible for a significant portion of the difference between models with and without thermo- elasticity, with a particularly important role in fractures where flow channeling is minimal. Thermo-elastic effects become stronger with increasing host rock bulk modulus, but the relative effect of the stress-ring remains the same. The work presented here chiefly investigates the thermal impact of natural convection in a fracture (chapter 2), the interactions between multiple convecting fractures (chapter 3), the effect of thermo-elasticity on convection patterns (chapter 4), and the impact of thermal strains on fracture flow patterns between a well doublet (chapter 5). We then conclude by summarizing our key findings and look forward to the ways future studies may build upon the work presented here (Chapter 6).