Different codes may have different definitions, but the following is probably typical.
For a 4 node shell element (quad), compute the mid-points of each side. Connect the mid-points on opposite sides with vectors. These 2 vectors are in general not perpendicular to each other, but they do define a plane which is taken to be the element x-y plane. The cross product of those 2 vectors defines the element normal (element z-axis). One of the in-plane vectors becomes the element x-axis. Finally, cross the x-vector with the z-vector to get the orthogonal y-axis. (The second mid-point vector was just an intermediate step).
Also note that for non-planar quad elements, the 4 corner nodes will not lie in the average mid-plane as defined above. 2 of the corner nodes will be above that plane by a certain distance, and 2 of them will be below by the same distance. This defines the "warping" of the element. Too much warp is bad. To minimize warping, it is best to define the mesh so that it "follows" the surface contours. For example, if you mesh a cylindrical surface with rows of nodes running axially and hoop-wise, the elements will be flat (no warp). But if the node lines run +/-45 degrees to the axial direction, the quads will be warped.