convertTrajectory enhancements Handle inclination-only surveys

Handle inclination-only surveys (stations have no azimuth) for "inputKind": "MD_Incl".

Context:

Directional surveys (as measured with MWD or Gyro) generally have Measured Depth, Inclination and Azimuth observed at all survey stations, as indicated by "inputKind": "MD_Incl_Azim". This enhancement introduces a new "inputKind": "MD_Incl" for inclination-only surveys.
Inclination-only data, also known as INC-ONLY or TOTCO, are survey stations that do not have azimuth observables. The survey program may flag the segments that use this as error model, or missing azimuth can be detected to identify such surveys. In this enhancement, all stations in the survey (segment) must be of the same type (i.e., MD_Incl).
An issue in some subsurface software packages is that only MD_INC_AZI can be handled. It happens that for inclination-only surveys the AZI is set to "0", which then results in a calculated path that tracks due north.
To avoid this, sometimes the "INC" is then also set to "0" (replaced) so that the path will go vertically down. But then of course the original measurements are removed, and TVD is slightly off. The TVD error is approximately:
TVD_error = tan(inc_avg)^2*TD, where inc_avg is the average inclination. For inc_avg=2 deg (which is a lot), this is about 0.1% of TD (only 10ft at 10,000ft depth).

The horizontal uncertainty due to not knowing the azimuth can be large (but presumably acceptable for the purpose of the wellbore). If the wellbore would have a constant azimuth (e.g., "0"), then the maximum horizontal error compared to the calculated vertical path for a well with TD=10,000ft is approximately:

max_hor_offset = tan(inc_avg)*TD=350ft (3.5%).

One could consider this acceptable since it is an assumed-vertical well anyway, and presumably the inc-only survey program was selected to satisfy the objective of the wellbore.

Even though these may be small enough for assumed vertical wells, the correct path TVD (vertically) and a max. horizontal error will be computed and returned (the computed path is straight vertical).

Request inputStation will be like the following:

discuss first if "method" is needed. Note that azimuthReference is not there.

  "trajectoryCRS": "{{NAMESPACE}}:reference-data--CoordinateReferenceSystem:BoundProjected:EPSG::32066_EPSG::15851:",
  "referencePoint": {
    "x": 2000000,
    "y": 10000000,
    "z": 100
  },
  "unitZ": "osdu:reference-data--UnitOfMeasure:m:",
  "unitMD": "osdu:reference-data--UnitOfMeasure:m:",
  "inputKind": "MD_Incl",
  "inputStations": [
    {
      "md": 0,
      "inclination": 0.0
    },
    {
      "md": 1000,
      "inclination": 1.5
    },    
    {
      "md": 2000,
      "inclination": 0.5
    },
    {
      "md": 3000,
      "inclination": 1.0
    }
...

The azimuth For "inputKind": "MD_Incl" calculate a path that goes straight vertical, as follows:

Temporarily add "Azimuth": 0 for all stations.
Make a "dummy call" to do the computations as for "inputKind": "MD_Incl_Azim" (with the given method or simply use "method": "AzimuthalEquidistant")
Get the error stats/messages to be returned as part of "operationsApplied":
- max_horizontal_error = n[last]
- TVD_correction = MD[last] - TVD[last]
Create a final path to be output by forcing it to be perfectly vertical below the first station, but keep the TVD computed in step 2.

Alternatively, call minimum curvature or use direct equations (see documentation via Bert Kampes, e.g., approximation Δn1,2 = n2 – n1 = ΔM * sin(Iavg)). This seems simple to do inline vs. making a "dummy call". However, it may be better to make the "dummy call", because then the response body contains almost all outputs, and special cases like a geographic CRS as trajectory CRS or different units for unitMD and unitZ are automatically dealt with.

Example Response (where azimuth is n/a. It is TBD if it is returned or can be blank or omitted):

    "stations": [
        {
            "md": 0.0,
            "inclination": 0.0,
            "azimuthTN": ,
            "azimuthGN": ,
            "dxTN": 0.0,
            "dyTN": 0.0,
            "dz": 0.0
            "point": {
                "x": 2000000,
                "y": 10000000,
                "z": 100.0
            },
            "wgs84Longitude": -85.88980921169528,
            "wgs84Latitude": 27.553258329196616,
            "dls": 0.0,
            "original": true,
        },
        {
            "md": 1000.0,
            "inclination": 0.0,
            "azimuthTN": ,
            "azimuthGN": ,
            "dxTN": 0.0,
            "dyTN": 0.0,
            "dz": 999.89
            "point": {
                "x": 2000000,
                "y": 10000000,
                "z": -899.89
            },
            "wgs84Longitude": -85.88980921169528,
            "wgs84Latitude": 27.553258329196616,
            "dls": 0.0,
            "original": true,
        },
        {
            "md": 2000.0,
            "inclination": 0.0,
            "azimuthTN": ,
            "azimuthGN": ,
            "dxTN": 0.0,
            "dyTN": 0.0,
            "dz": 1999.72
            "point": {
                "x": 2000000,
                "y": 10000000,
                "z": -1899.72
            },
            "wgs84Longitude": -85.88980921169528,
            "wgs84Latitude": 27.553258329196616,
            "dls": 0.0,
            "original": true,
        },
        {
            "md": 3000.0,
            "inclination": 0.0,
            "azimuthTN": ,
            "azimuthGN": ,
            "dxTN": 0.0,
            "dyTN": 0.0,
            "dz": 2999.63
            "point": {
                "x": 2000000,
                "y": 10000000,
                "z": 2899.63
            },
            "wgs84Longitude": -85.88980921169528,
            "wgs84Latitude": 27.553258329196616,
            "dls": 0.0,
            "original": true,
        },
    "operationsApplied": [
        ...
        "computed deflections via minimum curvature method",
        "computation method: AzimuthalEquidistant",
        "inputKind: MD_Incl (Inclination-only: a vertical path is calculated."
        "max_horizontal_error = 43.63 m <unitXY>",   // subst. unitXY (of projCRS)
        "TVD_correction applied = 0.37 m (unitZ)",   // use unitZ.

That's it.

Test cases:

different unitMD=ft and unitZ=m
geographicCRS as input with unitXY=ftUS (it is irrelevant though since output dxTN=dyTN=0 by definition)
ensure Request azimuthReference and method are not required on input.
Discuss exact OperationsApplied messages.

Click to expand

Following notes can be ignored for this implementation.

An improvement on the above is to detect missing or empty or zero "azimuth" for "inputKind": "MD_Incl_Azim" (i.e., this to handle inc-only surveys incorrectly flagged as having Azimuths. This happens sometimes for definitive surveys that are stitched together from an inclination-only survey in the top hole and an MWD in the deeper part).

The algorithm then for "inputKind": "MD_Incl_Azim" should detect if stations 1:N (where N>2 and may be all stations) are all missing azimuth values or all have "azimuth": 0. Then it should follow the algorithm as shown above, which would create a straight vertical path for stations 1:N. The path for stations N+1 to "last" is computed as normal and tied to station N.

For a survey that is deviated but is missing an azimuth value at a single station in the middle of a deviated part (for example) it is in general undesired to make it drop "straight down", because the inclination may be significant (i.e., this may be because the measurement is corrupted/not received with mud pump telemetry). It may be best to remove the station or to interpolate the azimuth at this MD, e.g., by copying the previous AZI. Therefore this algorithm will not check for this case.
(The better way of handling this is with a survey program and error model per station; or interpolating or even removing the station from the survey, but such sophisticated input is not available). In future if a survey object or witsml file is input then this could be done in OSDU).

Edited Jul 17, 2023 by Bert Kampes