The numerical modeling of radon concentrations in the fault zone of the underground excavations at Książ Castle was conducted using a stochastic Discrete Fracture Network (DFN) model. Due to the difficulties related with obtaining the exact fractures in a rock mass, the novel approach used in this study incorporates the stochastic model with known site data. The analysis utilized a dataset comprising long-term measurements of 222Rn activity concentration and geodetic measurements for twelve faults in the Książ unit. The parameters considered in the DFN model are: fracture length, Peclet number (Pe = 0.1 and 1.0, respectively), advection velocities (from 10–8 m/s to 10–6 m/s and from range from 10–7 m/s to 10–4 m/s, respectively), radon diffusion (D = 2.1 × 10–61/s), radon decay constant (λ = 1/s), and radon gas generation (q) along the fractures within the range of 1.5 × 10–3 Bq/m3·s to 3.5 × 10–3 Bq/m3·s. The calibration process obtained the best fit when the radon generation rate was uniformly distributed through the rock mass in addition to incorporating a higher value of radon generation rate (q = 3.0 × 10–3 Bq/m3·s) where elevated radon concentrations have been measured. The modeling results also confirmed that the radon generation rate should always be higher where elevated radon activity concentrations were measured regardless of the measurement period. For the indicated “area” the radon generation rate should be higher from 25% to 37.5% between May–October and 18.5% to 40% between November–April. The influence of fracture zones on the recorded radon activity concentrations was noticeable up to a depth of 15 m. Within this range, the highest values of 222Rn activity concentration, ranging from 1,600 Bq/m3 to 2,000 Bq/m3, were consistently observed regardless of the season. However, as the depth increased, the values of 222Rn activity concentration decreased from 800 Bq/m3 to 400 Bq/m3 and became more dispersed.