@MaximeJ thanks for pointing me towards those methods - I think I was able to get it to work. Python code for it is below.
@markrutten I also tried your method, and though it took longer to run than just propagating the covariance, it looked to yield more realistic results. Specifically the fact that covariance was largest in the In-Track direction with small components in the Radial and Cross-Track. I think the reason for this discrepancy is due to the equation used to propagate covariance. Specifically, the addition of the Process-Noise matrix tends to have an increasing effect as the difference between t and t0 increase.
I am using the same process noise matrix as described in Kalman Filter Covariance Propagation in RIC. Using the coefficients described in that paper (1.4e-9 radial, 2.4e-11 intrack, 5e-9 crosstrack) yield a Process-Noise matrix that has the largest growth in the radial direction which is not realistic. These coefficients can of course be tweaked to match Mark’s solution results, however it seems artificial and more like fitting the results to the data as opposed to the other way around.
I found a newer paper on covariance realism which uses similar equations, but also has a process-noise transition matrix that is added into the equation - is there a way to get that process-noise transition matrix out of Orekit? Paper is here https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/20160010501.pdf