First, we should clarify some points about the question,
ITRF and WGS84 are Geodetic Reference Frames
(a.k.a. "horizontal datums"), Which are calculated by different organizations, but are in general compatible.
The realization strategy of both systems differs. One key difference is that ITRF is recalculated from time to time, giving its "realizations" named ITRF, i.e.: ITRF97, ITFR2000, etc. Meanwhile, WGS84 is not usually renamed nor fully recalculated, instead it is periodically aligned to ITRF to ensure consistency with Geodetic mesurement. (see: this answer and its references in the comments).
It is probably not your use case, but there may be some differences in positions due to different epoch of materialization (the point in time for the coordinates which are published) of the reference frame, because tectonic dynamics must be considered when high precision is needed.
EGM96 is a geoid model, which is a model of the behavior of the gravity on the earth surface. As a derived product, equipotential surfaces of the earth can be computed from EGM96. More specifically, the equipotential surface of the earth which best fits to the mean level of the seas, is computed for its usage as a vertical reference surface, and orthometric heights are computed related to this surface. (see this question, and this other -disclaimer, accepted answer is mine-).
With this in mind, your needed conversion is probably from orthometric heights defined with respect to the EGM96 model to ellipsoidal heights with respect to the WGS84/ITRF97 reference frame(s) and its associated ellipsoid(s).
EGM96 is so common, that many libraries and online calculators have its data pre-loaded for computations. This answer should help for python
Edit: Epsg code and example
The EPSG code for EGM96 heights is EPSG:5773 (link here). If you look at the details, the definition includes a grid which contains the relevant conversion:
By example, the
cct tool, included on the PROJ tools can handle this (if you have the datum grids installed).
cct +ellps=WGS84 +proj=vgridshift +grids=egm96_15.gtx <<EOF
< coordinates here >
Here i have added a gist with an ipynb as example it uses the
cct commandline tool behind the scenes from the python notebook, you can open it on colab and make tests.
On pyproj (the python bindings for PROJ) the Transformer class should do the trick: See the docs
For information and data on Geoid models, a full reference is given on this link
For information on the geodetic reference frames, ITRF publications are usually a good and detailed reference, see the ITRF97 documentation