diff --git a/romancal/assign_wcs/assign_wcs_step.py b/romancal/assign_wcs/assign_wcs_step.py
index 44147e6cb..553fbeb2b 100644
--- a/romancal/assign_wcs/assign_wcs_step.py
+++ b/romancal/assign_wcs/assign_wcs_step.py
@@ -95,15 +95,25 @@ def load_wcs(input_model, reference_files=None):
         axes_names=("v2", "v3"),
         unit=(u.arcsec, u.arcsec),
     )
+    v2v3vacorr = cf.Frame2D(name='v2v3vacorr', axes_order=(0, 1),
+                            axes_names=('v2', 'v3'), unit=(u.arcsec, u.arcsec))
     world = cf.CelestialFrame(reference_frame=coord.ICRS(), name="world")
 
     # Transforms between frames
     distortion = wfi_distortion(output_model, reference_files)
     tel2sky = pointing.v23tosky(output_model)
 
+    # Compute differential velocity aberration (DVA) correction:
+    va_corr = pointing.dva_corr_model(
+        va_scale=input_model.meta.velocity_aberration.scale_factor,
+        v2_ref=input_model.meta.wcsinfo.v2_ref,
+        v3_ref=input_model.meta.wcsinfo.v3_ref
+    )
+
     pipeline = [
         Step(detector, distortion),
-        Step(v2v3, tel2sky),
+        Step(v2v3, va_corr),
+        Step(v2v3vacorr, tel2sky),
         Step(world, None),
     ]
     wcs = WCS(pipeline)
diff --git a/romancal/assign_wcs/pointing.py b/romancal/assign_wcs/pointing.py
index ce24a354d..edd967d06 100644
--- a/romancal/assign_wcs/pointing.py
+++ b/romancal/assign_wcs/pointing.py
@@ -1,5 +1,5 @@
 import numpy as np
-from astropy.modeling.models import RotationSequence3D, Scale
+from astropy.modeling.models import Identity, RotationSequence3D, Scale, Shift
 from gwcs.geometry import CartesianToSpherical, SphericalToCartesian
 
 
@@ -45,3 +45,58 @@ def v23tosky(input_model, wrap_v2_at=180, wrap_lon_at=360):
     )
     model.name = "v23tosky"
     return model
+
+
+def dva_corr_model(va_scale, v2_ref, v3_ref):
+    """
+    Create transformation that accounts for differential velocity aberration
+    (scale).
+
+    Parameters
+    ----------
+    va_scale : float, None
+        Ratio of the apparent plate scale to the true plate scale. When
+        ``va_scale`` is `None`, it is assumed to be identical to ``1`` and
+        an ``astropy.modeling.models.Identity`` model will be returned.
+
+    v2_ref : float, None
+        Telescope ``v2`` coordinate of the reference point in ``arcsec``. When
+        ``v2_ref`` is `None`, it is assumed to be identical to ``0``.
+
+    v3_ref : float, None
+        Telescope ``v3`` coordinate of the reference point in ``arcsec``. When
+        ``v3_ref`` is `None`, it is assumed to be identical to ``0``.
+
+    Returns
+    -------
+    va_corr : astropy.modeling.CompoundModel, astropy.modeling.models.Identity
+        A 2D compound model that corrects DVA. If ``va_scale`` is `None` or 1
+        then `astropy.modeling.models.Identity` will be returned.
+
+    """
+    if va_scale is None or va_scale == 1:
+        return Identity(2)
+
+    if va_scale <= 0:
+        raise ValueError("'Velocity aberration scale must be a positive number.")
+
+    va_corr = Scale(va_scale, name='dva_scale_v2') & Scale(va_scale, name='dva_scale_v3')
+
+    if v2_ref is None:
+        v2_ref = 0
+
+    if v3_ref is None:
+        v3_ref = 0
+
+    if v2_ref == 0 and v3_ref == 0:
+        return va_corr
+
+    # NOTE: it is assumed that v2, v3 angles and va scale are small enough
+    # so that for expected scale factors the issue of angle wrapping
+    # (180 degrees) can be neglected.
+    v2_shift = (1 - va_scale) * v2_ref
+    v3_shift = (1 - va_scale) * v3_ref
+
+    va_corr |= Shift(v2_shift, name='dva_v2_shift') & Shift(v3_shift, name='dva_v3_shift')
+    va_corr.name = 'DVA_Correction'
+    return va_corr