Properly account for energy deposition

[cjosey]

Current energy deposition in OpenDeplete is simply the "Total energy less neutrino energy" value in ENDF7.1 times the fission rate. The exact contents of this value include

• Kinetic energy of fragments after prompt neutron emission.
• Prompt / delayed neutron kinetic energies
• Total energy released prompt/delayed from gamma and beta

As such, if only fission occurred, the time-integral energy deposition would otherwise be exactly correct to the limit of the quality of the library.

However, there are three issues:

1. It assumes all fission energy deposition is instantaneous.
2. This does not include (n, *) non-fission reactions nor the decay heat of the daughters
3. It ignores the spectrum dependence.

If point 1 were only delayed neutrons, it would be of limited concern. The error for a 1 hour BOS time step under such circumstances is 4e-7. However, delayed gamma production is much slower and can introduce as large as a 3.4% error.

In point 2, the reactions themselves are trivial. The decay poses a much larger issue. With the current value of Q for fission, decay heat of fission fragments is already counted. To get the right value, I would need to know how much of each nuclide is fission fragment or irradiation. With the numerical integrators as-configured, this would require doubling the dimensions of the matrix exponential to track that and a huge increase in code complexity.

In point 3, ENDF-6 pages 46-47 indicate that these values are indeed spectrum-dependent, even if trivially so. It is to my understanding that OpenMC supports tallying these functions, but I do not know where to look.

The method below is my first attempt at a correct method. It's relatively simple, and, while it doesn't get the spatial distribution right (need neutron scatter energy deposition, photon transport or low order smearing), I think it does get the integral value right.

1. Tally EFR(E), ENP(E) and EGP(E), the three prompt energy-dependent energy deposition functions. Sum them up for all nuclides/cells and call it A.
2. Tally (n, *) reactions. Multiply these by the Q values of each reaction. Sum for all nuclides/cells, call it B.
3. Compute the Q_recover * decay constant * N, for all nuclides/cells, call it C.

The power normalization factor is then:

P_remain = P_true(t) - C # (if < 0, warn user that their power is nonphysical)
Factor = P_remain / (A + B)

The one remaining issue is "what is Q_recover" in point 3 above. My understanding is that the decay heat value we're using is the sum of the Egamma/Ebeta, etc. values. This is great as it implicitly ignores non-recoverable particles. It does ignore recoil energy as well. For alpha decay, this error is on the order of about a percent, so it should be considered.

[paulromano]

As you alluded to, OpenMC is capable of tallying incident energy-dependent fission energy release. The fission-q-prompt score would give you what you want for your A quantity. That being said, I don't know if the quality of the decay data is such that you can trust it to correctly capture the delayed energy release from gammas (or otherwise) and I believe that is one of the main reasons people resort to using the regular "recoverable" energy that includes delayed energy. @smharper looked into this pretty deeply and probably has a good opinion on the matter too.

[cjosey]

That is quite odd. I understand perhaps why it was done this way (fission heat is a much easier experiment), but then how do other depletion codes handle non-fission decay correctly? Are they seriously tracking the two populations separately?

For U235, the uncertainty for the prompt data is around a factor of 7 larger than for the sum (1 MeV vs 150 keV), which is a huge difference. However, it may be possible to get around that completely. We could still use the sum, except we form a decay chain with the yields and compute the total energy deposition of the daughters. We then subtract this from the sum to get an effective prompt value.

The advantage of this is that while the uncertainty will be technically higher, the decay chains covary in such a way that the integral error cancels. However, my suspicion is that this is what is done during the experiment anyways, and would be wasted effort.

[paulromano]

Historically, codes just assumed some value of total energy release originating from capture reactions -- a common value was 6.1 MeV/fission -- that would be added to the recoverable fission energy. CASMO-5 makes a nice improvement where they use the capture rates explicitly to account to get a better value of of the energy release from capture and, notably, the Q values that they use for capture reactions sometimes include the energy from subsequent decay of products (see Table 1 from the linked paper). Their approach probably involves some kind of cutoff on when to include decay energy, but the paper doesn't discuss how those choices were made.

[smharper]

I agree with Paul that the yield data is probably not quite good enough to use on it's own. It's probably best to use the MT=458 data as much as we can. It comes from a pretty complex analysis of many models and experiments, of which the independent yield tables are only one part. There's also a lot of attention to experimental and modeling uncertainties (and covariances) that you won't be able to recover from the ENDF data.

I would advocate for one of these two approaches for fission energy:

• If you don't care about the time dependence of fission energy deposition, just use the fission-Q-recoverable score and call it a day.
• If you do care about the time dependence then use EFR, EGP, ENP, and END (fission fragment, prompt gamma, prompt neutron, and delayed neutron) data from MT=458 and assume they are all deposited promptly (END is obviously not prompt, but I think it's prompt enough for depletion purposes) (we would have to add an OpenMC tally for END but that shouldn't be hard). Compute EB, EGD (delayed beta and gamma) as a function of time from yield tables and scale the integral over time to match MT=458.

In case you guys want to dive into more of the details, I'll share what I know about how MT 458 is computed.

tl;dr is that EB and EGD have been fudged away from the yield tables to be consistent with more precise datasets.

Probably the most important dataset here are the chain yields. Remember that beta decay is the dominant decay mode of fission products and that A (the mass number) doesn't change with beta decay. So experiments that compute chain yields lump together fission products that share the same A into a chain, and then they compare the relative yields of these chains. This experiment is easier than others because (after the delayed neutron precursors have died down) the yield of each chain doesn't vary with time after the fission event. Consequently, chain yields are known to a high degree of certainty.

ET (the total fission Q-value) is computed analytically from the mass defect of the chain yields. (In other words, look at the mass of the stable isotopes which terminate each A chain, then compute Q from the mass difference between these isotopes and the fission target.) This makes ET the most precisely known Q-value.

We also have sets of independent yields where we took the chain yields and applied models for yields of Z as a function of A. In theory we could use them to compute all the components of ET, but in practice they are not precise enough.

ENP and END (the energy carried by prompt and delayed neutrons) are computed directly from our data for nu values and spectra (e.g. MF=5, MT=18 and MF=1, MT=456 for prompt neutrons). These values are considered to be known extremely well. Yield tables are practically ignored for computing END.

The ratios between EB, EGD, and ENU (the energy carried by delayed betas, gammas, and neutrinos) are computed from the yield tables.

If ET, ENP, and END are well known, and the ratios between EB, EGD, and ENU are known, then the equation you need for consistency is

EFR + EGP + (EB + EGD + ENU) = ET - ENP - END

I think the basis for our modern 458 data is 1981 EPRI study that used that equation plus the independent yield tables, plus many (relatively uncertain) experiments of EFR, plus some experiments of EGP, plus the ANS decay heat standard in a big least-squares datafit. Since then, most of the terms have been updated individually, but they presumably don't deviate in a big way from this EPRI analysis.

[cjosey]

This is rather unfortunate. So, our options are:

1. Use prompt data only, compute delayed energy from decay chain.
   Suffer an integral accuracy penalty.
2. Use non-prompt data, track fission and non-fission products separately.
   Massive increase in code complexity, flops and memory requirements. Suffer temporal accuracy penalty but no fission integral accuracy penalty.
3. Generate our own data using some energy balance to match MT=458.
   You could make a whole PhD out of this.
4. Use non-prompt data, and use an approximation to integral energy per non-fission reaction.
   Might as well use CASMO. It's faster.

There is a possible shortcut for 3 in that we tally (ET - ENU) anyways, but then calculate a chain-dependent value given by:

With an initial N given by the independent yields. This can be approximately computed (as in t must not be infinity) using a relatively convoluted 7000x7000 matrix exponential. I subtract this value, treat that as prompt fission, and then use the decay chain as in 1. The result is that MT=458 is conserved at the limit of t to infinity. we get time dependence, and we don't have to modify the integrators. The downside is that if we add spectrum-dependent yields, this would have to be computed at every time step for every cell. Thoughts?

Edit:
On the note of spectrum-dependent yields, I am aware the data is roughly 50% nonsense. In addition, at worst, we'd be linearly interpolating and since the above equation is linear on initial condition, we can precompute the solution for all yield(E) values and interpolate on those.

[cjosey]

So, I implemented that matrix exponential. Turns out it's "merely" an n+1 x n+1 matrix where n is the number of nuclides.

According to it, a fission of U235 has a delayed emission of 11.750 MeV in the first day, 12.127 MeV in the first week, 12.598 MeV in the first year, 12.739 MeV in the first hundred years, and 13.136 MeV asymptotically. The non-prompt MT=458 value is 12.842 MeV. As such, it would be incorrect to use the non-prompt MT=458 data directly, as it would introduce a 0.3 MeV (+/- energy dependence) error with the current decay chain. Which, honestly, is still probably smaller than our uncertainty.

I also found an issue in the decay chain. I was trying to replicate the CASMO data when I noted that Cm244 decays to Cm244 via spontaneous fission. I will have to fix that eventually.

Regardless, I think I'll take a stab at implementing this. It should be pretty trivial.