Modelling physical processes such as flow and transport within a fracture network can be challenging, both conceptually and numerically. A typical approach is to upscale the properties of the network onto a regular grid of elements, which is then used to simulate the required physical processes. However, this method can be inaccurate if not carefully applied. For example, when most of the flow passes through a few choke points in the network, upscaling may introduce extra connectivity that is not present in the original network. The ConnectFlow software package has an alternative method for representing fracture networks, whereby the fractures are explicitly modelled as a network of intersecting two dimensional planes. The algorithms within ConnectFlow are very efficient, allowing millions of separate fractures to be simulated, each discretised using hundreds of finite elements. Here we present recent updates to this functionality: 1) to allow the advection-diffusion equation to be solved for multiple solute species (which is fully coupled to the pressure solution via the buoyancy term in Darcy’s equation); 2) to model the diffusion of solutes into the pore space of the surrounding rock (i.e. matrix diffusion); and 3) to carry out chemical reactions between solutes and minerals (which coat the fractures and/or the rock pores) using an interface to the IPhreeqc library. The ConnectFlow algorithms have also been parallelised to improve the tractability of this new functionality. This implementation represents a significant step forward in capability that allows groundwater flow, transport and hydrogeochemical reactions to be properly represented in the context of structurally constrained fractured bedrock. This has a wide range of potential applications, particularly for future safety assessments of nuclear waste facilities. Example calculations are presented for the Onkalo disposal facility in Olkiluoto, Finland, and the proposed Swedish repository for long-lived waste, SFL.