Modelling and Validation of CANDU Shim Operation Using Coupled TRACE/PARCS with Regulating System Response

,


Introduction
Te CANDU reactor has unique characteristics in its operation, primarily due to the use of natural uranium fuel and heavy water moderators. Most notably, CANDU reactors are fuelled online, with the capability to refuel individual fuel channels multiple times per day while at power. Tis "continuous refuelling" allows for excess reactivity to be kept at a minimum, improving the overall neutron economy to achieve a higher exit burnup out of the natural uranium fuel. CANDU reactors are also loosely spatially coupled, thus requiring spatial reactivity control to maintain a nominal fux distribution with no tilts. When modelling long-term reactor operation, it is important to model the core evolution, including depletion, saturating fssion products, and the actions of the reactor regulating system (RRS) [1].
During normal reactor operation, the power of the reactor is controlled primarily using 14 liquid zone compartments. Tese may be flled and drained with light water, which acts as an absorber in a heavy water-moderated reactor. Te RRS liquid zone control algorithm acts to control both the bulk reactor power as well as its spatial distribution and can counteract the reactivity imbalances that result from refuelling along with other minor reactivity perturbations. Adjuster rods, normally inserted, are used to fatten the fux profle in the core and may be withdrawn to provide additional reactivity when necessary. Mechanical control absorbers, normally withdrawn, are used to provide negative reactivity, and may also be partially or fully dropped to the core to rapidly reduce power without fully shutting down the reactor. Moderator poison may also be used to introduce negative reactivity.
Online refuelling is achieved using two fuelling machines, one on either end of the reactor. When refuelling a channel, one fuelling machine will push new fuel bundles into a channel, while the other machine will receive the used bundles. One possible operational occurrence is a fuelling machine outage, where a fuelling machine is unavailable due to maintenance and thus the reactor cannot be refuelled. In this scenario, the reactor may be operated in a "shim operation" mode, using the adjuster rods for reactivity shim. During the initial stage of shim operation, the reactivity loss due to depletion is compensated by draining the liquid zone compartments. When additional reactivity is required beyond what the liquid zone control can provide, banks of adjuster rods are sequentially withdrawn. Tis is continued until the reactor is either shut down or refuelling resumes.
A number of previous studies have been performed on adjuster rods, shim operation, and/or other adjuster-driven transients under various conditions. Some examples include (i) Evaluating reactivity devices using DUPIC fuel compared to natural uranium fuel [2]. Te study modelled CANDU-6 reactivity devices and shim operation using coupled RFSP and NUCIRC, fnding that the maximum shim duration was signifcantly shorter for DUPIC fuel. (ii) Evaluating reactivity devices using DUPIC fuel for xenon override in a CANDU-6 [3], using RFSP, fnding comparable xenon override capability to the natural uranium cycle. (iii) Optimisation of adjuster rods in thorium-based fuel cycles to achieve comparable xenon override capability and shim duration to natural uranium fuel cycles [4], using DRAGON and DONJON.
Unlike these studies, which used CANDU-specifc tools and evaluated theoretical shim operation performance for advanced fuel cycles when compared to natural uranium, this study focused on implementing the methodology using the PARCS and TRACE codes, which are typically used for the analysis of light water reactors rather than CANDU reactors, as a basis, and modifying them to support the necessary coupling capabilities. Tis study involved simulating the shim operation of a 900 MW class CANDU reactor using a standard natural uranium fuel cycle and comparing with operational data, to demonstrate that these tools could be adapted for modelling long CANDU transients and considered for analysis of either actual or theoretical reactor transients. In addition, the use of coupled transient simulations to supplement quasi-steady-state calculations is a novel contribution, as it allows for the reactor dynamics to be properly captured at key points in an otherwise quasistatic reactor evolution.
Te simulation and analysis of reactor operation may be performed using one of two methodologies based on the evolution of the reactor state. If the reactor fuel composition and transients evolve very slowly, then the state of the reactor at any given time can be modelled as a quasi-steady state. Under such a simulation, the burnup and saturating fssion product state of the core, which evolves slowly, is held constant, while more rapidly evolving variables such as the neutron fux and thermal-hydraulic conditions are solved for their equilibrium values. Te fux is then held constant and used to advance the burnup and saturating fssion product state by a given period of time. If the core state evolves rapidly, the reactor is simulated as a transient, with the neutron fux and thermal-hydraulic conditions evolving with time, and the efects of fuel composition changes are secondary.

Background.
Te CANDU reactor is a pressurised heavy water reactor (PHWR) design, difering signifcantly from the more popular light water reactor (LWR) designs. Some key traits of the CANDU design that difer from the PWR design include [5] (i) Heavy water is used as the coolant and moderator (ii) Natural (unenriched) uranium as fuel (iii) Horizontal fuel channel arrangement, with fuel channels serving as the pressure boundary. Coolant and moderator are physically separated, with the fuel channels surrounded by the moderator in an unpressurised calandria. (iv) Fuel channels have a signifcant degree of separation, to allow for moderation by heavy water, access by the fuelling machines, and the insertion of vertically oriented reactivity devices (v) Online refuelling, through the use of two fuelling machines which can latch onto either end of a fuel channel to insert and discharge fuel bundles while the reactor is operating at power (vi) Reactivity control primarily uses 14 liquid zone compartments, whose light water levels can be varied independently to control reactor power and its spatial distribution Online refuelling is important due to the use of natural uranium fuel. Tis fuel has far less excess reactivity when compared to enriched uranium fuel, leading to lower exit burnups and core residence times for CANDU fuel when compared to LWR fuel. A batch refuelling scheme for a CANDU reactor would require a refuelling outage every few months [5], which would greatly diminish the station's capacity factor. As well, the continuous refuelling allows for the excess reactivity on a core-wide basis to be kept at a consistent minimal value, improving the overall neutron economy and thus the exit burnup.
Online refuelling is achieved using a pair of fuelling machines, with one on each end of the reactor. Each fuelling machine is capable of latching onto one end of a fuel channel, removing the channel closure and shield plug, then either inserting new fuel bundles or receiving used fuel bundles, before replacing the shield plug and channel closure [5]. Typically, not all fuel bundles in a channel are replaced, but rather only 4 or 8 bundles (out of 12 or 13 bundles) are replaced, with the remaining bundles closer to the inlet end of the channel remaining in the core but being shifted towards the outlet end of the channel. Te horizontal orientation of the fuel channels allows for bidirectional refuelling, where adjacent channels are refuelled in opposite directions, along with coolant fow being in opposite directions, to achieve symmetry in the neutron fux distribution.
Tis introduces additional complexities in simulating reactor operation which is not present in the simulation of LWR operation. For an LWR, the fuel's physical location remains the same over the course of a fuel cycle, with the fuel being depleted over the course of the cycle and the core kept critical through the use of moderator poison, reactivity devices, and/or negative power feedbacks, depending on the type of reactor. In a CANDU reactor, the fuel distribution changes on a daily basis as individual channels are refuelled. Tus, high-fdelity simulations require a full accounting for both efects of burnup and fuel relocation due to fuelling prior to the event of interest, as well as the transient of interest. Diferent methods for performing the core follow are as follows: (i) Using actual station data, the operating history of the reactor is simulated, from the start of operations or from a given snapshot in time to a later snapshot in time using the fuelling and operating histories of that core. Tis may be performed using any difusion code, with the industry typically using RFSP, but studies have also been done using DONJON [6], to model the core state leading up to a transient. (ii) Starting with an initial core of fresh fuel, reactor operation is simulated, including the simulation of fuel channel selection for refuelling, until the core reaches an equilibrium. Tis typically takes 400 to 500 full-power days [5]. Such fuelling simulations have been previously performed using PARCS [7].
When simulating reactor operation, it is thus often important to have frequent snapshots which capture the salient core features at each discrete point in time. For a given snapshot, the slow-changing quantities, such as burnup, may be held constant, and the equilibrium state of the reactor is determined. Ten, a depletion step may be performed, where the slow-changing quantities are updated under the assumption that the equilibrium state of the reactor changes little between successive snapshots. Te burnups in one or more channels would then be updated to simulate refuelling. Depending on the model and particularly on the time between snapshots, the saturating fssion product densities may be treated as a slow-changing quantity to calculate transiently, or as a fast-changing quantity to calculate as an equilibrium. For Xe-135, the density reaches equilibrium on the order of a day for a reactor at full power, based on an approximately 6.7-hour half-life of I-135 [5]. Tus, for daily or less frequent snapshots, an equilibrium xenon model is reasonable except for recently refuelled channels.
In addition to determining the fux distribution in a core snapshot, it is also important to model the reactivity devices, since reactivity devices have an important efect on the fux distribution. Te following are used to control reactivity in the CANDU design: (1) 14 liquid zone compartments (LZCs) distributed throughout the core. Tese are flled with light water, which acts as a neutron poison due to its high absorption cross section compared to heavy water. Te unflled portion of each compartment is flled with helium gas. Te liquid levels are varied automatically by the reactor regulating system (RRS) running on a digital control computer (DCC). Both bulk control and spatial control are used. Bulk control adjusts the liquid levels equally in all 14 compartments to add positive or negative reactivity and control the reactor power. Spatial control adjusts the liquid levels independently to balance the spatial distribution of the fux [1]. (2) Adjuster absorbers, which are normally inserted in the inner part of the core to provide fux fattening. Flattening of the fux distribution is desirable as it allows for the total reactor power to be increased without exceeding power limits for individual fuel channels or fuel bundles. Adjuster absorbers are grouped into eight banks. Te RRS will withdraw one or more banks of adjusters if additional positive reactivity is required, based on either a low average liquid zone level (indicating that the LZCs cannot provide enough reactivity) or a large negative power error (indicating that positive reactivity is needed more rapidly) [1]. Withdrawing some of the adjusters requires the reactor power to be derated due to a loss of the fux fattening efect. (3) Mechanical control absorbers (MCAs), which are normally positioned out of the core. Te RRS will insert one or more banks of MCAs if additional negative reactivity is required, based on either a high average liquid zone level (indicating that the LZCs cannot provide enough reactivity) or a large positive power error (indicating that negative reactivity is needed more rapidly) [1]. (4) Moderator poison, either boron or gadolinium.
Boron is used on initial start-up to depress the excess core reactivity from a full load of fresh fuel during the frst weeks of operation. Gadolinium is used to compensate for transient excess positive reactivity beyond the worth of the LZCs and MCAs, particularly to manage the xenon transient resulting from an outage and subsequent restart [5].

Science and Technology of Nuclear Installations
For managing the normal operation of the core, the LZCs are most important, with the other reactivity device states typically not being actuated/changed. Te LZCs are used to manage small reactivity perturbations, particularly those that occur from depletion and refuelling. When a channel is refuelled, not only does the new fuel introduce positive reactivity, but it also locally increases the fux near the refuelled channel. For a large reactor such as the CANDU-9 reactor, the core radius is far greater than the neutron migration length, and the reactor is susceptible to certain fux harmonics, particularly the frst azimuthal (sideto-side or top-bottom) and frst axial (end-to-end) modes [8]. Terefore, spatial control, via the LZCs, is required to control these fux modes. Conceptually, the reactor is divided into 14 spatial control zones, with one LZC corresponding to each zone. Spatial control is used to maintain a design ratio of measured powers between the 14 zones, suppressing any tilts resulting from refuelling, xenon perturbations, or the movements of other reactivity devices. When calculating a snapshot, the 14 LZC water levels must be determined, both to achieve criticality as well as to eliminate any fux tilts.
It is also of interest to simulate operating conditions other than normal operation at full power. One such example is the simulation of a shim operation. CANDU online refuelling relies on a pair of fuelling machines. If the fuelling machines are not available, then the reactor will need to be operated without refuelling. In shim operation, the negative reactivity due to depletion is counteracted by the LZCs and the adjusters. Reactor power must be derated prior to withdrawing adjuster banks to counteract the resulting increase in channel and bundle power peaking factors. Tese power reductions also lead to xenon perturbations that must be managed by the LZC system. Once the fuelling machines are available, the reactor can gradually be refuelled until it reaches its normal level of reactivity. If the fuelling machines are unavailable for an extended period of time, the reactor may need to be shut down. In the case of planned maintenance, extra refuelling can be performed prior to shim operation, using LZCs or moderator poison to compensate for positive reactivity. Tis will extend the possible duration of the shim operation. Figure 1 provides the liquid zone designation for the CANDU 900 reactor used in this work.

Methodology for Modelling and Simulation.
Tis work uses the TRACE [9] and PARCS [10] computer codes, licensed through the NRC's Code Application and Maintenance Program (CAMP), to simulate system thermalhydraulics and core physics, respectively. More specifcally, this work uses modifed versions of the codes, hereafter named TRACE_Mac1.1 and PARCS_Mac1.1. It also uses the Exterior Communications Interface (ECI) [11] to couple an RRS model to PARCS through TRACE. Te RRS model was developed in the Python programming language and uses a Python ECI library built around the original Fortran-based ECI library. Tis coupled model has been previously evaluated against shorter transient events such as the fgure of eight fow oscillations (using a modifed system model without header interconnects) and a loss of fow transient initiated by a loss of Class IV power [12]. Tis section covers the additional work performed to model long-duration transients such as adjuster shim operation.
For modelling of CANDU adjuster shim operation, there are two simulation regimes to consider. During the majority of the operation, the core state evolves slowly, and the only signifcant efects to consider are the efects of depletion and saturating fssion product evolution. Te second regime corresponds to a time-dependent situation where device confguration (e.g., adjuster positions and liquid zone levels) change over a relatively short period of time due to a change in RRS control actions.
During the slowly evolving quasi-steady-state periods the reactor can be modelled as a steady-state simulation. During this simulation, the depletion, saturating fssion product densities, and reactor power are held at constant values, while the spatial fux distribution, thermal-hydraulic state, and control system state are computed. To model the core evolution during these steady-state steps, a PARCS depletion calculation is performed with PARCS_Mac1.1, with the reactivity device and thermal-hydraulic feedbacks held constant. Te depletion calculation is used to calculate an updated depletion and saturating fssion product density state for the next steady-state step. Tis methodology is very similar to a core follow, as it generates periodic core snapshots modelling the evolution of the core over a period of time. Core follow analysis has been previously performed using both PARCS [7] as well as other core physics tools, particularly RFSP.
For the steady-state model, certain simulation parameters have been adapted to support eigenvalue calculations, accelerate model convergence, and improve performance. a zone level error for steady-state analysis, to converge either the system k ef or average zone level to a target value. Te zone levels may also be frozen to a set of known values. (ii) Te PARCS eigenvalue calculation is performed only once every 20 TRACE time steps. Tis greatly accelerates the simulation as the PARCS eigenvalue calculation is signifcantly more computationally intensive than the TRACE time step. Sensitivity studies demonstrate that the transients predicted here are not sensitive to the exchange rate in this regard, as long as the three models (neutronics, thermal-hydraulics, RRS) are given sufcient time to converge and that the neutronics-RRS feedback loop is stable. (iii) Te RRS time step size is decoupled from the TRACE simulation and set to a fxed constant. Tis allows for the RRS to converge more quickly when the TRACE time step size is small. Te maximum step size is limited by stability due to the feedback between the RRS and PARCS models and is currently set to 4 seconds per PARCS step (0.2 seconds per TRACE time step). (iv) Convergence of thermal calibration of the ion chambers and in-core fux detectors (ICFDs) is greatly accelerated. Calibration is applied every two seconds instead of every 16 seconds, and the time constants of the flters in the thermal calibration routine are decreased to 1.875 seconds (instead of 30 seconds or 180 seconds). (v) Calculation of the delayed components in the ICFD response is performed using the fux from the steady-state calculation of the previous depletion step and the fux from the current calculation, under the assumption that the fux changed linearly over the course of the previous depletion step. Te ICFD detector response in this work is based on the response characteristics calculated in reference [13], with the four normalised cases averaged. Te resulting prompt fraction was 1.048 (4.8% overprompt), and the response also includes four different delayed components with time constants ranging from minutes to days.
Tis approach is very similar to RFSP's asymptotic zone control functionality [2], except that the liquid zone control module functions identically in the steady-state and transient modes, only manipulating its inputs (bulk power error and time step size) to implement the steady-state mode. Tis guarantees that the asymptotic zone levels are consistent with transient modelling, as long as the power error converges to zero.
Te second regime corresponds to time-dependent actuation of the device positions caused by RRS or operator actions. As the reactor is depleted, the reactor regulation system is expected to decrease the average liquid zone level in order to maintain criticality. Once the liquid zone level reaches a certain threshold, a bank of adjuster rods is to be withdrawn. Te dynamic nature of the adjuster rod withdrawal process is simulated using a coupled transient model with transient depletion and saturating fssion product evolution. In advance of withdrawing the adjuster bank, the reactor power will be derated by the operators to an appropriate level for withdrawing the adjuster bank consistent with station operating rules. Tis is done at a higher liquid zone level threshold, to allow the xenon-135 transient to trigger the adjuster withdrawal and provide enough reactivity margin so that only one adjuster bank is withdrawn. Te full transient procedure is thus as follows: (1) As the zone level nears a threshold where action to remove the adjuster may be required, the simulations will decrease reactor power to the setpoint for the next adjuster bank to be withdrawn to mimic operator actions. (2) Te simulations will then continue in time to allow the average LZC level to decrease to the adjuster withdraw threshold as the xenon transient progresses.
(3) Te simulations will then withdraw an adjuster bank at the time and rate dictated by the design parameters. Te average LZC level will increase to maintain criticality. (4) Average LZC level will continue to decrease to a minimum value corresponding to peak xenon, then subsequently increase as xenon density returns to a new equilibrium value. (5) As xenon approaches equilibrium, the average LZC level will begin to decrease again once depletion becomes the dominant efect.
For this work, the coupled transient is simulated from the start of the power derate until shortly after each adjuster bank is withdrawn. Te most rapid changes in fux distribution and reactivity device positions occur during adjuster bank movement, and it is desired to model the behaviour of the fux distribution as the bank moves out of the core. A signifcant transient fux tilt may degrade safety margins, and if sufcient fux tilt is detected, an RRS setback may be triggered, further reducing power below the intended setpoint. In this model, the setback threshold for fux tilt is set to 20%, based on zone power measurements excluding the centre zones 4/11. Once a setback is initiated, it continues until the fux tilt drops below 18%. Once a bank is fully out of the core, the spatial fux distribution is relatively stable, and the depletion and saturating fssion product evolution may be modelled using a steady-state depletion model once again.
Terefore, the overall procedure for the shim simulation is as follows: (1) Perform steady-state simulation and depletion steps while monitoring the average zone level (AZL). (2) As AZL approaches the threshold to derate power, determine the depletion step size to reach the AZL threshold. Iterate until the calculated AZL is at the threshold.
Science and Technology of Nuclear Installations For step 1, a 6-hour depletion interval has been selected. Once step 3 has been completed, step 4 is performed until 12 hours pass from the end of the previous adjuster bank withdrawal before returning to step 1. Te step size starts at 15 minutes and is doubled until it reaches the original 6hour interval again. Te 12-hour delay prevents the next adjuster bank withdrawal from being triggered by the AZL threshold due to the xenon transient, ensuring that each adjuster bank movement is triggered due to depletion and not xenon evolution. Figure 2 outlines a high-level methodology that can be used to implement this simulation. Te simulation alternates between performing a coupled steady-state calculation and a standalone depletion calculation, until the average zone level approaches the threshold to perform the adjuster bank withdrawal process. At this stage, a coupled transient run is performed to simulate the power manoeuvre and subsequent adjuster bank withdrawal.
A more detailed version of the methodology is provided in Figure 3. While similar to the high-level methodology presented in Figure 2, it includes more detailed logic for determining which action to perform after a steady-state simulation. First, it checks that there were no abnormalities in the simulation, such as a setback, stepback, or reactor trip. If such an abnormality is detected, then the simulation will either revert to before the previous depletion step and repeat it with a smaller depletion interval or will switch to coupled transient simulation. Tis allows for the simulation to pinpoint the occurrence of such an event and simulate it as a transient simulation. Secondly, the PARCS steady-state simulation is limited to a fxed reactor power, but the indicated power in RRS does not exactly match the power in the PARCS simulation. Terefore, at the end of a steady-state step, the indicated power in the simulated RRS is compared to its setpoint, and if their diference exceeds a certain tolerance, the PARCS power is adjusted, and the steady-state step is repeated.
Te simulation must then determine whether it has reached a part of the simulation that should be performed as a fully coupled transient. Tis occurs when either of the following conditions is met: (1) Te average zone level is at or below the threshold to perform the next manoeuvre, which consists of a power derate (if needed) followed by an adjuster pull. (2) Te average zone level is above the threshold, but the estimated time until the average zone level reaches the threshold is less than the minimum depletion step size (864 seconds, or 0.01 days). Tis estimate is calculated by extrapolating the average zone level from the previous two steady-state calculations. If the estimate is negative, due to the average zone level increasing over the previous depletion step, then it is ignored.
If neither condition is met, then the next step performed is a depletion step. A nominal depletion step size of 6 hours was selected in this work, but the size of a given depletion step is determined as follows: (1) If the estimated time from step 2 of the previous list is positive and no greater than the nominal depletion step size plus 15 minutes, then the depletion step size is set to this estimated value. Otherwise, the nominal step size is selected. (2) In addition, for a scheduled action that is not to be simulated using a coupled transient (e.g., refuelling), the depletion intervals are sequentially shortened as the action is approached, up until the action time is reached. (3) Te actual depletion step size is the lesser of this selected step size or twice the step size of the previous depletion step.
Tis creates the efect of performing fxed depletion intervals until the next transient portion of the simulation approaches, then performing a depletion interval which is estimated to bring the simulation to the starting time of the transient. After a transient is completed, the subsequent depletion step size is set to 1 hour, then allowed to double with each depletion step until it reaches the nominal step size. For actions not requiring transient simulation and that are based on a schedule, fner steps are taken around the time of the action, both before and after. For this transient, several refuelling events are scheduled. To simulate a refuelling, a script is used to read the binary depletion fle produced after a depletion step and modify the burnups in the selected channel (pushing 4 or 8 new bundles into the channel and pushing the existing bundles towards the end of the channel, discharging channel out the other end) before rewriting the fle.
Since exact controller parameters for the 900 MW CANDU design are not available in the literature, realistic approximate or representative values were selected. Te following values were used:  6 Science and Technology of Nuclear Installations (i) K P � 31.0 (LZC bulk controller gain) (ii) K T � 5.0 (LZC spatial controller gain for fux) (iii) K H � 0.6 (LZC spatial controller gain for level, when fux control is phased out) (iv) K L � 0 (LZC spatial controller gain for level, independent of fux control being phased in) (v) Adjuster movement speed � 10.5 cm/s (vi) LZC drain rate � 5.0 cm/s (vii) LZC max net fll rate � 5.0 cm/s (viii) LZC valve lift/bias for constant level � 60% (ix) Flux tilt setback threshold � 20% FP (diference from highest to lowest zone power, excludes centre zones 4 and 11) Te LZC drain/fll rate was selected based on taking approximately one minute to fll/drain a "long" zone compartment between 0% and 100% [1], the longest of which in this model is 285.8 cm, resulting in a rate of 4.76 cm/s, which was rounded to 5 cm/s. Te most sensitive parameters are the gains and detailed sensitivity results are performed to determine the efect of the gains on the transient. Outcomes related to reactor setback are also highly sensitive to the fux tilt setback threshold.

Station Description and Modelling.
Te station design being simulated is a 900 MW class CANDU reactor, based on those used in previous coupled simulations of TRACE_Mac1.0 and PARCS_Mac1.0 [12]. When comparing preliminary results to station data, it was found necessary to update the core physics model for better consistency with the specifc CANDU 900 under analysis. Tis included increasing the number of fuel bundles per channel from 12 to 13 both in the thermal-hydraulic and reactor physics models, as well as applying appropriate extrapolation distances as these are known to have a signifcant impact on full-core modelling outcomes [14] and was found to have a signifcant impact in this work.
In addition, adjuster rod ageing was accounted for by prorating adjuster incremental cross-sections, based on prior industry investigations, to account for their decrease in worth as a function of irradiation. An adjuster rod age of 182,500 efective full-power hours (EFPH) was selected for this work. Te efect of ageing was to reduce the efective thermal absorption cross-section of the adjuster by Other cross-sections were perturbed by varying amounts. Te RRS model was also updated, particularly to implement an approximate model for fully instrumented channels (FINCHes). Te FINCH powers are approximated as the PARCS computed channel powers with a 30-second decay flter. Te FINCH powers are compared against a nominal power distribution and are used to calibrate the zone power readings from the ICFDs as part of the RRS control features modelled in these simulations.
In the TRACE model, the 480 channels are reduced down to 40 representative channels, based on zone designation, then subdividing each zone into two regions, being the upper and lower half for the centre zone, and inner and outer regions based on average channel power for other zones. All channels in a given region that belong to the same fow loop and pass are grouped into a single channel in the TRACE model. Tis results in a total of 40 channels in the TRACE model to represent the fow through all 480 channels. Te coolant temperature and density, as well as fuel temperature, are shared for all channels in a given grouping.
Te 900 MW CANDU design uses fow-power matching to achieve a similar exit enthalpy in channels with varying average channel power, trimming the fow in lower-power channels [5]. To approximate this fowpower matching, a minor loss representing a fow orifce was added to the lower-power channel groups where the inlet feeder connects to the inlet end ftting such that the outlet enthalpy of all channels was approximately equivalent.
For most CANDU reactors, adjuster bank withdrawal typically requires power reductions to ensure thermal margins are maintained during these periods where radial core power peaking is increased. To simulate this, the simulations perform the following power changes with each adjuster bank:  To simulate this sequence of events, specifc conditions must be identifed for triggering the initiation of each event. Te simulated sequence is as follows: (i) Te simulation starts at the initial known burnup distribution, with 0.35 ppm of boron used to represent the equivalent moderator poison that was deemed to be present. Te moderator poison is gradually reduced to zero, following the known trend of equivalent moderator poison from station snapshots and linearly interpolating between them. (ii) Te seven-channel refuellings on days 5  drops to 10%, the reactor is shut down, as it would go subcritical anyway. Te simulation is ended when the actual reactor power has dropped below 2.0% FP or when the simulation time has exceeded 6 hours from the power reduction to 59% FP.

Figures of Merit.
Te following key fgures were identifed for evaluating the performance of the simulations compared to the station data: (i) Event timings (adjuster bank A, B, C withdrawals and subsequently reaching 25% average zone level) (ii) Trends in average liquid zone levels (iii) Trends in individual liquid zone levels (iv) Channel power distribution In addition, it is possible to calculate reactivity device worths in the simulated cases for comparative study.
For individual liquid zone level trends and channel power distribution, the average zone level of the two central zones (4 and 11) was selected as an indicator of the radial reactivity balance. More specifcally, it was checked 5 minutes prior to the end of the station data. Tis point was selected as the average zone level would be consistent (around 25%) and the level for the two central zones was generally within the range of 30% to 70% for both the station data and the simulations.
Te following phenomena have a signifcant efect on these fgures of merit as they afect the core reactivity: (i) Depletion-the rate of reactivity change due to depletion will govern the overall length of the transient. Science and Technology of Nuclear Installations 9 would be expected near the channel outlets for highpower channels. (iv) Absorption by reactivity devices-each reactivity device afects the core both through its absorption of neutrons as well as the resulting redistribution of fux in the core. (v) Saturating fssion products-when the power is decreased, there is an initial negative reactivity change as the concentration of Xenon-135 increases. After a couple of days, a new equilibrium is approached with a slight net positive change in reactivity. Te xenon-135 negative reactivity transient is also the driving force of the extended portion of the simulation when the power is reduced to 59% FP.
Te duration of each phase of the shim operation (one bank out, two banks out, and three banks out) is determined by the reactivity worth of the adjusters. For one and two banks out, there is also an efect from the power coefcient, void coefcient, and saturating fssion products, due to the power reduction immediately prior to the removal of the frst two adjuster banks. Te rate of reactivity change due to depletion also infuences the duration of the shim operation.
For the two-bank-out portion of the shim operation, four channels are refuelled, replacing 16 depleted bundles with 16 fresh bundles. Te resulting reactivity increases the delay in the withdrawal of bank C and subsequent events of the shim operation.
Te outcome of the 59% FP portion of the simulation is sensitive to the reactivity worth of the remaining fve adjuster banks as well as the xenon-135 transient, along with some contribution from the power coefcient of reactivity. Te power coefcient afects the initial part of the transient-the frst seconds to minutes after the power manoeuvre as the TH conditions settle to a new equilibrium. Beyond that, the core reactivity is dictated predominantly by the xenon-135 transient, and adjuster rod banks are withdrawn to compensate for this transient.
Te average liquid zone level may be used as an approximate analogue for the reactivity balance of the core. Te liquid zones contribute roughly 6 mk of reactivity. Other changes to core conditions that afect reactivity will introduce a compensating change to liquid zone levels from the RRS. However, two considerations must be kept in mind when considering the liquid zone levels as a reactivity measure.
(1) Te reactivity to zone level relationship depends on which zones are being flled and their levels. For example, the top of zones 3 and 10 are near the edge of the core and extend into the refector, thus the marginal reactivity worth of these zones is low when their fll level is high.
(2) Te model may overestimate or underestimate the reactivity worth of the liquid zones, thus underestimating or overestimating the change in zone level resulting from a given reactivity change.
Tus, the analysis of any discrepancies in the average zone level should consider the entire simulation to determine whether a discrepancy is due to the non-LZC reactivity change or due to the LZCs themselves. For example, if the LZC level change is consistently underestimated throughout the transient, it suggests that the worth of the LZCs themselves is being overestimated.
Te individual liquid zone levels may be used to evaluate the spatial distribution of reactivity, due to the distribution of fuel burnups, adjuster rods, and xenon-135. Te LZCs are used to spatially control the fux distribution to compensate for perturbations due to refuelling and xenon-135 and will also attempt to compensate for adjuster movements. Any diferences in individual zone levels from the data suggest a reactivity discrepancy for part of the core.

Sensitivity and Uncertainty Analysis.
For this study, the sensitivities of the shim operation to diferent uncertainties and parameters must be identifed. Tese may be compared to any diferences identifed between the model and station data to determine any potential sources of discrepancies. Te sensitivities may also be evaluated by comparing the different models to one another. Some of the key parameters being evaluated are as follows.

Burnup Distribution.
To evaluate the efect of the distribution of fuel burnups, a simulated core follow calculation was performed on the reactor. Using a random snapshot as a starting point, thousands of days of core operation and refuelling are simulated prior to the shim transient. Tis generates a set of burnup distributions that can be used as initial conditions for simulating the shim operation. Te results from diferent initial core snapshots may be used to identify whether the fuel burnup distribution can potentially afect the evolution of the transient. If different snapshots behave similarly during shim operation, it suggests that the actual station burnup distribution was unlikely to have a signifcant efect.

Nuclear Data.
Te neutronic calculation is highly reliant on the nuclear data used in the lattice physics code (SCALE/TRITON [15]) to calculate the homogenised crosssections for the difusion code (PARCS [10]). SCALE includes a stochastic uncertainty analysis tool, called Sampler, which runs a SCALE sequence using one or more perturbed sets of nuclear data. Tese data sets are preperturbed using the covariance matrix that accompanies the nuclear data and are distributed with SCALE.
For this work, lattice calculations were performed using 60 perturbed nuclear datasets from the SCALE 252-group library based on ENDF/B-VII.1. Te depletion and branch calculations are then performed using each dataset, for both the infnite lattice cell and the lattice edge cell, along with performing refector calculations. Te homogenised fuel and refector cross-sections are then combined with the device incremental from the supercell calculations, which are not perturbed by Sampler.
Tis creates a set of perturbed PARCS PMAXS crosssection datasets which may be used either to simulate the station transient from its actual burnup state or to perform simulated core follows to generate snapshots to simulate the station transient from diferent burnup states. A new core follow analysis is performed for each dataset to construct core snapshots consistent with the perturbed fuel properties, as described in [7].
It should be noted that the Shift module of SCALE [16] was not released and hence the nuclear data uncertainty impact on the device incrementals could not be directly calculated. Instead, incremental cross-sections were perturbed directly as described in the paper.

Incremental Modelling and Uncertainty.
While it is not possible to perform a proper uncertainty analysis for the reactivity devices in the current methodology, it is possible to perform a sensitivity study. Te device incrementals can be perturbed in order to increase or decrease their worth. Tere are several uncertainties in the modelling of reactivity devices.
(i) Uncertainties in the true device geometry due to manufacturing tolerances (assumed to be small). (ii) Uncertainties in the lattice supercell calculation, either from the microscopic cross-section data or from modelling approximations. (iii) Modelling uncertainties introduced by the application of the device incremental cross-sections to the PARCS core model, e.g., from the expansion of the device's efective volume to ft to the PARCS mesh. (iv) Adjuster rod irradiation.

Incremental Cross-Section Placement in the Lattice.
Te lattice supercell calculation uses a repeating lattice of 2 lattice pitches across by 1 lattice pitch vertically by 1 bundle length axially (2 × 1 × 1 region) and homogenises a 1 × 1 × 1 region centred on the reactivity device. Incremental crosssections are determined by comparing the macroscopic cross-sections with and without the device. It is assumed that this supercell geometry is representative of all possible device locations in the core, and the same incremental crosssections are applied for a given device for all PARCS nodes, even refector nodes, except that fssion incrementals are zeroed for the refector nodes.
While the homogenisation region is a 1 × 1 × 1 region centred on the reactivity device, since the device is vertically oriented (Y-direction in PARCS), running perpendicular to and between the fuel channels (Z-direction in PARCS), it is not possible to apply the incremental cross-sections directly to a 1 × 1 fuel channel region. Instead, the incremental crosssections are distributed to both adjacent fuel channels, with a weight of 0.5 assigned to each channel. Furthermore, as most reactivity devices are not perfectly centred on a single axial node, the incremental cross-sections must be further distributed between two axial planes, with the relative weights determined by the actual device location. Te distribution of device incrementals over a larger volume has a tendency to weaken the spatial self-shielding efect, increasing their predicted reactivity worth.
Te following three models are evaluated by modifying the axial (Z-direction, parallel to the horizontal fuel channels) nodalisation: (i) 12 axial plane model, with 12 fuel bundles per channel, and extrapolation distance applied axially. Te burnups from the station, which has 13 fuel bundles per channel, are averaged between adjacent bundles. In this model, the liquid zone compartments are centred on an axial plane, while the centre row of adjusters is distributed between the two centre axial planes. (ii) 13 axial plane model, with 13 fuel bundles per channel, with the lengths of the two end planes being equal to a ½ bundle plus an extrapolation length (while the actual length of the core is exactly 12 bundle lengths, the fux does not go exactly to zero at the ends of the core, with the extrapolated length being approximately 12 cm longer than the actual length (see page 22 of chapter 2 of [5]). For the 13 plane model, this extrapolation is applied by lengthening the nodes at either end). In this model, the centre row of adjusters is centred on an axial plane, while the liquid zone compartments are distributed between two axial planes. (iii) 26 axial plane model, with 13 fuel bundles per channel, with the length of the two end planes consisting solely of the extrapolation length used in the 13 plane model. Both the liquid zone compartments and the adjusters are distributed between two axial planes, but as the length of each axial plane is half of its value in the other cases, this is equivalent to one axial plane in the other cases.
Te reference case uses the 13 axial plane model. It should be noted that, while the modelling of the centre row of adjusters is signifcantly afected by the axial nodalisation, the efect on the other two rows is less as they do not fall neatly along a bundle or half-bundle position axially (they are ofset 80 cm from the centre row on either side or 1.62 bundle lengths). For these adjusters, the 12 axial plane model has these adjusters more heavily weighted towards one axial plane. Tus, the efect of the nodalisation of these adjusters is expected to be similar to the liquid zone compartments, except with a smaller magnitude.
2.6.5. Adjuster Irradiation/Ageing. As the adjuster rods are normally in-core, they are continuously exposed to the full in-core fux and are thus irradiated over their lifetime. Te neutron-absorbing materials in the rods will gradually be transmuted to less-strongly absorbing isotopes, over time decreasing the macroscopic absorption cross-section of the rods. Te efect can be signifcant over the lifetime of the adjuster rods. For the CANDU 900 data, the adjuster rods are signifcantly aged and therefore ageing factors must be Science and Technology of Nuclear Installations applied to achieve accurate results for k ef and fux distribution. Te incremental cross-sections for unirradiated adjusters, calculated using Serpent, are perturbed by an ageing factor dependent on the adjuster rod age, adjuster rod type, and cross-section type. Tese aging factors were derived from prior industry work using DRAGON. Perturbation factors are calculated separately for each reaction type and energy group.
For this work, an adjuster rod age of 182,500 efective full-power hours (EFPH) was selected. Adjuster rod ages of 195,000 EFPH and 154,008 EFPH were also selected as sensitivity cases, along with a "new adjusters" case.

Reference Case.
For the reference case, using a nodalisation of 13 axial planes (one bundle per plane, with the end bundles shortened to a ½ bundle plus extrapolation distance), using the initial CANDU 900 depletion snapshot and zone levels, unperturbed nuclear data, unperturbed device incremental cross-sections, and aged adjuster rods with irradiation of 182,500 EFPH, the coupled TRACE/PARCS model, running the modifed codes TRACE_Mac1.1 and PARCS_Mac1.1, it reproduces the overall behaviour of the transient at the CANDU 900 station. Figure 5 shows a comparison of PARCS_Mac1.1 channel powers against reference SORO channel powers obtained from independent simulations by station staf, given identical bundle depletions and equivalent boron moderator poison. Te channel power discrepancies for the vast majority of fuel channels do not exceed 2%, except for edge channels (channels with fewer than four neighbours) and some channels in the vicinity of the liquid zone compartments. In the PARCS model, the edge channels use a different set of cross-sections calculated by homogenising an edge channel lattice model. Table 1 summarises the initial zone levels calculated by the coupled RRS model compared to SORO predictions made by the operators. Most of the predicted zone levels in this simulation are very close to the SORO predictions, except for zones 3, 4, and 11. For zone 3, the discrepancy corresponds to a rather small reactivity diference as the zone is near the top of the core, and the change in fll from 70% to 84% corresponds to near the top of the core, where the fux is lower and thus the zone level change contributes less reactivity. Te other two zones, 4 and 11, are in the centre of the core, and the lower zone level may indicate either excess absorption near the centre of the core or excess refection from the periphery of the core or a combination of the two. Table 2, which summarises the liquid zone level results obtained throughout this work, shows that the centre zone level is particularly sensitive to nuclear data uncertainty. Tough the values in Table 2 are calculated at the end of the shim operation rather than at the start, they show that the standard deviation is roughly 15%, thus the zone level discrepancy is well within two standard deviations when considering nuclear data uncertainty, and thus not signifcant.
Te efect of the thermal-hydraulic channel grouping (40 channels) was evaluated by modelling 480 single channels using the conditions from the coupled model, then feeding the single-channel model conditions back into PARCS. Te efect was found to be insignifcant, with k ef changing by −0.03 mk and the vast majority of channel powers changing by less than 0.2%. Figure 6 presents the liquid zone level trends of the simulated reference case against the station data, with a rolling window used to reduce the noise for individual zone levels. Te two subplots are synchronised based on a known ofset between the initial snapshot used to start the simulation and the station data. Te trends are overall similar between the station data and simulation, but some diferences were noted.
(i) Te average liquid zone levels for the central zones (zones 4 and 11) are up to 15% lower in the simulation compared to the station data, for similar points in the transient. Tis indicates that there is slightly less reactivity in the centre of the core than the periphery in the simulation, compared to the station data, that the liquid zone controllers are compensating for. (ii) Te maximum liquid zone levels after each adjuster pull, as well as after the refuelling run, are slightly less in the simulation than in the station data (by less than 5%). (iii) Te withdrawal of adjuster bank A occurs approximately one day earlier in the simulation. Te withdrawal of adjuster bank B occurs approximately 3 days earlier in the simulation, giving a discrepancy of two days for the interval between adjuster banks A and B. Te remainder of the simulation shows little time discrepancy from the station data, besides the 3-day shift arising from the previous discrepancies. (iv) Te zone level tilt between the two zones in each zone pair is considerably larger in the simulated core than in the real core, particularly for the 3/10 and 4/11 zone pairs. Figure 7 shows the continuation of the simulation, where the power is reduced to 59% FP, showing the simulated reactor and zone powers along with liquid zone compartment levels. After roughly 2.5 hours, the average zone level reaches 10% with all adjusters withdrawn, and the reactor is shut down. Te zone levels behave as expected, with the zone levels in the central zones (4/11) increasing, in particular, to counteract the fux peaking from the removed adjusters. Te zone powers also follow expected trends, where fully drained zones decrease in power while the remaining zones increase in power. As well, towards the end of the simulation, the total power decreases by approximately 2% FP. Tis is due to the RRS response from the liquid zone control algorithm, from the following two efects: (1) With twelve liquid zones fully drained, the bulk controller's efective gain is reduced, as a given negative valve lift will drain only two zones instead of all 14. Tus, a larger power error is required for a given reactivity insertion rate (dictated by the xenon transient being ofset to maintain criticality). (2) Te spatial controller is applying positive valve lift to the 12 drained zones (to attempt to bring their fll levels up to the average fll level) as well as applying positive valve lift to the two central zones (to attempt to bring the zoning power down to average zone power). Te bulk controller must counteract this positive lift by applying a negative lift, which requires a negative power error.
Tis efect would not have been captured in a quasi-static calculation that assumes reactor power is held at the setpoint while controlling for k ef and spatial fux, thus demonstrating the advantage of using a dynamic (transient) approach at key points in a long-term simulation.
Te maximum power tilts immediately prior to shutdown, ignoring zones 4 and 11, is approximately 13%. Te fux tilt increases over the course of the transient as the zone   Science and Technology of Nuclear Installations 13 levels become fully drained or flled and cannot spatially control the fux, combined with the positive feedback efect from xenon. Table 3 gives the efect of applying the adjuster depletion on the reactivity worth of the adjusters. Te reactivity worth of each bank is calculated from the initial fuel depletion snapshot by calculating the reactivity change between the bank-in and bank-out states with all liquid zones flled to 50% and with all prior banks withdrawn. Overall, the reactivity worth of the adjusters is reduced by 10.92% when aged to 182,500 EFPH, with a greater efect on adjuster types 1, 2, and 4. Table 4 summarises the diferences in the timeline of the transient between the station data and the diferent simulations throughout this work. Overall, the length of the operation with zero and two banks out is slightly underpredicted, the length of the operation with three banks out is overpredicted, and the length of the operation with one bank out is greatly underpredicted. Te greater the adjuster age, the shorter each interval becomes, with the three-bank-out interval being the most affected. With unirradiated adjuster rods, the entire operation is 1.22 days longer than for the reference case.

Efect of Adjuster Rod Irradiation.
Te adjuster depletion also infuences the radial fux distribution in the core. Table 2 shows the magnitude of this efect, with the diference in zone levels in zones 4 and 11 between the simulations with irradiated versus unirradiated adjusters at 28.0% towards the end of the transient, with the unirradiated adjusters showing a greater underprediction against the station data when compared to the reference case.
Finally, Figure 8 shows the results of changing the effective adjuster rod age to the end of the simulated transient. Results are shown for fresh adjusters, as well as ages of 154,008 EFPH, 182,500 EFPH, and 195,000 EFPH. Te irradiation of the adjuster rods has a signifcant efect on the outcome of the transient. With fresh adjusters, the reactor survives xenon poison out with a minimum zone level of 35%, achieved with all adjusters out of the core. With an adjuster rod irradiation of 154,008 EFPH, the reactor initially survives xenon poison out with a minimum zone level of 12%. However, the spatial xenon transient continues until reaching a 20% power tilt, which triggers a setback in the simulation. Tis reduces the power setpoint to 55.8% FP, causing the xenon transient to poison out the reactor. With adjuster rod irradiations of 182,500 EFPH (reference case) and 195,000 EFPH, the reactor poisons out without a setback, with the poison out occurring 10 minutes sooner for the 195,000 EFPH case.
Te 154,008 EFPH case demonstrates that the reactor becomes unrecoverable once the average zone level is too low, as the lack of spatial control leads to a fux tilt setback. Tus, the cases that are shut down without a setback upon reaching the 10% average zone level would be expected to inevitably shut down, either from fully draining the liquid zones and becoming subcritical or from a fux tilt setback leading to the reactor becoming subcritical.

Efect of Liquid Zone Reactivity Worth.
For this sensitivity case, a perturbation was applied to the liquid zone compartment incremental cross-sections, decreasing their values by 10%. Tis resulted in decreasing their reactivity worth by approximately 8.6%, from 5.76 mk to 5.26 mk. Tis change only has a small impact on the poison-out trajectory compared to the reference case, with the adjuster movement timings essentially unchanged, and the poison out occurring roughly 10 minutes earlier. Te rest of the shim transient is slightly shorter. Tis efect is in line with expectations, as the reactivity gained from draining the liquid zones from 50% (initial state) to 25% (the threshold for adjuster movements) and then to 10% (the threshold for poison out) is slightly decreased. Tese results are quantifed in Tables 2 and 4.

Efect of Axial Nodalisation.
For this sensitivity case, the axial nodalisation is altered. Te efect on the fnal phase of the transient is shown in Figure 9. When a 26-plane nodalisation is used instead of the 13-plane nodalisation, the simulation predicts that the minimum zone level reaches 10%, then poisons out when the fux tilt triggers a setback. When a 12plane nodalisation is used, the simulation predicts that the reactor survives xenon poison out with a minimum zone level of approximately 18%. Tere is also a small efect on the duration of the shim operation; the efect is shown in Table 4. Table 5 shows how the nodalisation afects the predicted reactivity worth of each set of reactivity devices. Note that the 13-plane model predicts a higher liquid zone controller worth than the 12-plane model. Conversely, the 12-plane model predicts higher adjuster worths than the 13-plane model, with the efect coming from banks D, F, G, and H, which are the centre row banks. Tis is consistent with expectations, where a model that distributes the device incrementals over a larger volume will predict a larger reactivity worth due to the weaker spatial self-shielding. Tis also explains the results of the previous fgures. For Figure 9, the 12-plane model survives poison out as the adjuster reactivity worths are signifcantly greater. Te results of the 26plane model fall intermediately between the 12-plane and 13-plane models.  It should also be noted that the efect of the nodalisation on the liquid zone compartments is of a similar magnitude to the liquid zone sensitivity study. Terefore, the efect of the nodalisation of the liquid zone compartments, independent of the adjusters, will be similar to the results presented for the efect of the liquid zone reactivity worth. Table 4. Te station data fell within roughly two standard deviations of the average simulation results, except for the onebank-out interval. However, the durations of the intervals are highly correlated (>96%), while the diferences in station data for diferent intervals have opposite signs, thus nuclear data uncertainty alone cannot explain all of the diferences between simulations and station data for shim durations. Figure 10 plots the efect of nuclear data perturbation on the extended portion of the transient, with each sample, plotted as a coloured line, with time zero set as the time at which the reactor power setpoint is set to 59.0% FP. Each of the 60 samples follows a similar trend, though the timing varies. Te results of the samples are as follows:

Efect of Nuclear Data Uncertainty. Te comparison of the shim transient timings is included in
(i) 44 samples (73.3%) reach the 10% average zone level threshold, triggering the simulation to shut down the reactor. Tese can be considered to poison out from bulk xenon only, even without the efect of spatial xenon. (ii) 15 samples (25.0%) do not reach the 10% average zone level threshold at frst, but a reactor setback on fux tilt is subsequently triggered. Tese can be considered to not poison out directly from bulk xenon, but poison out due to spatial xenon increasing the fux tilt to an unacceptable level.
(iii) 1 sample (3.3%) does not reach the 10% average zone level threshold and does not experience a setback.
Due to the negative valve lift from the bulk controller along with spatial control being phased out below the 10% zone level, it is possible for the zoning power in drained zones to exceed the zoning power in spatially controlled zones while the zones remain drained. Tis was observed for zones 5 and 12. Tis continues until either a fux tilt setback occurs, or a crossover point is reached where the liquid zone compartment starts to refll, providing positive feedback as spatial control is phased in, at which point the liquid zone levels rapidly change to a new equilibrium.
Te sensitivity to liquid zone controller gains was evaluated, with a focus on the marginal cases (i.e., those closest to the threshold of whether a setback occurs) and it was found that realistic perturbations to the controller gains did not signifcantly alter the outcome of the cases evaluated in this work.
Tis indicates that, when only nuclear data uncertainty is considered, the simulation predicts a poison out with an approximately 73% confdence, or 0.62σ, if the setback is not considered. However, if setbacks within the simulation duration are also considered, then the confdence rises to 98%, or 2.1σ. In addition, the nuclear data uncertainty does not include uncertainty in the incremental cross-sections of the reactivity devices, as perturbed incremental crosssections were not available in the current methodology. However, the sensitivity of the outcome to the adjuster and liquid zone incremental cross-sections was determined in the prior subsections.
For comparison, 24 cases were run with unirradiated adjusters modelled and showed that none of these cases poison out with unirradiated adjusters and that the nuclear data uncertainty has far less of an impact than the adjuster rod irradiation efect. Figure 11 plots the efect of simulating the transient with diferent initial fuel burnup distributions. Unlike the simulations using the original CANDU 900 snapshot, these cases do not include any refuelling runs. Terefore, a case was also run using the  original CANDU 900 snapshot, but without any refuelling runs, to be the reference against which other burnup distributions are compared. All 14 snapshots lead to a poison out, with the reference case falling roughly in the middle of the 14 snapshots. Some reactor setbacks were observed earlier in the transient than in the above sections, as the generated burnup distributions were less optimal than the CANDU 900 snapshot.

Efect of Burnup Distribution.
When the transient was simulated with both nuclear data uncertainty and diferent initial fuel distributions, with each nuclear data sample having its own corresponding core snapshot, the distribution of shim operation trajectories was found to be similar to those observed for the CANDU 900 samples in Figure 10, except that earlier setbacks were observed as in Figure 11. All 20 of the perturbed cases led to a poison out. Tese results are quantifed in Tables 2 and 4.
For one of the nuclear data samples with its corresponding snapshot, the simulated reactor tripped on regional overpower while withdrawing adjuster bank A. Tis was due to the snapshot having a single channel at abnormally high power (>7.1 MW) which increased CPPF to an abnormal level (>1.15), reducing the trip margin (note: the reactor would never be operated in such a state in reality). For this sample, the snapshot was replaced with a diferent snapshot from the same simulated core follow.
A minority of the perturbed cases also experienced setbacks during the adjuster movements. In this case, the simulation was set up to restore the reactor power as soon as the setback condition cleared.

Discussion
Tis study demonstrates that, by coupling TRACE_Mac1.1 and PARCS_Mac1.1 along with an RRS model, shim operation can be simulated, and sensitive parameters can be evaluated. Tese sensitivities are most interesting to examine for the extended portion of the simulation, where the reactor power is reduced to 59% FP. Te importance of these sensitivities in this phase of the simulation, as identifed from this study, are presented by listing the sensitivities from most to least signifcant as follows, quantifed using the average zone level at 105 minutes.
(1) Adjuster rod irradiation (−19.5%, reference versus unirradiated) (2) Axial nodalisation (+6.9%, 12 planes versus 13 planes, due to efect on adjusters) (3) Nuclear data uncertainty (σ � 2.6%) (4) Burnup distribution (σ �1.0%) (5) Liquid zone compartment reactivity worth (−0.6% for 10% incremental change) Te adjuster rod irradiation was the greatest factor in defning the end of the transient, especially when comparing irradiated and unirradiated adjusters. When the transient is simulated with unirradiated adjusters, they had sufcient reactivity depth to maintain criticality and prevent the core from being poisoned out by xenon, with a reasonable margin. When the transient is simulated with irradiated adjusters, the core is poisoned out. Te reactivity worth of the fve banks of adjusters was reduced by 0.96 mk, or 11.6%, for the irradiated case compared to the unirradiated case. Varying the irradiation around the selected reference value had a proportionally smaller efect on the transient.
Te axial nodalisation had a signifcant efect on the end of the transient, as the centre row of adjusters, which dominate the adjuster contribution to the end of the transient, were signifcantly afected by the doubling of their efective volume. Te calculated worth of the fve banks of adjusters was 0.38 mk, or 5.3%, higher for the 12-plane model compared to the 13-plane model-nearly half of the adjuster irradiation efect! Since the adjuster worth, particularly for the centre row, dominates the outcome of the end of the transient, the 13-plane model should be more accurate in this regard, as the adjuster efective volume is closer to the original homogenisation volume. As well, it should be considered that the model used in this study also doubles the efective volume of the adjusters in the horizontal (x) direction, making this an area for future investigation.
Te nuclear data uncertainty was also found to have a signifcant efect. When the average zone level was taken 105 minutes after the reactor power setpoint is set to 59% FP, Table 2 showed that the efect of 28,500 EFPH of adjuster depletion was 2.8% zone level, while one standard deviation of the nuclear data uncertainty was 2.6% zone level, thus one standard deviation of nuclear data uncertainty is equivalent to approximately 28500 EFPH + 2.6% 2.8% � 26500 EFPH.
Tus, one standard deviation of nuclear data uncertainty, not including device incremental uncertainty, is roughly 15% of the total adjuster depletion efect. Te infuence of the fuel distribution in the core is equal to approximately 40% of the infuence of the nuclear data uncertainty, based on a standard deviation of 1.0% zone level, compared to 2.6% zone level for the nuclear data uncertainty. Te other signifcant fgure of merit was the timing of events from the initial snapshot up to the end of threebank-out operations at a 25% average zone level. Te reference case underpredicted the total length of the shim operation by 2.75 days (10%), primarily due to an underprediction of the one-bank-out duration by 1.73 days (32%), as well as an underprediction of the zero-bank-out duration by 1.14 days (9%). Te importance of the sensitivities, from most important to least important, is as follows: (1) Nuclear data uncertainty (σ � 1.31 days) ( Te discrepancy thus cannot be entirely explained by the uncertainties studied in this work, especially for the onebank-out duration. Te two-bank-out and three-bank-out durations are within two standard deviations (from nuclear data uncertainty) of the station data, but in opposite directions. Te discrepancies must be explained by phenomena that were not evaluated in this work, with potential phenomena as follows: (i) Nodalisation of the core physics model, outside of the evaluated axial nodalisation meshes. (ii) Diferences in depletions of diferent adjusters of the same type (up to a 5% variation in fux, and thus depletion, was identifed to exist, but not applied to the model). Tis efect is expected to be small (equivalent to a 5% variation in adjuster depletion, but for individual banks). (iii) Errors afecting the rate that fuel depletion changes the core reactivity, either from approximations in the lattice physics and core physics model, or parameters in the lattice physics model (e.g., fuel density). (iv) Errors afecting the power coefcient of reactivity, void coefcient of reactivity, or xenon reactivity, such as the modelling of fuel temperature or coolant voiding. A negative change to the power coefcient would add positive reactivity when the reactor power is decreased, lengthening the shim operation (results move closer to the station data). A negative change to the void coefcient would have a similar efect as the power coefcient, but only at high power. (v) Changes to thermal-hydraulic setpoints over the course of the transient that are neither captured in the model nor provided with the station data. No known changes to setpoints were identifed for this transient.
(vi) Te presence of moderator poisons and changes in their concentration during the shim operation. While the equivalent moderator poison was estimated periodically over the course of the real-world transient using SORO, no actual measurements were taken. Te simulation in this work used the equivalent moderator poison values, which went to zero prior to the frst adjuster bank being withdrawn. If these estimates were inaccurate, particularly for the initial state as well as for the time at which adjuster bank A is withdrawn, then this would create a discrepancy in the shim duration. More specifcally, if moderator poison was present when bank A was withdrawn, and subsequently removed, it would extend the duration of the onebank-out shim operation.
At the time of this writing, the most likely explanations for the discrepancies in shim duration are a combination of nuclear data uncertainty and errors in the estimated moderator poison concentrations.
Te discrepancy in radial power (as indicated by central zone compartment levels) was also measured. Te importance of the sensitivities, from most important to least important, is as follows: (1) Adjuster rod irradiation (+28%) (2) Nuclear data uncertainty (σ �15%) (3) Axial nodalisation (−5.0%, 12 planes versus 13 planes) (4) Liquid zone compartment reactivity worth (+5.1% for 10% incremental) Te LZC reactivity worth is ranked lower as the efect is predominantly due to an increase in the amount of water required for the same amount of spatial control. Tus, zone levels are not a good indication of radial power balance when the LZC worth changes signifcantly. Since the LZC worth also changes for the axial nodalisation sensitivity case, the LZC worth efect partially ofsets the adjuster worth efect. Note that the diference between the reference simulation and station data was −18.6%. Terefore, while other errors may be present, the discrepancy in the station data may be explained by the nuclear data uncertainty.

Conclusions
Tis work was able to demonstrate that a shim operation could be simulated using a coupled model of TRACE_-Mac1.1, PARCS_Mac1.1, and an ECI-coupled RRS model, and model depletion and adjuster movements comparably to real-world data. Te coupled model is able to simulate the efect of depletion, xenon-135, thermal-hydraulics, and reactor regulating system response on a CANDU reactor.
When the transient was extended to evaluate a possible continuation, with the power reduced to 59% FP, the reference case predicted that xenon poison out would occur 2.5 hours after the power reduction. Te model showed that Science and Technology of Nuclear Installations adjuster depletion had a major impact on the outcome of the simulation, with no poison out occurring with unaged adjusters. When including cases that poison out due to fux tilt setback, 59/60 nuclear data uncertainty cases led to poison out. Te efect that the burnup distribution had on the simulated transient is less signifcant than nuclear data uncertainty, with only roughly 40% as great of an efect on the fnal stage of the transient.
Overall, out of all parameters evaluated, the modelling of the adjusters had the greatest impact on the results and fgures of merit for the transient. Not only did the adjuster depletion have a signifcant impact on the results, but the nodalisation of the adjusters also played a signifcant role. Tis work demonstrated that nodalising the core such that an appropriate equivalent volume is achieved for the adjusters is important and that a larger equivalent volume increases the estimated adjuster's worth. While the modelling of the LZCs is similarly afected, the LZCs had far less impact on the results of the transient.
Te nuclear data uncertainty was shown to be able to signifcantly perturb the results, both in terms of the duration of the shim operation, the radial fux distribution, and the outcome at the end of the extended simulation. Terefore, the efect of nuclear data uncertainties cannot be ruled out as a contributing factor to the observations made in this work or for other similar studies.

Future Work
While this work was able to generally model the real-world shim operation and identify several important sensitivities in the modelling of a shim transient, some unaccounted-for discrepancies were still identifed, and some sensitivities were identifed as candidates for future analysis. Te results of this work present several areas for future investigation.
Tis work identifed nodalisation of the core physics model as a signifcant contributing factor, due to the efect on the modelling of the adjuster rods. However, this work only tested a few nodalisations in the axial direction, and these nodalisations do not completely solve the issue of adjuster incremental cross-sections being diluted over a larger volume. In particular, in the horizontal (x-axis) direction, the incremental cross-sections are distributed between two adjacent channels.
A horizontal subdivision of the model was not pursued in this study due to the extra efort required in PARCS when compared to axial renodalisation. Currently, there is no method in PARCS to specify independent meshes for fuel assemblies/channels and for the macroscopic cross-sections, the latter of which is needed for reactivity device incremental cross-sections. Tis contrasts with some CANDU-specifc core physics codes, such as RFSP. Terefore, PARCS would treat these subdivided channels as independent fuel assemblies/channels, requiring extra efort to map powers and fuxes to TRACE and particularly to map powers and fuxes to the ECI through signal variables as is done in the current version of PARCS_Mac1.1 and TRACE_Mac1.1. Te execution time of the simulation would also be greatly increased.
However, due to the identifed importance of nodalisation, it is recommended that future studies further quantify the efect of the nodalisation on the adjuster rod reactivity worth. Te current model was able to quantify the efect of changing adjuster worths on the simulation of a shim transient, thus a study that estimates the change in adjuster reactivity worth due to other modelling parameters (such as nodalisation in the x-axis) may be used to estimate the efect these parameters would have on the shim transient. Future work may also investigate other methods for correcting nodalisation efects, such as extending the rod cusping correction methodology already present in PARCS or applying similar methodologies to calculate correction factors. Monte Carlo models of a full core or mini-core may be created to calculate reactivity device worths that are unafected by nodalisation, to compare against PARCS models.
Secondly, there were some large discrepancies in the duration of diferent phases of the shim operation, particularly for the one-bank-out phase. One signifcant "unknown" that was identifed was the state of moderator poison in the core. While an "equivalent moderator poison" was calculated at the time of the operation using SORO, no true moderator poison measurements were made. Any error in the moderator poison estimation, either at the initial state or upon withdrawing an adjuster bank, would impact the timing of each phase of the shim operation, as the coupled simulation assumed these estimates to be accurate. Another possible area for future investigation is the thermalhydraulic feedback.
Several recommendations can thus be made for studies aiming to reproduce these results or further expand on this work. Te frst recommendation is to obtain data from other shim operations besides the one used in this work, ideally using a case where moderator poison can be defnitively ruled out as a contributing factor. Tese may include shim operations from either 900 MW class or 600 MW class CANDU reactors. Te second recommendation would be to obtain thermal-hydraulic data, as the dataset provided for this work did not include any thermal-hydraulic data against which to compare the simulated thermal-hydraulic model. Another related recommendation is to obtain real quantities for RRS parameters, such as liquid zone control gains, rather than estimating these parameters, to more accurately model RRS feedback, eliminating a source of uncertainty.
As an example, consider the hypothesis that the large underestimation of the duration of the one-bank-out phase of the operation was due to a discrepancy between the estimated moderator poison concentration (through equivalent boron used in SORO) and the actual moderator poison concentration. If similar discrepancies are not found when simulating diferent shim operations, it suggests that the hypothesis may still hold, and alternative hypotheses, such as a systematic modelling error, may be ruled out, or at least considered less likely. If the discrepancy does recur, then it suggests that either the hypothesis may be ruled out (and alternative hypotheses that may be common to both events considered), that a common discrepancy in estimating the moderator poison occurred in the events (the likelihood of which can be estimated by comparing the timelines of both events, particularly the use of moderator poison and the estimates of equivalent moderator poison), or that the discrepancy was co-incidental, with two diferent causes for the two events.
Another area to be investigated is the use of nuclear data and the estimation of the uncertainty of nuclear data. Since nuclear data uncertainty was found to have a signifcant impact, the use of newer nuclear datasets, such as ENDF/B-VIII, may impact both the best estimate as well as uncertainty results. Additionally, the methodology in this study was limited as it was unable to apply the stochastic nuclear data perturbations (from Sampler) to the reactivity device incremental cross-section calculations (using Serpent, which is outside of the SCALE framework supported by Sampler). An improved uncertainty analysis would integrate the reactivity device calculations into the nuclear data uncertainty analysis, such as by migrating the supercell calculation to CSAS-Shift in the upcoming SCALE 6.3 release [16]. Finally, the methodology described in reference [7] could be applied to the uncertainty analysis, performing a core follow from a much earlier state to the point of the shim event using the reactor's operational history, to achieve core depletion profles consistent with their respective nuclear data perturbations.
Te methodology used in this work may also be applied to the analysis of design changes, such as the adoption of advanced fuel cycles, which had been previously studied for DUPIC fuel [2,3] as well as thorium-based fuel [4], with PARCS_Mac serving as a viable alternative to RFSP and DONJON while allowing coupling to TRACE_Mac for thermal-hydraulic feedback as well as taking advantage of the uncertainty propagation capabilities present in SCALE using Sampler.

Data Availability
Te TRACE and PARCS model data used to support the fndings of this study have not been made available due to the confdentiality of the models. Te detailed data used to support this work is also unavailable as it was provided confdentially, with approval to present the data included in this article.