A computer implemented method for restoring a 3D model of the subsurface geology of the Earth from a present day geometry measured, using geological or seismic tomography of subsurface structures, at a present time to a predicted past geometry at a past restoration time, the method comprising: receiving the 3D model of the present day measured geometry comprising a network of faults, wherein a fault is a discontinuity that divides fault blocks that slide in opposite directions tangential to the surface of the fault as time approaches a modeled time; selecting a past restoration time τ that is prior to the present time and after a time when an oldest horizon surface in the 3D model was originally deposited; the method characterized by further comprising: dividing the network of faults into a subset of τ-active faults and a subset of τ -inactive faults, wherein a τ-active fault is a fault that is active at the past restoration time τ and a τ-inactive fault is a fault that is inactive at the past restoration time τ, wherein a fault is determined to be τ-active when the fault intersects a horizon H τ that was originally deposited at the past restoration time τ and a fault is determined to be τ-inactive when the fault does not intersect the horizon H τ that was originally deposited at the past restoration time τ; restoring the 3D model from the present day measured geometry to the predicted past geometry at the past restoration time τ by: modeling each τ-active fault to join end points of a horizon H τ separated on opposite sides of the fault in the present day model to merge into the same position in the restored model by sliding the end points towards each other in a direction tangential to the surface of the τ -active fault; and modeling each τ -inactive fault to keep collocated points on opposite sides of the fault together; and generating images of the interior subsurface of the Earth for the past restored geology based on said restored 3D model. receiving the 3D model of the present day measured geometry comprising a network of faults, wherein a fault is a discontinuity that divides fault blocks that slide in opposite directions tangential to the surface of the fault as time approaches a modeled time; selecting a past restoration time τ that is prior to the present time and after a time when an oldest horizon surface in the 3D model was originally deposited; the method characterized by further comprising: dividing the network of faults into a subset of τ-active faults and a subset of τ -inactive faults, wherein a τ-active fault is a fault that is active at the past restoration time τ and a τ-inactive fault is a fault that is inactive at the past restoration time τ, wherein a fault is determined to be τ-active when the fault intersects a horizon H τ that was originally deposited at the past restoration time τ and a fault is determined to be τ-inactive when the fault does not intersect the horizon H τ that was originally deposited at the past restoration time τ; restoring the 3D model from the present day measured geometry to the predicted past geometry at the past restoration time τ by: modeling each τ-active fault to join end points of a horizon H τ separated on opposite sides of the fault in the present day model to merge into the same position in the restored model by sliding the end points towards each other in a direction tangential to the surface of the τ -active fault; and modeling each τ -inactive fault to keep collocated points on opposite sides of the fault together; and generating images of the interior subsurface of the Earth for the past restored geology based on said restored 3D model. dividing the network of faults into a subset of τ-active faults and a subset of τ -inactive faults, wherein a τ-active fault is a fault that is active at the past restoration time τ and a τ-inactive fault is a fault that is inactive at the past restoration time τ, wherein a fault is determined to be τ-active when the fault intersects a horizon H τ that was originally deposited at the past restoration time τ and a fault is determined to be τ-inactive when the fault does not intersect the horizon H τ that was originally deposited at the past restoration time τ; restoring the 3D model from the present day measured geometry to the predicted past geometry at the past restoration time τ by: modeling each τ-active fault to join end points of a horizon H τ separated on opposite sides of the fault in the present day model to merge into the same position in the restored model by sliding the end points towards each other in a direction tangential to the surface of the τ -active fault; and modeling each τ -inactive fault to keep collocated points on opposite sides of the fault together; and generating images of the interior subsurface of the Earth for the past restored geology based on said restored 3D model. modeling each τ-active fault to join end points of a horizon H τ separated on opposite sides of the fault in the present day model to merge into the same position in the restored model by sliding the end points towards each other in a direction tangential to the surface of the τ -active fault; and modeling each τ -inactive fault to keep collocated points on opposite sides of the fault together; and generating images of the interior subsurface of the Earth for the past restored geology based on said restored 3D model.